1 : /******************************************************************************
2 : * $Id: envidataset.cpp 24764 2012-08-11 18:14:56Z 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 24764 2012-08-11 18:14:56Z 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 113 : ENVIDataset::ENVIDataset()
292 : {
293 113 : fpImage = NULL;
294 113 : fp = NULL;
295 113 : pszHDRFilename = NULL;
296 113 : pszProjection = CPLStrdup("");
297 :
298 113 : papszHeader = NULL;
299 :
300 113 : bFoundMapinfo = FALSE;
301 :
302 113 : bHeaderDirty = FALSE;
303 :
304 113 : adfGeoTransform[0] = 0.0;
305 113 : adfGeoTransform[1] = 1.0;
306 113 : adfGeoTransform[2] = 0.0;
307 113 : adfGeoTransform[3] = 0.0;
308 113 : adfGeoTransform[4] = 0.0;
309 113 : adfGeoTransform[5] = 1.0;
310 113 : }
311 :
312 : /************************************************************************/
313 : /* ~ENVIDataset() */
314 : /************************************************************************/
315 :
316 113 : ENVIDataset::~ENVIDataset()
317 :
318 : {
319 113 : FlushCache();
320 113 : if( fpImage )
321 111 : VSIFCloseL( fpImage );
322 113 : if( fp )
323 113 : VSIFCloseL( fp );
324 113 : CPLFree( pszProjection );
325 113 : CSLDestroy( papszHeader );
326 113 : CPLFree(pszHDRFilename);
327 113 : }
328 :
329 : /************************************************************************/
330 : /* FlushCache() */
331 : /************************************************************************/
332 :
333 113 : void ENVIDataset::FlushCache()
334 :
335 : {
336 113 : RawDataset::FlushCache();
337 :
338 113 : GDALRasterBand* band = (GetRasterCount() > 0) ? GetRasterBand(1) : NULL;
339 :
340 113 : if ( band == NULL || !(bHeaderDirty || band->GetCategoryNames() != NULL) )
341 58 : return;
342 :
343 55 : CPLLocaleC oLocaleEnforcer;
344 :
345 55 : 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 55 : iBigEndian = 0;
356 : #else
357 : iBigEndian = 1;
358 : #endif
359 :
360 55 : VSIFPrintfL( fp, "ENVI\n" );
361 55 : if ("" != sDescription)
362 55 : VSIFPrintfL( fp, "description = {\n%s}\n", sDescription.c_str());
363 : VSIFPrintfL( fp, "samples = %d\nlines = %d\nbands = %d\n",
364 55 : nRasterXSize, nRasterYSize, nBands );
365 :
366 55 : catNames = band->GetCategoryNames();
367 :
368 55 : VSIFPrintfL( fp, "header offset = 0\n");
369 55 : if (0 == catNames)
370 51 : VSIFPrintfL( fp, "file type = ENVI Standard\n" );
371 : else
372 4 : VSIFPrintfL( fp, "file type = ENVI Classification\n" );
373 :
374 55 : int iENVIType = GetEnviType(band->GetRasterDataType());
375 55 : VSIFPrintfL( fp, "data type = %d\n", iENVIType );
376 55 : 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 55 : pszInterleaving = "bsq"; // band sequental by default
386 55 : break;
387 : default:
388 0 : pszInterleaving = "bsq";
389 : break;
390 : }
391 55 : VSIFPrintfL( fp, "interleave = %s\n", pszInterleaving);
392 55 : VSIFPrintfL( fp, "byte order = %d\n", iBigEndian );
393 :
394 : /* -------------------------------------------------------------------- */
395 : /* Write class and color information */
396 : /* -------------------------------------------------------------------- */
397 55 : catNames = band->GetCategoryNames();
398 55 : if (0 != catNames)
399 : {
400 4 : int nrClasses = 0;
401 16 : while (*catNames++)
402 8 : ++nrClasses;
403 :
404 4 : if (nrClasses > 0)
405 : {
406 4 : VSIFPrintfL( fp, "classes = %d\n", nrClasses );
407 :
408 4 : GDALColorTable* colorTable = band->GetColorTable();
409 4 : if (0 != colorTable)
410 : {
411 4 : int nrColors = colorTable->GetColorEntryCount();
412 4 : if (nrColors > nrClasses)
413 0 : nrColors = nrClasses;
414 4 : VSIFPrintfL( fp, "class lookup = {\n");
415 12 : for (int i = 0; i < nrColors; ++i)
416 : {
417 8 : const GDALColorEntry* color = colorTable->GetColorEntry(i);
418 8 : VSIFPrintfL(fp, "%d, %d, %d", color->c1, color->c2, color->c3);
419 8 : if (i < nrColors - 1)
420 : {
421 4 : VSIFPrintfL(fp, ", ");
422 4 : if (0 == (i+1) % 5)
423 0 : VSIFPrintfL(fp, "\n");
424 : }
425 : }
426 4 : VSIFPrintfL(fp, "}\n");
427 : }
428 :
429 4 : catNames = band->GetCategoryNames();
430 4 : if (0 != *catNames)
431 : {
432 4 : VSIFPrintfL( fp, "class names = {\n%s", *catNames++);
433 4 : int i = 0;
434 12 : while (*catNames) {
435 4 : VSIFPrintfL( fp, ",");
436 4 : if (0 == (++i) % 5)
437 0 : VSIFPrintfL(fp, "\n");
438 4 : VSIFPrintfL( fp, " %s", *catNames++);
439 : }
440 4 : 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 55 : if ( !WriteRpcInfo() ) // are rpcs in the metadata
454 : {
455 51 : if ( !WritePseudoGcpInfo() ) // are gcps in the metadata
456 : {
457 51 : WriteProjectionInfo(); // standard - affine xform/coord sys str
458 : }
459 : }
460 :
461 :
462 55 : VSIFPrintfL( fp, "band names = {\n" );
463 150 : for ( int i = 1; i <= nBands; i++ )
464 : {
465 95 : CPLString sBandDesc = GetRasterBand( i )->GetDescription();
466 :
467 95 : if ( sBandDesc == "" )
468 84 : sBandDesc = CPLSPrintf( "Band %d", i );
469 95 : VSIFPrintfL( fp, "%s", sBandDesc.c_str() );
470 95 : if ( i != nBands )
471 40 : VSIFPrintfL( fp, ",\n" );
472 : }
473 55 : VSIFPrintfL( fp, "}\n" );
474 : }
475 :
476 : /************************************************************************/
477 : /* GetFileList() */
478 : /************************************************************************/
479 :
480 10 : char **ENVIDataset::GetFileList()
481 :
482 : {
483 10 : char **papszFileList = NULL;
484 :
485 : // Main data file, etc.
486 10 : papszFileList = GDALPamDataset::GetFileList();
487 :
488 : // Header file.
489 10 : papszFileList = CSLAddString( papszFileList, pszHDRFilename );
490 :
491 : // Statistics file
492 10 : if (osStaFilename.size() != 0)
493 0 : papszFileList = CSLAddString( papszFileList, osStaFilename );
494 :
495 10 : 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 46 : static int ENVIGetEPSGGeogCS( OGRSpatialReference *poThis )
509 :
510 : {
511 46 : const char *pszAuthName = poThis->GetAuthorityName( "GEOGCS" );
512 :
513 : /* -------------------------------------------------------------------- */
514 : /* Do we already have it? */
515 : /* -------------------------------------------------------------------- */
516 46 : if( pszAuthName != NULL && EQUAL(pszAuthName,"epsg") )
517 16 : return atoi(poThis->GetAuthorityCode( "GEOGCS" ));
518 :
519 : /* -------------------------------------------------------------------- */
520 : /* Get the datum and geogcs names. */
521 : /* -------------------------------------------------------------------- */
522 30 : const char *pszGEOGCS = poThis->GetAttrValue( "GEOGCS" );
523 30 : const char *pszDatum = poThis->GetAttrValue( "DATUM" );
524 :
525 : // We can only operate on coordinate systems with a geogcs.
526 30 : 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 30 : || 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 30 : || strstr(pszDatum, "North_American");
547 :
548 30 : if( bWGS && (strstr(pszGEOGCS,"84") || strstr(pszDatum,"84")) )
549 25 : return 4326;
550 :
551 5 : if( bWGS && (strstr(pszGEOGCS,"72") || strstr(pszDatum,"72")) )
552 0 : return 4322;
553 :
554 5 : if( bNAD && (strstr(pszGEOGCS,"83") || strstr(pszDatum,"83")) )
555 1 : return 4269;
556 :
557 4 : 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 4 : pszAuthName = poThis->GetAuthorityName( "GEOGCS|DATUM" );
565 :
566 4 : 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 4 : return -1;
577 : }
578 :
579 : /************************************************************************/
580 : /* WriteProjectionInfo() */
581 : /************************************************************************/
582 :
583 51 : void ENVIDataset::WriteProjectionInfo()
584 :
585 : {
586 : /* -------------------------------------------------------------------- */
587 : /* Format the location (geotransform) portion of the map info */
588 : /* line. */
589 : /* -------------------------------------------------------------------- */
590 51 : CPLString osLocation;
591 :
592 : osLocation.Printf( "1, 1, %.15g, %.15g, %.15g, %.15g",
593 : adfGeoTransform[0], adfGeoTransform[3],
594 51 : adfGeoTransform[1], fabs(adfGeoTransform[5]) );
595 :
596 : /* -------------------------------------------------------------------- */
597 : /* Minimal case - write out simple geotransform if we have a */
598 : /* non-default geotransform. */
599 : /* -------------------------------------------------------------------- */
600 51 : if( pszProjection == NULL || strlen(pszProjection) == 0 )
601 : {
602 21 : if( adfGeoTransform[0] != 0.0 || adfGeoTransform[1] != 1.0
603 8 : || adfGeoTransform[2] != 0.0 || adfGeoTransform[3] != 0.0
604 8 : || adfGeoTransform[4] != 0.0 || adfGeoTransform[5] != 1.0 )
605 : {
606 1 : const char* pszHemisphere = "North";
607 : VSIFPrintfL( fp, "map info = {Unknown, %s, %d, %s}\n",
608 1 : osLocation.c_str(), 0, pszHemisphere);
609 : }
610 : return;
611 : }
612 :
613 : /* -------------------------------------------------------------------- */
614 : /* Ingest WKT. */
615 : /* -------------------------------------------------------------------- */
616 46 : OGRSpatialReference oSRS;
617 :
618 46 : char *pszProj = pszProjection;
619 :
620 46 : 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 46 : int nEPSG_GCS = ENVIGetEPSGGeogCS( &oSRS );
628 46 : CPLString osDatum, osCommaDatum;
629 : double dfA, dfB;
630 :
631 46 : if( nEPSG_GCS == 4326 )
632 40 : osDatum = "WGS-84";
633 6 : else if( nEPSG_GCS == 4322 )
634 0 : osDatum = "WGS-72";
635 6 : else if( nEPSG_GCS == 4269 )
636 1 : osDatum = "North America 1983";
637 5 : else if( nEPSG_GCS == 4267 )
638 0 : osDatum = "North America 1927";
639 5 : else if( nEPSG_GCS == 4230 )
640 0 : osDatum = "European 1950";
641 5 : else if( nEPSG_GCS == 4277 )
642 1 : osDatum = "Ordnance Survey of Great Britain '36";
643 4 : else if( nEPSG_GCS == 4291 )
644 0 : osDatum = "SAD-69/Brazil";
645 4 : else if( nEPSG_GCS == 4283 )
646 0 : osDatum = "Geocentric Datum of Australia 1994";
647 4 : else if( nEPSG_GCS == 4275 )
648 0 : osDatum = "Nouvelle Triangulation Francaise IGN";
649 :
650 46 : if( osDatum != "" )
651 42 : osCommaDatum.Printf( ",%s", osDatum.c_str() );
652 :
653 46 : dfA = oSRS.GetSemiMajor();
654 46 : dfB = oSRS.GetSemiMinor();
655 :
656 : /* -------------------------------------------------------------------- */
657 : /* Do we have unusual linear units? */
658 : /* -------------------------------------------------------------------- */
659 46 : CPLString osOptionalUnits;
660 46 : 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 46 : const char *pszProjName = oSRS.GetAttrValue("PROJECTION");
668 : int bNorth;
669 : int iUTMZone;
670 :
671 46 : iUTMZone = oSRS.GetUTMZone( &bNorth );
672 46 : 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 46 : else if( oSRS.IsGeographic() )
684 : {
685 : VSIFPrintfL( fp, "map info = {Geographic Lat/Lon, %s%s}\n",
686 40 : osLocation.c_str(), osCommaDatum.c_str());
687 : }
688 6 : else if( pszProjName == NULL )
689 : {
690 : // what to do?
691 : }
692 6 : 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 6 : 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 1 : 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 1 : osCommaDatum.c_str() );
720 : }
721 6 : 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 1 : 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 1 : osCommaDatum.c_str() );
737 : }
738 4 : 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 4 : 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 4 : 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 4 : 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 3 : 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 3 : osCommaDatum.c_str() );
804 : }
805 1 : 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 1 : 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 1 : 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 1 : 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 46 : if ( oSRS.morphToESRI() == OGRERR_NONE )
869 : {
870 46 : char *pszProjESRI = NULL;
871 46 : if ( oSRS.exportToWkt(&pszProjESRI) == OGRERR_NONE )
872 : {
873 46 : if ( strlen(pszProjESRI) )
874 46 : VSIFPrintfL( fp, "coordinate system string = {%s}\n", pszProjESRI);
875 : }
876 46 : CPLFree(pszProjESRI);
877 46 : pszProjESRI = NULL;
878 46 : }
879 : }
880 :
881 : /************************************************************************/
882 : /* ParseRpcCoeffsMetaDataString() */
883 : /************************************************************************/
884 :
885 16 : 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 16 : const char *psz20Vals = GetMetadataItem(psName, "RPC");
890 16 : if (!psz20Vals)
891 0 : return FALSE;
892 :
893 16 : char** papszArr = CSLTokenizeString2(psz20Vals, " ", 0);
894 16 : if (!papszArr)
895 0 : return FALSE;
896 :
897 16 : int x = 0;
898 352 : while ((papszArr[x] != NULL) && (x < 20))
899 : {
900 320 : papszVal[idx++] = CPLStrdup(papszArr[x]);
901 320 : x++;
902 : }
903 :
904 16 : CSLDestroy(papszArr);
905 :
906 16 : return (x == 20);
907 : }
908 :
909 : #define CPLStrdupIfNotNull(x) ((x) ? CPLStrdup(x) : NULL)
910 :
911 : /************************************************************************/
912 : /* WriteRpcInfo() */
913 : /************************************************************************/
914 :
915 55 : 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 55 : int bRet = FALSE;
920 55 : int x, idx = 0;
921 : char* papszVal[93];
922 :
923 55 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("LINE_OFF", "RPC"));
924 55 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("SAMP_OFF", "RPC"));
925 55 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("LAT_OFF", "RPC"));
926 55 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("LONG_OFF", "RPC"));
927 55 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("HEIGHT_OFF", "RPC"));
928 55 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("LINE_SCALE", "RPC"));
929 55 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("SAMP_SCALE", "RPC"));
930 55 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("LAT_SCALE", "RPC"));
931 55 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("LONG_SCALE", "RPC"));
932 55 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("HEIGHT_SCALE", "RPC"));
933 :
934 95 : for (x=0; x<10; x++) // if we do not have 10 values we return 0
935 : {
936 91 : if (!papszVal[x])
937 51 : goto end;
938 : }
939 :
940 4 : if (!ParseRpcCoeffsMetaDataString("LINE_NUM_COEFF", papszVal, idx))
941 0 : goto end;
942 :
943 4 : if (!ParseRpcCoeffsMetaDataString("LINE_DEN_COEFF", papszVal, idx))
944 0 : goto end;
945 :
946 4 : if (!ParseRpcCoeffsMetaDataString("SAMP_NUM_COEFF", papszVal, idx))
947 0 : goto end;
948 :
949 4 : if (!ParseRpcCoeffsMetaDataString("SAMP_DEN_COEFF", papszVal, idx))
950 0 : goto end;
951 :
952 4 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("TILE_ROW_OFFSET", "RPC"));
953 4 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("TILE_COL_OFFSET", "RPC"));
954 4 : papszVal[idx++] = CPLStrdupIfNotNull(GetMetadataItem("ENVI_RPC_EMULATION", "RPC"));
955 4 : CPLAssert(idx == 93);
956 16 : for (x=90; x<93; x++)
957 : {
958 12 : 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 4 : x = 1;
964 4 : VSIFPrintfL(fp, "rpc info = {\n");
965 376 : for (int iR=0; iR<93; iR++)
966 : {
967 372 : if (papszVal[iR][0] == '-')
968 32 : VSIFPrintfL(fp, " %s", papszVal[iR]);
969 : else
970 340 : VSIFPrintfL(fp, " %s", papszVal[iR]);
971 :
972 372 : if (iR<92)
973 368 : VSIFPrintfL(fp, ",");
974 :
975 372 : if ((x % 4) == 0)
976 92 : VSIFPrintfL(fp, "\n");
977 :
978 372 : x++;
979 372 : if (x > 4)
980 92 : x = 1;
981 : }
982 :
983 4 : VSIFPrintfL(fp, "}\n" );
984 :
985 4 : bRet = TRUE;
986 :
987 : end:
988 937 : for (x=0;x<idx;x++)
989 882 : CPLFree(papszVal[x]);
990 :
991 55 : return bRet;
992 : }
993 :
994 : /************************************************************************/
995 : /* WritePseudoGcpInfo() */
996 : /************************************************************************/
997 :
998 51 : int ENVIDataset::WritePseudoGcpInfo()
999 : {
1000 : // write out gcps into the envi header
1001 : // returns 0 if the gcps are not present
1002 :
1003 51 : int iNum = GetGCPCount();
1004 51 : if (iNum == 0)
1005 51 : 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 20 : const char *ENVIDataset::GetProjectionRef()
1034 :
1035 : {
1036 20 : return pszProjection;
1037 : }
1038 :
1039 : /************************************************************************/
1040 : /* SetProjection() */
1041 : /************************************************************************/
1042 :
1043 46 : CPLErr ENVIDataset::SetProjection( const char *pszNewProjection )
1044 :
1045 : {
1046 46 : CPLFree( pszProjection );
1047 46 : pszProjection = CPLStrdup( pszNewProjection );
1048 :
1049 46 : bHeaderDirty = TRUE;
1050 :
1051 46 : return CE_None;
1052 : }
1053 :
1054 : /************************************************************************/
1055 : /* GetGeoTransform() */
1056 : /************************************************************************/
1057 :
1058 38 : CPLErr ENVIDataset::GetGeoTransform( double * padfTransform )
1059 :
1060 : {
1061 38 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
1062 :
1063 38 : if( bFoundMapinfo )
1064 36 : return CE_None;
1065 : else
1066 2 : return CE_Failure;
1067 : }
1068 :
1069 : /************************************************************************/
1070 : /* SetGeoTransform() */
1071 : /************************************************************************/
1072 :
1073 47 : CPLErr ENVIDataset::SetGeoTransform( double * padfTransform )
1074 : {
1075 47 : memcpy( adfGeoTransform, padfTransform, sizeof(double) * 6 );
1076 :
1077 47 : bHeaderDirty = TRUE;
1078 47 : bFoundMapinfo = TRUE;
1079 :
1080 47 : return CE_None;
1081 : }
1082 :
1083 : /************************************************************************/
1084 : /* SetMetadata() */
1085 : /************************************************************************/
1086 :
1087 21 : CPLErr ENVIDataset::SetMetadata( char ** papszMetadata,
1088 : const char * pszDomain )
1089 : {
1090 21 : if (pszDomain && EQUAL(pszDomain, "RPC"))
1091 1 : bHeaderDirty = TRUE;
1092 :
1093 21 : return RawDataset::SetMetadata(papszMetadata, pszDomain);
1094 : }
1095 :
1096 : /************************************************************************/
1097 : /* SetMetadataItem() */
1098 : /************************************************************************/
1099 :
1100 1387 : CPLErr ENVIDataset::SetMetadataItem( const char * pszName,
1101 : const char * pszValue,
1102 : const char * pszDomain )
1103 : {
1104 1387 : if (pszDomain && EQUAL(pszDomain, "RPC"))
1105 63 : bHeaderDirty = TRUE;
1106 :
1107 1387 : 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 196 : char **ENVIDataset::SplitList( const char *pszCleanInput )
1130 :
1131 : {
1132 196 : char **papszReturn = NULL;
1133 196 : char *pszInput = CPLStrdup(pszCleanInput);
1134 :
1135 196 : if( pszInput[0] != '{' )
1136 : {
1137 58 : CPLFree(pszInput);
1138 58 : return NULL;
1139 : }
1140 :
1141 138 : int iChar=1;
1142 :
1143 :
1144 1264 : while( pszInput[iChar] != '}' && pszInput[iChar] != '\0' )
1145 : {
1146 988 : int iFStart=-1, iFEnd=-1;
1147 :
1148 : // Find start of token.
1149 988 : iFStart = iChar;
1150 3028 : while( pszInput[iFStart] == ' ' )
1151 1052 : iFStart++;
1152 :
1153 988 : iFEnd = iFStart;
1154 19634 : while( pszInput[iFEnd] != ','
1155 5978 : && pszInput[iFEnd] != '}'
1156 5840 : && pszInput[iFEnd] != '\0' )
1157 5840 : iFEnd++;
1158 :
1159 988 : if( pszInput[iFEnd] == '\0' )
1160 0 : break;
1161 :
1162 988 : iChar = iFEnd + 1;
1163 988 : iFEnd = iFEnd - 1;
1164 :
1165 1976 : while( iFEnd > iFStart && pszInput[iFEnd] == ' ' )
1166 0 : iFEnd--;
1167 :
1168 988 : pszInput[iFEnd + 1] = '\0';
1169 988 : papszReturn = CSLAddString( papszReturn, pszInput + iFStart );
1170 : }
1171 :
1172 138 : CPLFree( pszInput );
1173 :
1174 138 : 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 7 : void ENVIDataset::SetENVIEllipse( OGRSpatialReference *poSRS,
1233 : char **papszPI_EI )
1234 :
1235 : {
1236 7 : double dfA = CPLAtofM(papszPI_EI[0]);
1237 7 : double dfB = CPLAtofM(papszPI_EI[1]);
1238 : double dfInvF;
1239 :
1240 7 : if( fabs(dfA-dfB) < 0.1 )
1241 0 : dfInvF = 0.0; // sphere
1242 : else
1243 7 : dfInvF = dfA / (dfA - dfB);
1244 :
1245 :
1246 : poSRS->SetGeogCS( "Ellipse Based", "Ellipse Based", "Unnamed",
1247 7 : dfA, dfInvF );
1248 7 : }
1249 :
1250 : /************************************************************************/
1251 : /* ProcessMapinfo() */
1252 : /* */
1253 : /* Extract projection, and geotransform from a mapinfo value in */
1254 : /* the header. */
1255 : /************************************************************************/
1256 :
1257 52 : int ENVIDataset::ProcessMapinfo( const char *pszMapinfo )
1258 :
1259 : {
1260 : char **papszFields;
1261 : int nCount;
1262 52 : OGRSpatialReference oSRS;
1263 :
1264 52 : papszFields = SplitList( pszMapinfo );
1265 52 : nCount = CSLCount(papszFields);
1266 :
1267 52 : 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 52 : char **papszCSS = NULL;
1277 52 : if( CSLFetchNameValue( papszHeader, "coordinate_system_string" ) != NULL )
1278 : {
1279 : papszCSS = CSLTokenizeString2(
1280 : CSLFetchNameValue( papszHeader, "coordinate_system_string" ),
1281 37 : "{}", CSLT_PRESERVEQUOTES );
1282 : }
1283 :
1284 : /* -------------------------------------------------------------------- */
1285 : /* Check if we have projection info, and if so parse it. */
1286 : /* -------------------------------------------------------------------- */
1287 52 : char **papszPI = NULL;
1288 52 : int nPICount = 0;
1289 52 : if( CSLFetchNameValue( papszHeader, "projection_info" ) != NULL )
1290 : {
1291 : papszPI = SplitList(
1292 19 : CSLFetchNameValue( papszHeader, "projection_info" ) );
1293 19 : nPICount = CSLCount(papszPI);
1294 : }
1295 :
1296 : /* -------------------------------------------------------------------- */
1297 : /* Capture geotransform. */
1298 : /* -------------------------------------------------------------------- */
1299 52 : adfGeoTransform[1] = atof(papszFields[5]); // Pixel width
1300 52 : adfGeoTransform[5] = -atof(papszFields[6]); // Pixel height
1301 : adfGeoTransform[0] = // Upper left X coordinate
1302 52 : atof(papszFields[3]) - (atof(papszFields[1]) - 1) * adfGeoTransform[1];
1303 : adfGeoTransform[3] = // Upper left Y coordinate
1304 52 : atof(papszFields[4]) - (atof(papszFields[2]) - 1) * adfGeoTransform[5];
1305 52 : adfGeoTransform[2] = 0.0;
1306 52 : adfGeoTransform[4] = 0.0;
1307 :
1308 : /* -------------------------------------------------------------------- */
1309 : /* Capture projection. */
1310 : /* -------------------------------------------------------------------- */
1311 52 : if ( oSRS.importFromESRI( papszCSS ) != OGRERR_NONE )
1312 : {
1313 15 : oSRS.Clear();
1314 :
1315 22 : if( EQUALN(papszFields[0],"UTM",3) && nCount >= 9 )
1316 : {
1317 7 : oSRS.SetUTM( atoi(papszFields[7]),
1318 14 : !EQUAL(papszFields[8],"South") );
1319 7 : if( nCount >= 10 && strstr(papszFields[9],"=") == NULL )
1320 0 : SetENVIDatum( &oSRS, papszFields[9] );
1321 : else
1322 7 : oSRS.SetWellKnownGeogCS( "NAD27" );
1323 : }
1324 8 : 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 8 : 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 8 : 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 8 : 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 8 : 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 8 : 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 8 : 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 8 : 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 15 : else if( nPICount > 8 && atoi(papszPI[0]) == 9 ) // Albers Equal Area
1376 : {
1377 14 : oSRS.SetACEA( CPLAtofM(papszPI[7]), CPLAtofM(papszPI[8]),
1378 14 : CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1379 35 : CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) );
1380 : }
1381 1 : 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 1 : 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 1 : 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 1 : 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 52 : 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 52 : if( oSRS.GetRoot() == NULL )
1412 1 : 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 52 : if( oSRS.IsProjected() && oSRS.GetAttrNode("GEOGCS") == NULL
1419 : && nPICount > 3 )
1420 : {
1421 : // Do we have a datum on the projection info line?
1422 7 : int iDatum = nPICount-1;
1423 :
1424 : // Ignore units= items.
1425 7 : if( strstr(papszPI[iDatum],"=") != NULL )
1426 0 : iDatum--;
1427 :
1428 : // Skip past the name.
1429 7 : iDatum--;
1430 :
1431 7 : CPLString osDatumName = papszPI[iDatum];
1432 7 : if( osDatumName.find_first_of("abcdefghijklmnopqrstuvwxyz"
1433 : "ABCDEFGHIJKLMNOPQRSTUVWXYZ")
1434 : != CPLString::npos )
1435 : {
1436 0 : SetENVIDatum( &oSRS, osDatumName );
1437 : }
1438 : else
1439 : {
1440 7 : SetENVIEllipse( &oSRS, papszPI + 1 );
1441 7 : }
1442 : }
1443 :
1444 :
1445 : /* -------------------------------------------------------------------- */
1446 : /* Try to process specialized units. */
1447 : /* -------------------------------------------------------------------- */
1448 52 : 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 52 : if( oSRS.GetRoot() != NULL )
1492 : {
1493 52 : oSRS.Fixup();
1494 52 : if ( pszProjection )
1495 : {
1496 52 : CPLFree( pszProjection );
1497 52 : pszProjection = NULL;
1498 : }
1499 52 : oSRS.exportToWkt( &pszProjection );
1500 : }
1501 :
1502 52 : CSLDestroy( papszFields );
1503 52 : CSLDestroy( papszPI );
1504 52 : 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 3 : void ENVIDataset::ProcessRPCinfo( const char *pszRPCinfo,
1515 : int numCols, int numRows)
1516 : {
1517 : char **papszFields;
1518 : char sVal[1280];
1519 : int nCount;
1520 :
1521 3 : papszFields = SplitList( pszRPCinfo );
1522 3 : nCount = CSLCount(papszFields);
1523 :
1524 3 : if( nCount < 90 )
1525 : {
1526 0 : CSLDestroy( papszFields );
1527 0 : return;
1528 : }
1529 :
1530 3 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[0]));
1531 3 : SetMetadataItem("LINE_OFF",sVal,"RPC");
1532 3 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[5]));
1533 3 : SetMetadataItem("LINE_SCALE",sVal,"RPC");
1534 3 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[1]));
1535 3 : SetMetadataItem("SAMP_OFF",sVal,"RPC");
1536 3 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[6]));
1537 3 : SetMetadataItem("SAMP_SCALE",sVal,"RPC");
1538 3 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[2]));
1539 3 : SetMetadataItem("LAT_OFF",sVal,"RPC");
1540 3 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[7]));
1541 3 : SetMetadataItem("LAT_SCALE",sVal,"RPC");
1542 3 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[3]));
1543 3 : SetMetadataItem("LONG_OFF",sVal,"RPC");
1544 3 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[8]));
1545 3 : SetMetadataItem("LONG_SCALE",sVal,"RPC");
1546 3 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[4]));
1547 3 : SetMetadataItem("HEIGHT_OFF",sVal,"RPC");
1548 3 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[9]));
1549 3 : SetMetadataItem("HEIGHT_SCALE",sVal,"RPC");
1550 :
1551 3 : sVal[0] = '\0';
1552 : int i;
1553 63 : for(i = 0; i < 20; i++ )
1554 : snprintf(sVal+strlen(sVal), sizeof(sVal), "%.16g ",
1555 60 : atof(papszFields[10+i]));
1556 3 : SetMetadataItem("LINE_NUM_COEFF",sVal,"RPC");
1557 :
1558 3 : sVal[0] = '\0';
1559 63 : for(i = 0; i < 20; i++ )
1560 : snprintf(sVal+strlen(sVal), sizeof(sVal), "%.16g ",
1561 60 : atof(papszFields[30+i]));
1562 3 : SetMetadataItem("LINE_DEN_COEFF",sVal,"RPC");
1563 :
1564 3 : sVal[0] = '\0';
1565 63 : for(i = 0; i < 20; i++ )
1566 : snprintf(sVal+strlen(sVal), sizeof(sVal), "%.16g ",
1567 60 : atof(papszFields[50+i]));
1568 3 : SetMetadataItem("SAMP_NUM_COEFF",sVal,"RPC");
1569 :
1570 3 : sVal[0] = '\0';
1571 63 : for(i = 0; i < 20; i++ )
1572 : snprintf(sVal+strlen(sVal), sizeof(sVal), "%.16g ",
1573 60 : atof(papszFields[70+i]));
1574 3 : SetMetadataItem("SAMP_DEN_COEFF",sVal,"RPC");
1575 :
1576 : snprintf(sVal, sizeof(sVal), "%.16g",
1577 3 : atof(papszFields[3]) - atof(papszFields[8]));
1578 3 : SetMetadataItem("MIN_LONG",sVal,"RPC");
1579 :
1580 : snprintf(sVal, sizeof(sVal), "%.16g",
1581 3 : atof(papszFields[3]) + atof(papszFields[8]) );
1582 3 : SetMetadataItem("MAX_LONG",sVal,"RPC");
1583 :
1584 : snprintf(sVal, sizeof(sVal), "%.16g",
1585 3 : atof(papszFields[2]) - atof(papszFields[7]));
1586 3 : SetMetadataItem("MIN_LAT",sVal,"RPC");
1587 :
1588 : snprintf(sVal, sizeof(sVal), "%.16g",
1589 3 : atof(papszFields[2]) + atof(papszFields[7]));
1590 3 : SetMetadataItem("MAX_LAT",sVal,"RPC");
1591 :
1592 3 : if (nCount == 93)
1593 : {
1594 3 : SetMetadataItem("TILE_ROW_OFFSET",papszFields[90],"RPC");
1595 3 : SetMetadataItem("TILE_COL_OFFSET",papszFields[91],"RPC");
1596 3 : 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 3 : rowOffset = (nCount == 93) ? atof(papszFields[90]) : 0;
1602 3 : colOffset = (nCount == 93) ? atof(papszFields[91]) : 0;
1603 3 : 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 3 : CSLDestroy(papszFields);
1634 : }
1635 :
1636 111 : void ENVIDataset::ProcessStatsFile()
1637 : {
1638 : VSILFILE *fpStaFile;
1639 :
1640 111 : osStaFilename = CPLResetExtension( pszHDRFilename, "sta" );
1641 111 : fpStaFile = VSIFOpenL( osStaFilename, "rb" );
1642 :
1643 111 : if (!fpStaFile)
1644 : {
1645 110 : osStaFilename = "";
1646 110 : return;
1647 : }
1648 :
1649 : int lTestHeader[10],lOffset;
1650 :
1651 1 : 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 1 : isFloat = (byteSwapInt(lTestHeader[0]) == 1111838282);
1660 :
1661 : int nb,i;
1662 : float * fStats;
1663 : double * dStats, dMin, dMax, dMean, dStd;
1664 :
1665 1 : nb=byteSwapInt(lTestHeader[3]);
1666 :
1667 1 : 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 1 : VSIFSeekL(fpStaFile,40+(nb+1)*4,SEEK_SET);
1675 :
1676 1 : if (VSIFReadL(&lOffset,sizeof(int),1,fpStaFile) == 1)
1677 : {
1678 1 : VSIFSeekL(fpStaFile,40+(nb+1)*8+byteSwapInt(lOffset)+nb,SEEK_SET);
1679 : // This should be the beginning of the statistics
1680 1 : 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 1 : dStats = (double*)CPLCalloc(nb*4,8);
1699 1 : if ((int)VSIFReadL(dStats,8,nb*4,fpStaFile) == nb*4)
1700 : {
1701 7 : for (i=0;i<nb;i++)
1702 : {
1703 6 : dMin = byteSwapDouble(dStats[i]);
1704 6 : dMax = byteSwapDouble(dStats[nb+i]);
1705 6 : dMean = byteSwapDouble(dStats[2*nb+i]);
1706 6 : dStd = byteSwapDouble(dStats[3*nb+i]);
1707 6 : if (dMin != dMax && dStd != 0)
1708 6 : GetRasterBand(i+1)->SetStatistics(dMin,dMax,dMean,dStd);
1709 : }
1710 : }
1711 1 : CPLFree(dStats);
1712 : }
1713 : }
1714 1 : VSIFCloseL( fpStaFile );
1715 : }
1716 :
1717 3 : int ENVIDataset::byteSwapInt(int swapMe)
1718 : {
1719 3 : CPL_MSBPTR32(&swapMe);
1720 3 : return swapMe;
1721 : }
1722 :
1723 0 : float ENVIDataset::byteSwapFloat(float swapMe)
1724 : {
1725 0 : CPL_MSBPTR32(&swapMe);
1726 0 : return swapMe;
1727 : }
1728 :
1729 24 : double ENVIDataset::byteSwapDouble(double swapMe)
1730 : {
1731 24 : CPL_MSBPTR64(&swapMe);
1732 24 : return swapMe;
1733 : }
1734 :
1735 :
1736 : /************************************************************************/
1737 : /* ReadHeader() */
1738 : /************************************************************************/
1739 :
1740 113 : int ENVIDataset::ReadHeader( VSILFILE * fpHdr )
1741 :
1742 : {
1743 :
1744 113 : CPLReadLineL( fpHdr );
1745 :
1746 : /* -------------------------------------------------------------------- */
1747 : /* Now start forming sets of name/value pairs. */
1748 : /* -------------------------------------------------------------------- */
1749 1141 : while( TRUE )
1750 : {
1751 : const char *pszNewLine;
1752 : char *pszWorkingLine;
1753 :
1754 1254 : pszNewLine = CPLReadLineL( fpHdr );
1755 1254 : if( pszNewLine == NULL )
1756 : break;
1757 :
1758 1141 : if( strstr(pszNewLine,"=") == NULL )
1759 0 : continue;
1760 :
1761 1141 : pszWorkingLine = CPLStrdup(pszNewLine);
1762 :
1763 : // Collect additional lines if we have open sqiggly bracket.
1764 1141 : if( strstr(pszWorkingLine,"{") != NULL
1765 : && strstr(pszWorkingLine,"}") == NULL )
1766 : {
1767 217 : do {
1768 217 : pszNewLine = CPLReadLineL( fpHdr );
1769 217 : if( pszNewLine )
1770 : {
1771 : pszWorkingLine = (char *)
1772 : CPLRealloc(pszWorkingLine,
1773 217 : strlen(pszWorkingLine)+strlen(pszNewLine)+1);
1774 217 : 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 27851 : for( iEqual = 0;
1784 26710 : pszWorkingLine[iEqual] != '\0' && pszWorkingLine[iEqual] != '=';
1785 : iEqual++ ) {}
1786 :
1787 1141 : if( pszWorkingLine[iEqual] == '=' )
1788 : {
1789 : int i;
1790 :
1791 1141 : pszValue = pszWorkingLine + iEqual + 1;
1792 3423 : while( *pszValue == ' ' || *pszValue == '\t' )
1793 1141 : pszValue++;
1794 :
1795 1141 : pszWorkingLine[iEqual--] = '\0';
1796 7750 : while( iEqual > 0
1797 2734 : && (pszWorkingLine[iEqual] == ' '
1798 1141 : || pszWorkingLine[iEqual] == '\t') )
1799 1593 : pszWorkingLine[iEqual--] = '\0';
1800 :
1801 : // Convert spaces in the name to underscores.
1802 11762 : for( i = 0; pszWorkingLine[i] != '\0'; i++ )
1803 : {
1804 10621 : if( pszWorkingLine[i] == ' ' )
1805 672 : pszWorkingLine[i] = '_';
1806 : }
1807 :
1808 : papszHeader = CSLSetNameValue( papszHeader,
1809 1141 : pszWorkingLine, pszValue );
1810 : }
1811 :
1812 1141 : CPLFree( pszWorkingLine );
1813 : }
1814 :
1815 113 : return TRUE;
1816 : }
1817 :
1818 : /************************************************************************/
1819 : /* Open() */
1820 : /************************************************************************/
1821 :
1822 12338 : 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 12338 : if( poOpenInfo->nHeaderBytes < 2 )
1831 11319 : 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 1019 : CPLString osHdrFilename;
1840 1019 : VSILFILE *fpHeader = NULL;
1841 :
1842 1019 : if( poOpenInfo->eAccess == GA_Update )
1843 237 : pszMode = "r+";
1844 : else
1845 782 : pszMode = "r";
1846 :
1847 1019 : if (poOpenInfo->papszSiblingFiles == NULL)
1848 : {
1849 4 : osHdrFilename = CPLResetExtension( poOpenInfo->pszFilename, "hdr" );
1850 4 : fpHeader = VSIFOpenL( osHdrFilename, pszMode );
1851 :
1852 4 : if( fpHeader == NULL && VSIIsCaseSensitiveFS(osHdrFilename) )
1853 : {
1854 4 : osHdrFilename = CPLResetExtension( poOpenInfo->pszFilename, "HDR" );
1855 4 : fpHeader = VSIFOpenL( osHdrFilename, pszMode );
1856 : }
1857 :
1858 4 : if( fpHeader == NULL )
1859 : {
1860 : osHdrFilename = CPLFormFilename( NULL, poOpenInfo->pszFilename,
1861 4 : "hdr" );
1862 4 : fpHeader = VSIFOpenL( osHdrFilename, pszMode );
1863 : }
1864 :
1865 4 : if( fpHeader == NULL && VSIIsCaseSensitiveFS(osHdrFilename) )
1866 : {
1867 : osHdrFilename = CPLFormFilename( NULL, poOpenInfo->pszFilename,
1868 4 : "HDR" );
1869 4 : 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 1015 : CPLString osPath = CPLGetPath( poOpenInfo->pszFilename );
1880 1015 : CPLString osName = CPLGetFilename( poOpenInfo->pszFilename );
1881 :
1882 : int iFile = CSLFindString(poOpenInfo->papszSiblingFiles,
1883 1015 : CPLResetExtension( osName, "hdr" ) );
1884 1015 : if( iFile >= 0 )
1885 : {
1886 220 : osHdrFilename = CPLFormFilename( osPath, poOpenInfo->papszSiblingFiles[iFile],
1887 440 : NULL );
1888 220 : fpHeader = VSIFOpenL( osHdrFilename, pszMode );
1889 : }
1890 : else
1891 : {
1892 : iFile = CSLFindString(poOpenInfo->papszSiblingFiles,
1893 795 : CPLFormFilename( NULL, osName, "hdr" ));
1894 795 : if( iFile >= 0 )
1895 : {
1896 0 : osHdrFilename = CPLFormFilename( osPath, poOpenInfo->papszSiblingFiles[iFile],
1897 0 : NULL );
1898 0 : fpHeader = VSIFOpenL( osHdrFilename, pszMode );
1899 : }
1900 1015 : }
1901 : }
1902 :
1903 1019 : if( fpHeader == NULL )
1904 799 : return NULL;
1905 :
1906 : /* -------------------------------------------------------------------- */
1907 : /* Check that the first line says "ENVI". */
1908 : /* -------------------------------------------------------------------- */
1909 : char szTestHdr[4];
1910 :
1911 220 : if( VSIFReadL( szTestHdr, 4, 1, fpHeader ) != 1 )
1912 : {
1913 0 : VSIFCloseL( fpHeader );
1914 0 : return NULL;
1915 : }
1916 220 : if( strncmp(szTestHdr,"ENVI",4) != 0 )
1917 : {
1918 107 : VSIFCloseL( fpHeader );
1919 107 : return NULL;
1920 : }
1921 :
1922 : /* -------------------------------------------------------------------- */
1923 : /* Create a corresponding GDALDataset. */
1924 : /* -------------------------------------------------------------------- */
1925 : ENVIDataset *poDS;
1926 :
1927 113 : poDS = new ENVIDataset();
1928 226 : poDS->pszHDRFilename = CPLStrdup(osHdrFilename);
1929 113 : poDS->fp = fpHeader;
1930 :
1931 : /* -------------------------------------------------------------------- */
1932 : /* Read the header. */
1933 : /* -------------------------------------------------------------------- */
1934 113 : 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 113 : 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 113 : 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 113 : int nLines = 0, nSamples = 0, nBands = 0, nHeaderSize = 0;
1976 113 : const char *pszInterleave = NULL;
1977 :
1978 113 : if( CSLFetchNameValue(poDS->papszHeader,"lines") )
1979 113 : nLines = atoi(CSLFetchNameValue(poDS->papszHeader,"lines"));
1980 :
1981 113 : if( CSLFetchNameValue(poDS->papszHeader,"samples") )
1982 113 : nSamples = atoi(CSLFetchNameValue(poDS->papszHeader,"samples"));
1983 :
1984 113 : if( CSLFetchNameValue(poDS->papszHeader,"bands") )
1985 113 : nBands = atoi(CSLFetchNameValue(poDS->papszHeader,"bands"));
1986 :
1987 113 : 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 113 : if( pszInterleave == NULL )
1992 : {
1993 0 : const char* pszExtension = CPLGetExtension(poOpenInfo->pszFilename);
1994 0 : pszInterleave = pszExtension;
1995 : }
1996 :
1997 113 : 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 113 : if (!GDALCheckDatasetDimensions(nSamples, nLines) || !GDALCheckBandCount(nBands, FALSE))
2006 : {
2007 2 : 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 2 : "keywords appears to be missing or invalid." );
2012 2 : return NULL;
2013 : }
2014 :
2015 111 : if( CSLFetchNameValue(poDS->papszHeader,"header_offset") )
2016 111 : nHeaderSize = atoi(CSLFetchNameValue(poDS->papszHeader,"header_offset"));
2017 :
2018 : /* -------------------------------------------------------------------- */
2019 : /* Translate the datatype. */
2020 : /* -------------------------------------------------------------------- */
2021 111 : GDALDataType eType = GDT_Byte;
2022 :
2023 111 : if( CSLFetchNameValue(poDS->papszHeader,"data_type" ) != NULL )
2024 : {
2025 111 : switch( atoi(CSLFetchNameValue(poDS->papszHeader,"data_type" )) )
2026 : {
2027 : case 1:
2028 65 : eType = GDT_Byte;
2029 65 : break;
2030 :
2031 : case 2:
2032 6 : eType = GDT_Int16;
2033 6 : break;
2034 :
2035 : case 3:
2036 6 : eType = GDT_Int32;
2037 6 : break;
2038 :
2039 : case 4:
2040 6 : eType = GDT_Float32;
2041 6 : break;
2042 :
2043 : case 5:
2044 6 : eType = GDT_Float64;
2045 6 : break;
2046 :
2047 : case 6:
2048 5 : eType = GDT_CFloat32;
2049 5 : break;
2050 :
2051 : case 9:
2052 5 : eType = GDT_CFloat64;
2053 5 : break;
2054 :
2055 : case 12:
2056 6 : eType = GDT_UInt16;
2057 6 : break;
2058 :
2059 : case 13:
2060 6 : eType = GDT_UInt32;
2061 6 : 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 111 : int bNativeOrder = TRUE;
2078 :
2079 111 : if( CSLFetchNameValue(poDS->papszHeader,"byte_order" ) != NULL )
2080 : {
2081 : #ifdef CPL_LSB
2082 : bNativeOrder = atoi(CSLFetchNameValue(poDS->papszHeader,
2083 111 : "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 111 : 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 111 : 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 111 : if(EQUAL(pszEnviFileType, "envi meta file") ||
2113 : EQUAL(pszEnviFileType, "envi virtual mosaic") ||
2114 : EQUAL(pszEnviFileType, "envi spectral library"))
2115 : {
2116 : CPLError( CE_Failure, CPLE_OpenFailed,
2117 : "File %s contains an invalid file type in the ENVI .hdr\n"
2118 : "GDAL does not support '%s' type files.",
2119 0 : poOpenInfo->pszFilename, pszEnviFileType );
2120 0 : delete poDS;
2121 0 : return NULL;
2122 : }
2123 : }
2124 :
2125 : /* -------------------------------------------------------------------- */
2126 : /* Detect (gzipped) compressed datasets. */
2127 : /* -------------------------------------------------------------------- */
2128 111 : int bIsCompressed = FALSE;
2129 111 : if( CSLFetchNameValue(poDS->papszHeader,"file_compression" ) != NULL )
2130 : {
2131 1 : if( atoi(CSLFetchNameValue(poDS->papszHeader,"file_compression" ))
2132 : != 0 )
2133 : {
2134 1 : bIsCompressed = TRUE;
2135 : }
2136 : }
2137 :
2138 : /* -------------------------------------------------------------------- */
2139 : /* Capture some information from the file that is of interest. */
2140 : /* -------------------------------------------------------------------- */
2141 111 : poDS->nRasterXSize = nSamples;
2142 111 : poDS->nRasterYSize = nLines;
2143 111 : poDS->eAccess = poOpenInfo->eAccess;
2144 :
2145 : /* -------------------------------------------------------------------- */
2146 : /* Reopen file in update mode if necessary. */
2147 : /* -------------------------------------------------------------------- */
2148 111 : CPLString osImageFilename(poOpenInfo->pszFilename);
2149 111 : if (bIsCompressed)
2150 1 : osImageFilename = "/vsigzip/" + osImageFilename;
2151 111 : if( poOpenInfo->eAccess == GA_Update )
2152 : {
2153 50 : if (bIsCompressed)
2154 : {
2155 0 : delete poDS;
2156 : CPLError( CE_Failure, CPLE_OpenFailed,
2157 0 : "Cannot open compressed file in update mode.\n");
2158 0 : return NULL;
2159 : }
2160 50 : poDS->fpImage = VSIFOpenL( osImageFilename, "rb+" );
2161 : }
2162 : else
2163 61 : poDS->fpImage = VSIFOpenL( osImageFilename, "rb" );
2164 :
2165 111 : if( poDS->fpImage == NULL )
2166 : {
2167 0 : delete poDS;
2168 : CPLError( CE_Failure, CPLE_OpenFailed,
2169 : "Failed to re-open %s within ENVI driver.\n",
2170 0 : poOpenInfo->pszFilename );
2171 0 : return NULL;
2172 : }
2173 :
2174 : /* -------------------------------------------------------------------- */
2175 : /* Compute the line offset. */
2176 : /* -------------------------------------------------------------------- */
2177 111 : int nDataSize = GDALGetDataTypeSize(eType)/8;
2178 : int nPixelOffset, nLineOffset;
2179 : vsi_l_offset nBandOffset;
2180 111 : int bIntOverflow = FALSE;
2181 :
2182 111 : if( EQUALN(pszInterleave, "bil", 3) )
2183 : {
2184 0 : poDS->interleave = BIL;
2185 0 : poDS->SetMetadataItem( "INTERLEAVE", "LINE", "IMAGE_STRUCTURE" );
2186 0 : if (nSamples > INT_MAX / (nDataSize * nBands)) bIntOverflow = TRUE;
2187 0 : nLineOffset = nDataSize * nSamples * nBands;
2188 0 : nPixelOffset = nDataSize;
2189 0 : nBandOffset = (vsi_l_offset)nDataSize * nSamples;
2190 : }
2191 111 : else if( EQUALN(pszInterleave, "bip", 3) )
2192 : {
2193 0 : poDS->interleave = BIP;
2194 0 : poDS->SetMetadataItem( "INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE" );
2195 0 : if (nSamples > INT_MAX / (nDataSize * nBands)) bIntOverflow = TRUE;
2196 0 : nLineOffset = nDataSize * nSamples * nBands;
2197 0 : nPixelOffset = nDataSize * nBands;
2198 0 : nBandOffset = nDataSize;
2199 : }
2200 : else /* bsq */
2201 : {
2202 111 : poDS->interleave = BSQ;
2203 111 : poDS->SetMetadataItem( "INTERLEAVE", "BAND", "IMAGE_STRUCTURE" );
2204 111 : if (nSamples > INT_MAX / nDataSize) bIntOverflow = TRUE;
2205 111 : nLineOffset = nDataSize * nSamples;
2206 111 : nPixelOffset = nDataSize;
2207 111 : nBandOffset = (vsi_l_offset)nLineOffset * nLines;
2208 : }
2209 :
2210 111 : if (bIntOverflow)
2211 : {
2212 0 : delete poDS;
2213 : CPLError( CE_Failure, CPLE_AppDefined,
2214 0 : "Int overflow occured.");
2215 0 : return NULL;
2216 : }
2217 :
2218 : /* -------------------------------------------------------------------- */
2219 : /* Create band information objects. */
2220 : /* -------------------------------------------------------------------- */
2221 111 : poDS->nBands = nBands;
2222 297 : for( i = 0; i < poDS->nBands; i++ )
2223 : {
2224 : poDS->SetBand( i + 1,
2225 : new RawRasterBand(poDS, i + 1, poDS->fpImage,
2226 : nHeaderSize + nBandOffset * i,
2227 : nPixelOffset, nLineOffset, eType,
2228 186 : bNativeOrder, TRUE) );
2229 : }
2230 :
2231 : /* -------------------------------------------------------------------- */
2232 : /* Apply band names if we have them. */
2233 : /* Use wavelength for more descriptive information if possible */
2234 : /* -------------------------------------------------------------------- */
2235 111 : if( CSLFetchNameValue( poDS->papszHeader, "band_names" ) != NULL ||
2236 : CSLFetchNameValue( poDS->papszHeader, "wavelength" ) != NULL)
2237 : {
2238 : char **papszBandNames =
2239 : poDS->SplitList( CSLFetchNameValue( poDS->papszHeader,
2240 58 : "band_names" ) );
2241 : char **papszWL =
2242 : poDS->SplitList( CSLFetchNameValue( poDS->papszHeader,
2243 58 : "wavelength" ) );
2244 :
2245 58 : const char *pszWLUnits = NULL;
2246 58 : if (papszWL)
2247 : {
2248 : /* If WL information is present, process wavelength units */
2249 : pszWLUnits = CSLFetchNameValue( poDS->papszHeader,
2250 0 : "wavelength_units" );
2251 0 : if (pszWLUnits)
2252 : {
2253 : /* Don't show unknown or index units */
2254 0 : if (EQUAL(pszWLUnits,"Unknown") ||
2255 : EQUAL(pszWLUnits,"Index") )
2256 0 : pszWLUnits=0;
2257 : }
2258 : }
2259 :
2260 146 : for( i = 0; i < nBands; i++ )
2261 : {
2262 88 : CPLString osBandId, osBandName, osWavelength;
2263 :
2264 : /* First set up the wavelength names and units if available */
2265 88 : if (papszWL && CSLCount(papszWL) > i)
2266 : {
2267 0 : osWavelength = papszWL[i];
2268 0 : if (pszWLUnits)
2269 : {
2270 0 : osWavelength += " ";
2271 0 : osWavelength += pszWLUnits;
2272 : }
2273 : }
2274 :
2275 : /* Build the final name for this band */
2276 88 : if (papszBandNames && CSLCount(papszBandNames) > i)
2277 : {
2278 88 : osBandName = papszBandNames[i];
2279 88 : if (strlen(osWavelength) > 0)
2280 : {
2281 0 : osBandName += " (";
2282 0 : osBandName += osWavelength;
2283 0 : osBandName += ")";
2284 : }
2285 : }
2286 : else /* WL but no band names */
2287 0 : osBandName = osWavelength;
2288 :
2289 : /* Description is for internal GDAL usage */
2290 88 : poDS->GetRasterBand(i + 1)->SetDescription( osBandName );
2291 :
2292 : /* Metadata field named Band_1, etc. needed for ArcGIS integration */
2293 88 : osBandId = CPLSPrintf("Band_%i", i+1);
2294 88 : poDS->SetMetadataItem(osBandId, osBandName);
2295 : }
2296 58 : CSLDestroy( papszWL );
2297 58 : CSLDestroy( papszBandNames );
2298 : }
2299 : /* -------------------------------------------------------------------- */
2300 : /* Apply class names if we have them. */
2301 : /* -------------------------------------------------------------------- */
2302 111 : if( CSLFetchNameValue( poDS->papszHeader, "class_names" ) != NULL )
2303 : {
2304 : char **papszClassNames =
2305 : poDS->SplitList( CSLFetchNameValue( poDS->papszHeader,
2306 3 : "class_names" ) );
2307 :
2308 3 : poDS->GetRasterBand(1)->SetCategoryNames( papszClassNames );
2309 3 : CSLDestroy( papszClassNames );
2310 : }
2311 :
2312 : /* -------------------------------------------------------------------- */
2313 : /* Apply colormap if we have one. */
2314 : /* -------------------------------------------------------------------- */
2315 111 : if( CSLFetchNameValue( poDS->papszHeader, "class_lookup" ) != NULL )
2316 : {
2317 : char **papszClassColors =
2318 : poDS->SplitList( CSLFetchNameValue( poDS->papszHeader,
2319 3 : "class_lookup" ) );
2320 3 : int nColorValueCount = CSLCount(papszClassColors);
2321 3 : GDALColorTable oCT;
2322 :
2323 9 : for( i = 0; i*3+2 < nColorValueCount; i++ )
2324 : {
2325 : GDALColorEntry sEntry;
2326 :
2327 6 : sEntry.c1 = (short) atoi(papszClassColors[i*3+0]);
2328 6 : sEntry.c2 = (short) atoi(papszClassColors[i*3+1]);
2329 6 : sEntry.c3 = (short) atoi(papszClassColors[i*3+2]);
2330 6 : sEntry.c4 = 255;
2331 6 : oCT.SetColorEntry( i, &sEntry );
2332 : }
2333 :
2334 3 : CSLDestroy( papszClassColors );
2335 :
2336 3 : poDS->GetRasterBand(1)->SetColorTable( &oCT );
2337 3 : poDS->GetRasterBand(1)->SetColorInterpretation( GCI_PaletteIndex );
2338 : }
2339 :
2340 : /* -------------------------------------------------------------------- */
2341 : /* Set the nodata value if it is present */
2342 : /* -------------------------------------------------------------------- */
2343 111 : if( CSLFetchNameValue(poDS->papszHeader,"data_ignore_value" ) != NULL )
2344 : {
2345 0 : for( i = 0; i < poDS->nBands; i++ )
2346 : ((RawRasterBand*)poDS->GetRasterBand(i+1))->SetNoDataValue(atof(
2347 0 : CSLFetchNameValue(poDS->papszHeader,"data_ignore_value")));
2348 : }
2349 :
2350 : /* -------------------------------------------------------------------- */
2351 : /* Set all the header metadata into the ENVI domain */
2352 : /* -------------------------------------------------------------------- */
2353 : {
2354 111 : char** pTmp = poDS->papszHeader;
2355 1347 : while (*pTmp != NULL)
2356 : {
2357 1125 : char** pTokens = CSLTokenizeString2(*pTmp, "=", CSLT_STRIPLEADSPACES | CSLT_STRIPENDSPACES);
2358 1125 : if (pTokens[0] != NULL && pTokens[1] != NULL && pTokens[2] == NULL)
2359 : {
2360 1125 : poDS->SetMetadataItem(pTokens[0], pTokens[1], "ENVI");
2361 : }
2362 1125 : CSLDestroy(pTokens);
2363 1125 : pTmp++;
2364 : }
2365 : }
2366 :
2367 : /* -------------------------------------------------------------------- */
2368 : /* Read the stats file if it is present */
2369 : /* -------------------------------------------------------------------- */
2370 111 : poDS->ProcessStatsFile();
2371 :
2372 : /* -------------------------------------------------------------------- */
2373 : /* Look for mapinfo */
2374 : /* -------------------------------------------------------------------- */
2375 111 : if( CSLFetchNameValue( poDS->papszHeader, "map_info" ) != NULL )
2376 : {
2377 : poDS->bFoundMapinfo =
2378 : poDS->ProcessMapinfo(
2379 52 : CSLFetchNameValue(poDS->papszHeader,"map_info") );
2380 : }
2381 :
2382 : /* -------------------------------------------------------------------- */
2383 : /* Look for RPC mapinfo */
2384 : /* -------------------------------------------------------------------- */
2385 111 : if( !poDS->bFoundMapinfo &&
2386 : CSLFetchNameValue( poDS->papszHeader, "rpc_info" ) != NULL )
2387 : {
2388 : poDS->ProcessRPCinfo(
2389 : CSLFetchNameValue(poDS->papszHeader,"rpc_info"),
2390 3 : poDS->nRasterXSize, poDS->nRasterYSize);
2391 : }
2392 :
2393 : /* -------------------------------------------------------------------- */
2394 : /* Initialize any PAM information. */
2395 : /* -------------------------------------------------------------------- */
2396 111 : poDS->SetDescription( poOpenInfo->pszFilename );
2397 111 : poDS->TryLoadXML();
2398 :
2399 : /* -------------------------------------------------------------------- */
2400 : /* Check for overviews. */
2401 : /* -------------------------------------------------------------------- */
2402 111 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
2403 :
2404 111 : return( poDS );
2405 : }
2406 :
2407 123 : int ENVIDataset::GetEnviType(GDALDataType eType)
2408 : {
2409 : int iENVIType;
2410 123 : switch( eType )
2411 : {
2412 : case GDT_Byte:
2413 59 : iENVIType = 1;
2414 59 : break;
2415 : case GDT_Int16:
2416 7 : iENVIType = 2;
2417 7 : break;
2418 : case GDT_Int32:
2419 7 : iENVIType = 3;
2420 7 : break;
2421 : case GDT_Float32:
2422 7 : iENVIType = 4;
2423 7 : break;
2424 : case GDT_Float64:
2425 7 : iENVIType = 5;
2426 7 : break;
2427 : case GDT_CFloat32:
2428 7 : iENVIType = 6;
2429 7 : break;
2430 : case GDT_CFloat64:
2431 7 : iENVIType = 9;
2432 7 : break;
2433 : case GDT_UInt16:
2434 7 : iENVIType = 12;
2435 7 : break;
2436 : case GDT_UInt32:
2437 7 : iENVIType = 13;
2438 7 : break;
2439 :
2440 : /* 14=Int64, 15=UInt64 */
2441 :
2442 : default:
2443 : CPLError( CE_Failure, CPLE_AppDefined,
2444 : "Attempt to create ENVI .hdr labelled dataset with an illegal\n"
2445 : "data type (%s).\n",
2446 8 : GDALGetDataTypeName(eType) );
2447 8 : return 1;
2448 : }
2449 115 : return iENVIType;
2450 : }
2451 :
2452 : /************************************************************************/
2453 : /* Create() */
2454 : /************************************************************************/
2455 :
2456 68 : GDALDataset *ENVIDataset::Create( const char * pszFilename,
2457 : int nXSize, int nYSize, int nBands,
2458 : GDALDataType eType,
2459 : char ** papszOptions )
2460 :
2461 : {
2462 : /* -------------------------------------------------------------------- */
2463 : /* Verify input options. */
2464 : /* -------------------------------------------------------------------- */
2465 68 : int iENVIType = GetEnviType(eType);
2466 68 : if (0 == iENVIType)
2467 0 : return 0;
2468 :
2469 : /* -------------------------------------------------------------------- */
2470 : /* Try to create the file. */
2471 : /* -------------------------------------------------------------------- */
2472 : VSILFILE *fp;
2473 :
2474 68 : fp = VSIFOpenL( pszFilename, "wb" );
2475 :
2476 68 : if( fp == NULL )
2477 : {
2478 : CPLError( CE_Failure, CPLE_OpenFailed,
2479 : "Attempt to create file `%s' failed.\n",
2480 17 : pszFilename );
2481 17 : return NULL;
2482 : }
2483 :
2484 : /* -------------------------------------------------------------------- */
2485 : /* Just write out a couple of bytes to establish the binary */
2486 : /* file, and then close it. */
2487 : /* -------------------------------------------------------------------- */
2488 51 : VSIFWriteL( (void *) "\0\0", 2, 1, fp );
2489 51 : VSIFCloseL( fp );
2490 :
2491 : /* -------------------------------------------------------------------- */
2492 : /* Create the .hdr filename. */
2493 : /* -------------------------------------------------------------------- */
2494 : const char *pszHDRFilename;
2495 : const char *pszSuffix;
2496 :
2497 51 : pszSuffix = CSLFetchNameValue( papszOptions, "SUFFIX" );
2498 51 : if ( pszSuffix && EQUALN( pszSuffix, "ADD", 3 ))
2499 0 : pszHDRFilename = CPLFormFilename( NULL, pszFilename, "hdr" );
2500 : else
2501 51 : pszHDRFilename = CPLResetExtension(pszFilename, "hdr" );
2502 :
2503 : /* -------------------------------------------------------------------- */
2504 : /* Open the file. */
2505 : /* -------------------------------------------------------------------- */
2506 51 : fp = VSIFOpenL( pszHDRFilename, "wt" );
2507 51 : if( fp == NULL )
2508 : {
2509 : CPLError( CE_Failure, CPLE_OpenFailed,
2510 : "Attempt to create file `%s' failed.\n",
2511 0 : pszHDRFilename );
2512 0 : return NULL;
2513 : }
2514 :
2515 : /* -------------------------------------------------------------------- */
2516 : /* Write out the header. */
2517 : /* -------------------------------------------------------------------- */
2518 : int iBigEndian;
2519 : const char *pszInterleaving;
2520 :
2521 : #ifdef CPL_LSB
2522 51 : iBigEndian = 0;
2523 : #else
2524 : iBigEndian = 1;
2525 : #endif
2526 :
2527 51 : VSIFPrintfL( fp, "ENVI\n" );
2528 : VSIFPrintfL( fp, "samples = %d\nlines = %d\nbands = %d\n",
2529 51 : nXSize, nYSize, nBands );
2530 51 : VSIFPrintfL( fp, "header offset = 0\nfile type = ENVI Standard\n" );
2531 51 : VSIFPrintfL( fp, "data type = %d\n", iENVIType );
2532 51 : pszInterleaving = CSLFetchNameValue( papszOptions, "INTERLEAVE" );
2533 51 : if ( pszInterleaving )
2534 : {
2535 0 : if ( EQUALN( pszInterleaving, "bip", 3 ) )
2536 0 : pszInterleaving = "bip"; // interleaved by pixel
2537 0 : else if ( EQUALN( pszInterleaving, "bil", 3 ) )
2538 0 : pszInterleaving = "bil"; // interleaved by line
2539 : else
2540 0 : pszInterleaving = "bsq"; // band sequental by default
2541 : }
2542 : else
2543 51 : pszInterleaving = "bsq";
2544 51 : VSIFPrintfL( fp, "interleave = %s\n", pszInterleaving);
2545 51 : VSIFPrintfL( fp, "byte order = %d\n", iBigEndian );
2546 :
2547 51 : VSIFCloseL( fp );
2548 :
2549 51 : return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
2550 : }
2551 :
2552 : /************************************************************************/
2553 : /* GDALRegister_ENVI() */
2554 : /************************************************************************/
2555 :
2556 582 : void GDALRegister_ENVI()
2557 : {
2558 : GDALDriver *poDriver;
2559 :
2560 582 : if( GDALGetDriverByName( "ENVI" ) == NULL )
2561 : {
2562 561 : poDriver = new GDALDriver();
2563 :
2564 561 : poDriver->SetDescription( "ENVI" );
2565 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
2566 561 : "ENVI .hdr Labelled" );
2567 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
2568 561 : "frmt_various.html#ENVI" );
2569 561 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "" );
2570 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
2571 : "Byte Int16 UInt16 Int32 UInt32 "
2572 561 : "Float32 Float64 CFloat32 CFloat64" );
2573 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
2574 : "<CreationOptionList>"
2575 : " <Option name='SUFFIX' type='string-select'>"
2576 : " <Value>ADD</Value>"
2577 : " </Option>"
2578 : " <Option name='INTERLEAVE' type='string-select'>"
2579 : " <Value>BIP</Value>"
2580 : " <Value>BIL</Value>"
2581 : " <Value>BSQ</Value>"
2582 : " </Option>"
2583 561 : "</CreationOptionList>" );
2584 :
2585 561 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
2586 561 : poDriver->pfnOpen = ENVIDataset::Open;
2587 561 : poDriver->pfnCreate = ENVIDataset::Create;
2588 :
2589 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
2590 : }
2591 582 : }
|
