1 : /******************************************************************************
2 : * $Id: idadataset.cpp 25311 2012-12-15 12:48:14Z rouault $
3 : *
4 : * Project: IDA Raster Driver
5 : * Purpose: Implemenents IDA driver/dataset/rasterband.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2004, Frank Warmerdam <warmerdam@pobox.com>
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "rawdataset.h"
31 : #include "ogr_spatialref.h"
32 : #include "gdal_rat.h"
33 :
34 : CPL_CVSID("$Id: idadataset.cpp 25311 2012-12-15 12:48:14Z rouault $");
35 :
36 : CPL_C_START
37 : void GDALRegister_IDA(void);
38 : CPL_C_END
39 :
40 : // convert a Turbo Pascal real into a double
41 : static double tp2c(GByte *r);
42 :
43 : // convert a double into a Turbo Pascal real
44 : static void c2tp(double n, GByte *r);
45 :
46 : /************************************************************************/
47 : /* ==================================================================== */
48 : /* IDADataset */
49 : /* ==================================================================== */
50 : /************************************************************************/
51 :
52 : class IDADataset : public RawDataset
53 : {
54 : friend class IDARasterBand;
55 :
56 : int nImageType;
57 : int nProjection;
58 : char szTitle[81];
59 : double dfLatCenter;
60 : double dfLongCenter;
61 : double dfXCenter;
62 : double dfYCenter;
63 : double dfDX;
64 : double dfDY;
65 : double dfParallel1;
66 : double dfParallel2;
67 : int nMissing;
68 : double dfM;
69 : double dfB;
70 :
71 : FILE *fpRaw;
72 :
73 : char *pszProjection;
74 : double adfGeoTransform[6];
75 :
76 : void ProcessGeoref();
77 :
78 : GByte abyHeader[512];
79 : int bHeaderDirty;
80 :
81 : void ReadColorTable();
82 :
83 : public:
84 : IDADataset();
85 : ~IDADataset();
86 :
87 : virtual void FlushCache();
88 : virtual const char *GetProjectionRef(void);
89 : virtual CPLErr SetProjection( const char * );
90 :
91 : virtual CPLErr GetGeoTransform( double * );
92 : virtual CPLErr SetGeoTransform( double * );
93 :
94 : static GDALDataset *Open( GDALOpenInfo * );
95 : static GDALDataset *Create( const char * pszFilename,
96 : int nXSize, int nYSize, int nBands,
97 : GDALDataType eType,
98 : char ** /* papszParmList */ );
99 :
100 : };
101 :
102 : /************************************************************************/
103 : /* ==================================================================== */
104 : /* IDARasterBand */
105 : /* ==================================================================== */
106 : /************************************************************************/
107 :
108 : class IDARasterBand : public RawRasterBand
109 : {
110 : friend class IDADataset;
111 :
112 : GDALRasterAttributeTable *poRAT;
113 : GDALColorTable *poColorTable;
114 :
115 : public:
116 : IDARasterBand( IDADataset *poDSIn, FILE *fpRaw, int nXSize );
117 : virtual ~IDARasterBand();
118 :
119 : virtual const GDALRasterAttributeTable *GetDefaultRAT();
120 : virtual GDALColorInterp GetColorInterpretation();
121 : virtual GDALColorTable *GetColorTable();
122 : virtual double GetOffset( int *pbSuccess = NULL );
123 : virtual CPLErr SetOffset( double dfNewValue );
124 : virtual double GetScale( int *pbSuccess = NULL );
125 : virtual CPLErr SetScale( double dfNewValue );
126 : virtual double GetNoDataValue( int *pbSuccess = NULL );
127 : };
128 :
129 : /************************************************************************/
130 : /* IDARasterBand */
131 : /************************************************************************/
132 :
133 25 : IDARasterBand::IDARasterBand( IDADataset *poDSIn,
134 : FILE *fpRaw, int nXSize )
135 : : RawRasterBand( poDSIn, 1, fpRaw, 512, 1, nXSize,
136 25 : GDT_Byte, FALSE, FALSE )
137 :
138 : {
139 25 : poColorTable = NULL;
140 25 : poRAT = NULL;
141 25 : }
142 :
143 : /************************************************************************/
144 : /* ~IDARasterBand() */
145 : /************************************************************************/
146 :
147 25 : IDARasterBand::~IDARasterBand()
148 :
149 : {
150 25 : delete poColorTable;
151 25 : delete poRAT;
152 25 : }
153 :
154 : /************************************************************************/
155 : /* GetNoDataValue() */
156 : /************************************************************************/
157 :
158 7 : double IDARasterBand::GetNoDataValue( int *pbSuccess )
159 :
160 : {
161 7 : if( pbSuccess != NULL )
162 6 : *pbSuccess = TRUE;
163 7 : return ((IDADataset *) poDS)->nMissing;
164 : }
165 :
166 : /************************************************************************/
167 : /* GetOffset() */
168 : /************************************************************************/
169 :
170 1 : double IDARasterBand::GetOffset( int *pbSuccess )
171 :
172 : {
173 1 : if( pbSuccess != NULL )
174 1 : *pbSuccess = TRUE;
175 1 : return ((IDADataset *) poDS)->dfB;
176 : }
177 :
178 : /************************************************************************/
179 : /* SetOffset() */
180 : /************************************************************************/
181 :
182 0 : CPLErr IDARasterBand::SetOffset( double dfNewValue )
183 :
184 : {
185 0 : IDADataset *poIDS = (IDADataset *) poDS;
186 :
187 0 : if( dfNewValue == poIDS->dfB )
188 0 : return CE_None;
189 :
190 0 : if( poIDS->nImageType != 200 )
191 : {
192 : CPLError( CE_Failure, CPLE_AppDefined,
193 0 : "Setting explicit offset only support for image type 200.");
194 0 : return CE_Failure;
195 : }
196 :
197 0 : poIDS->dfB = dfNewValue;
198 0 : c2tp( dfNewValue, poIDS->abyHeader + 177 );
199 0 : poIDS->bHeaderDirty = TRUE;
200 :
201 0 : return CE_None;
202 : }
203 :
204 : /************************************************************************/
205 : /* GetScale() */
206 : /************************************************************************/
207 :
208 1 : double IDARasterBand::GetScale( int *pbSuccess )
209 :
210 : {
211 1 : if( pbSuccess != NULL )
212 1 : *pbSuccess = TRUE;
213 1 : return ((IDADataset *) poDS)->dfM;
214 : }
215 :
216 : /************************************************************************/
217 : /* SetScale() */
218 : /************************************************************************/
219 :
220 1 : CPLErr IDARasterBand::SetScale( double dfNewValue )
221 :
222 : {
223 1 : IDADataset *poIDS = (IDADataset *) poDS;
224 :
225 1 : if( dfNewValue == poIDS->dfM )
226 0 : return CE_None;
227 :
228 1 : if( poIDS->nImageType != 200 )
229 : {
230 : CPLError( CE_Failure, CPLE_AppDefined,
231 0 : "Setting explicit scale only support for image type 200.");
232 0 : return CE_Failure;
233 : }
234 :
235 1 : poIDS->dfM = dfNewValue;
236 1 : c2tp( dfNewValue, poIDS->abyHeader + 171 );
237 1 : poIDS->bHeaderDirty = TRUE;
238 :
239 1 : return CE_None;
240 : }
241 :
242 : /************************************************************************/
243 : /* GetColorTable() */
244 : /************************************************************************/
245 :
246 1 : GDALColorTable *IDARasterBand::GetColorTable()
247 :
248 : {
249 1 : if( poColorTable )
250 0 : return poColorTable;
251 : else
252 1 : return RawRasterBand::GetColorTable();
253 : }
254 :
255 : /************************************************************************/
256 : /* GetColorInterpretation() */
257 : /************************************************************************/
258 :
259 2 : GDALColorInterp IDARasterBand::GetColorInterpretation()
260 :
261 : {
262 2 : if( poColorTable )
263 0 : return GCI_PaletteIndex;
264 : else
265 2 : return RawRasterBand::GetColorInterpretation();
266 : }
267 :
268 : /************************************************************************/
269 : /* GetDefaultRAT() */
270 : /************************************************************************/
271 :
272 0 : const GDALRasterAttributeTable *IDARasterBand::GetDefaultRAT()
273 :
274 : {
275 0 : if( poRAT )
276 0 : return poRAT;
277 : else
278 0 : return RawRasterBand::GetDefaultRAT();
279 : }
280 :
281 : /************************************************************************/
282 : /* ==================================================================== */
283 : /* IDADataset */
284 : /* ==================================================================== */
285 : /************************************************************************/
286 :
287 : /************************************************************************/
288 : /* IDADataset() */
289 : /************************************************************************/
290 :
291 25 : IDADataset::IDADataset()
292 : {
293 25 : fpRaw = NULL;
294 25 : pszProjection = NULL;
295 :
296 25 : adfGeoTransform[0] = 0.0;
297 25 : adfGeoTransform[1] = 1.0;
298 25 : adfGeoTransform[2] = 0.0;
299 25 : adfGeoTransform[3] = 0.0;
300 25 : adfGeoTransform[4] = 0.0;
301 25 : adfGeoTransform[5] = 1.0;
302 :
303 25 : bHeaderDirty = FALSE;
304 25 : }
305 :
306 : /************************************************************************/
307 : /* ~IDADataset() */
308 : /************************************************************************/
309 :
310 25 : IDADataset::~IDADataset()
311 :
312 : {
313 25 : FlushCache();
314 :
315 25 : if( fpRaw != NULL )
316 25 : VSIFClose( fpRaw );
317 25 : CPLFree( pszProjection );
318 25 : }
319 :
320 : /************************************************************************/
321 : /* ProcessGeoref() */
322 : /************************************************************************/
323 :
324 25 : void IDADataset::ProcessGeoref()
325 :
326 : {
327 25 : OGRSpatialReference oSRS;
328 :
329 25 : if( nProjection == 3 )
330 : {
331 10 : oSRS.SetWellKnownGeogCS( "WGS84" );
332 : }
333 15 : else if( nProjection == 4 )
334 : {
335 : oSRS.SetLCC( dfParallel1, dfParallel2,
336 : dfLatCenter, dfLongCenter,
337 2 : 0.0, 0.0 );
338 : oSRS.SetGeogCS( "Clarke 1866", "Clarke 1866", "Clarke 1866",
339 2 : 6378206.4, 293.97869821389662 );
340 : }
341 13 : else if( nProjection == 6 )
342 : {
343 2 : oSRS.SetLAEA( dfLatCenter, dfLongCenter, 0.0, 0.0 );
344 : oSRS.SetGeogCS( "Sphere", "Sphere", "Sphere",
345 2 : 6370997.0, 0.0 );
346 : }
347 11 : else if( nProjection == 8 )
348 : {
349 : oSRS.SetACEA( dfParallel1, dfParallel2,
350 : dfLatCenter, dfLongCenter,
351 2 : 0.0, 0.0 );
352 : oSRS.SetGeogCS( "Clarke 1866", "Clarke 1866", "Clarke 1866",
353 2 : 6378206.4, 293.97869821389662 );
354 : }
355 9 : else if( nProjection == 9 )
356 : {
357 2 : oSRS.SetGH( dfLongCenter, 0.0, 0.0 );
358 : oSRS.SetGeogCS( "Sphere", "Sphere", "Sphere",
359 2 : 6370997.0, 0.0 );
360 : }
361 :
362 25 : if( oSRS.GetRoot() != NULL )
363 : {
364 18 : CPLFree( pszProjection );
365 18 : pszProjection = NULL;
366 :
367 18 : oSRS.exportToWkt( &pszProjection );
368 : }
369 :
370 25 : adfGeoTransform[0] = 0 - dfDX * dfXCenter;
371 25 : adfGeoTransform[1] = dfDX;
372 25 : adfGeoTransform[2] = 0.0;
373 25 : adfGeoTransform[3] = dfDY * dfYCenter;
374 25 : adfGeoTransform[4] = 0.0;
375 25 : adfGeoTransform[5] = -dfDY;
376 :
377 25 : if( nProjection == 3 )
378 : {
379 10 : adfGeoTransform[0] += dfLongCenter;
380 10 : adfGeoTransform[3] += dfLatCenter;
381 25 : }
382 25 : }
383 :
384 : /************************************************************************/
385 : /* FlushCache() */
386 : /************************************************************************/
387 :
388 25 : void IDADataset::FlushCache()
389 :
390 : {
391 25 : RawDataset::FlushCache();
392 :
393 25 : if( bHeaderDirty )
394 : {
395 7 : VSIFSeek( fpRaw, 0, SEEK_SET );
396 7 : VSIFWrite( abyHeader, 512, 1, fpRaw );
397 7 : bHeaderDirty = FALSE;
398 : }
399 25 : }
400 :
401 : /************************************************************************/
402 : /* GetGeoTransform() */
403 : /************************************************************************/
404 :
405 4 : CPLErr IDADataset::GetGeoTransform( double *padfGeoTransform )
406 :
407 : {
408 4 : memcpy( padfGeoTransform, adfGeoTransform, sizeof(double) * 6 );
409 4 : return CE_None;
410 : }
411 :
412 : /************************************************************************/
413 : /* SetGeoTransform() */
414 : /************************************************************************/
415 :
416 7 : CPLErr IDADataset::SetGeoTransform( double *padfGeoTransform )
417 :
418 : {
419 7 : if( padfGeoTransform[2] != 0.0 || padfGeoTransform[4] != 0.0 )
420 0 : return GDALPamDataset::SetGeoTransform( padfGeoTransform );
421 :
422 7 : memcpy( adfGeoTransform, padfGeoTransform, sizeof(double) * 6 );
423 7 : bHeaderDirty = TRUE;
424 :
425 7 : dfDX = adfGeoTransform[1];
426 7 : dfDY = -adfGeoTransform[5];
427 7 : dfXCenter = -adfGeoTransform[0] / dfDX;
428 7 : dfYCenter = adfGeoTransform[3] / dfDY;
429 :
430 7 : c2tp( dfDX, abyHeader + 144 );
431 7 : c2tp( dfDY, abyHeader + 150 );
432 7 : c2tp( dfXCenter, abyHeader + 132 );
433 7 : c2tp( dfYCenter, abyHeader + 138 );
434 :
435 7 : return CE_None;
436 : }
437 :
438 : /************************************************************************/
439 : /* GetProjectionRef() */
440 : /************************************************************************/
441 :
442 10 : const char *IDADataset::GetProjectionRef()
443 :
444 : {
445 10 : if( pszProjection )
446 10 : return pszProjection;
447 : else
448 0 : return "";
449 : }
450 :
451 : /************************************************************************/
452 : /* SetProjection() */
453 : /************************************************************************/
454 :
455 7 : CPLErr IDADataset::SetProjection( const char *pszWKTIn )
456 :
457 : {
458 7 : OGRSpatialReference oSRS;
459 :
460 7 : oSRS.importFromWkt( (char **) &pszWKTIn );
461 :
462 7 : if( !oSRS.IsGeographic() && !oSRS.IsProjected() )
463 0 : GDALPamDataset::SetProjection( pszWKTIn );
464 :
465 : /* -------------------------------------------------------------------- */
466 : /* Clear projection parameters. */
467 : /* -------------------------------------------------------------------- */
468 7 : dfParallel1 = 0.0;
469 7 : dfParallel2 = 0.0;
470 7 : dfLatCenter = 0.0;
471 7 : dfLongCenter = 0.0;
472 :
473 : /* -------------------------------------------------------------------- */
474 : /* Geographic. */
475 : /* -------------------------------------------------------------------- */
476 7 : if( oSRS.IsGeographic() )
477 : {
478 : // If no change, just return.
479 3 : if( nProjection == 3 )
480 0 : return CE_None;
481 :
482 3 : nProjection = 3;
483 : }
484 :
485 : /* -------------------------------------------------------------------- */
486 : /* Verify we don't have a false easting or northing as these */
487 : /* will be ignored for the projections we do support. */
488 : /* -------------------------------------------------------------------- */
489 7 : if( oSRS.GetProjParm( SRS_PP_FALSE_EASTING ) != 0.0
490 : || oSRS.GetProjParm( SRS_PP_FALSE_NORTHING ) != 0.0 )
491 : {
492 : CPLError( CE_Failure, CPLE_AppDefined,
493 : "Attempt to set a projection on an IDA file with a non-zero\n"
494 0 : "false easting and/or northing. This is not supported." );
495 0 : return CE_Failure;
496 : }
497 :
498 : /* -------------------------------------------------------------------- */
499 : /* Lambert Conformal Conic. Note that we don't support false */
500 : /* eastings or nothings. */
501 : /* -------------------------------------------------------------------- */
502 7 : const char *pszProjection = oSRS.GetAttrValue( "PROJECTION" );
503 :
504 7 : if( pszProjection == NULL )
505 : {
506 : /* do nothing - presumably geographic */;
507 : }
508 4 : else if( EQUAL(pszProjection,SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP) )
509 : {
510 1 : nProjection = 4;
511 1 : dfParallel1 = oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0);
512 1 : dfParallel2 = oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2,0.0);
513 1 : dfLatCenter = oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0);
514 1 : dfLongCenter = oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0);
515 : }
516 3 : else if( EQUAL(pszProjection,SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA) )
517 : {
518 1 : nProjection = 6;
519 1 : dfLatCenter = oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0);
520 1 : dfLongCenter = oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0);
521 : }
522 2 : else if( EQUAL(pszProjection,SRS_PT_ALBERS_CONIC_EQUAL_AREA) )
523 : {
524 1 : nProjection = 8;
525 1 : dfParallel1 = oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0);
526 1 : dfParallel2 = oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2,0.0);
527 1 : dfLatCenter = oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0);
528 1 : dfLongCenter = oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0);
529 : }
530 1 : else if( EQUAL(pszProjection,SRS_PT_GOODE_HOMOLOSINE) )
531 : {
532 1 : nProjection = 9;
533 1 : dfLongCenter = oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0);
534 : }
535 : else
536 : {
537 0 : return GDALPamDataset::SetProjection( pszWKTIn );
538 : }
539 :
540 : /* -------------------------------------------------------------------- */
541 : /* Update header and mark it as dirty. */
542 : /* -------------------------------------------------------------------- */
543 7 : bHeaderDirty = TRUE;
544 :
545 7 : abyHeader[23] = (GByte) nProjection;
546 7 : c2tp( dfLatCenter, abyHeader + 120 );
547 7 : c2tp( dfLongCenter, abyHeader + 126 );
548 7 : c2tp( dfParallel1, abyHeader + 156 );
549 7 : c2tp( dfParallel2, abyHeader + 162 );
550 :
551 7 : return CE_None;
552 : }
553 :
554 : /************************************************************************/
555 : /* ReadColorTable() */
556 : /************************************************************************/
557 :
558 25 : void IDADataset::ReadColorTable()
559 :
560 : {
561 : /* -------------------------------------------------------------------- */
562 : /* Decide what .clr file to look for and try to open. */
563 : /* -------------------------------------------------------------------- */
564 25 : CPLString osCLRFilename;
565 :
566 25 : osCLRFilename = CPLGetConfigOption( "IDA_COLOR_FILE", "" );
567 25 : if( strlen(osCLRFilename) == 0 )
568 25 : osCLRFilename = CPLResetExtension(GetDescription(), "clr" );
569 :
570 :
571 25 : FILE *fp = VSIFOpen( osCLRFilename, "r" );
572 25 : if( fp == NULL )
573 : {
574 25 : osCLRFilename = CPLResetExtension(osCLRFilename, "CLR" );
575 25 : fp = VSIFOpen( osCLRFilename, "r" );
576 : }
577 :
578 25 : if( fp == NULL )
579 : return;
580 :
581 : /* -------------------------------------------------------------------- */
582 : /* Skip first line, with the column titles. */
583 : /* -------------------------------------------------------------------- */
584 0 : CPLReadLine( fp );
585 :
586 : /* -------------------------------------------------------------------- */
587 : /* Create a RAT to populate. */
588 : /* -------------------------------------------------------------------- */
589 0 : GDALRasterAttributeTable *poRAT = new GDALRasterAttributeTable();
590 :
591 0 : poRAT->CreateColumn( "FROM", GFT_Integer, GFU_Min );
592 0 : poRAT->CreateColumn( "TO", GFT_Integer, GFU_Max );
593 0 : poRAT->CreateColumn( "RED", GFT_Integer, GFU_Red );
594 0 : poRAT->CreateColumn( "GREEN", GFT_Integer, GFU_Green );
595 0 : poRAT->CreateColumn( "BLUE", GFT_Integer, GFU_Blue );
596 0 : poRAT->CreateColumn( "LEGEND", GFT_String, GFU_Name );
597 :
598 : /* -------------------------------------------------------------------- */
599 : /* Apply lines. */
600 : /* -------------------------------------------------------------------- */
601 0 : const char *pszLine = CPLReadLine( fp );
602 0 : int iRow = 0;
603 :
604 0 : while( pszLine != NULL )
605 : {
606 : char **papszTokens =
607 0 : CSLTokenizeStringComplex( pszLine, " \t", FALSE, FALSE );
608 :
609 0 : if( CSLCount( papszTokens ) >= 5 )
610 : {
611 0 : poRAT->SetValue( iRow, 0, atoi(papszTokens[0]) );
612 0 : poRAT->SetValue( iRow, 1, atoi(papszTokens[1]) );
613 0 : poRAT->SetValue( iRow, 2, atoi(papszTokens[2]) );
614 0 : poRAT->SetValue( iRow, 3, atoi(papszTokens[3]) );
615 0 : poRAT->SetValue( iRow, 4, atoi(papszTokens[4]) );
616 :
617 : // find name, first nonspace after 5th token.
618 0 : const char *pszName = pszLine;
619 :
620 : // skip from
621 0 : while( *pszName == ' ' || *pszName == '\t' )
622 0 : pszName++;
623 0 : while( *pszName != ' ' && *pszName != '\t' && *pszName != '\0' )
624 0 : pszName++;
625 :
626 : // skip to
627 0 : while( *pszName == ' ' || *pszName == '\t' )
628 0 : pszName++;
629 0 : while( *pszName != ' ' && *pszName != '\t' && *pszName != '\0' )
630 0 : pszName++;
631 :
632 : // skip red
633 0 : while( *pszName == ' ' || *pszName == '\t' )
634 0 : pszName++;
635 0 : while( *pszName != ' ' && *pszName != '\t' && *pszName != '\0' )
636 0 : pszName++;
637 :
638 : // skip green
639 0 : while( *pszName == ' ' || *pszName == '\t' )
640 0 : pszName++;
641 0 : while( *pszName != ' ' && *pszName != '\t' && *pszName != '\0' )
642 0 : pszName++;
643 :
644 : // skip blue
645 0 : while( *pszName == ' ' || *pszName == '\t' )
646 0 : pszName++;
647 0 : while( *pszName != ' ' && *pszName != '\t' && *pszName != '\0' )
648 0 : pszName++;
649 :
650 : // skip pre-name white space
651 0 : while( *pszName == ' ' || *pszName == '\t' )
652 0 : pszName++;
653 :
654 0 : poRAT->SetValue( iRow, 5, pszName );
655 :
656 0 : iRow++;
657 : }
658 :
659 0 : CSLDestroy( papszTokens );
660 0 : pszLine = CPLReadLine( fp );
661 : }
662 :
663 0 : VSIFClose( fp );
664 :
665 : /* -------------------------------------------------------------------- */
666 : /* Attach RAT to band. */
667 : /* -------------------------------------------------------------------- */
668 0 : ((IDARasterBand *) GetRasterBand( 1 ))->poRAT = poRAT;
669 :
670 : /* -------------------------------------------------------------------- */
671 : /* Build a conventional color table from this. */
672 : /* -------------------------------------------------------------------- */
673 : ((IDARasterBand *) GetRasterBand( 1 ))->poColorTable =
674 0 : poRAT->TranslateToColorTable();
675 : }
676 :
677 : /************************************************************************/
678 : /* Open() */
679 : /************************************************************************/
680 :
681 11981 : GDALDataset *IDADataset::Open( GDALOpenInfo * poOpenInfo )
682 :
683 : {
684 : /* -------------------------------------------------------------------- */
685 : /* Is this an IDA file? */
686 : /* -------------------------------------------------------------------- */
687 : int nXSize, nYSize;
688 : GIntBig nExpectedFileSize, nActualFileSize;
689 :
690 11981 : if( poOpenInfo->fp == NULL )
691 11308 : return NULL;
692 :
693 673 : if( poOpenInfo->nHeaderBytes < 512 )
694 180 : return NULL;
695 :
696 : // projection legal?
697 493 : if( poOpenInfo->pabyHeader[23] > 10 )
698 388 : return NULL;
699 :
700 : // imagetype legal?
701 278 : if( (poOpenInfo->pabyHeader[22] > 14
702 34 : && poOpenInfo->pabyHeader[22] < 100)
703 105 : || (poOpenInfo->pabyHeader[22] > 114
704 34 : && poOpenInfo->pabyHeader[22] != 200 ) )
705 8 : return NULL;
706 :
707 97 : nXSize = poOpenInfo->pabyHeader[30] + poOpenInfo->pabyHeader[31] * 256;
708 97 : nYSize = poOpenInfo->pabyHeader[32] + poOpenInfo->pabyHeader[33] * 256;
709 :
710 97 : if( nXSize == 0 || nYSize == 0 )
711 70 : return NULL;
712 :
713 : // The file just be exactly the image size + header size in length.
714 27 : nExpectedFileSize = nXSize * nYSize + 512;
715 :
716 27 : VSIFSeek( poOpenInfo->fp, 0, SEEK_END );
717 27 : nActualFileSize = VSIFTell( poOpenInfo->fp );
718 27 : VSIRewind( poOpenInfo->fp );
719 :
720 27 : if( nActualFileSize != nExpectedFileSize )
721 2 : return NULL;
722 :
723 : /* -------------------------------------------------------------------- */
724 : /* Create the dataset. */
725 : /* -------------------------------------------------------------------- */
726 25 : IDADataset *poDS = new IDADataset();
727 :
728 25 : memcpy( poDS->abyHeader, poOpenInfo->pabyHeader, 512 );
729 :
730 : /* -------------------------------------------------------------------- */
731 : /* Parse various values out of the header. */
732 : /* -------------------------------------------------------------------- */
733 25 : poDS->nImageType = poOpenInfo->pabyHeader[22];
734 25 : poDS->nProjection = poOpenInfo->pabyHeader[23];
735 :
736 25 : poDS->nRasterYSize = poOpenInfo->pabyHeader[30]
737 25 : + poOpenInfo->pabyHeader[31] * 256;
738 25 : poDS->nRasterXSize = poOpenInfo->pabyHeader[32]
739 25 : + poOpenInfo->pabyHeader[33] * 256;
740 :
741 25 : strncpy( poDS->szTitle, (const char *) poOpenInfo->pabyHeader+38, 80 );
742 25 : poDS->szTitle[80] = '\0';
743 :
744 25 : int nLastTitleChar = strlen(poDS->szTitle)-1;
745 1527 : while( nLastTitleChar > -1
746 371 : && (poDS->szTitle[nLastTitleChar] == 10
747 371 : || poDS->szTitle[nLastTitleChar] == 13
748 371 : || poDS->szTitle[nLastTitleChar] == ' ') )
749 364 : poDS->szTitle[nLastTitleChar--] = '\0';
750 :
751 25 : poDS->dfLatCenter = tp2c( poOpenInfo->pabyHeader + 120 );
752 25 : poDS->dfLongCenter = tp2c( poOpenInfo->pabyHeader + 126 );
753 25 : poDS->dfXCenter = tp2c( poOpenInfo->pabyHeader + 132 );
754 25 : poDS->dfYCenter = tp2c( poOpenInfo->pabyHeader + 138 );
755 25 : poDS->dfDX = tp2c( poOpenInfo->pabyHeader + 144 );
756 25 : poDS->dfDY = tp2c( poOpenInfo->pabyHeader + 150 );
757 25 : poDS->dfParallel1 = tp2c( poOpenInfo->pabyHeader + 156 );
758 25 : poDS->dfParallel2 = tp2c( poOpenInfo->pabyHeader + 162 );
759 :
760 25 : poDS->ProcessGeoref();
761 :
762 25 : poDS->SetMetadataItem( "TITLE", poDS->szTitle );
763 :
764 : /* -------------------------------------------------------------------- */
765 : /* Handle various image types. */
766 : /* -------------------------------------------------------------------- */
767 :
768 : /*
769 : GENERIC = 0
770 : FEW S NDVI = 1
771 : EROS NDVI = 6
772 : ARTEMIS CUTOFF = 10
773 : ARTEMIS RECODE = 11
774 : ARTEMIS NDVI = 12
775 : ARTEMIS FEWS = 13
776 : ARTEMIS NEWNASA = 14
777 : GENERIC DIFF = 100
778 : FEW S NDVI DIFF = 101
779 : EROS NDVI DIFF = 106
780 : ARTEMIS CUTOFF DIFF = 110
781 : ARTEMIS RECODE DIFF = 111
782 : ARTEMIS NDVI DIFF = 112
783 : ARTEMIS FEWS DIFF = 113
784 : ARTEMIS NEWNASA DIFF = 114
785 : CALCULATED =200
786 : */
787 :
788 25 : poDS->nMissing = 0;
789 :
790 25 : switch( poDS->nImageType )
791 : {
792 : case 1:
793 0 : poDS->SetMetadataItem( "IMAGETYPE", "1, FEWS NDVI" );
794 0 : poDS->dfM = 1/256.0;
795 0 : poDS->dfB = -82/256.0;
796 0 : break;
797 :
798 : case 6:
799 0 : poDS->SetMetadataItem( "IMAGETYPE", "6, EROS NDVI" );
800 0 : poDS->dfM = 1/100.0;
801 0 : poDS->dfB = -100/100.0;
802 0 : break;
803 :
804 : case 10:
805 0 : poDS->SetMetadataItem( "IMAGETYPE", "10, ARTEMIS CUTOFF" );
806 0 : poDS->dfM = 1.0;
807 0 : poDS->dfB = 0.0;
808 0 : poDS->nMissing = 254;
809 0 : break;
810 :
811 : case 11:
812 0 : poDS->SetMetadataItem( "IMAGETYPE", "11, ARTEMIS RECODE" );
813 0 : poDS->dfM = 4.0;
814 0 : poDS->dfB = 0.0;
815 0 : poDS->nMissing = 254;
816 0 : break;
817 :
818 : case 12: /* ANDVI */
819 0 : poDS->SetMetadataItem( "IMAGETYPE", "12, ARTEMIS NDVI" );
820 0 : poDS->dfM = 4/500.0;
821 0 : poDS->dfB = -3/500.0 - 1.0;
822 0 : poDS->nMissing = 254;
823 0 : break;
824 :
825 : case 13: /* AFEWS */
826 0 : poDS->SetMetadataItem( "IMAGETYPE", "13, ARTEMIS FEWS" );
827 0 : poDS->dfM = 1/256.0;
828 0 : poDS->dfB = -82/256.0;
829 0 : poDS->nMissing = 254;
830 0 : break;
831 :
832 : case 14: /* NEWNASA */
833 0 : poDS->SetMetadataItem( "IMAGETYPE", "13, ARTEMIS NEWNASA" );
834 0 : poDS->dfM = 0.75/250.0;
835 0 : poDS->dfB = 0.0;
836 0 : poDS->nMissing = 254;
837 0 : break;
838 :
839 : case 101: /* NDVI_DIFF (FEW S) */
840 0 : poDS->dfM = 1/128.0;
841 0 : poDS->dfB = -1.0;
842 0 : poDS->nMissing = 0;
843 0 : break;
844 :
845 : case 106: /* EROS_DIFF */
846 0 : poDS->dfM = 1/50.0;
847 0 : poDS->dfB = -128/50.0;
848 0 : poDS->nMissing = 0;
849 0 : break;
850 :
851 : case 110: /* CUTOFF_DIFF */
852 0 : poDS->dfM = 2.0;
853 0 : poDS->dfB = -128*2;
854 0 : poDS->nMissing = 254;
855 0 : break;
856 :
857 : case 111: /* RECODE_DIFF */
858 0 : poDS->dfM = 8;
859 0 : poDS->dfB = -128*8;
860 0 : poDS->nMissing = 254;
861 0 : break;
862 :
863 : case 112: /* ANDVI_DIFF */
864 0 : poDS->dfM = 8/1000.0;
865 0 : poDS->dfB = (-128*8)/1000.0;
866 0 : poDS->nMissing = 254;
867 0 : break;
868 :
869 : case 113: /* AFEWS_DIFF */
870 0 : poDS->dfM = 1/128.0;
871 0 : poDS->dfB = -1;
872 0 : poDS->nMissing = 254;
873 0 : break;
874 :
875 : case 114: /* NEWNASA_DIFF */
876 0 : poDS->dfM = 0.75/125.0;
877 0 : poDS->dfB = -128*poDS->dfM;
878 0 : poDS->nMissing = 254;
879 0 : break;
880 :
881 : case 200:
882 : /* we use the values from the header */
883 25 : poDS->dfM = tp2c( poOpenInfo->pabyHeader + 171 );
884 25 : poDS->dfB = tp2c( poOpenInfo->pabyHeader + 177 );
885 25 : poDS->nMissing = poOpenInfo->pabyHeader[170];
886 25 : break;
887 :
888 : default:
889 0 : poDS->dfM = 1.0;
890 0 : poDS->dfB = 0.0;
891 : break;
892 : }
893 :
894 : /* -------------------------------------------------------------------- */
895 : /* Create the band. */
896 : /* -------------------------------------------------------------------- */
897 25 : if( poOpenInfo->eAccess == GA_ReadOnly )
898 : {
899 18 : poDS->fpRaw = poOpenInfo->fp;
900 18 : poOpenInfo->fp = NULL;
901 : }
902 : else
903 : {
904 7 : poDS->fpRaw = VSIFOpen( poOpenInfo->pszFilename, "rb+" );
905 7 : poDS->eAccess = GA_Update;
906 7 : if( poDS->fpRaw == NULL )
907 : {
908 : CPLError( CE_Failure, CPLE_OpenFailed,
909 : "Failed to open %s for write access.",
910 0 : poOpenInfo->pszFilename );
911 0 : return NULL;
912 : }
913 : }
914 :
915 : poDS->SetBand( 1, new IDARasterBand( poDS, poDS->fpRaw,
916 25 : poDS->nRasterXSize ) );
917 :
918 : /* -------------------------------------------------------------------- */
919 : /* Check for a color table. */
920 : /* -------------------------------------------------------------------- */
921 25 : poDS->SetDescription( poOpenInfo->pszFilename );
922 25 : poDS->ReadColorTable();
923 :
924 : /* -------------------------------------------------------------------- */
925 : /* Initialize any PAM information. */
926 : /* -------------------------------------------------------------------- */
927 25 : poDS->TryLoadXML();
928 :
929 : /* -------------------------------------------------------------------- */
930 : /* Check for overviews. */
931 : /* -------------------------------------------------------------------- */
932 25 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
933 :
934 25 : return( poDS );
935 : }
936 :
937 : /************************************************************************/
938 : /* tp2c() */
939 : /* */
940 : /* convert a Turbo Pascal real into a double */
941 : /************************************************************************/
942 :
943 250 : static double tp2c(GByte *r)
944 : {
945 : double mant;
946 : int sign, exp, i;
947 :
948 : // handle 0 case
949 250 : if (r[0] == 0)
950 105 : return 0.0;
951 :
952 : // extract sign: bit 7 of byte 5
953 145 : sign = r[5] & 0x80 ? -1 : 1;
954 :
955 : // extract mantissa from first bit of byte 1 to bit 7 of byte 5
956 145 : mant = 0;
957 725 : for (i = 1; i < 5; i++)
958 580 : mant = (r[i] + mant) / 256;
959 145 : mant = (mant + (r[5] & 0x7F)) / 128 + 1;
960 :
961 : // extract exponent
962 145 : exp = r[0] - 129;
963 :
964 : // compute the damned number
965 145 : return sign * ldexp(mant, exp);
966 : }
967 :
968 : /************************************************************************/
969 : /* c2tp() */
970 : /* */
971 : /* convert a double into a Turbo Pascal real */
972 : /************************************************************************/
973 :
974 85 : static void c2tp(double x, GByte *r)
975 : {
976 : double mant, temp;
977 : int negative, exp, i;
978 :
979 : // handle 0 case
980 85 : if (x == 0.0)
981 : {
982 25 : for (i = 0; i < 6; r[i++] = 0);
983 25 : return;
984 : }
985 :
986 : // compute mantissa, sign and exponent
987 60 : mant = frexp(x, &exp) * 2 - 1;
988 60 : exp--;
989 60 : negative = 0;
990 60 : if (mant < 0)
991 : {
992 9 : mant = -mant;
993 9 : negative = 1;
994 : }
995 : // stuff mantissa into Turbo Pascal real
996 60 : mant = modf(mant * 128, &temp);
997 60 : r[5] = (unsigned char) temp;
998 300 : for (i = 4; i >= 1; i--)
999 : {
1000 240 : mant = modf(mant * 256, &temp);
1001 240 : r[i] = (unsigned char) temp;
1002 : }
1003 : // add sign
1004 60 : if (negative)
1005 9 : r[5] |= 0x80;
1006 :
1007 : // put exponent
1008 60 : r[0] = (GByte) (exp + 129);
1009 : }
1010 :
1011 : /************************************************************************/
1012 : /* Create() */
1013 : /************************************************************************/
1014 :
1015 49 : GDALDataset *IDADataset::Create( const char * pszFilename,
1016 : int nXSize, int nYSize, int nBands,
1017 : GDALDataType eType,
1018 : char ** /* papszParmList */ )
1019 :
1020 : {
1021 : /* -------------------------------------------------------------------- */
1022 : /* Verify input options. */
1023 : /* -------------------------------------------------------------------- */
1024 49 : if( eType != GDT_Byte || nBands != 1 )
1025 : {
1026 : CPLError( CE_Failure, CPLE_AppDefined,
1027 39 : "Only 1 band, Byte datasets supported for IDA format." );
1028 :
1029 39 : return NULL;
1030 : }
1031 :
1032 : /* -------------------------------------------------------------------- */
1033 : /* Try to create the file. */
1034 : /* -------------------------------------------------------------------- */
1035 : FILE *fp;
1036 :
1037 10 : fp = VSIFOpen( pszFilename, "wb" );
1038 :
1039 10 : if( fp == NULL )
1040 : {
1041 : CPLError( CE_Failure, CPLE_OpenFailed,
1042 : "Attempt to create file `%s' failed.\n",
1043 3 : pszFilename );
1044 3 : return NULL;
1045 : }
1046 :
1047 : /* -------------------------------------------------------------------- */
1048 : /* Prepare formatted header. */
1049 : /* -------------------------------------------------------------------- */
1050 : GByte abyHeader[512];
1051 :
1052 7 : memset( abyHeader, 0, sizeof(abyHeader) );
1053 :
1054 7 : abyHeader[22] = 200; /* image type - CALCULATED */
1055 7 : abyHeader[23] = 0; /* projection - NONE */
1056 7 : abyHeader[30] = nYSize % 256;
1057 7 : abyHeader[31] = (GByte) (nYSize / 256);
1058 7 : abyHeader[32] = nXSize % 256;
1059 7 : abyHeader[33] = (GByte) (nXSize / 256);
1060 :
1061 7 : abyHeader[170] = 255; /* missing = 255 */
1062 7 : c2tp( 1.0, abyHeader + 171 ); /* slope = 1.0 */
1063 7 : c2tp( 0.0, abyHeader + 177 ); /* offset = 0 */
1064 7 : abyHeader[168] = 0; // lower limit
1065 7 : abyHeader[169] = 254; // upper limit
1066 :
1067 : // pixel size = 1.0
1068 7 : c2tp( 1.0, abyHeader + 144 );
1069 7 : c2tp( 1.0, abyHeader + 150 );
1070 :
1071 7 : if( VSIFWrite( abyHeader, 1, 512, fp ) != 512 )
1072 : {
1073 : CPLError( CE_Failure, CPLE_AppDefined,
1074 : "IO error writing %s.\n%s",
1075 0 : pszFilename, VSIStrerror( errno ) );
1076 0 : VSIFClose( fp );
1077 0 : return NULL;
1078 : }
1079 :
1080 : /* -------------------------------------------------------------------- */
1081 : /* Now we need to extend the file to just the right number of */
1082 : /* bytes for the data we have to ensure it will open again */
1083 : /* properly. */
1084 : /* -------------------------------------------------------------------- */
1085 7 : if( VSIFSeek( fp, nXSize * nYSize - 1, SEEK_CUR ) != 0
1086 : || VSIFWrite( abyHeader, 1, 1, fp ) != 1 )
1087 : {
1088 : CPLError( CE_Failure, CPLE_AppDefined,
1089 : "IO error writing %s.\n%s",
1090 0 : pszFilename, VSIStrerror( errno ) );
1091 0 : VSIFClose( fp );
1092 0 : return NULL;
1093 : }
1094 :
1095 7 : VSIFClose( fp );
1096 :
1097 7 : return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
1098 : }
1099 :
1100 : /************************************************************************/
1101 : /* GDALRegister_IDA() */
1102 : /************************************************************************/
1103 :
1104 582 : void GDALRegister_IDA()
1105 :
1106 : {
1107 : GDALDriver *poDriver;
1108 :
1109 582 : if( GDALGetDriverByName( "IDA" ) == NULL )
1110 : {
1111 561 : poDriver = new GDALDriver();
1112 :
1113 561 : poDriver->SetDescription( "IDA" );
1114 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
1115 561 : "Image Data and Analysis" );
1116 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
1117 561 : "frmt_various.html#IDA" );
1118 561 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES, "Byte" );
1119 :
1120 561 : poDriver->pfnOpen = IDADataset::Open;
1121 561 : poDriver->pfnCreate = IDADataset::Create;
1122 :
1123 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
1124 : }
1125 582 : }
1126 :
|
