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