1 : /******************************************************************************
2 : * $Id: idadataset.cpp 16706 2009-04-02 03:44:07Z warmerdam $
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 16706 2009-04-02 03:44:07Z warmerdam $");
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 24 : IDARasterBand::IDARasterBand( IDADataset *poDSIn,
137 : FILE *fpRaw, int nXSize )
138 : : RawRasterBand( poDSIn, 1, fpRaw, 512, 1, nXSize,
139 24 : GDT_Byte, FALSE, FALSE )
140 :
141 : {
142 24 : poColorTable = NULL;
143 24 : poRAT = NULL;
144 24 : }
145 :
146 : /************************************************************************/
147 : /* ~IDARasterBand() */
148 : /************************************************************************/
149 :
150 48 : IDARasterBand::~IDARasterBand()
151 :
152 : {
153 24 : delete poColorTable;
154 24 : delete poRAT;
155 48 : }
156 :
157 : /************************************************************************/
158 : /* GetNoDataValue() */
159 : /************************************************************************/
160 :
161 7 : double IDARasterBand::GetNoDataValue( int *pbSuccess )
162 :
163 : {
164 7 : if( pbSuccess != NULL )
165 6 : *pbSuccess = TRUE;
166 7 : return ((IDADataset *) poDS)->nMissing;
167 : }
168 :
169 : /************************************************************************/
170 : /* GetOffset() */
171 : /************************************************************************/
172 :
173 1 : double IDARasterBand::GetOffset( int *pbSuccess )
174 :
175 : {
176 1 : if( pbSuccess != NULL )
177 1 : *pbSuccess = TRUE;
178 1 : 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 1 : double IDARasterBand::GetScale( int *pbSuccess )
212 :
213 : {
214 1 : if( pbSuccess != NULL )
215 1 : *pbSuccess = TRUE;
216 1 : return ((IDADataset *) poDS)->dfM;
217 : }
218 :
219 : /************************************************************************/
220 : /* SetScale() */
221 : /************************************************************************/
222 :
223 1 : CPLErr IDARasterBand::SetScale( double dfNewValue )
224 :
225 : {
226 1 : IDADataset *poIDS = (IDADataset *) poDS;
227 :
228 1 : if( dfNewValue == poIDS->dfM )
229 0 : return CE_None;
230 :
231 1 : 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 1 : poIDS->dfM = dfNewValue;
239 1 : c2tp( dfNewValue, poIDS->abyHeader + 171 );
240 1 : poIDS->bHeaderDirty = TRUE;
241 :
242 1 : return CE_None;
243 : }
244 :
245 : /************************************************************************/
246 : /* GetColorTable() */
247 : /************************************************************************/
248 :
249 1 : GDALColorTable *IDARasterBand::GetColorTable()
250 :
251 : {
252 1 : if( poColorTable )
253 0 : return poColorTable;
254 : else
255 1 : return RawRasterBand::GetColorTable();
256 : }
257 :
258 : /************************************************************************/
259 : /* GetColorInterpretation() */
260 : /************************************************************************/
261 :
262 2 : GDALColorInterp IDARasterBand::GetColorInterpretation()
263 :
264 : {
265 2 : if( poColorTable )
266 0 : return GCI_PaletteIndex;
267 : else
268 2 : 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 24 : IDADataset::IDADataset()
295 : {
296 24 : fpRaw = NULL;
297 24 : pszProjection = NULL;
298 :
299 24 : adfGeoTransform[0] = 0.0;
300 24 : adfGeoTransform[1] = 1.0;
301 24 : adfGeoTransform[2] = 0.0;
302 24 : adfGeoTransform[3] = 0.0;
303 24 : adfGeoTransform[4] = 0.0;
304 24 : adfGeoTransform[5] = 1.0;
305 :
306 24 : bHeaderDirty = FALSE;
307 24 : }
308 :
309 : /************************************************************************/
310 : /* ~IDADataset() */
311 : /************************************************************************/
312 :
313 48 : IDADataset::~IDADataset()
314 :
315 : {
316 24 : FlushCache();
317 :
318 24 : if( fpRaw != NULL )
319 24 : VSIFClose( fpRaw );
320 24 : CPLFree( pszProjection );
321 48 : }
322 :
323 : /************************************************************************/
324 : /* ProcessGeoref() */
325 : /************************************************************************/
326 :
327 24 : void IDADataset::ProcessGeoref()
328 :
329 : {
330 24 : OGRSpatialReference oSRS;
331 :
332 24 : if( nProjection == 3 )
333 : {
334 9 : oSRS.SetWellKnownGeogCS( "WGS84" );
335 : }
336 15 : else if( nProjection == 4 )
337 : {
338 : oSRS.SetLCC( dfParallel1, dfParallel2,
339 : dfLatCenter, dfLongCenter,
340 2 : 0.0, 0.0 );
341 : oSRS.SetGeogCS( "Clarke 1866", "Clarke 1866", "Clarke 1866",
342 2 : 6378206.4, 293.97869821389662 );
343 : }
344 13 : else if( nProjection == 6 )
345 : {
346 2 : oSRS.SetLAEA( dfLatCenter, dfLongCenter, 0.0, 0.0 );
347 : oSRS.SetGeogCS( "Sphere", "Sphere", "Sphere",
348 2 : 6370997.0, 0.0 );
349 : }
350 11 : else if( nProjection == 8 )
351 : {
352 : oSRS.SetACEA( dfParallel1, dfParallel2,
353 : dfLatCenter, dfLongCenter,
354 2 : 0.0, 0.0 );
355 : oSRS.SetGeogCS( "Clarke 1866", "Clarke 1866", "Clarke 1866",
356 2 : 6378206.4, 293.97869821389662 );
357 : }
358 9 : else if( nProjection == 9 )
359 : {
360 2 : oSRS.SetGH( dfLongCenter, 0.0, 0.0 );
361 : oSRS.SetGeogCS( "Sphere", "Sphere", "Sphere",
362 2 : 6370997.0, 0.0 );
363 : }
364 :
365 24 : if( oSRS.GetRoot() != NULL )
366 : {
367 17 : CPLFree( pszProjection );
368 17 : pszProjection = NULL;
369 :
370 17 : oSRS.exportToWkt( &pszProjection );
371 : }
372 :
373 24 : adfGeoTransform[0] = 0 - dfDX * dfXCenter;
374 24 : adfGeoTransform[1] = dfDX;
375 24 : adfGeoTransform[2] = 0.0;
376 24 : adfGeoTransform[3] = dfDY * dfYCenter;
377 24 : adfGeoTransform[4] = 0.0;
378 24 : adfGeoTransform[5] = -dfDY;
379 :
380 24 : if( nProjection == 3 )
381 : {
382 9 : adfGeoTransform[0] += dfLongCenter;
383 9 : adfGeoTransform[3] += dfLatCenter;
384 24 : }
385 24 : }
386 :
387 : /************************************************************************/
388 : /* FlushCache() */
389 : /************************************************************************/
390 :
391 24 : void IDADataset::FlushCache()
392 :
393 : {
394 24 : RawDataset::FlushCache();
395 :
396 24 : if( bHeaderDirty )
397 : {
398 6 : VSIFSeek( fpRaw, 0, SEEK_SET );
399 6 : VSIFWrite( abyHeader, 512, 1, fpRaw );
400 6 : bHeaderDirty = FALSE;
401 : }
402 24 : }
403 :
404 : /************************************************************************/
405 : /* GetGeoTransform() */
406 : /************************************************************************/
407 :
408 3 : CPLErr IDADataset::GetGeoTransform( double *padfGeoTransform )
409 :
410 : {
411 3 : memcpy( padfGeoTransform, adfGeoTransform, sizeof(double) * 6 );
412 3 : return CE_None;
413 : }
414 :
415 : /************************************************************************/
416 : /* SetGeoTransform() */
417 : /************************************************************************/
418 :
419 6 : CPLErr IDADataset::SetGeoTransform( double *padfGeoTransform )
420 :
421 : {
422 6 : if( padfGeoTransform[2] != 0.0 || padfGeoTransform[4] != 0.0 )
423 0 : return GDALPamDataset::SetGeoTransform( padfGeoTransform );
424 :
425 6 : memcpy( adfGeoTransform, padfGeoTransform, sizeof(double) * 6 );
426 6 : bHeaderDirty = TRUE;
427 :
428 6 : dfDX = adfGeoTransform[1];
429 6 : dfDY = -adfGeoTransform[5];
430 6 : dfXCenter = -adfGeoTransform[0] / dfDX;
431 6 : dfYCenter = adfGeoTransform[3] / dfDY;
432 :
433 6 : c2tp( dfDX, abyHeader + 144 );
434 6 : c2tp( dfDY, abyHeader + 150 );
435 6 : c2tp( dfXCenter, abyHeader + 132 );
436 6 : c2tp( dfYCenter, abyHeader + 138 );
437 :
438 6 : return CE_None;
439 : }
440 :
441 : /************************************************************************/
442 : /* GetProjectionRef() */
443 : /************************************************************************/
444 :
445 10 : const char *IDADataset::GetProjectionRef()
446 :
447 : {
448 10 : if( pszProjection )
449 10 : return pszProjection;
450 : else
451 0 : return "";
452 : }
453 :
454 : /************************************************************************/
455 : /* SetProjection() */
456 : /************************************************************************/
457 :
458 6 : CPLErr IDADataset::SetProjection( const char *pszWKTIn )
459 :
460 : {
461 6 : OGRSpatialReference oSRS;
462 :
463 6 : oSRS.importFromWkt( (char **) &pszWKTIn );
464 :
465 6 : if( !oSRS.IsGeographic() && !oSRS.IsProjected() )
466 0 : GDALPamDataset::SetProjection( pszWKTIn );
467 :
468 : /* -------------------------------------------------------------------- */
469 : /* Clear projection parameters. */
470 : /* -------------------------------------------------------------------- */
471 6 : dfParallel1 = 0.0;
472 6 : dfParallel2 = 0.0;
473 6 : dfLatCenter = 0.0;
474 6 : dfLongCenter = 0.0;
475 :
476 : /* -------------------------------------------------------------------- */
477 : /* Geographic. */
478 : /* -------------------------------------------------------------------- */
479 6 : if( oSRS.IsGeographic() )
480 : {
481 : // If no change, just return.
482 2 : if( nProjection == 3 )
483 0 : return CE_None;
484 :
485 2 : 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 6 : 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 6 : const char *pszProjection = oSRS.GetAttrValue( "PROJECTION" );
506 :
507 6 : if( pszProjection == NULL )
508 : {
509 : /* do nothing - presumably geographic */;
510 : }
511 4 : else if( EQUAL(pszProjection,SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP) )
512 : {
513 1 : nProjection = 4;
514 1 : dfParallel1 = oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0);
515 1 : dfParallel2 = oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2,0.0);
516 1 : dfLatCenter = oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0);
517 1 : dfLongCenter = oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0);
518 : }
519 3 : else if( EQUAL(pszProjection,SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA) )
520 : {
521 1 : nProjection = 6;
522 1 : dfLatCenter = oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0);
523 1 : dfLongCenter = oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0);
524 : }
525 2 : else if( EQUAL(pszProjection,SRS_PT_ALBERS_CONIC_EQUAL_AREA) )
526 : {
527 1 : nProjection = 8;
528 1 : dfParallel1 = oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0);
529 1 : dfParallel2 = oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2,0.0);
530 1 : dfLatCenter = oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0);
531 1 : dfLongCenter = oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0);
532 : }
533 1 : else if( EQUAL(pszProjection,SRS_PT_GOODE_HOMOLOSINE) )
534 : {
535 1 : nProjection = 9;
536 1 : 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 6 : bHeaderDirty = TRUE;
547 :
548 6 : abyHeader[23] = nProjection;
549 6 : c2tp( dfLatCenter, abyHeader + 120 );
550 6 : c2tp( dfLongCenter, abyHeader + 126 );
551 6 : c2tp( dfParallel1, abyHeader + 156 );
552 6 : c2tp( dfParallel2, abyHeader + 162 );
553 :
554 6 : return CE_None;
555 : }
556 :
557 : /************************************************************************/
558 : /* ReadColorTable() */
559 : /************************************************************************/
560 :
561 24 : void IDADataset::ReadColorTable()
562 :
563 : {
564 : /* -------------------------------------------------------------------- */
565 : /* Decide what .clr file to look for and try to open. */
566 : /* -------------------------------------------------------------------- */
567 24 : CPLString osCLRFilename;
568 :
569 24 : osCLRFilename = CPLGetConfigOption( "IDA_COLOR_FILE", "" );
570 24 : if( strlen(osCLRFilename) == 0 )
571 24 : osCLRFilename = CPLResetExtension(GetDescription(), "clr" );
572 :
573 :
574 24 : FILE *fp = VSIFOpen( osCLRFilename, "r" );
575 24 : if( fp == NULL )
576 : {
577 24 : osCLRFilename = CPLResetExtension(osCLRFilename, "CLR" );
578 24 : fp = VSIFOpen( osCLRFilename, "r" );
579 : }
580 :
581 24 : 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 : /* -------------------------------------------------------------------- */
667 : /* Attach RAT to band. */
668 : /* -------------------------------------------------------------------- */
669 0 : ((IDARasterBand *) GetRasterBand( 1 ))->poRAT = poRAT;
670 :
671 : /* -------------------------------------------------------------------- */
672 : /* Build a conventional color table from this. */
673 : /* -------------------------------------------------------------------- */
674 : ((IDARasterBand *) GetRasterBand( 1 ))->poColorTable =
675 0 : poRAT->TranslateToColorTable();
676 : }
677 :
678 : /************************************************************************/
679 : /* Open() */
680 : /************************************************************************/
681 :
682 8520 : GDALDataset *IDADataset::Open( GDALOpenInfo * poOpenInfo )
683 :
684 : {
685 : /* -------------------------------------------------------------------- */
686 : /* Is this an IDA file? */
687 : /* -------------------------------------------------------------------- */
688 : int nXSize, nYSize;
689 : GIntBig nExpectedFileSize, nActualFileSize;
690 :
691 8520 : if( poOpenInfo->fp == NULL )
692 8276 : return NULL;
693 :
694 244 : if( poOpenInfo->nHeaderBytes < 512 )
695 83 : return NULL;
696 :
697 : // projection legal?
698 161 : if( poOpenInfo->pabyHeader[23] > 10 )
699 98 : return NULL;
700 :
701 : // imagetype legal?
702 190 : if( (poOpenInfo->pabyHeader[22] > 14
703 32 : && poOpenInfo->pabyHeader[22] < 100)
704 63 : || (poOpenInfo->pabyHeader[22] > 114
705 32 : && poOpenInfo->pabyHeader[22] != 200 ) )
706 8 : return NULL;
707 :
708 55 : nXSize = poOpenInfo->pabyHeader[30] + poOpenInfo->pabyHeader[31] * 256;
709 55 : nYSize = poOpenInfo->pabyHeader[32] + poOpenInfo->pabyHeader[33] * 256;
710 :
711 55 : if( nXSize == 0 || nYSize == 0 )
712 27 : return NULL;
713 :
714 : // The file just be exactly the image size + header size in length.
715 28 : nExpectedFileSize = nXSize * nYSize + 512;
716 :
717 28 : VSIFSeek( poOpenInfo->fp, 0, SEEK_END );
718 28 : nActualFileSize = VSIFTell( poOpenInfo->fp );
719 28 : VSIRewind( poOpenInfo->fp );
720 :
721 28 : if( nActualFileSize != nExpectedFileSize )
722 4 : return NULL;
723 :
724 : /* -------------------------------------------------------------------- */
725 : /* Create the dataset. */
726 : /* -------------------------------------------------------------------- */
727 24 : IDADataset *poDS = new IDADataset();
728 :
729 24 : memcpy( poDS->abyHeader, poOpenInfo->pabyHeader, 512 );
730 :
731 : /* -------------------------------------------------------------------- */
732 : /* Parse various values out of the header. */
733 : /* -------------------------------------------------------------------- */
734 24 : poDS->nImageType = poOpenInfo->pabyHeader[22];
735 24 : poDS->nProjection = poOpenInfo->pabyHeader[23];
736 :
737 24 : poDS->nRasterYSize = poOpenInfo->pabyHeader[30]
738 24 : + poOpenInfo->pabyHeader[31] * 256;
739 24 : poDS->nRasterXSize = poOpenInfo->pabyHeader[32]
740 24 : + poOpenInfo->pabyHeader[33] * 256;
741 :
742 24 : strncpy( poDS->szTitle, (const char *) poOpenInfo->pabyHeader+38, 80 );
743 24 : poDS->szTitle[80] = '\0';
744 :
745 24 : int nLastTitleChar = strlen(poDS->szTitle)-1;
746 1525 : while( nLastTitleChar > -1
747 371 : && (poDS->szTitle[nLastTitleChar] == 10
748 371 : || poDS->szTitle[nLastTitleChar] == 13
749 371 : || poDS->szTitle[nLastTitleChar] == ' ') )
750 364 : poDS->szTitle[nLastTitleChar--] = '\0';
751 :
752 24 : poDS->dfLatCenter = tp2c( poOpenInfo->pabyHeader + 120 );
753 24 : poDS->dfLongCenter = tp2c( poOpenInfo->pabyHeader + 126 );
754 24 : poDS->dfXCenter = tp2c( poOpenInfo->pabyHeader + 132 );
755 24 : poDS->dfYCenter = tp2c( poOpenInfo->pabyHeader + 138 );
756 24 : poDS->dfDX = tp2c( poOpenInfo->pabyHeader + 144 );
757 24 : poDS->dfDY = tp2c( poOpenInfo->pabyHeader + 150 );
758 24 : poDS->dfParallel1 = tp2c( poOpenInfo->pabyHeader + 156 );
759 24 : poDS->dfParallel2 = tp2c( poOpenInfo->pabyHeader + 162 );
760 :
761 24 : poDS->ProcessGeoref();
762 :
763 24 : poDS->SetMetadataItem( "TITLE", poDS->szTitle );
764 :
765 : /* -------------------------------------------------------------------- */
766 : /* Handle various image types. */
767 : /* -------------------------------------------------------------------- */
768 :
769 : /*
770 : GENERIC = 0
771 : FEW S NDVI = 1
772 : EROS NDVI = 6
773 : ARTEMIS CUTOFF = 10
774 : ARTEMIS RECODE = 11
775 : ARTEMIS NDVI = 12
776 : ARTEMIS FEWS = 13
777 : ARTEMIS NEWNASA = 14
778 : GENERIC DIFF = 100
779 : FEW S NDVI DIFF = 101
780 : EROS NDVI DIFF = 106
781 : ARTEMIS CUTOFF DIFF = 110
782 : ARTEMIS RECODE DIFF = 111
783 : ARTEMIS NDVI DIFF = 112
784 : ARTEMIS FEWS DIFF = 113
785 : ARTEMIS NEWNASA DIFF = 114
786 : CALCULATED =200
787 : */
788 :
789 24 : poDS->nMissing = 0;
790 :
791 24 : switch( poDS->nImageType )
792 : {
793 : case 1:
794 0 : poDS->SetMetadataItem( "IMAGETYPE", "1, FEWS NDVI" );
795 0 : poDS->dfM = 1/256.0;
796 0 : poDS->dfB = -82/256.0;
797 0 : break;
798 :
799 : case 6:
800 0 : poDS->SetMetadataItem( "IMAGETYPE", "6, EROS NDVI" );
801 0 : poDS->dfM = 1/100.0;
802 0 : poDS->dfB = -100/100.0;
803 0 : break;
804 :
805 : case 10:
806 0 : poDS->SetMetadataItem( "IMAGETYPE", "10, ARTEMIS CUTOFF" );
807 0 : poDS->dfM = 1.0;
808 0 : poDS->dfB = 0.0;
809 0 : poDS->nMissing = 254;
810 0 : break;
811 :
812 : case 11:
813 0 : poDS->SetMetadataItem( "IMAGETYPE", "11, ARTEMIS RECODE" );
814 0 : poDS->dfM = 4.0;
815 0 : poDS->dfB = 0.0;
816 0 : poDS->nMissing = 254;
817 0 : break;
818 :
819 : case 12: /* ANDVI */
820 0 : poDS->SetMetadataItem( "IMAGETYPE", "12, ARTEMIS NDVI" );
821 0 : poDS->dfM = 4/500.0;
822 0 : poDS->dfB = -3/500.0 - 1.0;
823 0 : poDS->nMissing = 254;
824 0 : break;
825 :
826 : case 13: /* AFEWS */
827 0 : poDS->SetMetadataItem( "IMAGETYPE", "13, ARTEMIS FEWS" );
828 0 : poDS->dfM = 1/256.0;
829 0 : poDS->dfB = -82/256.0;
830 0 : poDS->nMissing = 254;
831 0 : break;
832 :
833 : case 14: /* NEWNASA */
834 0 : poDS->SetMetadataItem( "IMAGETYPE", "13, ARTEMIS NEWNASA" );
835 0 : poDS->dfM = 0.75/250.0;
836 0 : poDS->dfB = 0.0;
837 0 : poDS->nMissing = 254;
838 0 : break;
839 :
840 : case 101: /* NDVI_DIFF (FEW S) */
841 0 : poDS->dfM = 1/128.0;
842 0 : poDS->dfB = -1.0;
843 0 : poDS->nMissing = 0;
844 0 : break;
845 :
846 : case 106: /* EROS_DIFF */
847 0 : poDS->dfM = 1/50.0;
848 0 : poDS->dfB = -128/50.0;
849 0 : poDS->nMissing = 0;
850 0 : break;
851 :
852 : case 110: /* CUTOFF_DIFF */
853 0 : poDS->dfM = 2.0;
854 0 : poDS->dfB = -128*2;
855 0 : poDS->nMissing = 254;
856 0 : break;
857 :
858 : case 111: /* RECODE_DIFF */
859 0 : poDS->dfM = 8;
860 0 : poDS->dfB = -128*8;
861 0 : poDS->nMissing = 254;
862 0 : break;
863 :
864 : case 112: /* ANDVI_DIFF */
865 0 : poDS->dfM = 8/1000.0;
866 0 : poDS->dfB = (-128*8)/1000.0;
867 0 : poDS->nMissing = 254;
868 0 : break;
869 :
870 : case 113: /* AFEWS_DIFF */
871 0 : poDS->dfM = 1/128.0;
872 0 : poDS->dfB = -1;
873 0 : poDS->nMissing = 254;
874 0 : break;
875 :
876 : case 114: /* NEWNASA_DIFF */
877 0 : poDS->dfM = 0.75/125.0;
878 0 : poDS->dfB = -128*poDS->dfM;
879 0 : poDS->nMissing = 254;
880 0 : break;
881 :
882 : case 200:
883 : /* we use the values from the header */
884 24 : poDS->dfM = tp2c( poOpenInfo->pabyHeader + 171 );
885 24 : poDS->dfB = tp2c( poOpenInfo->pabyHeader + 177 );
886 24 : poDS->nMissing = poOpenInfo->pabyHeader[170];
887 24 : break;
888 :
889 : default:
890 0 : poDS->dfM = 1.0;
891 0 : poDS->dfB = 0.0;
892 : break;
893 : }
894 :
895 : /* -------------------------------------------------------------------- */
896 : /* Create the band. */
897 : /* -------------------------------------------------------------------- */
898 24 : if( poOpenInfo->eAccess == GA_ReadOnly )
899 : {
900 17 : poDS->fpRaw = poOpenInfo->fp;
901 17 : poOpenInfo->fp = NULL;
902 : }
903 : else
904 : {
905 7 : poDS->fpRaw = VSIFOpen( poOpenInfo->pszFilename, "rb+" );
906 7 : poDS->eAccess = GA_Update;
907 7 : if( poDS->fpRaw == NULL )
908 : {
909 : CPLError( CE_Failure, CPLE_OpenFailed,
910 : "Failed to open %s for write access.",
911 0 : poOpenInfo->pszFilename );
912 0 : return NULL;
913 : }
914 : }
915 :
916 : poDS->SetBand( 1, new IDARasterBand( poDS, poDS->fpRaw,
917 24 : poDS->nRasterXSize ) );
918 :
919 : /* -------------------------------------------------------------------- */
920 : /* Check for a color table. */
921 : /* -------------------------------------------------------------------- */
922 24 : poDS->SetDescription( poOpenInfo->pszFilename );
923 24 : poDS->ReadColorTable();
924 :
925 : /* -------------------------------------------------------------------- */
926 : /* Initialize any PAM information. */
927 : /* -------------------------------------------------------------------- */
928 24 : poDS->TryLoadXML();
929 :
930 : /* -------------------------------------------------------------------- */
931 : /* Check for overviews. */
932 : /* -------------------------------------------------------------------- */
933 24 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
934 :
935 24 : return( poDS );
936 : }
937 :
938 : /************************************************************************/
939 : /* tp2c() */
940 : /* */
941 : /* convert a Turbo Pascal real into a double */
942 : /************************************************************************/
943 :
944 240 : static double tp2c(GByte *r)
945 : {
946 : double mant;
947 : int sign, exp, i;
948 :
949 : // handle 0 case
950 240 : if (r[0] == 0)
951 100 : return 0.0;
952 :
953 : // extract sign: bit 7 of byte 5
954 140 : sign = r[5] & 0x80 ? -1 : 1;
955 :
956 : // extract mantissa from first bit of byte 1 to bit 7 of byte 5
957 140 : mant = 0;
958 700 : for (i = 1; i < 5; i++)
959 560 : mant = (r[i] + mant) / 256;
960 140 : mant = (mant + (r[5] & 0x7F)) / 128 + 1;
961 :
962 : // extract exponent
963 140 : exp = r[0] - 129;
964 :
965 : // compute the damned number
966 140 : return sign * ldexp(mant, exp);
967 : }
968 :
969 : /************************************************************************/
970 : /* c2tp() */
971 : /* */
972 : /* convert a double into a Turbo Pascal real */
973 : /************************************************************************/
974 :
975 77 : static void c2tp(double x, GByte *r)
976 : {
977 : double mant, temp;
978 : int negative, exp, i;
979 :
980 : // handle 0 case
981 77 : if (x == 0.0)
982 : {
983 21 : for (i = 0; i < 6; r[i++] = 0);
984 21 : return;
985 : }
986 :
987 : // compute mantissa, sign and exponent
988 56 : mant = frexp(x, &exp) * 2 - 1;
989 56 : exp--;
990 56 : negative = 0;
991 56 : if (mant < 0)
992 : {
993 8 : mant = -mant;
994 8 : negative = 1;
995 : }
996 : // stuff mantissa into Turbo Pascal real
997 56 : mant = modf(mant * 128, &temp);
998 56 : r[5] = (unsigned char) temp;
999 280 : for (i = 4; i >= 1; i--)
1000 : {
1001 224 : mant = modf(mant * 256, &temp);
1002 224 : r[i] = (unsigned char) temp;
1003 : }
1004 : // add sign
1005 56 : if (negative)
1006 8 : r[5] |= 0x80;
1007 :
1008 : // put exponent
1009 56 : r[0] = exp + 129;
1010 : }
1011 :
1012 : /************************************************************************/
1013 : /* Create() */
1014 : /************************************************************************/
1015 :
1016 36 : GDALDataset *IDADataset::Create( const char * pszFilename,
1017 : int nXSize, int nYSize, int nBands,
1018 : GDALDataType eType,
1019 : char ** /* papszParmList */ )
1020 :
1021 : {
1022 : /* -------------------------------------------------------------------- */
1023 : /* Verify input options. */
1024 : /* -------------------------------------------------------------------- */
1025 36 : if( eType != GDT_Byte || nBands != 1 )
1026 : {
1027 : CPLError( CE_Failure, CPLE_AppDefined,
1028 29 : "Only 1 band, Byte datasets supported for IDA format." );
1029 :
1030 29 : return NULL;
1031 : }
1032 :
1033 : /* -------------------------------------------------------------------- */
1034 : /* Try to create the file. */
1035 : /* -------------------------------------------------------------------- */
1036 : FILE *fp;
1037 :
1038 7 : fp = VSIFOpen( pszFilename, "wb" );
1039 :
1040 7 : if( fp == NULL )
1041 : {
1042 : CPLError( CE_Failure, CPLE_OpenFailed,
1043 : "Attempt to create file `%s' failed.\n",
1044 0 : pszFilename );
1045 0 : return NULL;
1046 : }
1047 :
1048 : /* -------------------------------------------------------------------- */
1049 : /* Prepare formatted header. */
1050 : /* -------------------------------------------------------------------- */
1051 : GByte abyHeader[512];
1052 :
1053 7 : memset( abyHeader, 0, sizeof(abyHeader) );
1054 :
1055 7 : abyHeader[22] = 200; /* image type - CALCULATED */
1056 7 : abyHeader[23] = 0; /* projection - NONE */
1057 7 : abyHeader[30] = nYSize % 256;
1058 7 : abyHeader[31] = nYSize / 256;
1059 7 : abyHeader[32] = nXSize % 256;
1060 7 : abyHeader[33] = nXSize / 256;
1061 :
1062 7 : abyHeader[170] = 255; /* missing = 255 */
1063 7 : c2tp( 1.0, abyHeader + 171 ); /* slope = 1.0 */
1064 7 : c2tp( 0.0, abyHeader + 177 ); /* offset = 0 */
1065 7 : abyHeader[168] = 0; // lower limit
1066 7 : abyHeader[169] = 254; // upper limit
1067 :
1068 : // pixel size = 1.0
1069 7 : c2tp( 1.0, abyHeader + 144 );
1070 7 : c2tp( 1.0, abyHeader + 150 );
1071 :
1072 7 : if( VSIFWrite( abyHeader, 1, 512, fp ) != 512 )
1073 : {
1074 : CPLError( CE_Failure, CPLE_AppDefined,
1075 : "IO error writing %s.\n%s",
1076 0 : pszFilename, VSIStrerror( errno ) );
1077 0 : VSIFClose( fp );
1078 0 : return NULL;
1079 : }
1080 :
1081 : /* -------------------------------------------------------------------- */
1082 : /* Now we need to extend the file to just the right number of */
1083 : /* bytes for the data we have to ensure it will open again */
1084 : /* properly. */
1085 : /* -------------------------------------------------------------------- */
1086 7 : if( VSIFSeek( fp, nXSize * nYSize - 1, SEEK_CUR ) != 0
1087 : || VSIFWrite( abyHeader, 1, 1, fp ) != 1 )
1088 : {
1089 : CPLError( CE_Failure, CPLE_AppDefined,
1090 : "IO error writing %s.\n%s",
1091 0 : pszFilename, VSIStrerror( errno ) );
1092 0 : VSIFClose( fp );
1093 0 : return NULL;
1094 : }
1095 :
1096 7 : VSIFClose( fp );
1097 :
1098 7 : return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
1099 : }
1100 :
1101 : /************************************************************************/
1102 : /* GDALRegister_IDA() */
1103 : /************************************************************************/
1104 :
1105 338 : void GDALRegister_IDA()
1106 :
1107 : {
1108 : GDALDriver *poDriver;
1109 :
1110 338 : if( GDALGetDriverByName( "IDA" ) == NULL )
1111 : {
1112 336 : poDriver = new GDALDriver();
1113 :
1114 336 : poDriver->SetDescription( "IDA" );
1115 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
1116 336 : "Image Data and Analysis" );
1117 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
1118 336 : "frmt_various.html#IDA" );
1119 336 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES, "Byte" );
1120 :
1121 336 : poDriver->pfnOpen = IDADataset::Open;
1122 336 : poDriver->pfnCreate = IDADataset::Create;
1123 :
1124 336 : GetGDALDriverManager()->RegisterDriver( poDriver );
1125 : }
1126 338 : }
1127 :
|
