1 : /******************************************************************************
2 : * $Id: envidataset.cpp 23913 2012-02-06 09:04:41Z rouault $
3 : *
4 : * Project: ENVI .hdr Driver
5 : * Purpose: Implementation of ENVI .hdr labelled raw raster support.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : * Maintainer: Chris Padwick (cpadwick at ittvis.com)
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 2002, Frank Warmerdam
11 : *
12 : * Permission is hereby granted, free of charge, to any person obtaining a
13 : * copy of this software and associated documentation files (the "Software"),
14 : * to deal in the Software without restriction, including without limitation
15 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : * and/or sell copies of the Software, and to permit persons to whom the
17 : * Software is furnished to do so, subject to the following conditions:
18 : *
19 : * The above copyright notice and this permission notice shall be included
20 : * in all copies or substantial portions of the Software.
21 : *
22 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 : * DEALINGS IN THE SOFTWARE.
29 : ****************************************************************************/
30 :
31 : #include "rawdataset.h"
32 : #include "ogr_spatialref.h"
33 : #include "cpl_string.h"
34 :
35 : CPL_CVSID("$Id: envidataset.cpp 23913 2012-02-06 09:04:41Z rouault $");
36 :
37 : CPL_C_START
38 : void GDALRegister_ENVI(void);
39 : CPL_C_END
40 :
41 : static const int anUsgsEsriZones[] =
42 : {
43 : 101, 3101,
44 : 102, 3126,
45 : 201, 3151,
46 : 202, 3176,
47 : 203, 3201,
48 : 301, 3226,
49 : 302, 3251,
50 : 401, 3276,
51 : 402, 3301,
52 : 403, 3326,
53 : 404, 3351,
54 : 405, 3376,
55 : 406, 3401,
56 : 407, 3426,
57 : 501, 3451,
58 : 502, 3476,
59 : 503, 3501,
60 : 600, 3526,
61 : 700, 3551,
62 : 901, 3601,
63 : 902, 3626,
64 : 903, 3576,
65 : 1001, 3651,
66 : 1002, 3676,
67 : 1101, 3701,
68 : 1102, 3726,
69 : 1103, 3751,
70 : 1201, 3776,
71 : 1202, 3801,
72 : 1301, 3826,
73 : 1302, 3851,
74 : 1401, 3876,
75 : 1402, 3901,
76 : 1501, 3926,
77 : 1502, 3951,
78 : 1601, 3976,
79 : 1602, 4001,
80 : 1701, 4026,
81 : 1702, 4051,
82 : 1703, 6426,
83 : 1801, 4076,
84 : 1802, 4101,
85 : 1900, 4126,
86 : 2001, 4151,
87 : 2002, 4176,
88 : 2101, 4201,
89 : 2102, 4226,
90 : 2103, 4251,
91 : 2111, 6351,
92 : 2112, 6376,
93 : 2113, 6401,
94 : 2201, 4276,
95 : 2202, 4301,
96 : 2203, 4326,
97 : 2301, 4351,
98 : 2302, 4376,
99 : 2401, 4401,
100 : 2402, 4426,
101 : 2403, 4451,
102 : 2500, 0,
103 : 2501, 4476,
104 : 2502, 4501,
105 : 2503, 4526,
106 : 2600, 0,
107 : 2601, 4551,
108 : 2602, 4576,
109 : 2701, 4601,
110 : 2702, 4626,
111 : 2703, 4651,
112 : 2800, 4676,
113 : 2900, 4701,
114 : 3001, 4726,
115 : 3002, 4751,
116 : 3003, 4776,
117 : 3101, 4801,
118 : 3102, 4826,
119 : 3103, 4851,
120 : 3104, 4876,
121 : 3200, 4901,
122 : 3301, 4926,
123 : 3302, 4951,
124 : 3401, 4976,
125 : 3402, 5001,
126 : 3501, 5026,
127 : 3502, 5051,
128 : 3601, 5076,
129 : 3602, 5101,
130 : 3701, 5126,
131 : 3702, 5151,
132 : 3800, 5176,
133 : 3900, 0,
134 : 3901, 5201,
135 : 3902, 5226,
136 : 4001, 5251,
137 : 4002, 5276,
138 : 4100, 5301,
139 : 4201, 5326,
140 : 4202, 5351,
141 : 4203, 5376,
142 : 4204, 5401,
143 : 4205, 5426,
144 : 4301, 5451,
145 : 4302, 5476,
146 : 4303, 5501,
147 : 4400, 5526,
148 : 4501, 5551,
149 : 4502, 5576,
150 : 4601, 5601,
151 : 4602, 5626,
152 : 4701, 5651,
153 : 4702, 5676,
154 : 4801, 5701,
155 : 4802, 5726,
156 : 4803, 5751,
157 : 4901, 5776,
158 : 4902, 5801,
159 : 4903, 5826,
160 : 4904, 5851,
161 : 5001, 6101,
162 : 5002, 6126,
163 : 5003, 6151,
164 : 5004, 6176,
165 : 5005, 6201,
166 : 5006, 6226,
167 : 5007, 6251,
168 : 5008, 6276,
169 : 5009, 6301,
170 : 5010, 6326,
171 : 5101, 5876,
172 : 5102, 5901,
173 : 5103, 5926,
174 : 5104, 5951,
175 : 5105, 5976,
176 : 5201, 6001,
177 : 5200, 6026,
178 : 5200, 6076,
179 : 5201, 6051,
180 : 5202, 6051,
181 : 5300, 0,
182 : 5400, 0
183 : };
184 :
185 : /************************************************************************/
186 : /* ITTVISToUSGSZone() */
187 : /* */
188 : /* Convert ITTVIS style state plane zones to NOS style state */
189 : /* plane zones. The ENVI default is to use the new NOS zones, */
190 : /* but the old state plane zones can be used. Handle this. */
191 : /************************************************************************/
192 :
193 0 : static int ITTVISToUSGSZone( int nITTVISZone )
194 :
195 : {
196 0 : int nPairs = sizeof(anUsgsEsriZones) / (2*sizeof(int));
197 : int i;
198 :
199 : // Default is to use the zone as-is, as long as it is in the
200 : // available list
201 0 : for( i = 0; i < nPairs; i++ )
202 : {
203 0 : if( anUsgsEsriZones[i*2] == nITTVISZone )
204 0 : return anUsgsEsriZones[i*2];
205 : }
206 :
207 : // If not found in the new style, see if it is present in the
208 : // old style list and convert it. We don't expect to see this
209 : // often, but older files allowed it and may still exist.
210 0 : for( i = 0; i < nPairs; i++ )
211 : {
212 0 : if( anUsgsEsriZones[i*2+1] == nITTVISZone )
213 0 : return anUsgsEsriZones[i*2];
214 : }
215 :
216 0 : return nITTVISZone; // perhaps it *is* the USGS zone?
217 : }
218 :
219 : /************************************************************************/
220 : /* ==================================================================== */
221 : /* ENVIDataset */
222 : /* ==================================================================== */
223 : /************************************************************************/
224 :
225 : class ENVIDataset : public RawDataset
226 : {
227 : VSILFILE *fpImage; // image data file.
228 : VSILFILE *fp; // header file
229 : char *pszHDRFilename;
230 :
231 : int bFoundMapinfo;
232 :
233 : int bHeaderDirty;
234 :
235 : double adfGeoTransform[6];
236 :
237 : char *pszProjection;
238 :
239 : char **papszHeader;
240 :
241 : CPLString osStaFilename;
242 :
243 : int ReadHeader( VSILFILE * );
244 : int ProcessMapinfo( const char * );
245 : void ProcessRPCinfo( const char * ,int ,int);
246 : void ProcessStatsFile();
247 : int byteSwapInt(int);
248 : float byteSwapFloat(float);
249 : double byteSwapDouble(double);
250 : void SetENVIDatum( OGRSpatialReference *, const char * );
251 : void SetENVIEllipse( OGRSpatialReference *, char ** );
252 : void WriteProjectionInfo();
253 : int ParseRpcCoeffsMetaDataString(const char *psName, char *papszVal[], int& idx);
254 : int WriteRpcInfo();
255 : int WritePseudoGcpInfo();
256 :
257 : char **SplitList( const char * );
258 :
259 : enum Interleave { BSQ, BIL, BIP } interleave;
260 : static int GetEnviType(GDALDataType eType);
261 :
262 : public:
263 : ENVIDataset();
264 : ~ENVIDataset();
265 :
266 : virtual void FlushCache( void );
267 : virtual CPLErr GetGeoTransform( double * padfTransform );
268 : virtual CPLErr SetGeoTransform( double * );
269 : virtual const char *GetProjectionRef(void);
270 : virtual CPLErr SetProjection( const char * );
271 : virtual char **GetFileList(void);
272 :
273 : virtual CPLErr SetMetadata( char ** papszMetadata,
274 : const char * pszDomain = "" );
275 : virtual CPLErr SetMetadataItem( const char * pszName,
276 : const char * pszValue,
277 : const char * pszDomain = "" );
278 : virtual CPLErr SetGCPs( int nGCPCount, const GDAL_GCP *pasGCPList,
279 : const char *pszGCPProjection );
280 :
281 : static GDALDataset *Open( GDALOpenInfo * );
282 : static GDALDataset *Create( const char * pszFilename,
283 : int nXSize, int nYSize, int nBands,
284 : GDALDataType eType, char ** papszOptions );
285 : };
286 :
287 : /************************************************************************/
288 : /* ENVIDataset() */
289 : /************************************************************************/
290 :
291 226 : ENVIDataset::ENVIDataset()
292 : {
293 226 : fpImage = NULL;
294 226 : fp = NULL;
295 226 : pszHDRFilename = NULL;
296 226 : pszProjection = CPLStrdup("");
297 :
298 226 : papszHeader = NULL;
299 :
300 226 : bFoundMapinfo = FALSE;
301 :
302 226 : bHeaderDirty = FALSE;
303 :
304 226 : adfGeoTransform[0] = 0.0;
305 226 : adfGeoTransform[1] = 1.0;
306 226 : adfGeoTransform[2] = 0.0;
307 226 : adfGeoTransform[3] = 0.0;
308 226 : adfGeoTransform[4] = 0.0;
309 226 : adfGeoTransform[5] = 1.0;
310 226 : }
311 :
312 : /************************************************************************/
313 : /* ~ENVIDataset() */
314 : /************************************************************************/
315 :
316 226 : ENVIDataset::~ENVIDataset()
317 :
318 : {
319 226 : FlushCache();
320 226 : if( fpImage )
321 222 : VSIFCloseL( fpImage );
322 226 : if( fp )
323 226 : VSIFCloseL( fp );
324 226 : CPLFree( pszProjection );
325 226 : CSLDestroy( papszHeader );
326 226 : CPLFree(pszHDRFilename);
327 226 : }
328 :
329 : /************************************************************************/
330 : /* FlushCache() */
331 : /************************************************************************/
332 :
333 226 : void ENVIDataset::FlushCache()
334 :
335 : {
336 226 : RawDataset::FlushCache();
337 :
338 226 : GDALRasterBand* band = (GetRasterCount() > 0) ? GetRasterBand(1) : NULL;
339 :
340 226 : if ( band == NULL || !(bHeaderDirty || band->GetCategoryNames() != NULL) )
341 116 : return;
342 :
343 110 : CPLLocaleC oLocaleEnforcer;
344 :
345 110 : VSIFSeekL( fp, 0, SEEK_SET );
346 : /* -------------------------------------------------------------------- */
347 : /* Rewrite out the header. */
348 : /* -------------------------------------------------------------------- */
349 : int iBigEndian;
350 :
351 : const char *pszInterleaving;
352 : char** catNames;
353 :
354 : #ifdef CPL_LSB
355 110 : iBigEndian = 0;
356 : #else
357 : iBigEndian = 1;
358 : #endif
359 :
360 110 : VSIFPrintfL( fp, "ENVI\n" );
361 110 : if ("" != sDescription)
362 110 : VSIFPrintfL( fp, "description = {\n%s}\n", sDescription.c_str());
363 : VSIFPrintfL( fp, "samples = %d\nlines = %d\nbands = %d\n",
364 110 : nRasterXSize, nRasterYSize, nBands );
365 :
366 110 : catNames = band->GetCategoryNames();
367 :
368 110 : VSIFPrintfL( fp, "header offset = 0\n");
369 110 : if (0 == catNames)
370 102 : VSIFPrintfL( fp, "file type = ENVI Standard\n" );
371 : else
372 8 : VSIFPrintfL( fp, "file type = ENVI Classification\n" );
373 :
374 110 : int iENVIType = GetEnviType(band->GetRasterDataType());
375 110 : VSIFPrintfL( fp, "data type = %d\n", iENVIType );
376 110 : switch (interleave)
377 : {
378 : case BIP:
379 0 : pszInterleaving = "bip"; // interleaved by pixel
380 0 : break;
381 : case BIL:
382 0 : pszInterleaving = "bil"; // interleaved by line
383 0 : break;
384 : case BSQ:
385 110 : pszInterleaving = "bsq"; // band sequental by default
386 110 : break;
387 : default:
388 0 : pszInterleaving = "bsq";
389 : break;
390 : }
391 110 : VSIFPrintfL( fp, "interleave = %s\n", pszInterleaving);
392 110 : VSIFPrintfL( fp, "byte order = %d\n", iBigEndian );
393 :
394 : /* -------------------------------------------------------------------- */
395 : /* Write class and color information */
396 : /* -------------------------------------------------------------------- */
397 110 : catNames = band->GetCategoryNames();
398 110 : if (0 != catNames)
399 : {
400 8 : int nrClasses = 0;
401 32 : while (*catNames++)
402 16 : ++nrClasses;
403 :
404 8 : if (nrClasses > 0)
405 : {
406 8 : VSIFPrintfL( fp, "classes = %d\n", nrClasses );
407 :
408 8 : GDALColorTable* colorTable = band->GetColorTable();
409 8 : if (0 != colorTable)
410 : {
411 8 : int nrColors = colorTable->GetColorEntryCount();
412 8 : if (nrColors > nrClasses)
413 0 : nrColors = nrClasses;
414 8 : VSIFPrintfL( fp, "class lookup = {\n");
415 24 : for (int i = 0; i < nrColors; ++i)
416 : {
417 16 : const GDALColorEntry* color = colorTable->GetColorEntry(i);
418 16 : VSIFPrintfL(fp, "%d, %d, %d", color->c1, color->c2, color->c3);
419 16 : if (i < nrColors - 1)
420 : {
421 8 : VSIFPrintfL(fp, ", ");
422 8 : if (0 == (i+1) % 5)
423 0 : VSIFPrintfL(fp, "\n");
424 : }
425 : }
426 8 : VSIFPrintfL(fp, "}\n");
427 : }
428 :
429 8 : catNames = band->GetCategoryNames();
430 8 : if (0 != *catNames)
431 : {
432 8 : VSIFPrintfL( fp, "class names = {\n%s", *catNames++);
433 8 : int i = 0;
434 24 : while (*catNames) {
435 8 : VSIFPrintfL( fp, ",");
436 8 : if (0 == (++i) % 5)
437 0 : VSIFPrintfL(fp, "\n");
438 8 : VSIFPrintfL( fp, " %s", *catNames++);
439 : }
440 8 : VSIFPrintfL( fp, "}\n");
441 : }
442 : }
443 : }
444 :
445 : /* -------------------------------------------------------------------- */
446 : /* Write the rest of header. */
447 : /* -------------------------------------------------------------------- */
448 :
449 : // only one map info type should be set
450 : // - rpc
451 : // - pseudo/gcp
452 : // - standard
453 110 : if ( !WriteRpcInfo() ) // are rpcs in the metadata
454 : {
455 102 : if ( !WritePseudoGcpInfo() ) // are gcps in the metadata
456 : {
457 102 : WriteProjectionInfo(); // standard - affine xform/coord sys str
458 : }
459 : }
460 :
461 :
462 110 : VSIFPrintfL( fp, "band names = {\n" );
463 300 : for ( int i = 1; i <= nBands; i++ )
464 : {
465 190 : CPLString sBandDesc = GetRasterBand( i )->GetDescription();
466 :
467 190 : if ( sBandDesc == "" )
468 168 : sBandDesc = CPLSPrintf( "Band %d", i );
469 190 : VSIFPrintfL( fp, "%s", sBandDesc.c_str() );
470 190 : if ( i != nBands )
471 80 : VSIFPrintfL( fp, ",\n" );
472 : }
473 110 : VSIFPrintfL( fp, "}\n" );
474 : }
475 :
476 : /************************************************************************/
477 : /* GetFileList() */
478 : /************************************************************************/
479 :
480 20 : char **ENVIDataset::GetFileList()
481 :
482 : {
483 20 : char **papszFileList = NULL;
484 :
485 : // Main data file, etc.
486 20 : papszFileList = GDALPamDataset::GetFileList();
487 :
488 : // Header file.
489 20 : papszFileList = CSLAddString( papszFileList, pszHDRFilename );
490 :
491 : // Statistics file
492 20 : if (osStaFilename.size() != 0)
493 0 : papszFileList = CSLAddString( papszFileList, osStaFilename );
494 :
495 20 : return papszFileList;
496 : }
497 :
498 : /************************************************************************/
499 : /* GetEPSGGeogCS() */
500 : /* */
501 : /* Try to establish what the EPSG code for this coordinate */
502 : /* systems GEOGCS might be. Returns -1 if no reasonable guess */
503 : /* can be made. */
504 : /* */
505 : /* TODO: We really need to do some name lookups. */
506 : /************************************************************************/
507 :
508 92 : static int ENVIGetEPSGGeogCS( OGRSpatialReference *poThis )
509 :
510 : {
511 92 : const char *pszAuthName = poThis->GetAuthorityName( "GEOGCS" );
512 :
513 : /* -------------------------------------------------------------------- */
514 : /* Do we already have it? */
515 : /* -------------------------------------------------------------------- */
516 92 : if( pszAuthName != NULL && EQUAL(pszAuthName,"epsg") )
517 32 : return atoi(poThis->GetAuthorityCode( "GEOGCS" ));
518 :
519 : /* -------------------------------------------------------------------- */
520 : /* Get the datum and geogcs names. */
521 : /* -------------------------------------------------------------------- */
522 60 : const char *pszGEOGCS = poThis->GetAttrValue( "GEOGCS" );
523 60 : const char *pszDatum = poThis->GetAttrValue( "DATUM" );
524 :
525 : // We can only operate on coordinate systems with a geogcs.
526 60 : if( pszGEOGCS == NULL || pszDatum == NULL )
527 0 : return -1;
528 :
529 : /* -------------------------------------------------------------------- */
530 : /* Is this a "well known" geographic coordinate system? */
531 : /* -------------------------------------------------------------------- */
532 : int bWGS, bNAD;
533 :
534 : bWGS = strstr(pszGEOGCS,"WGS") != NULL
535 : || strstr(pszDatum, "WGS")
536 : || strstr(pszGEOGCS,"World Geodetic System")
537 : || strstr(pszGEOGCS,"World_Geodetic_System")
538 : || strstr(pszDatum, "World Geodetic System")
539 60 : || strstr(pszDatum, "World_Geodetic_System");
540 :
541 : bNAD = strstr(pszGEOGCS,"NAD") != NULL
542 : || strstr(pszDatum, "NAD")
543 : || strstr(pszGEOGCS,"North American")
544 : || strstr(pszGEOGCS,"North_American")
545 : || strstr(pszDatum, "North American")
546 60 : || strstr(pszDatum, "North_American");
547 :
548 60 : if( bWGS && (strstr(pszGEOGCS,"84") || strstr(pszDatum,"84")) )
549 50 : return 4326;
550 :
551 10 : if( bWGS && (strstr(pszGEOGCS,"72") || strstr(pszDatum,"72")) )
552 0 : return 4322;
553 :
554 10 : if( bNAD && (strstr(pszGEOGCS,"83") || strstr(pszDatum,"83")) )
555 2 : return 4269;
556 :
557 8 : if( bNAD && (strstr(pszGEOGCS,"27") || strstr(pszDatum,"27")) )
558 0 : return 4267;
559 :
560 : /* -------------------------------------------------------------------- */
561 : /* If we know the datum, associate the most likely GCS with */
562 : /* it. */
563 : /* -------------------------------------------------------------------- */
564 8 : pszAuthName = poThis->GetAuthorityName( "GEOGCS|DATUM" );
565 :
566 8 : if( pszAuthName != NULL
567 : && EQUAL(pszAuthName,"epsg")
568 : && poThis->GetPrimeMeridian() == 0.0 )
569 : {
570 0 : int nDatum = atoi(poThis->GetAuthorityCode("GEOGCS|DATUM"));
571 :
572 0 : if( nDatum >= 6000 && nDatum <= 6999 )
573 0 : return nDatum - 2000;
574 : }
575 :
576 8 : return -1;
577 : }
578 :
579 : /************************************************************************/
580 : /* WriteProjectionInfo() */
581 : /************************************************************************/
582 :
583 102 : void ENVIDataset::WriteProjectionInfo()
584 :
585 : {
586 : /* -------------------------------------------------------------------- */
587 : /* Format the location (geotransform) portion of the map info */
588 : /* line. */
589 : /* -------------------------------------------------------------------- */
590 102 : CPLString osLocation;
591 :
592 : osLocation.Printf( "1, 1, %.15g, %.15g, %.15g, %.15g",
593 : adfGeoTransform[0], adfGeoTransform[3],
594 102 : adfGeoTransform[1], fabs(adfGeoTransform[5]) );
595 :
596 : /* -------------------------------------------------------------------- */
597 : /* Minimal case - write out simple geotransform if we have a */
598 : /* non-default geotransform. */
599 : /* -------------------------------------------------------------------- */
600 102 : if( pszProjection == NULL || strlen(pszProjection) == 0 )
601 : {
602 42 : if( adfGeoTransform[0] != 0.0 || adfGeoTransform[1] != 1.0
603 16 : || adfGeoTransform[2] != 0.0 || adfGeoTransform[3] != 0.0
604 16 : || adfGeoTransform[4] != 0.0 || adfGeoTransform[5] != 1.0 )
605 : {
606 2 : const char* pszHemisphere = "North";
607 : VSIFPrintfL( fp, "map info = {Unknown, %s, %d, %s}\n",
608 2 : osLocation.c_str(), 0, pszHemisphere);
609 : }
610 : return;
611 : }
612 :
613 : /* -------------------------------------------------------------------- */
614 : /* Ingest WKT. */
615 : /* -------------------------------------------------------------------- */
616 92 : OGRSpatialReference oSRS;
617 :
618 92 : char *pszProj = pszProjection;
619 :
620 92 : if( oSRS.importFromWkt( &pszProj ) != OGRERR_NONE )
621 : return;
622 :
623 : /* -------------------------------------------------------------------- */
624 : /* Try to translate the datum and get major/minor ellipsoid */
625 : /* values. */
626 : /* -------------------------------------------------------------------- */
627 92 : int nEPSG_GCS = ENVIGetEPSGGeogCS( &oSRS );
628 92 : CPLString osDatum, osCommaDatum;
629 : double dfA, dfB;
630 :
631 92 : if( nEPSG_GCS == 4326 )
632 80 : osDatum = "WGS-84";
633 12 : else if( nEPSG_GCS == 4322 )
634 0 : osDatum = "WGS-72";
635 12 : else if( nEPSG_GCS == 4269 )
636 2 : osDatum = "North America 1983";
637 10 : else if( nEPSG_GCS == 4267 )
638 0 : osDatum = "North America 1927";
639 10 : else if( nEPSG_GCS == 4230 )
640 0 : osDatum = "European 1950";
641 10 : else if( nEPSG_GCS == 4277 )
642 2 : osDatum = "Ordnance Survey of Great Britain '36";
643 8 : else if( nEPSG_GCS == 4291 )
644 0 : osDatum = "SAD-69/Brazil";
645 8 : else if( nEPSG_GCS == 4283 )
646 0 : osDatum = "Geocentric Datum of Australia 1994";
647 8 : else if( nEPSG_GCS == 4275 )
648 0 : osDatum = "Nouvelle Triangulation Francaise IGN";
649 :
650 92 : if( osDatum != "" )
651 84 : osCommaDatum.Printf( ",%s", osDatum.c_str() );
652 :
653 92 : dfA = oSRS.GetSemiMajor();
654 92 : dfB = oSRS.GetSemiMinor();
655 :
656 : /* -------------------------------------------------------------------- */
657 : /* Do we have unusual linear units? */
658 : /* -------------------------------------------------------------------- */
659 92 : CPLString osOptionalUnits;
660 92 : if( fabs(oSRS.GetLinearUnits()-0.3048) < 0.0001 )
661 0 : osOptionalUnits = ", units=Feet";
662 :
663 : /* -------------------------------------------------------------------- */
664 : /* Handle UTM case. */
665 : /* -------------------------------------------------------------------- */
666 : const char *pszHemisphere;
667 92 : const char *pszProjName = oSRS.GetAttrValue("PROJECTION");
668 : int bNorth;
669 : int iUTMZone;
670 :
671 92 : iUTMZone = oSRS.GetUTMZone( &bNorth );
672 92 : if ( iUTMZone )
673 : {
674 0 : if ( bNorth )
675 0 : pszHemisphere = "North";
676 : else
677 0 : pszHemisphere = "South";
678 :
679 : VSIFPrintfL( fp, "map info = {UTM, %s, %d, %s%s%s}\n",
680 : osLocation.c_str(), iUTMZone, pszHemisphere,
681 0 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
682 : }
683 92 : else if( oSRS.IsGeographic() )
684 : {
685 : VSIFPrintfL( fp, "map info = {Geographic Lat/Lon, %s%s}\n",
686 80 : osLocation.c_str(), osCommaDatum.c_str());
687 : }
688 12 : else if( pszProjName == NULL )
689 : {
690 : // what to do?
691 : }
692 12 : else if( EQUAL(pszProjName,SRS_PT_NEW_ZEALAND_MAP_GRID) )
693 : {
694 : VSIFPrintfL( fp, "map info = {New Zealand Map Grid, %s%s%s}\n",
695 : osLocation.c_str(),
696 0 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
697 :
698 : VSIFPrintfL( fp, "projection info = {39, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, New Zealand Map Grid}\n",
699 : dfA, dfB,
700 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
701 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
702 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
703 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
704 0 : osCommaDatum.c_str() );
705 : }
706 12 : else if( EQUAL(pszProjName,SRS_PT_TRANSVERSE_MERCATOR) )
707 : {
708 : VSIFPrintfL( fp, "map info = {Transverse Mercator, %s%s%s}\n",
709 : osLocation.c_str(),
710 2 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
711 :
712 : VSIFPrintfL( fp, "projection info = {3, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Transverse Mercator}\n",
713 : dfA, dfB,
714 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
715 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
716 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
717 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
718 : oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR,1.0),
719 2 : osCommaDatum.c_str() );
720 : }
721 12 : else if( EQUAL(pszProjName,SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP)
722 : || EQUAL(pszProjName,SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP_BELGIUM) )
723 : {
724 : VSIFPrintfL( fp, "map info = {Lambert Conformal Conic, %s%s%s}\n",
725 : osLocation.c_str(),
726 2 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
727 :
728 : VSIFPrintfL( fp, "projection info = {4, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Lambert Conformal Conic}\n",
729 : dfA, dfB,
730 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
731 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
732 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
733 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
734 : oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0),
735 : oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2,0.0),
736 2 : osCommaDatum.c_str() );
737 : }
738 8 : else if( EQUAL(pszProjName,
739 : SRS_PT_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN) )
740 : {
741 : VSIFPrintfL( fp, "map info = {Hotine Oblique Mercator A, %s%s%s}\n",
742 : osLocation.c_str(),
743 0 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
744 :
745 : VSIFPrintfL( fp, "projection info = {5, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Hotine Oblique Mercator A}\n",
746 : dfA, dfB,
747 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
748 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_POINT_1,0.0),
749 : oSRS.GetNormProjParm(SRS_PP_LONGITUDE_OF_POINT_1,0.0),
750 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_POINT_2,0.0),
751 : oSRS.GetNormProjParm(SRS_PP_LONGITUDE_OF_POINT_2,0.0),
752 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
753 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
754 : oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR,1.0),
755 0 : osCommaDatum.c_str() );
756 : }
757 8 : else if( EQUAL(pszProjName,SRS_PT_HOTINE_OBLIQUE_MERCATOR) )
758 : {
759 : VSIFPrintfL( fp, "map info = {Hotine Oblique Mercator B, %s%s%s}\n",
760 : osLocation.c_str(),
761 0 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
762 :
763 : VSIFPrintfL( fp, "projection info = {6, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Hotine Oblique Mercator B}\n",
764 : dfA, dfB,
765 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
766 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
767 : oSRS.GetNormProjParm(SRS_PP_AZIMUTH,0.0),
768 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
769 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
770 : oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR,1.0),
771 0 : osCommaDatum.c_str() );
772 : }
773 8 : else if( EQUAL(pszProjName,SRS_PT_STEREOGRAPHIC)
774 : || EQUAL(pszProjName,SRS_PT_OBLIQUE_STEREOGRAPHIC) )
775 : {
776 : VSIFPrintfL( fp, "map info = {Stereographic (ellipsoid), %s%s%s}\n",
777 : osLocation.c_str(),
778 0 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
779 :
780 : VSIFPrintfL( fp, "projection info = {7, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %s, Stereographic (ellipsoid)}\n",
781 : dfA, dfB,
782 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
783 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
784 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
785 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
786 : oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR,1.0),
787 0 : osCommaDatum.c_str() );
788 : }
789 8 : else if( EQUAL(pszProjName,SRS_PT_ALBERS_CONIC_EQUAL_AREA) )
790 : {
791 : VSIFPrintfL( fp, "map info = {Albers Conical Equal Area, %s%s%s}\n",
792 : osLocation.c_str(),
793 6 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
794 :
795 : VSIFPrintfL( fp, "projection info = {9, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Albers Conical Equal Area}\n",
796 : dfA, dfB,
797 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
798 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
799 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
800 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
801 : oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0),
802 : oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2,0.0),
803 6 : osCommaDatum.c_str() );
804 : }
805 2 : else if( EQUAL(pszProjName,SRS_PT_POLYCONIC) )
806 : {
807 : VSIFPrintfL( fp, "map info = {Polyconic, %s%s%s}\n",
808 : osLocation.c_str(),
809 0 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
810 :
811 : VSIFPrintfL( fp, "projection info = {10, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Polyconic}\n",
812 : dfA, dfB,
813 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
814 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
815 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
816 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
817 0 : osCommaDatum.c_str() );
818 : }
819 2 : else if( EQUAL(pszProjName,SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA) )
820 : {
821 : VSIFPrintfL( fp, "map info = {Lambert Azimuthal Equal Area, %s%s%s}\n",
822 : osLocation.c_str(),
823 2 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
824 :
825 : VSIFPrintfL( fp, "projection info = {11, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Lambert Azimuthal Equal Area}\n",
826 : dfA, dfB,
827 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
828 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
829 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
830 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
831 2 : osCommaDatum.c_str() );
832 : }
833 0 : else if( EQUAL(pszProjName,SRS_PT_AZIMUTHAL_EQUIDISTANT) )
834 : {
835 : VSIFPrintfL( fp, "map info = {Azimuthal Equadistant, %s%s%s}\n",
836 : osLocation.c_str(),
837 0 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
838 :
839 : VSIFPrintfL( fp, "projection info = {12, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Azimuthal Equadistant}\n",
840 : dfA, dfB,
841 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
842 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
843 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
844 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
845 0 : osCommaDatum.c_str() );
846 : }
847 0 : else if( EQUAL(pszProjName,SRS_PT_POLAR_STEREOGRAPHIC) )
848 : {
849 : VSIFPrintfL( fp, "map info = {Polar Stereographic, %s%s%s}\n",
850 : osLocation.c_str(),
851 0 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
852 :
853 : VSIFPrintfL( fp, "projection info = {31, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Polar Stereographic}\n",
854 : dfA, dfB,
855 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,90.0),
856 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
857 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
858 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
859 0 : osCommaDatum.c_str() );
860 : }
861 : else
862 : {
863 : VSIFPrintfL( fp, "map info = {%s, %s}\n",
864 0 : pszProjName, osLocation.c_str());
865 : }
866 :
867 : // write out coordinate system string
868 92 : if ( oSRS.morphToESRI() == OGRERR_NONE )
869 : {
870 92 : char *pszProjESRI = NULL;
871 92 : if ( oSRS.exportToWkt(&pszProjESRI) == OGRERR_NONE )
872 : {
873 92 : if ( strlen(pszProjESRI) )
874 92 : VSIFPrintfL( fp, "coordinate system string = {%s}\n", pszProjESRI);
875 : }
876 92 : CPLFree(pszProjESRI);
877 92 : pszProjESRI = NULL;
878 92 : }
879 : }
880 :
881 : /************************************************************************/
882 : /* ParseRpcCoeffsMetaDataString() */
883 : /************************************************************************/
884 :
885 32 : int ENVIDataset::ParseRpcCoeffsMetaDataString(const char *psName, char **papszVal,
886 : int& idx)
887 : {
888 : // separate one string with 20 coefficients into an array of 20 strings.
889 32 : const char *psz20Vals = GetMetadataItem(psName, "RPC");
890 32 : if (!psz20Vals)
891 0 : return FALSE;
892 :
893 32 : char** papszArr = CSLTokenizeString2(psz20Vals, " ", 0);
894 32 : if (!papszArr)
895 0 : return FALSE;
896 :
897 32 : int x = 0;
898 704 : while ((papszArr[x] != NULL) && (x < 20))
899 : {
900 640 : papszVal[idx++] = CPLStrdup(papszArr[x]);
901 640 : x++;
902 : }
903 :
904 32 : CSLDestroy(papszArr);
905 :
906 32 : return (x == 20);
907 : }
908 :
909 : #define CPLStrdupIfNotNull(x) ((x) ? CPLStrdup(x) : NULL)
910 :
911 : /************************************************************************/
912 : /* WriteRpcInfo() */
913 : /************************************************************************/
914 :
915 110 : int ENVIDataset::WriteRpcInfo()
916 : {
917 : // write out 90 rpc coeffs into the envi header plus 3 envi specific rpc values
918 : // returns 0 if the coeffs are not present or not valid
919 110 : int bRet = FALSE;
920 110 : int x, idx = 0;
921 : char* papszVal[93];
922 :
923 110 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("LINE_OFF", "RPC"));
924 110 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("SAMP_OFF", "RPC"));
925 110 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("LAT_OFF", "RPC"));
926 110 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("LONG_OFF", "RPC"));
927 110 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("HEIGHT_OFF", "RPC"));
928 110 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("LINE_SCALE", "RPC"));
929 110 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("SAMP_SCALE", "RPC"));
930 110 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("LAT_SCALE", "RPC"));
931 110 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("LONG_SCALE", "RPC"));
932 110 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("HEIGHT_SCALE", "RPC"));
933 :
934 190 : for (x=0; x<10; x++) // if we do not have 10 values we return 0
935 : {
936 182 : if (!papszVal[x])
937 102 : goto end;
938 : }
939 :
940 8 : if (!ParseRpcCoeffsMetaDataString("LINE_NUM_COEFF", papszVal, idx))
941 0 : goto end;
942 :
943 8 : if (!ParseRpcCoeffsMetaDataString("LINE_DEN_COEFF", papszVal, idx))
944 0 : goto end;
945 :
946 8 : if (!ParseRpcCoeffsMetaDataString("SAMP_NUM_COEFF", papszVal, idx))
947 0 : goto end;
948 :
949 8 : if (!ParseRpcCoeffsMetaDataString("SAMP_DEN_COEFF", papszVal, idx))
950 0 : goto end;
951 :
952 8 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("TILE_ROW_OFFSET", "RPC"));
953 8 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("TILE_COL_OFFSET", "RPC"));
954 8 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("ENVI_RPC_EMULATION", "RPC"));
955 8 : CPLAssert(idx == 93);
956 32 : for (x=90; x<93; x++)
957 : {
958 24 : if (!papszVal[x])
959 0 : goto end;
960 : }
961 :
962 : // ok all the needed 93 values are present so write the rpcs into the envi header
963 8 : x = 1;
964 8 : VSIFPrintfL(fp, "rpc info = {\n");
965 752 : for (int iR=0; iR<93; iR++)
966 : {
967 744 : if (papszVal[iR][0] == '-')
968 64 : VSIFPrintfL(fp, " %s", papszVal[iR]);
969 : else
970 680 : VSIFPrintfL(fp, " %s", papszVal[iR]);
971 :
972 744 : if (iR<92)
973 736 : VSIFPrintfL(fp, ",");
974 :
975 744 : if ((x % 4) == 0)
976 184 : VSIFPrintfL(fp, "\n");
977 :
978 744 : x++;
979 744 : if (x > 4)
980 184 : x = 1;
981 : }
982 :
983 8 : VSIFPrintfL(fp, "}\n" );
984 :
985 8 : bRet = TRUE;
986 :
987 : end:
988 1874 : for (x=0;x<idx;x++)
989 1764 : CPLFree(papszVal[x]);
990 :
991 110 : return bRet;
992 : }
993 :
994 : /************************************************************************/
995 : /* WritePseudoGcpInfo() */
996 : /************************************************************************/
997 :
998 102 : int ENVIDataset::WritePseudoGcpInfo()
999 : {
1000 : // write out gcps into the envi header
1001 : // returns 0 if the gcps are not present
1002 :
1003 102 : int iNum = GetGCPCount();
1004 102 : if (iNum == 0)
1005 102 : return FALSE;
1006 :
1007 0 : const GDAL_GCP *pGcpStructs = GetGCPs();
1008 :
1009 : // double dfGCPPixel; /** Pixel (x) location of GCP on raster */
1010 : // double dfGCPLine; /** Line (y) location of GCP on raster */
1011 : // double dfGCPX; /** X position of GCP in georeferenced space */
1012 : // double dfGCPY; /** Y position of GCP in georeferenced space */
1013 :
1014 0 : VSIFPrintfL(fp, "geo points = {\n");
1015 0 : for (int iR=0; iR<iNum; iR++)
1016 : {
1017 : VSIFPrintfL(fp, " %#0.4f, %#0.4f, %#0.8f, %#0.8f",
1018 0 : pGcpStructs[iR].dfGCPPixel, pGcpStructs[iR].dfGCPLine,
1019 0 : pGcpStructs[iR].dfGCPY, pGcpStructs[iR].dfGCPX);
1020 0 : if (iR<iNum-1)
1021 0 : VSIFPrintfL(fp, ",\n");
1022 : }
1023 :
1024 0 : VSIFPrintfL(fp, "}\n" );
1025 :
1026 0 : return TRUE;
1027 : }
1028 :
1029 : /************************************************************************/
1030 : /* GetProjectionRef() */
1031 : /************************************************************************/
1032 :
1033 40 : const char *ENVIDataset::GetProjectionRef()
1034 :
1035 : {
1036 40 : return pszProjection;
1037 : }
1038 :
1039 : /************************************************************************/
1040 : /* SetProjection() */
1041 : /************************************************************************/
1042 :
1043 92 : CPLErr ENVIDataset::SetProjection( const char *pszNewProjection )
1044 :
1045 : {
1046 92 : CPLFree( pszProjection );
1047 92 : pszProjection = CPLStrdup( pszNewProjection );
1048 :
1049 92 : bHeaderDirty = TRUE;
1050 :
1051 92 : return CE_None;
1052 : }
1053 :
1054 : /************************************************************************/
1055 : /* GetGeoTransform() */
1056 : /************************************************************************/
1057 :
1058 76 : CPLErr ENVIDataset::GetGeoTransform( double * padfTransform )
1059 :
1060 : {
1061 76 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
1062 :
1063 76 : if( bFoundMapinfo )
1064 72 : return CE_None;
1065 : else
1066 4 : return CE_Failure;
1067 : }
1068 :
1069 : /************************************************************************/
1070 : /* SetGeoTransform() */
1071 : /************************************************************************/
1072 :
1073 94 : CPLErr ENVIDataset::SetGeoTransform( double * padfTransform )
1074 : {
1075 94 : memcpy( adfGeoTransform, padfTransform, sizeof(double) * 6 );
1076 :
1077 94 : bHeaderDirty = TRUE;
1078 94 : bFoundMapinfo = TRUE;
1079 :
1080 94 : return CE_None;
1081 : }
1082 :
1083 : /************************************************************************/
1084 : /* SetMetadata() */
1085 : /************************************************************************/
1086 :
1087 42 : CPLErr ENVIDataset::SetMetadata( char ** papszMetadata,
1088 : const char * pszDomain )
1089 : {
1090 42 : if (pszDomain && EQUAL(pszDomain, "RPC"))
1091 2 : bHeaderDirty = TRUE;
1092 :
1093 42 : return RawDataset::SetMetadata(papszMetadata, pszDomain);
1094 : }
1095 :
1096 : /************************************************************************/
1097 : /* SetMetadataItem() */
1098 : /************************************************************************/
1099 :
1100 524 : CPLErr ENVIDataset::SetMetadataItem( const char * pszName,
1101 : const char * pszValue,
1102 : const char * pszDomain )
1103 : {
1104 524 : if (pszDomain && EQUAL(pszDomain, "RPC"))
1105 126 : bHeaderDirty = TRUE;
1106 :
1107 524 : return RawDataset::SetMetadataItem(pszName, pszValue, pszDomain);
1108 : }
1109 :
1110 : /************************************************************************/
1111 : /* SetGCPs() */
1112 : /************************************************************************/
1113 :
1114 0 : CPLErr ENVIDataset::SetGCPs( int nGCPCount, const GDAL_GCP *pasGCPList,
1115 : const char *pszGCPProjection )
1116 : {
1117 0 : bHeaderDirty = TRUE;
1118 :
1119 0 : return RawDataset::SetGCPs(nGCPCount, pasGCPList, pszGCPProjection);
1120 : }
1121 :
1122 : /************************************************************************/
1123 : /* SplitList() */
1124 : /* */
1125 : /* Split an ENVI value list into component fields, and strip */
1126 : /* white space. */
1127 : /************************************************************************/
1128 :
1129 392 : char **ENVIDataset::SplitList( const char *pszCleanInput )
1130 :
1131 : {
1132 392 : char **papszReturn = NULL;
1133 392 : char *pszInput = CPLStrdup(pszCleanInput);
1134 :
1135 392 : if( pszInput[0] != '{' )
1136 : {
1137 116 : CPLFree(pszInput);
1138 116 : return NULL;
1139 : }
1140 :
1141 276 : int iChar=1;
1142 :
1143 :
1144 2528 : while( pszInput[iChar] != '}' && pszInput[iChar] != '\0' )
1145 : {
1146 1976 : int iFStart=-1, iFEnd=-1;
1147 :
1148 : // Find start of token.
1149 1976 : iFStart = iChar;
1150 6056 : while( pszInput[iFStart] == ' ' )
1151 2104 : iFStart++;
1152 :
1153 1976 : iFEnd = iFStart;
1154 39268 : while( pszInput[iFEnd] != ','
1155 11956 : && pszInput[iFEnd] != '}'
1156 11680 : && pszInput[iFEnd] != '\0' )
1157 11680 : iFEnd++;
1158 :
1159 1976 : if( pszInput[iFEnd] == '\0' )
1160 0 : break;
1161 :
1162 1976 : iChar = iFEnd + 1;
1163 1976 : iFEnd = iFEnd - 1;
1164 :
1165 3952 : while( iFEnd > iFStart && pszInput[iFEnd] == ' ' )
1166 0 : iFEnd--;
1167 :
1168 1976 : pszInput[iFEnd + 1] = '\0';
1169 1976 : papszReturn = CSLAddString( papszReturn, pszInput + iFStart );
1170 : }
1171 :
1172 276 : CPLFree( pszInput );
1173 :
1174 276 : return papszReturn;
1175 : }
1176 :
1177 : /************************************************************************/
1178 : /* SetENVIDatum() */
1179 : /************************************************************************/
1180 :
1181 0 : void ENVIDataset::SetENVIDatum( OGRSpatialReference *poSRS,
1182 : const char *pszENVIDatumName )
1183 :
1184 : {
1185 : // datums
1186 0 : if( EQUAL(pszENVIDatumName, "WGS-84") )
1187 0 : poSRS->SetWellKnownGeogCS( "WGS84" );
1188 0 : else if( EQUAL(pszENVIDatumName, "WGS-72") )
1189 0 : poSRS->SetWellKnownGeogCS( "WGS72" );
1190 0 : else if( EQUAL(pszENVIDatumName, "North America 1983") )
1191 0 : poSRS->SetWellKnownGeogCS( "NAD83" );
1192 0 : else if( EQUAL(pszENVIDatumName, "North America 1927")
1193 : || strstr(pszENVIDatumName,"NAD27")
1194 : || strstr(pszENVIDatumName,"NAD-27") )
1195 0 : poSRS->SetWellKnownGeogCS( "NAD27" );
1196 0 : else if( EQUALN(pszENVIDatumName, "European 1950",13) )
1197 0 : poSRS->SetWellKnownGeogCS( "EPSG:4230" );
1198 0 : else if( EQUAL(pszENVIDatumName, "Ordnance Survey of Great Britain '36") )
1199 0 : poSRS->SetWellKnownGeogCS( "EPSG:4277" );
1200 0 : else if( EQUAL(pszENVIDatumName, "SAD-69/Brazil") )
1201 0 : poSRS->SetWellKnownGeogCS( "EPSG:4291" );
1202 0 : else if( EQUAL(pszENVIDatumName, "Geocentric Datum of Australia 1994") )
1203 0 : poSRS->SetWellKnownGeogCS( "EPSG:4283" );
1204 0 : else if( EQUAL(pszENVIDatumName, "Australian Geodetic 1984") )
1205 0 : poSRS->SetWellKnownGeogCS( "EPSG:4203" );
1206 0 : else if( EQUAL(pszENVIDatumName, "Nouvelle Triangulation Francaise IGN") )
1207 0 : poSRS->SetWellKnownGeogCS( "EPSG:4275" );
1208 :
1209 : // Ellipsoids
1210 0 : else if( EQUAL(pszENVIDatumName, "GRS 80") )
1211 0 : poSRS->SetWellKnownGeogCS( "NAD83" );
1212 0 : else if( EQUAL(pszENVIDatumName, "Airy") )
1213 0 : poSRS->SetWellKnownGeogCS( "EPSG:4001" );
1214 0 : else if( EQUAL(pszENVIDatumName, "Australian National") )
1215 0 : poSRS->SetWellKnownGeogCS( "EPSG:4003" );
1216 0 : else if( EQUAL(pszENVIDatumName, "Bessel 1841") )
1217 0 : poSRS->SetWellKnownGeogCS( "EPSG:4004" );
1218 0 : else if( EQUAL(pszENVIDatumName, "Clark 1866") )
1219 0 : poSRS->SetWellKnownGeogCS( "EPSG:4008" );
1220 : else
1221 : {
1222 : CPLError( CE_Warning, CPLE_AppDefined,
1223 0 : "Unrecognised datum '%s', defaulting to WGS84.", pszENVIDatumName);
1224 0 : poSRS->SetWellKnownGeogCS( "WGS84" );
1225 : }
1226 0 : }
1227 :
1228 : /************************************************************************/
1229 : /* SetENVIEllipse() */
1230 : /************************************************************************/
1231 :
1232 14 : void ENVIDataset::SetENVIEllipse( OGRSpatialReference *poSRS,
1233 : char **papszPI_EI )
1234 :
1235 : {
1236 14 : double dfA = CPLAtofM(papszPI_EI[0]);
1237 14 : double dfB = CPLAtofM(papszPI_EI[1]);
1238 : double dfInvF;
1239 :
1240 14 : if( fabs(dfA-dfB) < 0.1 )
1241 0 : dfInvF = 0.0; // sphere
1242 : else
1243 14 : dfInvF = dfA / (dfA - dfB);
1244 :
1245 :
1246 : poSRS->SetGeogCS( "Ellipse Based", "Ellipse Based", "Unnamed",
1247 14 : dfA, dfInvF );
1248 14 : }
1249 :
1250 : /************************************************************************/
1251 : /* ProcessMapinfo() */
1252 : /* */
1253 : /* Extract projection, and geotransform from a mapinfo value in */
1254 : /* the header. */
1255 : /************************************************************************/
1256 :
1257 104 : int ENVIDataset::ProcessMapinfo( const char *pszMapinfo )
1258 :
1259 : {
1260 : char **papszFields;
1261 : int nCount;
1262 104 : OGRSpatialReference oSRS;
1263 :
1264 104 : papszFields = SplitList( pszMapinfo );
1265 104 : nCount = CSLCount(papszFields);
1266 :
1267 104 : if( nCount < 7 )
1268 : {
1269 0 : CSLDestroy( papszFields );
1270 0 : return FALSE;
1271 : }
1272 :
1273 : /* -------------------------------------------------------------------- */
1274 : /* Check if we have coordinate system string, and if so parse it. */
1275 : /* -------------------------------------------------------------------- */
1276 104 : char **papszCSS = NULL;
1277 104 : if( CSLFetchNameValue( papszHeader, "coordinate_system_string" ) != NULL )
1278 : {
1279 : papszCSS = CSLTokenizeString2(
1280 : CSLFetchNameValue( papszHeader, "coordinate_system_string" ),
1281 74 : "{}", CSLT_PRESERVEQUOTES );
1282 : }
1283 :
1284 : /* -------------------------------------------------------------------- */
1285 : /* Check if we have projection info, and if so parse it. */
1286 : /* -------------------------------------------------------------------- */
1287 104 : char **papszPI = NULL;
1288 104 : int nPICount = 0;
1289 104 : if( CSLFetchNameValue( papszHeader, "projection_info" ) != NULL )
1290 : {
1291 : papszPI = SplitList(
1292 38 : CSLFetchNameValue( papszHeader, "projection_info" ) );
1293 38 : nPICount = CSLCount(papszPI);
1294 : }
1295 :
1296 : /* -------------------------------------------------------------------- */
1297 : /* Capture geotransform. */
1298 : /* -------------------------------------------------------------------- */
1299 104 : adfGeoTransform[1] = atof(papszFields[5]); // Pixel width
1300 104 : adfGeoTransform[5] = -atof(papszFields[6]); // Pixel height
1301 : adfGeoTransform[0] = // Upper left X coordinate
1302 104 : atof(papszFields[3]) - (atof(papszFields[1]) - 1) * adfGeoTransform[1];
1303 : adfGeoTransform[3] = // Upper left Y coordinate
1304 104 : atof(papszFields[4]) - (atof(papszFields[2]) - 1) * adfGeoTransform[5];
1305 104 : adfGeoTransform[2] = 0.0;
1306 104 : adfGeoTransform[4] = 0.0;
1307 :
1308 : /* -------------------------------------------------------------------- */
1309 : /* Capture projection. */
1310 : /* -------------------------------------------------------------------- */
1311 104 : if ( oSRS.importFromESRI( papszCSS ) != OGRERR_NONE )
1312 : {
1313 30 : oSRS.Clear();
1314 :
1315 44 : if( EQUALN(papszFields[0],"UTM",3) && nCount >= 9 )
1316 : {
1317 14 : oSRS.SetUTM( atoi(papszFields[7]),
1318 28 : !EQUAL(papszFields[8],"South") );
1319 14 : if( nCount >= 10 && strstr(papszFields[9],"=") == NULL )
1320 0 : SetENVIDatum( &oSRS, papszFields[9] );
1321 : else
1322 14 : oSRS.SetWellKnownGeogCS( "NAD27" );
1323 : }
1324 16 : else if( EQUALN(papszFields[0],"State Plane (NAD 27)",19)
1325 : && nCount >= 7 )
1326 : {
1327 0 : oSRS.SetStatePlane( ITTVISToUSGSZone(atoi(papszFields[7])), FALSE );
1328 : }
1329 16 : else if( EQUALN(papszFields[0],"State Plane (NAD 83)",19)
1330 : && nCount >= 7 )
1331 : {
1332 0 : oSRS.SetStatePlane( ITTVISToUSGSZone(atoi(papszFields[7])), TRUE );
1333 : }
1334 16 : else if( EQUALN(papszFields[0],"Geographic Lat",14)
1335 : && nCount >= 8 )
1336 : {
1337 0 : if( nCount >= 8 && strstr(papszFields[7],"=") == NULL )
1338 0 : SetENVIDatum( &oSRS, papszFields[7] );
1339 : else
1340 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
1341 : }
1342 16 : else if( nPICount > 8 && atoi(papszPI[0]) == 3 ) // TM
1343 : {
1344 0 : oSRS.SetTM( CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1345 0 : CPLAtofM(papszPI[7]),
1346 0 : CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) );
1347 : }
1348 16 : else if( nPICount > 8 && atoi(papszPI[0]) == 4 ) // Lambert Conformal Conic
1349 : {
1350 0 : oSRS.SetLCC( CPLAtofM(papszPI[7]), CPLAtofM(papszPI[8]),
1351 0 : CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1352 0 : CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) );
1353 : }
1354 16 : else if( nPICount > 10 && atoi(papszPI[0]) == 5 ) // Oblique Merc (2 point)
1355 : {
1356 0 : oSRS.SetHOM2PNO( CPLAtofM(papszPI[3]),
1357 0 : CPLAtofM(papszPI[4]), CPLAtofM(papszPI[5]),
1358 0 : CPLAtofM(papszPI[6]), CPLAtofM(papszPI[7]),
1359 0 : CPLAtofM(papszPI[10]),
1360 0 : CPLAtofM(papszPI[8]), CPLAtofM(papszPI[9]) );
1361 : }
1362 16 : else if( nPICount > 8 && atoi(papszPI[0]) == 6 ) // Oblique Merc
1363 : {
1364 0 : oSRS.SetHOM(CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1365 0 : CPLAtofM(papszPI[5]), 0.0,
1366 0 : CPLAtofM(papszPI[8]),
1367 0 : CPLAtofM(papszPI[6]), CPLAtofM(papszPI[7]) );
1368 : }
1369 16 : else if( nPICount > 8 && atoi(papszPI[0]) == 7 ) // Stereographic
1370 : {
1371 0 : oSRS.SetStereographic( CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1372 0 : CPLAtofM(papszPI[7]),
1373 0 : CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) );
1374 : }
1375 30 : else if( nPICount > 8 && atoi(papszPI[0]) == 9 ) // Albers Equal Area
1376 : {
1377 28 : oSRS.SetACEA( CPLAtofM(papszPI[7]), CPLAtofM(papszPI[8]),
1378 28 : CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1379 70 : CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) );
1380 : }
1381 2 : else if( nPICount > 6 && atoi(papszPI[0]) == 10 ) // Polyconic
1382 : {
1383 0 : oSRS.SetPolyconic(CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1384 0 : CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) );
1385 : }
1386 2 : else if( nPICount > 6 && atoi(papszPI[0]) == 11 ) // LAEA
1387 : {
1388 0 : oSRS.SetLAEA(CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1389 0 : CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) );
1390 : }
1391 2 : else if( nPICount > 6 && atoi(papszPI[0]) == 12 ) // Azimuthal Equid.
1392 : {
1393 0 : oSRS.SetAE(CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1394 0 : CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) );
1395 : }
1396 2 : else if( nPICount > 6 && atoi(papszPI[0]) == 31 ) // Polar Stereographic
1397 : {
1398 0 : oSRS.SetPS(CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1399 : 1.0,
1400 0 : CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) );
1401 : }
1402 : }
1403 :
1404 104 : CSLDestroy( papszCSS );
1405 :
1406 : // Still lots more that could be added for someone with the patience.
1407 :
1408 : /* -------------------------------------------------------------------- */
1409 : /* fallback to localcs if we don't recognise things. */
1410 : /* -------------------------------------------------------------------- */
1411 104 : if( oSRS.GetRoot() == NULL )
1412 2 : oSRS.SetLocalCS( papszFields[0] );
1413 :
1414 : /* -------------------------------------------------------------------- */
1415 : /* Try to set datum from projection info line if we have a */
1416 : /* projected coordinate system without a GEOGCS. */
1417 : /* -------------------------------------------------------------------- */
1418 104 : if( oSRS.IsProjected() && oSRS.GetAttrNode("GEOGCS") == NULL
1419 : && nPICount > 3 )
1420 : {
1421 : // Do we have a datum on the projection info line?
1422 14 : int iDatum = nPICount-1;
1423 :
1424 : // Ignore units= items.
1425 14 : if( strstr(papszPI[iDatum],"=") != NULL )
1426 0 : iDatum--;
1427 :
1428 : // Skip past the name.
1429 14 : iDatum--;
1430 :
1431 14 : CPLString osDatumName = papszPI[iDatum];
1432 14 : if( osDatumName.find_first_of("abcdefghijklmnopqrstuvwxyz"
1433 : "ABCDEFGHIJKLMNOPQRSTUVWXYZ")
1434 : != CPLString::npos )
1435 : {
1436 0 : SetENVIDatum( &oSRS, osDatumName );
1437 : }
1438 : else
1439 : {
1440 14 : SetENVIEllipse( &oSRS, papszPI + 1 );
1441 14 : }
1442 : }
1443 :
1444 :
1445 : /* -------------------------------------------------------------------- */
1446 : /* Try to process specialized units. */
1447 : /* -------------------------------------------------------------------- */
1448 104 : if( EQUALN( papszFields[nCount-1],"units",5))
1449 : {
1450 : /* Handle linear units first. */
1451 0 : if (EQUAL(papszFields[nCount-1],"units=Feet") )
1452 0 : oSRS.SetLinearUnitsAndUpdateParameters( SRS_UL_FOOT, atof(SRS_UL_FOOT_CONV) );
1453 0 : else if (EQUAL(papszFields[nCount-1],"units=Meters") )
1454 0 : oSRS.SetLinearUnitsAndUpdateParameters( SRS_UL_METER, 1. );
1455 0 : else if (EQUAL(papszFields[nCount-1],"units=Km") )
1456 0 : oSRS.SetLinearUnitsAndUpdateParameters( "Kilometer", 1000. );
1457 0 : else if (EQUAL(papszFields[nCount-1],"units=Yards") )
1458 0 : oSRS.SetLinearUnitsAndUpdateParameters( "Yard", .9144 );
1459 0 : else if (EQUAL(papszFields[nCount-1],"units=Miles") )
1460 0 : oSRS.SetLinearUnitsAndUpdateParameters( "Mile", 1609.344 );
1461 0 : else if (EQUAL(papszFields[nCount-1],"units=Nautical Miles") )
1462 0 : oSRS.SetLinearUnitsAndUpdateParameters( SRS_UL_NAUTICAL_MILE, atof(SRS_UL_NAUTICAL_MILE_CONV) );
1463 :
1464 : /* Only handle angular units if we know the projection is geographic. */
1465 0 : if (oSRS.IsGeographic())
1466 : {
1467 0 : if (EQUAL(papszFields[nCount-1],"units=Radians") )
1468 0 : oSRS.SetAngularUnits( SRS_UA_RADIAN, 1. );
1469 : else
1470 : {
1471 : /* Degrees, minutes and seconds will all be represented as degrees. */
1472 0 : oSRS.SetAngularUnits( SRS_UA_DEGREE, atof(SRS_UA_DEGREE_CONV));
1473 :
1474 0 : double conversionFactor = 1.;
1475 0 : if (EQUAL(papszFields[nCount-1],"units=Minutes") )
1476 0 : conversionFactor = 60.;
1477 0 : else if( EQUAL(papszFields[nCount-1],"units=Seconds") )
1478 0 : conversionFactor = 3600.;
1479 0 : adfGeoTransform[0] /= conversionFactor;
1480 0 : adfGeoTransform[1] /= conversionFactor;
1481 0 : adfGeoTransform[2] /= conversionFactor;
1482 0 : adfGeoTransform[3] /= conversionFactor;
1483 0 : adfGeoTransform[4] /= conversionFactor;
1484 0 : adfGeoTransform[5] /= conversionFactor;
1485 : }
1486 : }
1487 : }
1488 : /* -------------------------------------------------------------------- */
1489 : /* Turn back into WKT. */
1490 : /* -------------------------------------------------------------------- */
1491 104 : if( oSRS.GetRoot() != NULL )
1492 : {
1493 104 : oSRS.Fixup();
1494 104 : if ( pszProjection )
1495 : {
1496 104 : CPLFree( pszProjection );
1497 104 : pszProjection = NULL;
1498 : }
1499 104 : oSRS.exportToWkt( &pszProjection );
1500 : }
1501 :
1502 104 : CSLDestroy( papszFields );
1503 104 : CSLDestroy( papszPI );
1504 104 : return TRUE;
1505 : }
1506 :
1507 : /************************************************************************/
1508 : /* ProcessRPCinfo() */
1509 : /* */
1510 : /* Extract RPC transformation coefficients if they are present */
1511 : /* and sets into the standard metadata fields for RPC. */
1512 : /************************************************************************/
1513 :
1514 6 : void ENVIDataset::ProcessRPCinfo( const char *pszRPCinfo,
1515 : int numCols, int numRows)
1516 : {
1517 : char **papszFields;
1518 : char sVal[1280];
1519 : int nCount;
1520 :
1521 6 : papszFields = SplitList( pszRPCinfo );
1522 6 : nCount = CSLCount(papszFields);
1523 :
1524 6 : if( nCount < 90 )
1525 : {
1526 0 : CSLDestroy( papszFields );
1527 0 : return;
1528 : }
1529 :
1530 6 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[0]));
1531 6 : SetMetadataItem("LINE_OFF",sVal,"RPC");
1532 6 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[5]));
1533 6 : SetMetadataItem("LINE_SCALE",sVal,"RPC");
1534 6 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[1]));
1535 6 : SetMetadataItem("SAMP_OFF",sVal,"RPC");
1536 6 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[6]));
1537 6 : SetMetadataItem("SAMP_SCALE",sVal,"RPC");
1538 6 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[2]));
1539 6 : SetMetadataItem("LAT_OFF",sVal,"RPC");
1540 6 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[7]));
1541 6 : SetMetadataItem("LAT_SCALE",sVal,"RPC");
1542 6 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[3]));
1543 6 : SetMetadataItem("LONG_OFF",sVal,"RPC");
1544 6 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[8]));
1545 6 : SetMetadataItem("LONG_SCALE",sVal,"RPC");
1546 6 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[4]));
1547 6 : SetMetadataItem("HEIGHT_OFF",sVal,"RPC");
1548 6 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[9]));
1549 6 : SetMetadataItem("HEIGHT_SCALE",sVal,"RPC");
1550 :
1551 6 : sVal[0] = '\0';
1552 : int i;
1553 126 : for(i = 0; i < 20; i++ )
1554 : snprintf(sVal+strlen(sVal), sizeof(sVal), "%.16g ",
1555 120 : atof(papszFields[10+i]));
1556 6 : SetMetadataItem("LINE_NUM_COEFF",sVal,"RPC");
1557 :
1558 6 : sVal[0] = '\0';
1559 126 : for(i = 0; i < 20; i++ )
1560 : snprintf(sVal+strlen(sVal), sizeof(sVal), "%.16g ",
1561 120 : atof(papszFields[30+i]));
1562 6 : SetMetadataItem("LINE_DEN_COEFF",sVal,"RPC");
1563 :
1564 6 : sVal[0] = '\0';
1565 126 : for(i = 0; i < 20; i++ )
1566 : snprintf(sVal+strlen(sVal), sizeof(sVal), "%.16g ",
1567 120 : atof(papszFields[50+i]));
1568 6 : SetMetadataItem("SAMP_NUM_COEFF",sVal,"RPC");
1569 :
1570 6 : sVal[0] = '\0';
1571 126 : for(i = 0; i < 20; i++ )
1572 : snprintf(sVal+strlen(sVal), sizeof(sVal), "%.16g ",
1573 120 : atof(papszFields[70+i]));
1574 6 : SetMetadataItem("SAMP_DEN_COEFF",sVal,"RPC");
1575 :
1576 : snprintf(sVal, sizeof(sVal), "%.16g",
1577 6 : atof(papszFields[3]) - atof(papszFields[8]));
1578 6 : SetMetadataItem("MIN_LONG",sVal,"RPC");
1579 :
1580 : snprintf(sVal, sizeof(sVal), "%.16g",
1581 6 : atof(papszFields[3]) + atof(papszFields[8]) );
1582 6 : SetMetadataItem("MAX_LONG",sVal,"RPC");
1583 :
1584 : snprintf(sVal, sizeof(sVal), "%.16g",
1585 6 : atof(papszFields[2]) - atof(papszFields[7]));
1586 6 : SetMetadataItem("MIN_LAT",sVal,"RPC");
1587 :
1588 : snprintf(sVal, sizeof(sVal), "%.16g",
1589 6 : atof(papszFields[2]) + atof(papszFields[7]));
1590 6 : SetMetadataItem("MAX_LAT",sVal,"RPC");
1591 :
1592 6 : if (nCount == 93)
1593 : {
1594 6 : SetMetadataItem("TILE_ROW_OFFSET",papszFields[90],"RPC");
1595 6 : SetMetadataItem("TILE_COL_OFFSET",papszFields[91],"RPC");
1596 6 : SetMetadataItem("ENVI_RPC_EMULATION",papszFields[92],"RPC");
1597 : }
1598 :
1599 : /* Handle the chipping case where the image is a subset. */
1600 : double rowOffset, colOffset;
1601 6 : rowOffset = (nCount == 93) ? atof(papszFields[90]) : 0;
1602 6 : colOffset = (nCount == 93) ? atof(papszFields[91]) : 0;
1603 6 : if (rowOffset || colOffset)
1604 : {
1605 0 : SetMetadataItem("ICHIP_SCALE_FACTOR", "1");
1606 0 : SetMetadataItem("ICHIP_ANAMORPH_CORR", "0");
1607 0 : SetMetadataItem("ICHIP_SCANBLK_NUM", "0");
1608 :
1609 0 : SetMetadataItem("ICHIP_OP_ROW_11", "0.5");
1610 0 : SetMetadataItem("ICHIP_OP_COL_11", "0.5");
1611 0 : SetMetadataItem("ICHIP_OP_ROW_12", "0.5");
1612 0 : SetMetadataItem("ICHIP_OP_COL_21", "0.5");
1613 0 : snprintf(sVal, sizeof(sVal), "%.16g", numCols - 0.5);
1614 0 : SetMetadataItem("ICHIP_OP_COL_12", sVal);
1615 0 : SetMetadataItem("ICHIP_OP_COL_22", sVal);
1616 0 : snprintf(sVal, sizeof(sVal), "%.16g", numRows - 0.5);
1617 0 : SetMetadataItem("ICHIP_OP_ROW_21", sVal);
1618 0 : SetMetadataItem("ICHIP_OP_ROW_22", sVal);
1619 :
1620 0 : snprintf(sVal, sizeof(sVal), "%.16g", rowOffset + 0.5);
1621 0 : SetMetadataItem("ICHIP_FI_ROW_11", sVal);
1622 0 : SetMetadataItem("ICHIP_FI_ROW_12", sVal);
1623 0 : snprintf(sVal, sizeof(sVal), "%.16g", colOffset + 0.5);
1624 0 : SetMetadataItem("ICHIP_FI_COL_11", sVal);
1625 0 : SetMetadataItem("ICHIP_FI_COL_21", sVal);
1626 0 : snprintf(sVal, sizeof(sVal), "%.16g", colOffset + numCols - 0.5);
1627 0 : SetMetadataItem("ICHIP_FI_COL_12", sVal);
1628 0 : SetMetadataItem("ICHIP_FI_COL_22", sVal);
1629 0 : snprintf(sVal, sizeof(sVal), "%.16g", rowOffset + numRows - 0.5);
1630 0 : SetMetadataItem("ICHIP_FI_ROW_21", sVal);
1631 0 : SetMetadataItem("ICHIP_FI_ROW_22", sVal);
1632 : }
1633 6 : CSLDestroy(papszFields);
1634 : }
1635 :
1636 222 : void ENVIDataset::ProcessStatsFile()
1637 : {
1638 : VSILFILE *fpStaFile;
1639 :
1640 222 : osStaFilename = CPLResetExtension( pszHDRFilename, "sta" );
1641 222 : fpStaFile = VSIFOpenL( osStaFilename, "rb" );
1642 :
1643 222 : if (!fpStaFile)
1644 : {
1645 220 : osStaFilename = "";
1646 220 : return;
1647 : }
1648 :
1649 : int lTestHeader[10],lOffset;
1650 :
1651 2 : if( VSIFReadL( lTestHeader, sizeof(int), 10, fpStaFile ) != 10 )
1652 : {
1653 0 : VSIFCloseL( fpStaFile );
1654 0 : osStaFilename = "";
1655 0 : return;
1656 : }
1657 :
1658 : int isFloat;
1659 2 : isFloat = (byteSwapInt(lTestHeader[0]) == 1111838282);
1660 :
1661 : int nb,i;
1662 : float * fStats;
1663 : double * dStats, dMin, dMax, dMean, dStd;
1664 :
1665 2 : nb=byteSwapInt(lTestHeader[3]);
1666 :
1667 2 : if (nb < 0 || nb > nBands)
1668 : {
1669 : CPLDebug("ENVI", ".sta file has statistics for %d bands, "
1670 0 : "whereas the dataset has only %d bands", nb, nBands);
1671 0 : nb = nBands;
1672 : }
1673 :
1674 2 : VSIFSeekL(fpStaFile,40+(nb+1)*4,SEEK_SET);
1675 :
1676 2 : if (VSIFReadL(&lOffset,sizeof(int),1,fpStaFile) == 1)
1677 : {
1678 2 : VSIFSeekL(fpStaFile,40+(nb+1)*8+byteSwapInt(lOffset)+nb,SEEK_SET);
1679 : // This should be the beginning of the statistics
1680 2 : if (isFloat)
1681 : {
1682 0 : fStats = (float*)CPLCalloc(nb*4,4);
1683 0 : if ((int)VSIFReadL(fStats,4,nb*4,fpStaFile) == nb*4)
1684 : {
1685 0 : for (i=0;i<nb;i++)
1686 : {
1687 : GetRasterBand(i+1)->SetStatistics(
1688 : byteSwapFloat(fStats[i]),
1689 : byteSwapFloat(fStats[nb+i]),
1690 : byteSwapFloat(fStats[2*nb+i]),
1691 0 : byteSwapFloat(fStats[3*nb+i]));
1692 : }
1693 : }
1694 0 : CPLFree(fStats);
1695 : }
1696 : else
1697 : {
1698 2 : dStats = (double*)CPLCalloc(nb*4,8);
1699 2 : if ((int)VSIFReadL(dStats,8,nb*4,fpStaFile) == nb*4)
1700 : {
1701 14 : for (i=0;i<nb;i++)
1702 : {
1703 12 : dMin = byteSwapDouble(dStats[i]);
1704 12 : dMax = byteSwapDouble(dStats[nb+i]);
1705 12 : dMean = byteSwapDouble(dStats[2*nb+i]);
1706 12 : dStd = byteSwapDouble(dStats[3*nb+i]);
1707 12 : if (dMin != dMax && dStd != 0)
1708 12 : GetRasterBand(i+1)->SetStatistics(dMin,dMax,dMean,dStd);
1709 : }
1710 : }
1711 2 : CPLFree(dStats);
1712 : }
1713 : }
1714 2 : VSIFCloseL( fpStaFile );
1715 : }
1716 :
1717 6 : int ENVIDataset::byteSwapInt(int swapMe)
1718 : {
1719 6 : CPL_MSBPTR32(&swapMe);
1720 6 : return swapMe;
1721 : }
1722 :
1723 0 : float ENVIDataset::byteSwapFloat(float swapMe)
1724 : {
1725 0 : CPL_MSBPTR32(&swapMe);
1726 0 : return swapMe;
1727 : }
1728 :
1729 48 : double ENVIDataset::byteSwapDouble(double swapMe)
1730 : {
1731 48 : CPL_MSBPTR64(&swapMe);
1732 48 : return swapMe;
1733 : }
1734 :
1735 :
1736 : /************************************************************************/
1737 : /* ReadHeader() */
1738 : /************************************************************************/
1739 :
1740 226 : int ENVIDataset::ReadHeader( VSILFILE * fpHdr )
1741 :
1742 : {
1743 :
1744 226 : CPLReadLineL( fpHdr );
1745 :
1746 : /* -------------------------------------------------------------------- */
1747 : /* Now start forming sets of name/value pairs. */
1748 : /* -------------------------------------------------------------------- */
1749 2282 : while( TRUE )
1750 : {
1751 : const char *pszNewLine;
1752 : char *pszWorkingLine;
1753 :
1754 2508 : pszNewLine = CPLReadLineL( fpHdr );
1755 2508 : if( pszNewLine == NULL )
1756 : break;
1757 :
1758 2282 : if( strstr(pszNewLine,"=") == NULL )
1759 0 : continue;
1760 :
1761 2282 : pszWorkingLine = CPLStrdup(pszNewLine);
1762 :
1763 : // Collect additional lines if we have open sqiggly bracket.
1764 2282 : if( strstr(pszWorkingLine,"{") != NULL
1765 : && strstr(pszWorkingLine,"}") == NULL )
1766 : {
1767 434 : do {
1768 434 : pszNewLine = CPLReadLineL( fpHdr );
1769 434 : if( pszNewLine )
1770 : {
1771 : pszWorkingLine = (char *)
1772 : CPLRealloc(pszWorkingLine,
1773 434 : strlen(pszWorkingLine)+strlen(pszNewLine)+1);
1774 434 : strcat( pszWorkingLine, pszNewLine );
1775 : }
1776 : } while( pszNewLine != NULL && strstr(pszNewLine,"}") == NULL );
1777 : }
1778 :
1779 : // Try to break input into name and value portions. Trim whitespace.
1780 : const char *pszValue;
1781 : int iEqual;
1782 :
1783 55702 : for( iEqual = 0;
1784 53420 : pszWorkingLine[iEqual] != '\0' && pszWorkingLine[iEqual] != '=';
1785 : iEqual++ ) {}
1786 :
1787 2282 : if( pszWorkingLine[iEqual] == '=' )
1788 : {
1789 : int i;
1790 :
1791 2282 : pszValue = pszWorkingLine + iEqual + 1;
1792 6846 : while( *pszValue == ' ' || *pszValue == '\t' )
1793 2282 : pszValue++;
1794 :
1795 2282 : pszWorkingLine[iEqual--] = '\0';
1796 15500 : while( iEqual > 0
1797 5468 : && (pszWorkingLine[iEqual] == ' '
1798 2282 : || pszWorkingLine[iEqual] == '\t') )
1799 3186 : pszWorkingLine[iEqual--] = '\0';
1800 :
1801 : // Convert spaces in the name to underscores.
1802 23524 : for( i = 0; pszWorkingLine[i] != '\0'; i++ )
1803 : {
1804 21242 : if( pszWorkingLine[i] == ' ' )
1805 1344 : pszWorkingLine[i] = '_';
1806 : }
1807 :
1808 : papszHeader = CSLSetNameValue( papszHeader,
1809 2282 : pszWorkingLine, pszValue );
1810 : }
1811 :
1812 2282 : CPLFree( pszWorkingLine );
1813 : }
1814 :
1815 226 : return TRUE;
1816 : }
1817 :
1818 : /************************************************************************/
1819 : /* Open() */
1820 : /************************************************************************/
1821 :
1822 23256 : GDALDataset *ENVIDataset::Open( GDALOpenInfo * poOpenInfo )
1823 :
1824 : {
1825 : int i;
1826 :
1827 : /* -------------------------------------------------------------------- */
1828 : /* We assume the user is pointing to the binary (ie. .bil) file. */
1829 : /* -------------------------------------------------------------------- */
1830 23256 : if( poOpenInfo->nHeaderBytes < 2 )
1831 21302 : return NULL;
1832 :
1833 : /* -------------------------------------------------------------------- */
1834 : /* Do we have a .hdr file? Try upper and lower case, and */
1835 : /* replacing the extension as well as appending the extension */
1836 : /* to whatever we currently have. */
1837 : /* -------------------------------------------------------------------- */
1838 : const char *pszMode;
1839 1954 : CPLString osHdrFilename;
1840 1954 : VSILFILE *fpHeader = NULL;
1841 :
1842 1954 : if( poOpenInfo->eAccess == GA_Update )
1843 472 : pszMode = "r+";
1844 : else
1845 1482 : pszMode = "r";
1846 :
1847 1954 : if (poOpenInfo->papszSiblingFiles == NULL)
1848 : {
1849 8 : osHdrFilename = CPLResetExtension( poOpenInfo->pszFilename, "hdr" );
1850 8 : fpHeader = VSIFOpenL( osHdrFilename, pszMode );
1851 :
1852 8 : if( fpHeader == NULL && VSIIsCaseSensitiveFS(osHdrFilename) )
1853 : {
1854 8 : osHdrFilename = CPLResetExtension( poOpenInfo->pszFilename, "HDR" );
1855 8 : fpHeader = VSIFOpenL( osHdrFilename, pszMode );
1856 : }
1857 :
1858 8 : if( fpHeader == NULL && VSIIsCaseSensitiveFS(osHdrFilename) )
1859 : {
1860 : osHdrFilename = CPLFormFilename( NULL, poOpenInfo->pszFilename,
1861 8 : "hdr" );
1862 8 : fpHeader = VSIFOpenL( osHdrFilename, pszMode );
1863 : }
1864 :
1865 8 : if( fpHeader == NULL && VSIIsCaseSensitiveFS(osHdrFilename) )
1866 : {
1867 : osHdrFilename = CPLFormFilename( NULL, poOpenInfo->pszFilename,
1868 8 : "HDR" );
1869 8 : fpHeader = VSIFOpenL( osHdrFilename, pszMode );
1870 : }
1871 :
1872 : }
1873 : else
1874 : {
1875 : /* -------------------------------------------------------------------- */
1876 : /* Now we need to tear apart the filename to form a .HDR */
1877 : /* filename. */
1878 : /* -------------------------------------------------------------------- */
1879 1946 : CPLString osPath = CPLGetPath( poOpenInfo->pszFilename );
1880 1946 : CPLString osName = CPLGetFilename( poOpenInfo->pszFilename );
1881 :
1882 : int iFile = CSLFindString(poOpenInfo->papszSiblingFiles,
1883 1946 : CPLResetExtension( osName, "hdr" ) );
1884 1946 : if( iFile >= 0 )
1885 : {
1886 440 : osHdrFilename = CPLFormFilename( osPath, poOpenInfo->papszSiblingFiles[iFile],
1887 880 : NULL );
1888 440 : fpHeader = VSIFOpenL( osHdrFilename, pszMode );
1889 : }
1890 : else
1891 : {
1892 : iFile = CSLFindString(poOpenInfo->papszSiblingFiles,
1893 1506 : CPLFormFilename( NULL, osName, "hdr" ));
1894 1506 : if( iFile >= 0 )
1895 : {
1896 0 : osHdrFilename = CPLFormFilename( osPath, poOpenInfo->papszSiblingFiles[iFile],
1897 0 : NULL );
1898 0 : fpHeader = VSIFOpenL( osHdrFilename, pszMode );
1899 : }
1900 1946 : }
1901 : }
1902 :
1903 1954 : if( fpHeader == NULL )
1904 1514 : return NULL;
1905 :
1906 : /* -------------------------------------------------------------------- */
1907 : /* Check that the first line says "ENVI". */
1908 : /* -------------------------------------------------------------------- */
1909 : char szTestHdr[4];
1910 :
1911 440 : if( VSIFReadL( szTestHdr, 4, 1, fpHeader ) != 1 )
1912 : {
1913 0 : VSIFCloseL( fpHeader );
1914 0 : return NULL;
1915 : }
1916 440 : if( strncmp(szTestHdr,"ENVI",4) != 0 )
1917 : {
1918 214 : VSIFCloseL( fpHeader );
1919 214 : return NULL;
1920 : }
1921 :
1922 : /* -------------------------------------------------------------------- */
1923 : /* Create a corresponding GDALDataset. */
1924 : /* -------------------------------------------------------------------- */
1925 : ENVIDataset *poDS;
1926 :
1927 226 : poDS = new ENVIDataset();
1928 452 : poDS->pszHDRFilename = CPLStrdup(osHdrFilename);
1929 226 : poDS->fp = fpHeader;
1930 :
1931 : /* -------------------------------------------------------------------- */
1932 : /* Read the header. */
1933 : /* -------------------------------------------------------------------- */
1934 226 : if( !poDS->ReadHeader( fpHeader ) )
1935 : {
1936 0 : delete poDS;
1937 0 : return NULL;
1938 : }
1939 :
1940 : /* -------------------------------------------------------------------- */
1941 : /* Has the user selected the .hdr file to open? */
1942 : /* -------------------------------------------------------------------- */
1943 226 : if( EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "hdr") )
1944 : {
1945 0 : delete poDS;
1946 : CPLError( CE_Failure, CPLE_AppDefined,
1947 : "The selected file is an ENVI header file, but to\n"
1948 : "open ENVI datasets, the data file should be selected\n"
1949 : "instead of the .hdr file. Please try again selecting\n"
1950 : "the data file corresponding to the header file:\n"
1951 : " %s\n",
1952 0 : poOpenInfo->pszFilename );
1953 0 : return NULL;
1954 : }
1955 :
1956 : /* -------------------------------------------------------------------- */
1957 : /* Has the user selected the .sta (stats) file to open? */
1958 : /* -------------------------------------------------------------------- */
1959 226 : if( EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "sta") )
1960 : {
1961 0 : delete poDS;
1962 : CPLError( CE_Failure, CPLE_AppDefined,
1963 : "The selected file is an ENVI statistics file.\n"
1964 : "To open ENVI datasets, the data file should be selected\n"
1965 : "instead of the .sta file. Please try again selecting\n"
1966 : "the data file corresponding to the statistics file:\n"
1967 : " %s\n",
1968 0 : poOpenInfo->pszFilename );
1969 0 : return NULL;
1970 : }
1971 :
1972 : /* -------------------------------------------------------------------- */
1973 : /* Extract required values from the .hdr. */
1974 : /* -------------------------------------------------------------------- */
1975 226 : int nLines = 0, nSamples = 0, nBands = 0, nHeaderSize = 0;
1976 226 : const char *pszInterleave = NULL;
1977 :
1978 226 : if( CSLFetchNameValue(poDS->papszHeader,"lines") )
1979 226 : nLines = atoi(CSLFetchNameValue(poDS->papszHeader,"lines"));
1980 :
1981 226 : if( CSLFetchNameValue(poDS->papszHeader,"samples") )
1982 226 : nSamples = atoi(CSLFetchNameValue(poDS->papszHeader,"samples"));
1983 :
1984 226 : if( CSLFetchNameValue(poDS->papszHeader,"bands") )
1985 226 : nBands = atoi(CSLFetchNameValue(poDS->papszHeader,"bands"));
1986 :
1987 226 : pszInterleave = CSLFetchNameValue(poDS->papszHeader,"interleave");
1988 :
1989 : /* In case, there is no interleave keyword, we try to derive it from the */
1990 : /* file extension. */
1991 226 : if( pszInterleave == NULL )
1992 : {
1993 0 : const char* pszExtension = CPLGetExtension(poOpenInfo->pszFilename);
1994 0 : pszInterleave = pszExtension;
1995 : }
1996 :
1997 226 : if ( !EQUALN(pszInterleave, "BSQ",3) &&
1998 : !EQUALN(pszInterleave, "BIP",3) &&
1999 : !EQUALN(pszInterleave, "BIL",3) )
2000 : {
2001 0 : CPLDebug("ENVI", "Unset or unknown value for 'interleave' keyword --> assuming BSQ interleaving");
2002 0 : pszInterleave = "bsq";
2003 : }
2004 :
2005 226 : if (!GDALCheckDatasetDimensions(nSamples, nLines) || !GDALCheckBandCount(nBands, FALSE))
2006 : {
2007 4 : delete poDS;
2008 : CPLError( CE_Failure, CPLE_AppDefined,
2009 : "The file appears to have an associated ENVI header, but\n"
2010 : "one or more of the samples, lines and bands\n"
2011 4 : "keywords appears to be missing or invalid." );
2012 4 : return NULL;
2013 : }
2014 :
2015 222 : if( CSLFetchNameValue(poDS->papszHeader,"header_offset") )
2016 222 : nHeaderSize = atoi(CSLFetchNameValue(poDS->papszHeader,"header_offset"));
2017 :
2018 : /* -------------------------------------------------------------------- */
2019 : /* Translate the datatype. */
2020 : /* -------------------------------------------------------------------- */
2021 222 : GDALDataType eType = GDT_Byte;
2022 :
2023 222 : if( CSLFetchNameValue(poDS->papszHeader,"data_type" ) != NULL )
2024 : {
2025 222 : switch( atoi(CSLFetchNameValue(poDS->papszHeader,"data_type" )) )
2026 : {
2027 : case 1:
2028 130 : eType = GDT_Byte;
2029 130 : break;
2030 :
2031 : case 2:
2032 12 : eType = GDT_Int16;
2033 12 : break;
2034 :
2035 : case 3:
2036 12 : eType = GDT_Int32;
2037 12 : break;
2038 :
2039 : case 4:
2040 12 : eType = GDT_Float32;
2041 12 : break;
2042 :
2043 : case 5:
2044 12 : eType = GDT_Float64;
2045 12 : break;
2046 :
2047 : case 6:
2048 10 : eType = GDT_CFloat32;
2049 10 : break;
2050 :
2051 : case 9:
2052 10 : eType = GDT_CFloat64;
2053 10 : break;
2054 :
2055 : case 12:
2056 12 : eType = GDT_UInt16;
2057 12 : break;
2058 :
2059 : case 13:
2060 12 : eType = GDT_UInt32;
2061 12 : break;
2062 :
2063 : /* 14=Int64, 15=UInt64 */
2064 :
2065 : default:
2066 0 : delete poDS;
2067 : CPLError( CE_Failure, CPLE_AppDefined,
2068 : "The file does not have a value for the data_type\n"
2069 0 : "that is recognised by the GDAL ENVI driver.");
2070 0 : return NULL;
2071 : }
2072 : }
2073 :
2074 : /* -------------------------------------------------------------------- */
2075 : /* Translate the byte order. */
2076 : /* -------------------------------------------------------------------- */
2077 222 : int bNativeOrder = TRUE;
2078 :
2079 222 : if( CSLFetchNameValue(poDS->papszHeader,"byte_order" ) != NULL )
2080 : {
2081 : #ifdef CPL_LSB
2082 : bNativeOrder = atoi(CSLFetchNameValue(poDS->papszHeader,
2083 222 : "byte_order" )) == 0;
2084 : #else
2085 : bNativeOrder = atoi(CSLFetchNameValue(poDS->papszHeader,
2086 : "byte_order" )) != 0;
2087 : #endif
2088 : }
2089 :
2090 : /* -------------------------------------------------------------------- */
2091 : /* Warn about unsupported file types virtual mosaic and meta file.*/
2092 : /* -------------------------------------------------------------------- */
2093 222 : if( CSLFetchNameValue(poDS->papszHeader,"file_type" ) != NULL )
2094 : {
2095 : // when the file type is one of these we return an invalid file type err
2096 : //'envi meta file'
2097 : //'envi virtual mosaic'
2098 : //'envi spectral library'
2099 : //'envi fft result'
2100 :
2101 : // when the file type is one of these we open it
2102 : //'envi standard'
2103 : //'envi classification'
2104 :
2105 : // when the file type is anything else we attempt to open it as a raster.
2106 :
2107 : const char * pszEnviFileType;
2108 222 : pszEnviFileType = CSLFetchNameValue(poDS->papszHeader,"file_type");
2109 :
2110 : // envi gdal does not support any of these
2111 : // all others we will attempt to open
2112 222 : if(EQUAL(pszEnviFileType, "envi meta file") ||
2113 : EQUAL(pszEnviFileType, "envi virtual mosaic") ||
2114 : EQUAL(pszEnviFileType, "envi spectral library") ||
2115 : EQUAL(pszEnviFileType, "envi fft result"))
2116 : {
2117 : CPLError( CE_Failure, CPLE_OpenFailed,
2118 : "File %s contains an invalid file type in the ENVI .hdr\n"
2119 : "GDAL does not support '%s' type files.",
2120 0 : poOpenInfo->pszFilename, pszEnviFileType );
2121 0 : delete poDS;
2122 0 : return NULL;
2123 : }
2124 : }
2125 :
2126 : /* -------------------------------------------------------------------- */
2127 : /* Detect (gzipped) compressed datasets. */
2128 : /* -------------------------------------------------------------------- */
2129 222 : int bIsCompressed = FALSE;
2130 222 : if( CSLFetchNameValue(poDS->papszHeader,"file_compression" ) != NULL )
2131 : {
2132 2 : if( atoi(CSLFetchNameValue(poDS->papszHeader,"file_compression" ))
2133 : != 0 )
2134 : {
2135 2 : bIsCompressed = TRUE;
2136 : }
2137 : }
2138 :
2139 : /* -------------------------------------------------------------------- */
2140 : /* Capture some information from the file that is of interest. */
2141 : /* -------------------------------------------------------------------- */
2142 222 : poDS->nRasterXSize = nSamples;
2143 222 : poDS->nRasterYSize = nLines;
2144 222 : poDS->eAccess = poOpenInfo->eAccess;
2145 :
2146 : /* -------------------------------------------------------------------- */
2147 : /* Reopen file in update mode if necessary. */
2148 : /* -------------------------------------------------------------------- */
2149 222 : CPLString osImageFilename(poOpenInfo->pszFilename);
2150 222 : if (bIsCompressed)
2151 2 : osImageFilename = "/vsigzip/" + osImageFilename;
2152 222 : if( poOpenInfo->eAccess == GA_Update )
2153 : {
2154 100 : if (bIsCompressed)
2155 : {
2156 0 : delete poDS;
2157 : CPLError( CE_Failure, CPLE_OpenFailed,
2158 0 : "Cannot open compressed file in update mode.\n");
2159 0 : return NULL;
2160 : }
2161 100 : poDS->fpImage = VSIFOpenL( osImageFilename, "rb+" );
2162 : }
2163 : else
2164 122 : poDS->fpImage = VSIFOpenL( osImageFilename, "rb" );
2165 :
2166 222 : if( poDS->fpImage == NULL )
2167 : {
2168 0 : delete poDS;
2169 : CPLError( CE_Failure, CPLE_OpenFailed,
2170 : "Failed to re-open %s within ENVI driver.\n",
2171 0 : poOpenInfo->pszFilename );
2172 0 : return NULL;
2173 : }
2174 :
2175 : /* -------------------------------------------------------------------- */
2176 : /* Compute the line offset. */
2177 : /* -------------------------------------------------------------------- */
2178 222 : int nDataSize = GDALGetDataTypeSize(eType)/8;
2179 : int nPixelOffset, nLineOffset;
2180 : vsi_l_offset nBandOffset;
2181 222 : int bIntOverflow = FALSE;
2182 :
2183 222 : if( EQUALN(pszInterleave, "bil", 3) )
2184 : {
2185 0 : poDS->interleave = BIL;
2186 0 : poDS->SetMetadataItem( "INTERLEAVE", "LINE", "IMAGE_STRUCTURE" );
2187 0 : if (nSamples > INT_MAX / (nDataSize * nBands)) bIntOverflow = TRUE;
2188 0 : nLineOffset = nDataSize * nSamples * nBands;
2189 0 : nPixelOffset = nDataSize;
2190 0 : nBandOffset = (vsi_l_offset)nDataSize * nSamples;
2191 : }
2192 222 : else if( EQUALN(pszInterleave, "bip", 3) )
2193 : {
2194 0 : poDS->interleave = BIP;
2195 0 : poDS->SetMetadataItem( "INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE" );
2196 0 : if (nSamples > INT_MAX / (nDataSize * nBands)) bIntOverflow = TRUE;
2197 0 : nLineOffset = nDataSize * nSamples * nBands;
2198 0 : nPixelOffset = nDataSize * nBands;
2199 0 : nBandOffset = nDataSize;
2200 : }
2201 : else /* bsq */
2202 : {
2203 222 : poDS->interleave = BSQ;
2204 222 : poDS->SetMetadataItem( "INTERLEAVE", "BAND", "IMAGE_STRUCTURE" );
2205 222 : if (nSamples > INT_MAX / nDataSize) bIntOverflow = TRUE;
2206 222 : nLineOffset = nDataSize * nSamples;
2207 222 : nPixelOffset = nDataSize;
2208 222 : nBandOffset = (vsi_l_offset)nLineOffset * nLines;
2209 : }
2210 :
2211 222 : if (bIntOverflow)
2212 : {
2213 0 : delete poDS;
2214 : CPLError( CE_Failure, CPLE_AppDefined,
2215 0 : "Int overflow occured.");
2216 0 : return NULL;
2217 : }
2218 :
2219 : /* -------------------------------------------------------------------- */
2220 : /* Create band information objects. */
2221 : /* -------------------------------------------------------------------- */
2222 222 : poDS->nBands = nBands;
2223 594 : for( i = 0; i < poDS->nBands; i++ )
2224 : {
2225 : poDS->SetBand( i + 1,
2226 : new RawRasterBand(poDS, i + 1, poDS->fpImage,
2227 : nHeaderSize + nBandOffset * i,
2228 : nPixelOffset, nLineOffset, eType,
2229 372 : bNativeOrder, TRUE) );
2230 : }
2231 :
2232 : /* -------------------------------------------------------------------- */
2233 : /* Apply band names if we have them. */
2234 : /* Use wavelength for more descriptive information if possible */
2235 : /* -------------------------------------------------------------------- */
2236 222 : if( CSLFetchNameValue( poDS->papszHeader, "band_names" ) != NULL ||
2237 : CSLFetchNameValue( poDS->papszHeader, "wavelength" ) != NULL)
2238 : {
2239 : char **papszBandNames =
2240 : poDS->SplitList( CSLFetchNameValue( poDS->papszHeader,
2241 116 : "band_names" ) );
2242 : char **papszWL =
2243 : poDS->SplitList( CSLFetchNameValue( poDS->papszHeader,
2244 116 : "wavelength" ) );
2245 :
2246 116 : const char *pszWLUnits = NULL;
2247 116 : if (papszWL)
2248 : {
2249 : /* If WL information is present, process wavelength units */
2250 : pszWLUnits = CSLFetchNameValue( poDS->papszHeader,
2251 0 : "wavelength_units" );
2252 0 : if (pszWLUnits)
2253 : {
2254 : /* Don't show unknown or index units */
2255 0 : if (EQUAL(pszWLUnits,"Unknown") ||
2256 : EQUAL(pszWLUnits,"Index") )
2257 0 : pszWLUnits=0;
2258 : }
2259 : }
2260 :
2261 292 : for( i = 0; i < nBands; i++ )
2262 : {
2263 176 : CPLString osBandId, osBandName, osWavelength;
2264 :
2265 : /* First set up the wavelength names and units if available */
2266 176 : if (papszWL && CSLCount(papszWL) > i)
2267 : {
2268 0 : osWavelength = papszWL[i];
2269 0 : if (pszWLUnits)
2270 : {
2271 0 : osWavelength += " ";
2272 0 : osWavelength += pszWLUnits;
2273 : }
2274 : }
2275 :
2276 : /* Build the final name for this band */
2277 176 : if (papszBandNames && CSLCount(papszBandNames) > i)
2278 : {
2279 176 : osBandName = papszBandNames[i];
2280 176 : if (strlen(osWavelength) > 0)
2281 : {
2282 0 : osBandName += " (";
2283 0 : osBandName += osWavelength;
2284 0 : osBandName += ")";
2285 : }
2286 : }
2287 : else /* WL but no band names */
2288 0 : osBandName = osWavelength;
2289 :
2290 : /* Description is for internal GDAL usage */
2291 176 : poDS->GetRasterBand(i + 1)->SetDescription( osBandName );
2292 :
2293 : /* Metadata field named Band_1, etc. needed for ArcGIS integration */
2294 176 : osBandId = CPLSPrintf("Band_%i", i+1);
2295 176 : poDS->SetMetadataItem(osBandId, osBandName);
2296 : }
2297 116 : CSLDestroy( papszWL );
2298 116 : CSLDestroy( papszBandNames );
2299 : }
2300 : /* -------------------------------------------------------------------- */
2301 : /* Apply class names if we have them. */
2302 : /* -------------------------------------------------------------------- */
2303 222 : if( CSLFetchNameValue( poDS->papszHeader, "class_names" ) != NULL )
2304 : {
2305 : char **papszClassNames =
2306 : poDS->SplitList( CSLFetchNameValue( poDS->papszHeader,
2307 6 : "class_names" ) );
2308 :
2309 6 : poDS->GetRasterBand(1)->SetCategoryNames( papszClassNames );
2310 6 : CSLDestroy( papszClassNames );
2311 : }
2312 :
2313 : /* -------------------------------------------------------------------- */
2314 : /* Apply colormap if we have one. */
2315 : /* -------------------------------------------------------------------- */
2316 222 : if( CSLFetchNameValue( poDS->papszHeader, "class_lookup" ) != NULL )
2317 : {
2318 : char **papszClassColors =
2319 : poDS->SplitList( CSLFetchNameValue( poDS->papszHeader,
2320 6 : "class_lookup" ) );
2321 6 : int nColorValueCount = CSLCount(papszClassColors);
2322 6 : GDALColorTable oCT;
2323 :
2324 18 : for( i = 0; i*3+2 < nColorValueCount; i++ )
2325 : {
2326 : GDALColorEntry sEntry;
2327 :
2328 12 : sEntry.c1 = (short) atoi(papszClassColors[i*3+0]);
2329 12 : sEntry.c2 = (short) atoi(papszClassColors[i*3+1]);
2330 12 : sEntry.c3 = (short) atoi(papszClassColors[i*3+2]);
2331 12 : sEntry.c4 = 255;
2332 12 : oCT.SetColorEntry( i, &sEntry );
2333 : }
2334 :
2335 6 : CSLDestroy( papszClassColors );
2336 :
2337 6 : poDS->GetRasterBand(1)->SetColorTable( &oCT );
2338 6 : poDS->GetRasterBand(1)->SetColorInterpretation( GCI_PaletteIndex );
2339 : }
2340 :
2341 : /* -------------------------------------------------------------------- */
2342 : /* Set the nodata value if it is present */
2343 : /* -------------------------------------------------------------------- */
2344 222 : if( CSLFetchNameValue(poDS->papszHeader,"data_ignore_value" ) != NULL )
2345 : {
2346 0 : for( i = 0; i < poDS->nBands; i++ )
2347 : ((RawRasterBand*)poDS->GetRasterBand(i+1))->SetNoDataValue(atof(
2348 0 : CSLFetchNameValue(poDS->papszHeader,"data_ignore_value")));
2349 : }
2350 :
2351 : /* -------------------------------------------------------------------- */
2352 : /* Read the stats file if it is present */
2353 : /* -------------------------------------------------------------------- */
2354 222 : poDS->ProcessStatsFile();
2355 :
2356 : /* -------------------------------------------------------------------- */
2357 : /* Look for mapinfo */
2358 : /* -------------------------------------------------------------------- */
2359 222 : if( CSLFetchNameValue( poDS->papszHeader, "map_info" ) != NULL )
2360 : {
2361 : poDS->bFoundMapinfo =
2362 : poDS->ProcessMapinfo(
2363 104 : CSLFetchNameValue(poDS->papszHeader,"map_info") );
2364 : }
2365 :
2366 : /* -------------------------------------------------------------------- */
2367 : /* Look for RPC mapinfo */
2368 : /* -------------------------------------------------------------------- */
2369 222 : if( !poDS->bFoundMapinfo &&
2370 : CSLFetchNameValue( poDS->papszHeader, "rpc_info" ) != NULL )
2371 : {
2372 : poDS->ProcessRPCinfo(
2373 : CSLFetchNameValue(poDS->papszHeader,"rpc_info"),
2374 6 : poDS->nRasterXSize, poDS->nRasterYSize);
2375 : }
2376 :
2377 : /* -------------------------------------------------------------------- */
2378 : /* Initialize any PAM information. */
2379 : /* -------------------------------------------------------------------- */
2380 222 : poDS->SetDescription( poOpenInfo->pszFilename );
2381 222 : poDS->TryLoadXML();
2382 :
2383 : /* -------------------------------------------------------------------- */
2384 : /* Check for overviews. */
2385 : /* -------------------------------------------------------------------- */
2386 222 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
2387 :
2388 222 : return( poDS );
2389 : }
2390 :
2391 216 : int ENVIDataset::GetEnviType(GDALDataType eType)
2392 : {
2393 : int iENVIType;
2394 216 : switch( eType )
2395 : {
2396 : case GDT_Byte:
2397 108 : iENVIType = 1;
2398 108 : break;
2399 : case GDT_Int16:
2400 12 : iENVIType = 2;
2401 12 : break;
2402 : case GDT_Int32:
2403 12 : iENVIType = 3;
2404 12 : break;
2405 : case GDT_Float32:
2406 12 : iENVIType = 4;
2407 12 : break;
2408 : case GDT_Float64:
2409 12 : iENVIType = 5;
2410 12 : break;
2411 : case GDT_CFloat32:
2412 12 : iENVIType = 6;
2413 12 : break;
2414 : case GDT_CFloat64:
2415 12 : iENVIType = 9;
2416 12 : break;
2417 : case GDT_UInt16:
2418 12 : iENVIType = 12;
2419 12 : break;
2420 : case GDT_UInt32:
2421 12 : iENVIType = 13;
2422 12 : break;
2423 :
2424 : /* 14=Int64, 15=UInt64 */
2425 :
2426 : default:
2427 : CPLError( CE_Failure, CPLE_AppDefined,
2428 : "Attempt to create ENVI .hdr labelled dataset with an illegal\n"
2429 : "data type (%s).\n",
2430 12 : GDALGetDataTypeName(eType) );
2431 12 : return 1;
2432 : }
2433 204 : return iENVIType;
2434 : }
2435 :
2436 : /************************************************************************/
2437 : /* Create() */
2438 : /************************************************************************/
2439 :
2440 106 : GDALDataset *ENVIDataset::Create( const char * pszFilename,
2441 : int nXSize, int nYSize, int nBands,
2442 : GDALDataType eType,
2443 : char ** papszOptions )
2444 :
2445 : {
2446 : /* -------------------------------------------------------------------- */
2447 : /* Verify input options. */
2448 : /* -------------------------------------------------------------------- */
2449 106 : int iENVIType = GetEnviType(eType);
2450 106 : if (0 == iENVIType)
2451 0 : return 0;
2452 :
2453 : /* -------------------------------------------------------------------- */
2454 : /* Try to create the file. */
2455 : /* -------------------------------------------------------------------- */
2456 : VSILFILE *fp;
2457 :
2458 106 : fp = VSIFOpenL( pszFilename, "wb" );
2459 :
2460 106 : if( fp == NULL )
2461 : {
2462 : CPLError( CE_Failure, CPLE_OpenFailed,
2463 : "Attempt to create file `%s' failed.\n",
2464 4 : pszFilename );
2465 4 : return NULL;
2466 : }
2467 :
2468 : /* -------------------------------------------------------------------- */
2469 : /* Just write out a couple of bytes to establish the binary */
2470 : /* file, and then close it. */
2471 : /* -------------------------------------------------------------------- */
2472 102 : VSIFWriteL( (void *) "\0\0", 2, 1, fp );
2473 102 : VSIFCloseL( fp );
2474 :
2475 : /* -------------------------------------------------------------------- */
2476 : /* Create the .hdr filename. */
2477 : /* -------------------------------------------------------------------- */
2478 : const char *pszHDRFilename;
2479 : const char *pszSuffix;
2480 :
2481 102 : pszSuffix = CSLFetchNameValue( papszOptions, "SUFFIX" );
2482 102 : if ( pszSuffix && EQUALN( pszSuffix, "ADD", 3 ))
2483 0 : pszHDRFilename = CPLFormFilename( NULL, pszFilename, "hdr" );
2484 : else
2485 102 : pszHDRFilename = CPLResetExtension(pszFilename, "hdr" );
2486 :
2487 : /* -------------------------------------------------------------------- */
2488 : /* Open the file. */
2489 : /* -------------------------------------------------------------------- */
2490 102 : fp = VSIFOpenL( pszHDRFilename, "wt" );
2491 102 : if( fp == NULL )
2492 : {
2493 : CPLError( CE_Failure, CPLE_OpenFailed,
2494 : "Attempt to create file `%s' failed.\n",
2495 0 : pszHDRFilename );
2496 0 : return NULL;
2497 : }
2498 :
2499 : /* -------------------------------------------------------------------- */
2500 : /* Write out the header. */
2501 : /* -------------------------------------------------------------------- */
2502 : int iBigEndian;
2503 : const char *pszInterleaving;
2504 :
2505 : #ifdef CPL_LSB
2506 102 : iBigEndian = 0;
2507 : #else
2508 : iBigEndian = 1;
2509 : #endif
2510 :
2511 102 : VSIFPrintfL( fp, "ENVI\n" );
2512 : VSIFPrintfL( fp, "samples = %d\nlines = %d\nbands = %d\n",
2513 102 : nXSize, nYSize, nBands );
2514 102 : VSIFPrintfL( fp, "header offset = 0\nfile type = ENVI Standard\n" );
2515 102 : VSIFPrintfL( fp, "data type = %d\n", iENVIType );
2516 102 : pszInterleaving = CSLFetchNameValue( papszOptions, "INTERLEAVE" );
2517 102 : if ( pszInterleaving )
2518 : {
2519 0 : if ( EQUALN( pszInterleaving, "bip", 3 ) )
2520 0 : pszInterleaving = "bip"; // interleaved by pixel
2521 0 : else if ( EQUALN( pszInterleaving, "bil", 3 ) )
2522 0 : pszInterleaving = "bil"; // interleaved by line
2523 : else
2524 0 : pszInterleaving = "bsq"; // band sequental by default
2525 : }
2526 : else
2527 102 : pszInterleaving = "bsq";
2528 102 : VSIFPrintfL( fp, "interleave = %s\n", pszInterleaving);
2529 102 : VSIFPrintfL( fp, "byte order = %d\n", iBigEndian );
2530 :
2531 102 : VSIFCloseL( fp );
2532 :
2533 102 : return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
2534 : }
2535 :
2536 : /************************************************************************/
2537 : /* GDALRegister_ENVI() */
2538 : /************************************************************************/
2539 :
2540 1135 : void GDALRegister_ENVI()
2541 : {
2542 : GDALDriver *poDriver;
2543 :
2544 1135 : if( GDALGetDriverByName( "ENVI" ) == NULL )
2545 : {
2546 1093 : poDriver = new GDALDriver();
2547 :
2548 1093 : poDriver->SetDescription( "ENVI" );
2549 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
2550 1093 : "ENVI .hdr Labelled" );
2551 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
2552 1093 : "frmt_various.html#ENVI" );
2553 1093 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "" );
2554 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
2555 : "Byte Int16 UInt16 Int32 UInt32 "
2556 1093 : "Float32 Float64 CFloat32 CFloat64" );
2557 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
2558 : "<CreationOptionList>"
2559 : " <Option name='SUFFIX' type='string-select'>"
2560 : " <Value>ADD</Value>"
2561 : " </Option>"
2562 : " <Option name='INTERLEAVE' type='string-select'>"
2563 : " <Value>BIP</Value>"
2564 : " <Value>BIL</Value>"
2565 : " <Value>BSQ</Value>"
2566 : " </Option>"
2567 1093 : "</CreationOptionList>" );
2568 :
2569 1093 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
2570 1093 : poDriver->pfnOpen = ENVIDataset::Open;
2571 1093 : poDriver->pfnCreate = ENVIDataset::Create;
2572 :
2573 1093 : GetGDALDriverManager()->RegisterDriver( poDriver );
2574 : }
2575 1135 : }
|
