1 : /******************************************************************************
2 : * $Id: envidataset.cpp 23170 2011-10-03 09:35:21Z rprinceley $
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 23170 2011-10-03 09:35:21Z rprinceley $");
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 : int ReadHeader( VSILFILE * );
242 : int ProcessMapinfo( const char * );
243 : void ProcessRPCinfo( const char * ,int ,int);
244 : void ProcessStatsFile();
245 : long byteSwapLong(long);
246 : float byteSwapFloat(float);
247 : double byteSwapDouble(double);
248 : void SetENVIDatum( OGRSpatialReference *, const char * );
249 : void SetENVIEllipse( OGRSpatialReference *, char ** );
250 : void WriteProjectionInfo();
251 :
252 : char **SplitList( const char * );
253 :
254 : enum Interleave { BSQ, BIL, BIP } interleave;
255 : static int GetEnviType(GDALDataType eType);
256 :
257 : public:
258 : ENVIDataset();
259 : ~ENVIDataset();
260 :
261 : virtual void FlushCache( void );
262 : virtual CPLErr GetGeoTransform( double * padfTransform );
263 : virtual CPLErr SetGeoTransform( double * );
264 : virtual const char *GetProjectionRef(void);
265 : virtual CPLErr SetProjection( const char * );
266 : virtual char **GetFileList(void);
267 :
268 : static GDALDataset *Open( GDALOpenInfo * );
269 : static GDALDataset *Create( const char * pszFilename,
270 : int nXSize, int nYSize, int nBands,
271 : GDALDataType eType, char ** papszOptions );
272 : };
273 :
274 : /************************************************************************/
275 : /* ENVIDataset() */
276 : /************************************************************************/
277 :
278 104 : ENVIDataset::ENVIDataset()
279 : {
280 104 : fpImage = NULL;
281 104 : fp = NULL;
282 104 : pszHDRFilename = NULL;
283 104 : pszProjection = CPLStrdup("");
284 :
285 104 : papszHeader = NULL;
286 :
287 104 : bFoundMapinfo = FALSE;
288 :
289 104 : bHeaderDirty = FALSE;
290 :
291 104 : adfGeoTransform[0] = 0.0;
292 104 : adfGeoTransform[1] = 1.0;
293 104 : adfGeoTransform[2] = 0.0;
294 104 : adfGeoTransform[3] = 0.0;
295 104 : adfGeoTransform[4] = 0.0;
296 104 : adfGeoTransform[5] = 1.0;
297 104 : }
298 :
299 : /************************************************************************/
300 : /* ~ENVIDataset() */
301 : /************************************************************************/
302 :
303 104 : ENVIDataset::~ENVIDataset()
304 :
305 : {
306 104 : FlushCache();
307 104 : if( fpImage )
308 102 : VSIFCloseL( fpImage );
309 104 : if( fp )
310 104 : VSIFCloseL( fp );
311 104 : CPLFree( pszProjection );
312 104 : CSLDestroy( papszHeader );
313 104 : CPLFree(pszHDRFilename);
314 104 : }
315 :
316 : /************************************************************************/
317 : /* FlushCache() */
318 : /************************************************************************/
319 :
320 104 : void ENVIDataset::FlushCache()
321 :
322 : {
323 104 : RawDataset::FlushCache();
324 :
325 104 : if ( !bHeaderDirty)
326 57 : return;
327 :
328 47 : CPLLocaleC oLocaleEnforcer;
329 :
330 47 : VSIFSeekL( fp, 0, SEEK_SET );
331 : /* -------------------------------------------------------------------- */
332 : /* Rewrite out the header. */
333 : /* -------------------------------------------------------------------- */
334 : int iBigEndian;
335 :
336 : const char *pszInterleaving;
337 : char** catNames;
338 :
339 : #ifdef CPL_LSB
340 47 : iBigEndian = 0;
341 : #else
342 : iBigEndian = 1;
343 : #endif
344 :
345 47 : VSIFPrintfL( fp, "ENVI\n" );
346 47 : if ("" != sDescription)
347 47 : VSIFPrintfL( fp, "description = {\n%s}\n", sDescription.c_str());
348 : VSIFPrintfL( fp, "samples = %d\nlines = %d\nbands = %d\n",
349 47 : nRasterXSize, nRasterYSize, nBands );
350 :
351 47 : GDALRasterBand* band = GetRasterBand(1);
352 47 : catNames = band->GetCategoryNames();
353 :
354 47 : VSIFPrintfL( fp, "header offset = 0\n");
355 47 : if (0 == catNames)
356 47 : VSIFPrintfL( fp, "file type = ENVI Standard\n" );
357 : else
358 0 : VSIFPrintfL( fp, "file type = ENVI Classification\n" );
359 :
360 47 : int iENVIType = GetEnviType(band->GetRasterDataType());
361 47 : VSIFPrintfL( fp, "data type = %d\n", iENVIType );
362 47 : switch (interleave)
363 : {
364 : case BIP:
365 0 : pszInterleaving = "bip"; // interleaved by pixel
366 0 : break;
367 : case BIL:
368 0 : pszInterleaving = "bil"; // interleaved by line
369 0 : break;
370 : case BSQ:
371 47 : pszInterleaving = "bsq"; // band sequental by default
372 47 : break;
373 : default:
374 0 : pszInterleaving = "bsq";
375 : break;
376 : }
377 47 : VSIFPrintfL( fp, "interleave = %s\n", pszInterleaving);
378 47 : VSIFPrintfL( fp, "byte order = %d\n", iBigEndian );
379 :
380 : /* -------------------------------------------------------------------- */
381 : /* Write class and color information */
382 : /* -------------------------------------------------------------------- */
383 47 : catNames = band->GetCategoryNames();
384 47 : if (0 != catNames)
385 : {
386 0 : int nrClasses = 0;
387 0 : while (*catNames++)
388 0 : ++nrClasses;
389 :
390 0 : if (nrClasses > 0)
391 : {
392 0 : VSIFPrintfL( fp, "classes = %d\n", nrClasses );
393 :
394 0 : GDALColorTable* colorTable = band->GetColorTable();
395 0 : if (0 != colorTable)
396 : {
397 0 : int nrColors = colorTable->GetColorEntryCount();
398 0 : if (nrColors > nrClasses)
399 0 : nrColors = nrClasses;
400 0 : VSIFPrintfL( fp, "class lookup = {\n");
401 0 : for (int i = 0; i < nrColors; ++i)
402 : {
403 0 : const GDALColorEntry* color = colorTable->GetColorEntry(i);
404 0 : VSIFPrintfL(fp, "%d, %d, %d", color->c1, color->c2, color->c3);
405 0 : if (i < nrColors - 1)
406 : {
407 0 : VSIFPrintfL(fp, ", ");
408 0 : if (0 == (i+1) % 5)
409 0 : VSIFPrintfL(fp, "\n");
410 : }
411 : }
412 0 : VSIFPrintfL(fp, "}\n");
413 : }
414 :
415 0 : catNames = band->GetCategoryNames();
416 0 : if (0 != *catNames)
417 : {
418 0 : VSIFPrintfL( fp, "class names = {\n%s", *catNames++);
419 0 : int i = 0;
420 0 : while (*catNames) {
421 0 : VSIFPrintfL( fp, ",");
422 0 : if (0 == (++i) % 5)
423 0 : VSIFPrintfL(fp, "\n");
424 0 : VSIFPrintfL( fp, " %s", *catNames++);
425 : }
426 0 : VSIFPrintfL( fp, "}\n");
427 : }
428 : }
429 : }
430 :
431 : /* -------------------------------------------------------------------- */
432 : /* Write the rest of header. */
433 : /* -------------------------------------------------------------------- */
434 47 : WriteProjectionInfo();
435 :
436 :
437 47 : VSIFPrintfL( fp, "band names = {\n" );
438 134 : for ( int i = 1; i <= nBands; i++ )
439 : {
440 87 : CPLString sBandDesc = GetRasterBand( i )->GetDescription();
441 :
442 87 : if ( sBandDesc == "" )
443 84 : sBandDesc = CPLSPrintf( "Band %d", i );
444 87 : VSIFPrintfL( fp, "%s", sBandDesc.c_str() );
445 87 : if ( i != nBands )
446 40 : VSIFPrintfL( fp, ",\n" );
447 : }
448 47 : VSIFPrintfL( fp, "}\n" );
449 : }
450 :
451 : /************************************************************************/
452 : /* GetFileList() */
453 : /************************************************************************/
454 :
455 8 : char **ENVIDataset::GetFileList()
456 :
457 : {
458 8 : char **papszFileList = NULL;
459 :
460 : // Main data file, etc.
461 8 : papszFileList = GDALPamDataset::GetFileList();
462 :
463 : // Header file.
464 8 : papszFileList = CSLAddString( papszFileList, pszHDRFilename );
465 :
466 8 : return papszFileList;
467 : }
468 :
469 : /************************************************************************/
470 : /* GetEPSGGeogCS() */
471 : /* */
472 : /* Try to establish what the EPSG code for this coordinate */
473 : /* systems GEOGCS might be. Returns -1 if no reasonable guess */
474 : /* can be made. */
475 : /* */
476 : /* TODO: We really need to do some name lookups. */
477 : /************************************************************************/
478 :
479 46 : static int ENVIGetEPSGGeogCS( OGRSpatialReference *poThis )
480 :
481 : {
482 46 : const char *pszAuthName = poThis->GetAuthorityName( "GEOGCS" );
483 :
484 : /* -------------------------------------------------------------------- */
485 : /* Do we already have it? */
486 : /* -------------------------------------------------------------------- */
487 46 : if( pszAuthName != NULL && EQUAL(pszAuthName,"epsg") )
488 16 : return atoi(poThis->GetAuthorityCode( "GEOGCS" ));
489 :
490 : /* -------------------------------------------------------------------- */
491 : /* Get the datum and geogcs names. */
492 : /* -------------------------------------------------------------------- */
493 30 : const char *pszGEOGCS = poThis->GetAttrValue( "GEOGCS" );
494 30 : const char *pszDatum = poThis->GetAttrValue( "DATUM" );
495 :
496 : // We can only operate on coordinate systems with a geogcs.
497 30 : if( pszGEOGCS == NULL || pszDatum == NULL )
498 0 : return -1;
499 :
500 : /* -------------------------------------------------------------------- */
501 : /* Is this a "well known" geographic coordinate system? */
502 : /* -------------------------------------------------------------------- */
503 : int bWGS, bNAD;
504 :
505 : bWGS = strstr(pszGEOGCS,"WGS") != NULL
506 : || strstr(pszDatum, "WGS")
507 : || strstr(pszGEOGCS,"World Geodetic System")
508 : || strstr(pszGEOGCS,"World_Geodetic_System")
509 : || strstr(pszDatum, "World Geodetic System")
510 30 : || strstr(pszDatum, "World_Geodetic_System");
511 :
512 : bNAD = strstr(pszGEOGCS,"NAD") != NULL
513 : || strstr(pszDatum, "NAD")
514 : || strstr(pszGEOGCS,"North American")
515 : || strstr(pszGEOGCS,"North_American")
516 : || strstr(pszDatum, "North American")
517 30 : || strstr(pszDatum, "North_American");
518 :
519 30 : if( bWGS && (strstr(pszGEOGCS,"84") || strstr(pszDatum,"84")) )
520 25 : return 4326;
521 :
522 5 : if( bWGS && (strstr(pszGEOGCS,"72") || strstr(pszDatum,"72")) )
523 0 : return 4322;
524 :
525 5 : if( bNAD && (strstr(pszGEOGCS,"83") || strstr(pszDatum,"83")) )
526 1 : return 4269;
527 :
528 4 : if( bNAD && (strstr(pszGEOGCS,"27") || strstr(pszDatum,"27")) )
529 0 : return 4267;
530 :
531 : /* -------------------------------------------------------------------- */
532 : /* If we know the datum, associate the most likely GCS with */
533 : /* it. */
534 : /* -------------------------------------------------------------------- */
535 4 : pszAuthName = poThis->GetAuthorityName( "GEOGCS|DATUM" );
536 :
537 4 : if( pszAuthName != NULL
538 : && EQUAL(pszAuthName,"epsg")
539 : && poThis->GetPrimeMeridian() == 0.0 )
540 : {
541 0 : int nDatum = atoi(poThis->GetAuthorityCode("GEOGCS|DATUM"));
542 :
543 0 : if( nDatum >= 6000 && nDatum <= 6999 )
544 0 : return nDatum - 2000;
545 : }
546 :
547 4 : return -1;
548 : }
549 :
550 : /************************************************************************/
551 : /* WriteProjectionInfo() */
552 : /************************************************************************/
553 :
554 47 : void ENVIDataset::WriteProjectionInfo()
555 :
556 : {
557 : /* -------------------------------------------------------------------- */
558 : /* Format the location (geotransform) portion of the map info */
559 : /* line. */
560 : /* -------------------------------------------------------------------- */
561 47 : CPLString osLocation;
562 :
563 : osLocation.Printf( "1, 1, %.15g, %.15g, %.15g, %.15g",
564 : adfGeoTransform[0], adfGeoTransform[3],
565 47 : adfGeoTransform[1], fabs(adfGeoTransform[5]) );
566 :
567 : /* -------------------------------------------------------------------- */
568 : /* Minimal case - write out simple geotransform if we have a */
569 : /* non-default geotransform. */
570 : /* -------------------------------------------------------------------- */
571 47 : if( pszProjection == NULL || strlen(pszProjection) == 0 )
572 : {
573 1 : if( adfGeoTransform[0] != 0.0 || adfGeoTransform[1] != 1.0
574 0 : || adfGeoTransform[2] != 0.0 || adfGeoTransform[3] != 0.0
575 0 : || adfGeoTransform[4] != 0.0 || adfGeoTransform[5] != 1.0 )
576 : {
577 1 : const char* pszHemisphere = "North";
578 : VSIFPrintfL( fp, "map info = {Unknown, %s, %d, %s}\n",
579 1 : osLocation.c_str(), 0, pszHemisphere);
580 : }
581 : return;
582 : }
583 :
584 : /* -------------------------------------------------------------------- */
585 : /* Ingest WKT. */
586 : /* -------------------------------------------------------------------- */
587 46 : OGRSpatialReference oSRS;
588 :
589 46 : char *pszProj = pszProjection;
590 :
591 46 : if( oSRS.importFromWkt( &pszProj ) != OGRERR_NONE )
592 : return;
593 :
594 : /* -------------------------------------------------------------------- */
595 : /* Try to translate the datum and get major/minor ellipsoid */
596 : /* values. */
597 : /* -------------------------------------------------------------------- */
598 46 : int nEPSG_GCS = ENVIGetEPSGGeogCS( &oSRS );
599 46 : CPLString osDatum, osCommaDatum;
600 : double dfA, dfB;
601 :
602 46 : if( nEPSG_GCS == 4326 )
603 40 : osDatum = "WGS-84";
604 6 : else if( nEPSG_GCS == 4322 )
605 0 : osDatum = "WGS-72";
606 6 : else if( nEPSG_GCS == 4269 )
607 1 : osDatum = "North America 1983";
608 5 : else if( nEPSG_GCS == 4267 )
609 0 : osDatum = "North America 1927";
610 5 : else if( nEPSG_GCS == 4230 )
611 0 : osDatum = "European 1950";
612 5 : else if( nEPSG_GCS == 4277 )
613 1 : osDatum = "Ordnance Survey of Great Britain '36";
614 4 : else if( nEPSG_GCS == 4291 )
615 0 : osDatum = "SAD-69/Brazil";
616 4 : else if( nEPSG_GCS == 4283 )
617 0 : osDatum = "Geocentric Datum of Australia 1994";
618 4 : else if( nEPSG_GCS == 4275 )
619 0 : osDatum = "Nouvelle Triangulation Francaise IGN";
620 :
621 46 : if( osDatum != "" )
622 42 : osCommaDatum.Printf( ",%s", osDatum.c_str() );
623 :
624 46 : dfA = oSRS.GetSemiMajor();
625 46 : dfB = oSRS.GetSemiMinor();
626 :
627 : /* -------------------------------------------------------------------- */
628 : /* Do we have unusual linear units? */
629 : /* -------------------------------------------------------------------- */
630 46 : CPLString osOptionalUnits;
631 46 : if( fabs(oSRS.GetLinearUnits()-0.3048) < 0.0001 )
632 0 : osOptionalUnits = ", units=Feet";
633 :
634 : /* -------------------------------------------------------------------- */
635 : /* Handle UTM case. */
636 : /* -------------------------------------------------------------------- */
637 : const char *pszHemisphere;
638 46 : const char *pszProjName = oSRS.GetAttrValue("PROJECTION");
639 : int bNorth;
640 : int iUTMZone;
641 :
642 46 : iUTMZone = oSRS.GetUTMZone( &bNorth );
643 46 : if ( iUTMZone )
644 : {
645 0 : if ( bNorth )
646 0 : pszHemisphere = "North";
647 : else
648 0 : pszHemisphere = "South";
649 :
650 : VSIFPrintfL( fp, "map info = {UTM, %s, %d, %s%s%s}\n",
651 : osLocation.c_str(), iUTMZone, pszHemisphere,
652 0 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
653 : }
654 46 : else if( oSRS.IsGeographic() )
655 : {
656 : VSIFPrintfL( fp, "map info = {Geographic Lat/Lon, %s%s}\n",
657 40 : osLocation.c_str(), osCommaDatum.c_str());
658 : }
659 6 : else if( pszProjName == NULL )
660 : {
661 : // what to do?
662 : }
663 6 : else if( EQUAL(pszProjName,SRS_PT_NEW_ZEALAND_MAP_GRID) )
664 : {
665 : VSIFPrintfL( fp, "map info = {New Zealand Map Grid, %s%s%s}\n",
666 : osLocation.c_str(),
667 0 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
668 :
669 : VSIFPrintfL( fp, "projection info = {39, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, New Zealand Map Grid}\n",
670 : dfA, dfB,
671 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
672 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
673 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
674 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
675 0 : osCommaDatum.c_str() );
676 : }
677 6 : else if( EQUAL(pszProjName,SRS_PT_TRANSVERSE_MERCATOR) )
678 : {
679 : VSIFPrintfL( fp, "map info = {Transverse Mercator, %s%s%s}\n",
680 : osLocation.c_str(),
681 1 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
682 :
683 : VSIFPrintfL( fp, "projection info = {3, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Transverse Mercator}\n",
684 : dfA, dfB,
685 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
686 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
687 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
688 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
689 : oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR,1.0),
690 1 : osCommaDatum.c_str() );
691 : }
692 6 : else if( EQUAL(pszProjName,SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP)
693 : || EQUAL(pszProjName,SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP_BELGIUM) )
694 : {
695 : VSIFPrintfL( fp, "map info = {Lambert Conformal Conic, %s%s%s}\n",
696 : osLocation.c_str(),
697 1 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
698 :
699 : VSIFPrintfL( fp, "projection info = {4, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Lambert Conformal Conic}\n",
700 : dfA, dfB,
701 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
702 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
703 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
704 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
705 : oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0),
706 : oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2,0.0),
707 1 : osCommaDatum.c_str() );
708 : }
709 4 : else if( EQUAL(pszProjName,
710 : SRS_PT_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN) )
711 : {
712 : VSIFPrintfL( fp, "map info = {Hotine Oblique Mercator A, %s%s%s}\n",
713 : osLocation.c_str(),
714 0 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
715 :
716 : VSIFPrintfL( fp, "projection info = {5, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Hotine Oblique Mercator A}\n",
717 : dfA, dfB,
718 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
719 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_POINT_1,0.0),
720 : oSRS.GetNormProjParm(SRS_PP_LONGITUDE_OF_POINT_1,0.0),
721 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_POINT_2,0.0),
722 : oSRS.GetNormProjParm(SRS_PP_LONGITUDE_OF_POINT_2,0.0),
723 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
724 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
725 : oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR,1.0),
726 0 : osCommaDatum.c_str() );
727 : }
728 4 : else if( EQUAL(pszProjName,SRS_PT_HOTINE_OBLIQUE_MERCATOR) )
729 : {
730 : VSIFPrintfL( fp, "map info = {Hotine Oblique Mercator B, %s%s%s}\n",
731 : osLocation.c_str(),
732 0 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
733 :
734 : VSIFPrintfL( fp, "projection info = {6, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Hotine Oblique Mercator B}\n",
735 : dfA, dfB,
736 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
737 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
738 : oSRS.GetNormProjParm(SRS_PP_AZIMUTH,0.0),
739 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
740 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
741 : oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR,1.0),
742 0 : osCommaDatum.c_str() );
743 : }
744 4 : else if( EQUAL(pszProjName,SRS_PT_STEREOGRAPHIC)
745 : || EQUAL(pszProjName,SRS_PT_OBLIQUE_STEREOGRAPHIC) )
746 : {
747 : VSIFPrintfL( fp, "map info = {Stereographic (ellipsoid), %s%s%s}\n",
748 : osLocation.c_str(),
749 0 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
750 :
751 : VSIFPrintfL( fp, "projection info = {7, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %s, Stereographic (ellipsoid)}\n",
752 : dfA, dfB,
753 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
754 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
755 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
756 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
757 : oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR,1.0),
758 0 : osCommaDatum.c_str() );
759 : }
760 4 : else if( EQUAL(pszProjName,SRS_PT_ALBERS_CONIC_EQUAL_AREA) )
761 : {
762 : VSIFPrintfL( fp, "map info = {Albers Conical Equal Area, %s%s%s}\n",
763 : osLocation.c_str(),
764 3 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
765 :
766 : VSIFPrintfL( fp, "projection info = {9, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Albers Conical Equal Area}\n",
767 : dfA, dfB,
768 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
769 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
770 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
771 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
772 : oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1,0.0),
773 : oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2,0.0),
774 3 : osCommaDatum.c_str() );
775 : }
776 1 : else if( EQUAL(pszProjName,SRS_PT_POLYCONIC) )
777 : {
778 : VSIFPrintfL( fp, "map info = {Polyconic, %s%s%s}\n",
779 : osLocation.c_str(),
780 0 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
781 :
782 : VSIFPrintfL( fp, "projection info = {10, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Polyconic}\n",
783 : dfA, dfB,
784 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
785 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
786 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
787 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
788 0 : osCommaDatum.c_str() );
789 : }
790 1 : else if( EQUAL(pszProjName,SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA) )
791 : {
792 : VSIFPrintfL( fp, "map info = {Lambert Azimuthal Equal Area, %s%s%s}\n",
793 : osLocation.c_str(),
794 1 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
795 :
796 : VSIFPrintfL( fp, "projection info = {11, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Lambert Azimuthal Equal Area}\n",
797 : dfA, dfB,
798 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
799 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
800 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
801 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
802 1 : osCommaDatum.c_str() );
803 : }
804 0 : else if( EQUAL(pszProjName,SRS_PT_AZIMUTHAL_EQUIDISTANT) )
805 : {
806 : VSIFPrintfL( fp, "map info = {Azimuthal Equadistant, %s%s%s}\n",
807 : osLocation.c_str(),
808 0 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
809 :
810 : VSIFPrintfL( fp, "projection info = {12, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Azimuthal Equadistant}\n",
811 : dfA, dfB,
812 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,0.0),
813 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
814 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
815 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
816 0 : osCommaDatum.c_str() );
817 : }
818 0 : else if( EQUAL(pszProjName,SRS_PT_POLAR_STEREOGRAPHIC) )
819 : {
820 : VSIFPrintfL( fp, "map info = {Polar Stereographic, %s%s%s}\n",
821 : osLocation.c_str(),
822 0 : osCommaDatum.c_str(), osOptionalUnits.c_str() );
823 :
824 : VSIFPrintfL( fp, "projection info = {31, %.16g, %.16g, %.16g, %.16g, %.16g, %.16g%s, Polar Stereographic}\n",
825 : dfA, dfB,
826 : oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN,90.0),
827 : oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN,0.0),
828 : oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING,0.0),
829 : oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING,0.0),
830 0 : osCommaDatum.c_str() );
831 : }
832 : else
833 : {
834 : VSIFPrintfL( fp, "map info = {%s, %s}\n",
835 0 : pszProjName, osLocation.c_str());
836 : }
837 :
838 : // write out coordinate system string
839 46 : if ( oSRS.morphToESRI() == OGRERR_NONE )
840 : {
841 46 : char *pszProjESRI = NULL;
842 46 : if ( oSRS.exportToWkt(&pszProjESRI) == OGRERR_NONE )
843 : {
844 46 : if ( strlen(pszProjESRI) )
845 46 : VSIFPrintfL( fp, "coordinate system string = {%s}\n", pszProjESRI);
846 : }
847 46 : CPLFree(pszProjESRI);
848 46 : pszProjESRI = NULL;
849 46 : }
850 : }
851 :
852 :
853 : /************************************************************************/
854 : /* GetProjectionRef() */
855 : /************************************************************************/
856 :
857 16 : const char *ENVIDataset::GetProjectionRef()
858 :
859 : {
860 16 : return pszProjection;
861 : }
862 :
863 : /************************************************************************/
864 : /* SetProjection() */
865 : /************************************************************************/
866 :
867 46 : CPLErr ENVIDataset::SetProjection( const char *pszNewProjection )
868 :
869 : {
870 46 : CPLFree( pszProjection );
871 46 : pszProjection = CPLStrdup( pszNewProjection );
872 :
873 46 : bHeaderDirty = TRUE;
874 :
875 46 : return CE_None;
876 : }
877 :
878 : /************************************************************************/
879 : /* GetGeoTransform() */
880 : /************************************************************************/
881 :
882 36 : CPLErr ENVIDataset::GetGeoTransform( double * padfTransform )
883 :
884 : {
885 36 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
886 :
887 36 : if( bFoundMapinfo )
888 36 : return CE_None;
889 : else
890 0 : return CE_Failure;
891 : }
892 :
893 : /************************************************************************/
894 : /* SetGeoTransform() */
895 : /************************************************************************/
896 :
897 47 : CPLErr ENVIDataset::SetGeoTransform( double * padfTransform )
898 : {
899 47 : memcpy( adfGeoTransform, padfTransform, sizeof(double) * 6 );
900 :
901 47 : bHeaderDirty = TRUE;
902 47 : bFoundMapinfo = TRUE;
903 :
904 47 : return CE_None;
905 : }
906 :
907 : /************************************************************************/
908 : /* SplitList() */
909 : /* */
910 : /* Split an ENVI value list into component fields, and strip */
911 : /* white space. */
912 : /************************************************************************/
913 :
914 175 : char **ENVIDataset::SplitList( const char *pszCleanInput )
915 :
916 : {
917 175 : char **papszReturn = NULL;
918 175 : char *pszInput = CPLStrdup(pszCleanInput);
919 :
920 175 : if( pszInput[0] != '{' )
921 : {
922 52 : CPLFree(pszInput);
923 52 : return NULL;
924 : }
925 :
926 123 : int iChar=1;
927 :
928 :
929 925 : while( pszInput[iChar] != '}' && pszInput[iChar] != '\0' )
930 : {
931 679 : int iFStart=-1, iFEnd=-1;
932 :
933 : // Find start of token.
934 679 : iFStart = iChar;
935 1858 : while( pszInput[iFStart] == ' ' )
936 500 : iFStart++;
937 :
938 679 : iFEnd = iFStart;
939 15329 : while( pszInput[iFEnd] != ','
940 4739 : && pszInput[iFEnd] != '}'
941 4616 : && pszInput[iFEnd] != '\0' )
942 4616 : iFEnd++;
943 :
944 679 : if( pszInput[iFEnd] == '\0' )
945 0 : break;
946 :
947 679 : iChar = iFEnd + 1;
948 679 : iFEnd = iFEnd - 1;
949 :
950 1358 : while( iFEnd > iFStart && pszInput[iFEnd] == ' ' )
951 0 : iFEnd--;
952 :
953 679 : pszInput[iFEnd + 1] = '\0';
954 679 : papszReturn = CSLAddString( papszReturn, pszInput + iFStart );
955 : }
956 :
957 123 : CPLFree( pszInput );
958 :
959 123 : return papszReturn;
960 : }
961 :
962 : /************************************************************************/
963 : /* SetENVIDatum() */
964 : /************************************************************************/
965 :
966 0 : void ENVIDataset::SetENVIDatum( OGRSpatialReference *poSRS,
967 : const char *pszENVIDatumName )
968 :
969 : {
970 : // datums
971 0 : if( EQUAL(pszENVIDatumName, "WGS-84") )
972 0 : poSRS->SetWellKnownGeogCS( "WGS84" );
973 0 : else if( EQUAL(pszENVIDatumName, "WGS-72") )
974 0 : poSRS->SetWellKnownGeogCS( "WGS72" );
975 0 : else if( EQUAL(pszENVIDatumName, "North America 1983") )
976 0 : poSRS->SetWellKnownGeogCS( "NAD83" );
977 0 : else if( EQUAL(pszENVIDatumName, "North America 1927")
978 : || strstr(pszENVIDatumName,"NAD27")
979 : || strstr(pszENVIDatumName,"NAD-27") )
980 0 : poSRS->SetWellKnownGeogCS( "NAD27" );
981 0 : else if( EQUALN(pszENVIDatumName, "European 1950",13) )
982 0 : poSRS->SetWellKnownGeogCS( "EPSG:4230" );
983 0 : else if( EQUAL(pszENVIDatumName, "Ordnance Survey of Great Britain '36") )
984 0 : poSRS->SetWellKnownGeogCS( "EPSG:4277" );
985 0 : else if( EQUAL(pszENVIDatumName, "SAD-69/Brazil") )
986 0 : poSRS->SetWellKnownGeogCS( "EPSG:4291" );
987 0 : else if( EQUAL(pszENVIDatumName, "Geocentric Datum of Australia 1994") )
988 0 : poSRS->SetWellKnownGeogCS( "EPSG:4283" );
989 0 : else if( EQUAL(pszENVIDatumName, "Australian Geodetic 1984") )
990 0 : poSRS->SetWellKnownGeogCS( "EPSG:4203" );
991 0 : else if( EQUAL(pszENVIDatumName, "Nouvelle Triangulation Francaise IGN") )
992 0 : poSRS->SetWellKnownGeogCS( "EPSG:4275" );
993 :
994 : // Ellipsoids
995 0 : else if( EQUAL(pszENVIDatumName, "GRS 80") )
996 0 : poSRS->SetWellKnownGeogCS( "NAD83" );
997 0 : else if( EQUAL(pszENVIDatumName, "Airy") )
998 0 : poSRS->SetWellKnownGeogCS( "EPSG:4001" );
999 0 : else if( EQUAL(pszENVIDatumName, "Australian National") )
1000 0 : poSRS->SetWellKnownGeogCS( "EPSG:4003" );
1001 0 : else if( EQUAL(pszENVIDatumName, "Bessel 1841") )
1002 0 : poSRS->SetWellKnownGeogCS( "EPSG:4004" );
1003 0 : else if( EQUAL(pszENVIDatumName, "Clark 1866") )
1004 0 : poSRS->SetWellKnownGeogCS( "EPSG:4008" );
1005 : else
1006 : {
1007 : CPLError( CE_Warning, CPLE_AppDefined,
1008 0 : "Unrecognised datum '%s', defaulting to WGS84.", pszENVIDatumName);
1009 0 : poSRS->SetWellKnownGeogCS( "WGS84" );
1010 : }
1011 0 : }
1012 :
1013 : /************************************************************************/
1014 : /* SetENVIEllipse() */
1015 : /************************************************************************/
1016 :
1017 7 : void ENVIDataset::SetENVIEllipse( OGRSpatialReference *poSRS,
1018 : char **papszPI_EI )
1019 :
1020 : {
1021 7 : double dfA = CPLAtofM(papszPI_EI[0]);
1022 7 : double dfB = CPLAtofM(papszPI_EI[1]);
1023 : double dfInvF;
1024 :
1025 7 : if( fabs(dfA-dfB) < 0.1 )
1026 0 : dfInvF = 0.0; // sphere
1027 : else
1028 7 : dfInvF = dfA / (dfA - dfB);
1029 :
1030 :
1031 : poSRS->SetGeogCS( "Ellipse Based", "Ellipse Based", "Unnamed",
1032 7 : dfA, dfInvF );
1033 7 : }
1034 :
1035 : /************************************************************************/
1036 : /* ProcessMapinfo() */
1037 : /* */
1038 : /* Extract projection, and geotransform from a mapinfo value in */
1039 : /* the header. */
1040 : /************************************************************************/
1041 :
1042 52 : int ENVIDataset::ProcessMapinfo( const char *pszMapinfo )
1043 :
1044 : {
1045 : char **papszFields;
1046 : int nCount;
1047 52 : OGRSpatialReference oSRS;
1048 :
1049 52 : papszFields = SplitList( pszMapinfo );
1050 52 : nCount = CSLCount(papszFields);
1051 :
1052 52 : if( nCount < 7 )
1053 : {
1054 0 : CSLDestroy( papszFields );
1055 0 : return FALSE;
1056 : }
1057 :
1058 : /* -------------------------------------------------------------------- */
1059 : /* Check if we have coordinate system string, and if so parse it. */
1060 : /* -------------------------------------------------------------------- */
1061 52 : char **papszCSS = NULL;
1062 52 : if( CSLFetchNameValue( papszHeader, "coordinate_system_string" ) != NULL )
1063 : {
1064 : papszCSS = CSLTokenizeString2(
1065 : CSLFetchNameValue( papszHeader, "coordinate_system_string" ),
1066 37 : "{}", CSLT_PRESERVEQUOTES );
1067 : }
1068 :
1069 : /* -------------------------------------------------------------------- */
1070 : /* Check if we have projection info, and if so parse it. */
1071 : /* -------------------------------------------------------------------- */
1072 52 : char **papszPI = NULL;
1073 52 : int nPICount = 0;
1074 52 : if( CSLFetchNameValue( papszHeader, "projection_info" ) != NULL )
1075 : {
1076 : papszPI = SplitList(
1077 19 : CSLFetchNameValue( papszHeader, "projection_info" ) );
1078 19 : nPICount = CSLCount(papszPI);
1079 : }
1080 :
1081 : /* -------------------------------------------------------------------- */
1082 : /* Capture geotransform. */
1083 : /* -------------------------------------------------------------------- */
1084 52 : adfGeoTransform[1] = atof(papszFields[5]); // Pixel width
1085 52 : adfGeoTransform[5] = -atof(papszFields[6]); // Pixel height
1086 : adfGeoTransform[0] = // Upper left X coordinate
1087 52 : atof(papszFields[3]) - (atof(papszFields[1]) - 1) * adfGeoTransform[1];
1088 : adfGeoTransform[3] = // Upper left Y coordinate
1089 52 : atof(papszFields[4]) - (atof(papszFields[2]) - 1) * adfGeoTransform[5];
1090 52 : adfGeoTransform[2] = 0.0;
1091 52 : adfGeoTransform[4] = 0.0;
1092 :
1093 : /* -------------------------------------------------------------------- */
1094 : /* Capture projection. */
1095 : /* -------------------------------------------------------------------- */
1096 52 : if ( oSRS.importFromESRI( papszCSS ) != OGRERR_NONE )
1097 : {
1098 15 : oSRS.Clear();
1099 :
1100 22 : if( EQUALN(papszFields[0],"UTM",3) && nCount >= 9 )
1101 : {
1102 7 : oSRS.SetUTM( atoi(papszFields[7]),
1103 14 : !EQUAL(papszFields[8],"South") );
1104 7 : if( nCount >= 10 && strstr(papszFields[9],"=") == NULL )
1105 0 : SetENVIDatum( &oSRS, papszFields[9] );
1106 : else
1107 7 : oSRS.SetWellKnownGeogCS( "NAD27" );
1108 : }
1109 8 : else if( EQUALN(papszFields[0],"State Plane (NAD 27)",19)
1110 : && nCount >= 7 )
1111 : {
1112 0 : oSRS.SetStatePlane( ITTVISToUSGSZone(atoi(papszFields[7])), FALSE );
1113 : }
1114 8 : else if( EQUALN(papszFields[0],"State Plane (NAD 83)",19)
1115 : && nCount >= 7 )
1116 : {
1117 0 : oSRS.SetStatePlane( ITTVISToUSGSZone(atoi(papszFields[7])), TRUE );
1118 : }
1119 8 : else if( EQUALN(papszFields[0],"Geographic Lat",14)
1120 : && nCount >= 8 )
1121 : {
1122 0 : if( nCount >= 8 && strstr(papszFields[7],"=") == NULL )
1123 0 : SetENVIDatum( &oSRS, papszFields[7] );
1124 : else
1125 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
1126 : }
1127 8 : else if( nPICount > 8 && atoi(papszPI[0]) == 3 ) // TM
1128 : {
1129 0 : oSRS.SetTM( CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1130 0 : CPLAtofM(papszPI[7]),
1131 0 : CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) );
1132 : }
1133 8 : else if( nPICount > 8 && atoi(papszPI[0]) == 4 ) // Lambert Conformal Conic
1134 : {
1135 0 : oSRS.SetLCC( CPLAtofM(papszPI[7]), CPLAtofM(papszPI[8]),
1136 0 : CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1137 0 : CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) );
1138 : }
1139 8 : else if( nPICount > 10 && atoi(papszPI[0]) == 5 ) // Oblique Merc (2 point)
1140 : {
1141 0 : oSRS.SetHOM2PNO( CPLAtofM(papszPI[3]),
1142 0 : CPLAtofM(papszPI[4]), CPLAtofM(papszPI[5]),
1143 0 : CPLAtofM(papszPI[6]), CPLAtofM(papszPI[7]),
1144 0 : CPLAtofM(papszPI[10]),
1145 0 : CPLAtofM(papszPI[8]), CPLAtofM(papszPI[9]) );
1146 : }
1147 8 : else if( nPICount > 8 && atoi(papszPI[0]) == 6 ) // Oblique Merc
1148 : {
1149 0 : oSRS.SetHOM(CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1150 0 : CPLAtofM(papszPI[5]), 0.0,
1151 0 : CPLAtofM(papszPI[8]),
1152 0 : CPLAtofM(papszPI[6]), CPLAtofM(papszPI[7]) );
1153 : }
1154 8 : else if( nPICount > 8 && atoi(papszPI[0]) == 7 ) // Stereographic
1155 : {
1156 0 : oSRS.SetStereographic( CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1157 0 : CPLAtofM(papszPI[7]),
1158 0 : CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) );
1159 : }
1160 15 : else if( nPICount > 8 && atoi(papszPI[0]) == 9 ) // Albers Equal Area
1161 : {
1162 14 : oSRS.SetACEA( CPLAtofM(papszPI[7]), CPLAtofM(papszPI[8]),
1163 14 : CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1164 35 : CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) );
1165 : }
1166 1 : else if( nPICount > 6 && atoi(papszPI[0]) == 10 ) // Polyconic
1167 : {
1168 0 : oSRS.SetPolyconic(CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1169 0 : CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) );
1170 : }
1171 1 : else if( nPICount > 6 && atoi(papszPI[0]) == 11 ) // LAEA
1172 : {
1173 0 : oSRS.SetLAEA(CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1174 0 : CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) );
1175 : }
1176 1 : else if( nPICount > 6 && atoi(papszPI[0]) == 12 ) // Azimuthal Equid.
1177 : {
1178 0 : oSRS.SetAE(CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1179 0 : CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) );
1180 : }
1181 1 : else if( nPICount > 6 && atoi(papszPI[0]) == 31 ) // Polar Stereographic
1182 : {
1183 0 : oSRS.SetPS(CPLAtofM(papszPI[3]), CPLAtofM(papszPI[4]),
1184 : 1.0,
1185 0 : CPLAtofM(papszPI[5]), CPLAtofM(papszPI[6]) );
1186 : }
1187 : }
1188 :
1189 52 : CSLDestroy( papszCSS );
1190 :
1191 : // Still lots more that could be added for someone with the patience.
1192 :
1193 : /* -------------------------------------------------------------------- */
1194 : /* fallback to localcs if we don't recognise things. */
1195 : /* -------------------------------------------------------------------- */
1196 52 : if( oSRS.GetRoot() == NULL )
1197 1 : oSRS.SetLocalCS( papszFields[0] );
1198 :
1199 : /* -------------------------------------------------------------------- */
1200 : /* Try to set datum from projection info line if we have a */
1201 : /* projected coordinate system without a GEOGCS. */
1202 : /* -------------------------------------------------------------------- */
1203 52 : if( oSRS.IsProjected() && oSRS.GetAttrNode("GEOGCS") == NULL
1204 : && nPICount > 3 )
1205 : {
1206 : // Do we have a datum on the projection info line?
1207 7 : int iDatum = nPICount-1;
1208 :
1209 : // Ignore units= items.
1210 7 : if( strstr(papszPI[iDatum],"=") != NULL )
1211 0 : iDatum--;
1212 :
1213 : // Skip past the name.
1214 7 : iDatum--;
1215 :
1216 7 : CPLString osDatumName = papszPI[iDatum];
1217 7 : if( osDatumName.find_first_of("abcdefghijklmnopqrstuvwxyz"
1218 : "ABCDEFGHIJKLMNOPQRSTUVWXYZ")
1219 : != CPLString::npos )
1220 : {
1221 0 : SetENVIDatum( &oSRS, osDatumName );
1222 : }
1223 : else
1224 : {
1225 7 : SetENVIEllipse( &oSRS, papszPI + 1 );
1226 7 : }
1227 : }
1228 :
1229 :
1230 : /* -------------------------------------------------------------------- */
1231 : /* Try to process specialized units. */
1232 : /* -------------------------------------------------------------------- */
1233 52 : if( EQUALN( papszFields[nCount-1],"units",5))
1234 : {
1235 : /* Handle linear units first. */
1236 0 : if (EQUAL(papszFields[nCount-1],"units=Feet") )
1237 0 : oSRS.SetLinearUnitsAndUpdateParameters( SRS_UL_FOOT, atof(SRS_UL_FOOT_CONV) );
1238 0 : else if (EQUAL(papszFields[nCount-1],"units=Meters") )
1239 0 : oSRS.SetLinearUnitsAndUpdateParameters( SRS_UL_METER, 1. );
1240 0 : else if (EQUAL(papszFields[nCount-1],"units=Km") )
1241 0 : oSRS.SetLinearUnitsAndUpdateParameters( "Kilometer", 1000. );
1242 0 : else if (EQUAL(papszFields[nCount-1],"units=Yards") )
1243 0 : oSRS.SetLinearUnitsAndUpdateParameters( "Yard", .9144 );
1244 0 : else if (EQUAL(papszFields[nCount-1],"units=Miles") )
1245 0 : oSRS.SetLinearUnitsAndUpdateParameters( "Mile", 1609.344 );
1246 0 : else if (EQUAL(papszFields[nCount-1],"units=Nautical Miles") )
1247 0 : oSRS.SetLinearUnitsAndUpdateParameters( SRS_UL_NAUTICAL_MILE, atof(SRS_UL_NAUTICAL_MILE_CONV) );
1248 :
1249 : /* Only handle angular units if we know the projection is geographic. */
1250 0 : if (oSRS.IsGeographic())
1251 : {
1252 0 : if (EQUAL(papszFields[nCount-1],"units=Radians") )
1253 0 : oSRS.SetAngularUnits( SRS_UA_RADIAN, 1. );
1254 : else
1255 : {
1256 : /* Degrees, minutes and seconds will all be represented as degrees. */
1257 0 : oSRS.SetAngularUnits( SRS_UA_DEGREE, atof(SRS_UA_DEGREE_CONV));
1258 :
1259 0 : double conversionFactor = 1.;
1260 0 : if (EQUAL(papszFields[nCount-1],"units=Minutes") )
1261 0 : conversionFactor = 60.;
1262 0 : else if( EQUAL(papszFields[nCount-1],"units=Seconds") )
1263 0 : conversionFactor = 3600.;
1264 0 : adfGeoTransform[0] /= conversionFactor;
1265 0 : adfGeoTransform[1] /= conversionFactor;
1266 0 : adfGeoTransform[2] /= conversionFactor;
1267 0 : adfGeoTransform[3] /= conversionFactor;
1268 0 : adfGeoTransform[4] /= conversionFactor;
1269 0 : adfGeoTransform[5] /= conversionFactor;
1270 : }
1271 : }
1272 : }
1273 : /* -------------------------------------------------------------------- */
1274 : /* Turn back into WKT. */
1275 : /* -------------------------------------------------------------------- */
1276 52 : if( oSRS.GetRoot() != NULL )
1277 : {
1278 52 : oSRS.Fixup();
1279 52 : if ( pszProjection )
1280 : {
1281 52 : CPLFree( pszProjection );
1282 52 : pszProjection = NULL;
1283 : }
1284 52 : oSRS.exportToWkt( &pszProjection );
1285 : }
1286 :
1287 52 : CSLDestroy( papszFields );
1288 52 : CSLDestroy( papszPI );
1289 52 : return TRUE;
1290 : }
1291 :
1292 : /************************************************************************/
1293 : /* ProcessRPCinfo() */
1294 : /* */
1295 : /* Extract RPC transformation coefficients if they are present */
1296 : /* and sets into the standard metadata fields for RPC. */
1297 : /************************************************************************/
1298 :
1299 0 : void ENVIDataset::ProcessRPCinfo( const char *pszRPCinfo,
1300 : int numCols, int numRows)
1301 : {
1302 : char **papszFields;
1303 : char sVal[1280];
1304 : int nCount;
1305 :
1306 0 : papszFields = SplitList( pszRPCinfo );
1307 0 : nCount = CSLCount(papszFields);
1308 :
1309 0 : if( nCount < 90 )
1310 : {
1311 0 : CSLDestroy( papszFields );
1312 0 : return;
1313 : }
1314 :
1315 0 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[0]));
1316 0 : SetMetadataItem("LINE_OFF",sVal,"RPC");
1317 0 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[5]));
1318 0 : SetMetadataItem("LINE_SCALE",sVal,"RPC");
1319 0 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[1]));
1320 0 : SetMetadataItem("SAMP_OFF",sVal,"RPC");
1321 0 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[6]));
1322 0 : SetMetadataItem("SAMP_SCALE",sVal,"RPC");
1323 0 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[2]));
1324 0 : SetMetadataItem("LAT_OFF",sVal,"RPC");
1325 0 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[7]));
1326 0 : SetMetadataItem("LAT_SCALE",sVal,"RPC");
1327 0 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[3]));
1328 0 : SetMetadataItem("LONG_OFF",sVal,"RPC");
1329 0 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[8]));
1330 0 : SetMetadataItem("LONG_SCALE",sVal,"RPC");
1331 0 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[4]));
1332 0 : SetMetadataItem("HEIGHT_OFF",sVal,"RPC");
1333 0 : snprintf(sVal, sizeof(sVal), "%.16g",atof(papszFields[9]));
1334 0 : SetMetadataItem("HEIGHT_SCALE",sVal,"RPC");
1335 :
1336 0 : sVal[0] = '\0';
1337 : int i;
1338 0 : for(i = 0; i < 20; i++ )
1339 : snprintf(sVal+strlen(sVal), sizeof(sVal), "%.16g ",
1340 0 : atof(papszFields[10+i]));
1341 0 : SetMetadataItem("LINE_NUM_COEFF",sVal,"RPC");
1342 :
1343 0 : sVal[0] = '\0';
1344 0 : for(i = 0; i < 20; i++ )
1345 : snprintf(sVal+strlen(sVal), sizeof(sVal), "%.16g ",
1346 0 : atof(papszFields[30+i]));
1347 0 : SetMetadataItem("LINE_DEN_COEFF",sVal,"RPC");
1348 :
1349 0 : sVal[0] = '\0';
1350 0 : for(i = 0; i < 20; i++ )
1351 : snprintf(sVal+strlen(sVal), sizeof(sVal), "%.16g ",
1352 0 : atof(papszFields[50+i]));
1353 0 : SetMetadataItem("SAMP_NUM_COEFF",sVal,"RPC");
1354 :
1355 0 : sVal[0] = '\0';
1356 0 : for(i = 0; i < 20; i++ )
1357 : snprintf(sVal+strlen(sVal), sizeof(sVal), "%.16g ",
1358 0 : atof(papszFields[70+i]));
1359 0 : SetMetadataItem("SAMP_DEN_COEFF",sVal,"RPC");
1360 :
1361 : snprintf(sVal, sizeof(sVal), "%.16g",
1362 0 : atof(papszFields[3]) - atof(papszFields[8]));
1363 0 : SetMetadataItem("MIN_LONG",sVal,"RPC");
1364 :
1365 : snprintf(sVal, sizeof(sVal), "%.16g",
1366 0 : atof(papszFields[3]) + atof(papszFields[8]) );
1367 0 : SetMetadataItem("MAX_LONG",sVal,"RPC");
1368 :
1369 : snprintf(sVal, sizeof(sVal), "%.16g",
1370 0 : atof(papszFields[2]) - atof(papszFields[7]));
1371 0 : SetMetadataItem("MIN_LAT",sVal,"RPC");
1372 :
1373 : snprintf(sVal, sizeof(sVal), "%.16g",
1374 0 : atof(papszFields[2]) + atof(papszFields[7]));
1375 0 : SetMetadataItem("MAX_LAT",sVal,"RPC");
1376 :
1377 : /* Handle the chipping case where the image is a subset. */
1378 : double rowOffset, colOffset;
1379 0 : rowOffset = atof(papszFields[90]);
1380 0 : colOffset = atof(papszFields[91]);
1381 0 : if (rowOffset || colOffset)
1382 : {
1383 0 : SetMetadataItem("ICHIP_SCALE_FACTOR", "1");
1384 0 : SetMetadataItem("ICHIP_ANAMORPH_CORR", "0");
1385 0 : SetMetadataItem("ICHIP_SCANBLK_NUM", "0");
1386 :
1387 0 : SetMetadataItem("ICHIP_OP_ROW_11", "0.5");
1388 0 : SetMetadataItem("ICHIP_OP_COL_11", "0.5");
1389 0 : SetMetadataItem("ICHIP_OP_ROW_12", "0.5");
1390 0 : SetMetadataItem("ICHIP_OP_COL_21", "0.5");
1391 0 : snprintf(sVal, sizeof(sVal), "%.16g", numCols - 0.5);
1392 0 : SetMetadataItem("ICHIP_OP_COL_12", sVal);
1393 0 : SetMetadataItem("ICHIP_OP_COL_22", sVal);
1394 0 : snprintf(sVal, sizeof(sVal), "%.16g", numRows - 0.5);
1395 0 : SetMetadataItem("ICHIP_OP_ROW_21", sVal);
1396 0 : SetMetadataItem("ICHIP_OP_ROW_22", sVal);
1397 :
1398 0 : snprintf(sVal, sizeof(sVal), "%.16g", rowOffset + 0.5);
1399 0 : SetMetadataItem("ICHIP_FI_ROW_11", sVal);
1400 0 : SetMetadataItem("ICHIP_FI_ROW_12", sVal);
1401 0 : snprintf(sVal, sizeof(sVal), "%.16g", colOffset + 0.5);
1402 0 : SetMetadataItem("ICHIP_FI_COL_11", sVal);
1403 0 : SetMetadataItem("ICHIP_FI_COL_21", sVal);
1404 0 : snprintf(sVal, sizeof(sVal), "%.16g", colOffset + numCols - 0.5);
1405 0 : SetMetadataItem("ICHIP_FI_COL_12", sVal);
1406 0 : SetMetadataItem("ICHIP_FI_COL_22", sVal);
1407 0 : snprintf(sVal, sizeof(sVal), "%.16g", rowOffset + numRows - 0.5);
1408 0 : SetMetadataItem("ICHIP_FI_ROW_21", sVal);
1409 0 : SetMetadataItem("ICHIP_FI_ROW_22", sVal);
1410 : }
1411 0 : CSLDestroy(papszFields);
1412 : }
1413 :
1414 102 : void ENVIDataset::ProcessStatsFile()
1415 : {
1416 102 : CPLString osStaFilename;
1417 : VSILFILE *fpStaFile;
1418 :
1419 102 : osStaFilename = CPLResetExtension( pszHDRFilename, "sta" );
1420 102 : fpStaFile = VSIFOpenL( osStaFilename, "rb" );
1421 :
1422 102 : if (!fpStaFile)
1423 : return;
1424 :
1425 : long lTestHeader[10],lOffset;
1426 :
1427 0 : if( VSIFReadL( lTestHeader, sizeof(long), 10, fpStaFile ) != 10 )
1428 : {
1429 0 : VSIFCloseL( fpStaFile );
1430 : return;
1431 : }
1432 :
1433 : int isFloat;
1434 0 : isFloat = (byteSwapLong(lTestHeader[0]) == 1111838282);
1435 :
1436 : unsigned long nb,i;
1437 : float * fStats;
1438 : double * dStats, dMin, dMax, dMean, dStd;
1439 :
1440 0 : nb=byteSwapLong(lTestHeader[3]);
1441 :
1442 0 : if (nb > (unsigned long)nBands)
1443 : {
1444 : CPLDebug("ENVI", ".sta file has statistics for %ld bands, "
1445 0 : "whereas the dataset has only %d bands", nb, nBands);
1446 0 : nb = nBands;
1447 : }
1448 :
1449 0 : VSIFSeekL(fpStaFile,40+(nb+1)*4,SEEK_SET);
1450 :
1451 0 : if (VSIFReadL(&lOffset,sizeof(long),1,fpStaFile) == 1)
1452 : {
1453 0 : VSIFSeekL(fpStaFile,40+(nb+1)*8+byteSwapLong(lOffset)+nb,SEEK_SET);
1454 : // This should be the beginning of the statistics
1455 0 : if (isFloat)
1456 : {
1457 0 : fStats = (float*)CPLCalloc(nb*4,4);
1458 0 : if (VSIFReadL(fStats,4,nb*4,fpStaFile) == nb*4)
1459 : {
1460 0 : for (i=0;i<nb;i++)
1461 : {
1462 : GetRasterBand(i+1)->SetStatistics(
1463 : byteSwapFloat(fStats[i]),
1464 : byteSwapFloat(fStats[nb+i]),
1465 : byteSwapFloat(fStats[2*nb+i]),
1466 0 : byteSwapFloat(fStats[3*nb+i]));
1467 : }
1468 : }
1469 0 : CPLFree(fStats);
1470 : }
1471 : else
1472 : {
1473 0 : dStats = (double*)CPLCalloc(nb*4,8);
1474 0 : if (VSIFReadL(dStats,8,nb*4,fpStaFile) == nb*4)
1475 : {
1476 0 : for (i=0;i<nb;i++)
1477 : {
1478 0 : dMin = byteSwapDouble(dStats[i]);
1479 0 : dMax = byteSwapDouble(dStats[nb+i]);
1480 0 : dMean = byteSwapDouble(dStats[2*nb+i]);
1481 0 : dStd = byteSwapDouble(dStats[3*nb+i]);
1482 0 : if (dMin != dMax && dStd != 0)
1483 0 : GetRasterBand(i+1)->SetStatistics(dMin,dMax,dMean,dStd);
1484 : }
1485 : }
1486 0 : CPLFree(dStats);
1487 : }
1488 : }
1489 0 : VSIFCloseL( fpStaFile );
1490 : }
1491 0 : long ENVIDataset::byteSwapLong(long swapMe)
1492 : {
1493 : #ifndef CPL_LSB
1494 : return swapMe;
1495 : #endif
1496 : return ((swapMe&0xFF)<<24)+((swapMe&(0xFF00))<<8)+
1497 0 : ((swapMe&(0xFF0000))>>8)+((swapMe&(0xFF000000))>>24);
1498 : }
1499 :
1500 0 : float ENVIDataset::byteSwapFloat(float swapMe)
1501 : {
1502 : #ifndef CPL_LSB
1503 : return swapMe;
1504 : #endif
1505 : long foo,bar;
1506 : float retval;
1507 0 : memmove(&foo, &swapMe,4);
1508 0 : bar = byteSwapLong(foo);
1509 0 : memmove(&retval, &bar,4);
1510 :
1511 0 : return retval;
1512 : }
1513 :
1514 0 : double ENVIDataset::byteSwapDouble(double swapMe)
1515 : {
1516 : #ifndef CPL_LSB
1517 : return swapMe;
1518 : #endif
1519 :
1520 : char dblArrSrc[8],dblArrDest[8];
1521 : double retval;
1522 : int i;
1523 :
1524 0 : memmove(dblArrSrc, &swapMe, 8);
1525 :
1526 0 : for (i=0;i<8;i++)
1527 0 : dblArrDest[i]=dblArrSrc[7-i];
1528 :
1529 0 : memmove(&retval,dblArrDest, 8);
1530 0 : return retval;
1531 : }
1532 :
1533 :
1534 : /************************************************************************/
1535 : /* ReadHeader() */
1536 : /************************************************************************/
1537 :
1538 104 : int ENVIDataset::ReadHeader( VSILFILE * fpHdr )
1539 :
1540 : {
1541 :
1542 104 : CPLReadLineL( fpHdr );
1543 :
1544 : /* -------------------------------------------------------------------- */
1545 : /* Now start forming sets of name/value pairs. */
1546 : /* -------------------------------------------------------------------- */
1547 1045 : while( TRUE )
1548 : {
1549 : const char *pszNewLine;
1550 : char *pszWorkingLine;
1551 :
1552 1149 : pszNewLine = CPLReadLineL( fpHdr );
1553 1149 : if( pszNewLine == NULL )
1554 : break;
1555 :
1556 1045 : if( strstr(pszNewLine,"=") == NULL )
1557 0 : continue;
1558 :
1559 1045 : pszWorkingLine = CPLStrdup(pszNewLine);
1560 :
1561 : // Collect additional lines if we have open sqiggly bracket.
1562 1045 : if( strstr(pszWorkingLine,"{") != NULL
1563 : && strstr(pszWorkingLine,"}") == NULL )
1564 : {
1565 127 : do {
1566 127 : pszNewLine = CPLReadLineL( fpHdr );
1567 127 : if( pszNewLine )
1568 : {
1569 : pszWorkingLine = (char *)
1570 : CPLRealloc(pszWorkingLine,
1571 127 : strlen(pszWorkingLine)+strlen(pszNewLine)+1);
1572 127 : strcat( pszWorkingLine, pszNewLine );
1573 : }
1574 : } while( pszNewLine != NULL && strstr(pszNewLine,"}") == NULL );
1575 : }
1576 :
1577 : // Try to break input into name and value portions. Trim whitespace.
1578 : const char *pszValue;
1579 : int iEqual;
1580 :
1581 25595 : for( iEqual = 0;
1582 24550 : pszWorkingLine[iEqual] != '\0' && pszWorkingLine[iEqual] != '=';
1583 : iEqual++ ) {}
1584 :
1585 1045 : if( pszWorkingLine[iEqual] == '=' )
1586 : {
1587 : int i;
1588 :
1589 1045 : pszValue = pszWorkingLine + iEqual + 1;
1590 3135 : while( *pszValue == ' ' || *pszValue == '\t' )
1591 1045 : pszValue++;
1592 :
1593 1045 : pszWorkingLine[iEqual--] = '\0';
1594 7102 : while( iEqual > 0
1595 2506 : && (pszWorkingLine[iEqual] == ' '
1596 1045 : || pszWorkingLine[iEqual] == '\t') )
1597 1461 : pszWorkingLine[iEqual--] = '\0';
1598 :
1599 : // Convert spaces in the name to underscores.
1600 10814 : for( i = 0; pszWorkingLine[i] != '\0'; i++ )
1601 : {
1602 9769 : if( pszWorkingLine[i] == ' ' )
1603 621 : pszWorkingLine[i] = '_';
1604 : }
1605 :
1606 : papszHeader = CSLSetNameValue( papszHeader,
1607 1045 : pszWorkingLine, pszValue );
1608 : }
1609 :
1610 1045 : CPLFree( pszWorkingLine );
1611 : }
1612 :
1613 104 : return TRUE;
1614 : }
1615 :
1616 : /************************************************************************/
1617 : /* Open() */
1618 : /************************************************************************/
1619 :
1620 11158 : GDALDataset *ENVIDataset::Open( GDALOpenInfo * poOpenInfo )
1621 :
1622 : {
1623 : int i;
1624 :
1625 : /* -------------------------------------------------------------------- */
1626 : /* We assume the user is pointing to the binary (ie. .bil) file. */
1627 : /* -------------------------------------------------------------------- */
1628 11158 : if( poOpenInfo->nHeaderBytes < 2 )
1629 10405 : return NULL;
1630 :
1631 : /* -------------------------------------------------------------------- */
1632 : /* Do we have a .hdr file? Try upper and lower case, and */
1633 : /* replacing the extension as well as appending the extension */
1634 : /* to whatever we currently have. */
1635 : /* -------------------------------------------------------------------- */
1636 : const char *pszMode;
1637 753 : CPLString osHdrFilename;
1638 753 : VSILFILE *fpHeader = NULL;
1639 :
1640 753 : if( poOpenInfo->eAccess == GA_Update )
1641 209 : pszMode = "r+";
1642 : else
1643 544 : pszMode = "r";
1644 :
1645 753 : if (poOpenInfo->papszSiblingFiles == NULL)
1646 : {
1647 4 : osHdrFilename = CPLResetExtension( poOpenInfo->pszFilename, "hdr" );
1648 4 : fpHeader = VSIFOpenL( osHdrFilename, pszMode );
1649 :
1650 4 : if( fpHeader == NULL && VSIIsCaseSensitiveFS(osHdrFilename) )
1651 : {
1652 4 : osHdrFilename = CPLResetExtension( poOpenInfo->pszFilename, "HDR" );
1653 4 : fpHeader = VSIFOpenL( osHdrFilename, pszMode );
1654 : }
1655 :
1656 4 : if( fpHeader == NULL && VSIIsCaseSensitiveFS(osHdrFilename) )
1657 : {
1658 : osHdrFilename = CPLFormFilename( NULL, poOpenInfo->pszFilename,
1659 4 : "hdr" );
1660 4 : fpHeader = VSIFOpenL( osHdrFilename, pszMode );
1661 : }
1662 :
1663 4 : if( fpHeader == NULL && VSIIsCaseSensitiveFS(osHdrFilename) )
1664 : {
1665 : osHdrFilename = CPLFormFilename( NULL, poOpenInfo->pszFilename,
1666 4 : "HDR" );
1667 4 : fpHeader = VSIFOpenL( osHdrFilename, pszMode );
1668 : }
1669 :
1670 : }
1671 : else
1672 : {
1673 : /* -------------------------------------------------------------------- */
1674 : /* Now we need to tear apart the filename to form a .HDR */
1675 : /* filename. */
1676 : /* -------------------------------------------------------------------- */
1677 749 : CPLString osPath = CPLGetPath( poOpenInfo->pszFilename );
1678 749 : CPLString osName = CPLGetFilename( poOpenInfo->pszFilename );
1679 :
1680 : int iFile = CSLFindString(poOpenInfo->papszSiblingFiles,
1681 749 : CPLResetExtension( osName, "hdr" ) );
1682 749 : if( iFile >= 0 )
1683 : {
1684 211 : osHdrFilename = CPLFormFilename( osPath, poOpenInfo->papszSiblingFiles[iFile],
1685 422 : NULL );
1686 211 : fpHeader = VSIFOpenL( osHdrFilename, pszMode );
1687 : }
1688 : else
1689 : {
1690 : iFile = CSLFindString(poOpenInfo->papszSiblingFiles,
1691 538 : CPLFormFilename( NULL, osName, "hdr" ));
1692 538 : if( iFile >= 0 )
1693 : {
1694 0 : osHdrFilename = CPLFormFilename( osPath, poOpenInfo->papszSiblingFiles[iFile],
1695 0 : NULL );
1696 0 : fpHeader = VSIFOpenL( osHdrFilename, pszMode );
1697 : }
1698 749 : }
1699 : }
1700 :
1701 753 : if( fpHeader == NULL )
1702 542 : return NULL;
1703 :
1704 : /* -------------------------------------------------------------------- */
1705 : /* Check that the first line says "ENVI". */
1706 : /* -------------------------------------------------------------------- */
1707 : char szTestHdr[4];
1708 :
1709 211 : if( VSIFReadL( szTestHdr, 4, 1, fpHeader ) != 1 )
1710 : {
1711 0 : VSIFCloseL( fpHeader );
1712 0 : return NULL;
1713 : }
1714 211 : if( strncmp(szTestHdr,"ENVI",4) != 0 )
1715 : {
1716 107 : VSIFCloseL( fpHeader );
1717 107 : return NULL;
1718 : }
1719 :
1720 : /* -------------------------------------------------------------------- */
1721 : /* Create a corresponding GDALDataset. */
1722 : /* -------------------------------------------------------------------- */
1723 : ENVIDataset *poDS;
1724 :
1725 104 : poDS = new ENVIDataset();
1726 208 : poDS->pszHDRFilename = CPLStrdup(osHdrFilename);
1727 104 : poDS->fp = fpHeader;
1728 :
1729 : /* -------------------------------------------------------------------- */
1730 : /* Read the header. */
1731 : /* -------------------------------------------------------------------- */
1732 104 : if( !poDS->ReadHeader( fpHeader ) )
1733 : {
1734 0 : delete poDS;
1735 0 : return NULL;
1736 : }
1737 :
1738 : /* -------------------------------------------------------------------- */
1739 : /* Has the user selected the .hdr file to open? */
1740 : /* -------------------------------------------------------------------- */
1741 104 : if( EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "hdr") )
1742 : {
1743 0 : delete poDS;
1744 : CPLError( CE_Failure, CPLE_AppDefined,
1745 : "The selected file is an ENVI header file, but to\n"
1746 : "open ENVI datasets, the data file should be selected\n"
1747 : "instead of the .hdr file. Please try again selecting\n"
1748 : "the data file corresponding to the header file:\n"
1749 : " %s\n",
1750 0 : poOpenInfo->pszFilename );
1751 0 : return NULL;
1752 : }
1753 :
1754 : /* -------------------------------------------------------------------- */
1755 : /* Has the user selected the .sta (stats) file to open? */
1756 : /* -------------------------------------------------------------------- */
1757 104 : if( EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "sta") )
1758 : {
1759 0 : delete poDS;
1760 : CPLError( CE_Failure, CPLE_AppDefined,
1761 : "The selected file is an ENVI statistics file.\n"
1762 : "To open ENVI datasets, the data file should be selected\n"
1763 : "instead of the .sta file. Please try again selecting\n"
1764 : "the data file corresponding to the statistics file:\n"
1765 : " %s\n",
1766 0 : poOpenInfo->pszFilename );
1767 0 : return NULL;
1768 : }
1769 :
1770 : /* -------------------------------------------------------------------- */
1771 : /* Extract required values from the .hdr. */
1772 : /* -------------------------------------------------------------------- */
1773 104 : int nLines = 0, nSamples = 0, nBands = 0, nHeaderSize = 0;
1774 104 : const char *pszInterleave = NULL;
1775 :
1776 104 : if( CSLFetchNameValue(poDS->papszHeader,"lines") )
1777 104 : nLines = atoi(CSLFetchNameValue(poDS->papszHeader,"lines"));
1778 :
1779 104 : if( CSLFetchNameValue(poDS->papszHeader,"samples") )
1780 104 : nSamples = atoi(CSLFetchNameValue(poDS->papszHeader,"samples"));
1781 :
1782 104 : if( CSLFetchNameValue(poDS->papszHeader,"bands") )
1783 104 : nBands = atoi(CSLFetchNameValue(poDS->papszHeader,"bands"));
1784 :
1785 104 : pszInterleave = CSLFetchNameValue(poDS->papszHeader,"interleave");
1786 :
1787 : /* In case, there is no interleave keyword, we try to derive it from the */
1788 : /* file extension, if it matches one of the expected interleaving mode */
1789 104 : if( pszInterleave == NULL )
1790 : {
1791 0 : const char* pszExtension = CPLGetExtension(poOpenInfo->pszFilename);
1792 0 : if ( EQUAL(pszExtension, "BSQ") ||
1793 : EQUAL(pszExtension, "BIP") ||
1794 : EQUAL(pszExtension, "BIL") )
1795 : {
1796 0 : pszInterleave = pszExtension;
1797 : }
1798 : }
1799 :
1800 104 : if (!GDALCheckDatasetDimensions(nSamples, nLines) || !GDALCheckBandCount(nBands, FALSE) ||
1801 : pszInterleave == NULL )
1802 : {
1803 2 : delete poDS;
1804 : CPLError( CE_Failure, CPLE_AppDefined,
1805 : "The file appears to have an associated ENVI header, but\n"
1806 : "one or more of the samples, lines, bands and interleave\n"
1807 2 : "keywords appears to be missing or invalid." );
1808 2 : return NULL;
1809 : }
1810 :
1811 102 : if( CSLFetchNameValue(poDS->papszHeader,"header_offset") )
1812 102 : nHeaderSize = atoi(CSLFetchNameValue(poDS->papszHeader,"header_offset"));
1813 :
1814 : /* -------------------------------------------------------------------- */
1815 : /* Translate the datatype. */
1816 : /* -------------------------------------------------------------------- */
1817 102 : GDALDataType eType = GDT_Byte;
1818 :
1819 102 : if( CSLFetchNameValue(poDS->papszHeader,"data_type" ) != NULL )
1820 : {
1821 102 : switch( atoi(CSLFetchNameValue(poDS->papszHeader,"data_type" )) )
1822 : {
1823 : case 1:
1824 56 : eType = GDT_Byte;
1825 56 : break;
1826 :
1827 : case 2:
1828 6 : eType = GDT_Int16;
1829 6 : break;
1830 :
1831 : case 3:
1832 6 : eType = GDT_Int32;
1833 6 : break;
1834 :
1835 : case 4:
1836 6 : eType = GDT_Float32;
1837 6 : break;
1838 :
1839 : case 5:
1840 6 : eType = GDT_Float64;
1841 6 : break;
1842 :
1843 : case 6:
1844 5 : eType = GDT_CFloat32;
1845 5 : break;
1846 :
1847 : case 9:
1848 5 : eType = GDT_CFloat64;
1849 5 : break;
1850 :
1851 : case 12:
1852 6 : eType = GDT_UInt16;
1853 6 : break;
1854 :
1855 : case 13:
1856 6 : eType = GDT_UInt32;
1857 6 : break;
1858 :
1859 : /* 14=Int64, 15=UInt64 */
1860 :
1861 : default:
1862 0 : delete poDS;
1863 : CPLError( CE_Failure, CPLE_AppDefined,
1864 : "The file does not have a value for the data_type\n"
1865 0 : "that is recognised by the GDAL ENVI driver.");
1866 0 : return NULL;
1867 : }
1868 : }
1869 :
1870 : /* -------------------------------------------------------------------- */
1871 : /* Translate the byte order. */
1872 : /* -------------------------------------------------------------------- */
1873 102 : int bNativeOrder = TRUE;
1874 :
1875 102 : if( CSLFetchNameValue(poDS->papszHeader,"byte_order" ) != NULL )
1876 : {
1877 : #ifdef CPL_LSB
1878 : bNativeOrder = atoi(CSLFetchNameValue(poDS->papszHeader,
1879 102 : "byte_order" )) == 0;
1880 : #else
1881 : bNativeOrder = atoi(CSLFetchNameValue(poDS->papszHeader,
1882 : "byte_order" )) != 0;
1883 : #endif
1884 : }
1885 :
1886 : /* -------------------------------------------------------------------- */
1887 : /* Warn about unsupported file types virtual mosaic and meta file.*/
1888 : /* -------------------------------------------------------------------- */
1889 102 : if( CSLFetchNameValue(poDS->papszHeader,"file_type" ) != NULL )
1890 : {
1891 : // when the file type is one of these we return an invalid file type err
1892 : //'envi meta file'
1893 : //'envi virtual mosaic'
1894 : //'envi spectral library'
1895 : //'envi fft result'
1896 :
1897 : // when the file type is one of these we open it
1898 : //'envi standard'
1899 : //'envi classification'
1900 :
1901 : // when the file type is anything else we attempt to open it as a raster.
1902 :
1903 : const char * pszEnviFileType;
1904 102 : pszEnviFileType = CSLFetchNameValue(poDS->papszHeader,"file_type");
1905 :
1906 : // envi gdal does not support any of these
1907 : // all others we will attempt to open
1908 102 : if(EQUAL(pszEnviFileType, "envi meta file") ||
1909 : EQUAL(pszEnviFileType, "envi virtual mosaic") ||
1910 : EQUAL(pszEnviFileType, "envi spectral library") ||
1911 : EQUAL(pszEnviFileType, "envi fft result"))
1912 : {
1913 : CPLError( CE_Failure, CPLE_OpenFailed,
1914 : "File %s contains an invalid file type in the ENVI .hdr\n"
1915 : "GDAL does not support '%s' type files.",
1916 0 : poOpenInfo->pszFilename, pszEnviFileType );
1917 0 : delete poDS;
1918 0 : return NULL;
1919 : }
1920 : }
1921 :
1922 : /* -------------------------------------------------------------------- */
1923 : /* Detect (gzipped) compressed datasets. */
1924 : /* -------------------------------------------------------------------- */
1925 102 : int bIsCompressed = FALSE;
1926 102 : if( CSLFetchNameValue(poDS->papszHeader,"file_compression" ) != NULL )
1927 : {
1928 1 : if( atoi(CSLFetchNameValue(poDS->papszHeader,"file_compression" ))
1929 : != 0 )
1930 : {
1931 1 : bIsCompressed = TRUE;
1932 : }
1933 : }
1934 :
1935 : /* -------------------------------------------------------------------- */
1936 : /* Capture some information from the file that is of interest. */
1937 : /* -------------------------------------------------------------------- */
1938 102 : poDS->nRasterXSize = nSamples;
1939 102 : poDS->nRasterYSize = nLines;
1940 102 : poDS->eAccess = poOpenInfo->eAccess;
1941 :
1942 : /* -------------------------------------------------------------------- */
1943 : /* Reopen file in update mode if necessary. */
1944 : /* -------------------------------------------------------------------- */
1945 102 : CPLString osImageFilename(poOpenInfo->pszFilename);
1946 102 : if (bIsCompressed)
1947 1 : osImageFilename = "/vsigzip/" + osImageFilename;
1948 102 : if( poOpenInfo->eAccess == GA_Update )
1949 : {
1950 48 : if (bIsCompressed)
1951 : {
1952 0 : delete poDS;
1953 : CPLError( CE_Failure, CPLE_OpenFailed,
1954 0 : "Cannot open compressed file in update mode.\n");
1955 0 : return NULL;
1956 : }
1957 48 : poDS->fpImage = VSIFOpenL( osImageFilename, "rb+" );
1958 : }
1959 : else
1960 54 : poDS->fpImage = VSIFOpenL( osImageFilename, "rb" );
1961 :
1962 102 : if( poDS->fpImage == NULL )
1963 : {
1964 0 : delete poDS;
1965 : CPLError( CE_Failure, CPLE_OpenFailed,
1966 : "Failed to re-open %s within ENVI driver.\n",
1967 0 : poOpenInfo->pszFilename );
1968 0 : return NULL;
1969 : }
1970 :
1971 : /* -------------------------------------------------------------------- */
1972 : /* Compute the line offset. */
1973 : /* -------------------------------------------------------------------- */
1974 102 : int nDataSize = GDALGetDataTypeSize(eType)/8;
1975 : int nPixelOffset, nLineOffset;
1976 : vsi_l_offset nBandOffset;
1977 102 : int bIntOverflow = FALSE;
1978 :
1979 102 : if( EQUALN(pszInterleave, "bsq", 3) )
1980 : {
1981 102 : poDS->interleave = BSQ;
1982 102 : poDS->SetMetadataItem( "INTERLEAVE", "BAND", "IMAGE_STRUCTURE" );
1983 102 : if (nSamples > INT_MAX / nDataSize) bIntOverflow = TRUE;
1984 102 : nLineOffset = nDataSize * nSamples;
1985 102 : nPixelOffset = nDataSize;
1986 102 : nBandOffset = (vsi_l_offset)nLineOffset * nLines;
1987 : }
1988 0 : else if( EQUALN(pszInterleave, "bil", 3) )
1989 : {
1990 0 : poDS->interleave = BIL;
1991 0 : poDS->SetMetadataItem( "INTERLEAVE", "LINE", "IMAGE_STRUCTURE" );
1992 0 : if (nSamples > INT_MAX / (nDataSize * nBands)) bIntOverflow = TRUE;
1993 0 : nLineOffset = nDataSize * nSamples * nBands;
1994 0 : nPixelOffset = nDataSize;
1995 0 : nBandOffset = (vsi_l_offset)nDataSize * nSamples;
1996 : }
1997 0 : else if( EQUALN(pszInterleave, "bip", 3) )
1998 : {
1999 0 : poDS->interleave = BIP;
2000 0 : poDS->SetMetadataItem( "INTERLEAVE", "PIXEL", "IMAGE_STRUCTURE" );
2001 0 : if (nSamples > INT_MAX / (nDataSize * nBands)) bIntOverflow = TRUE;
2002 0 : nLineOffset = nDataSize * nSamples * nBands;
2003 0 : nPixelOffset = nDataSize * nBands;
2004 0 : nBandOffset = nDataSize;
2005 : }
2006 : else
2007 : {
2008 0 : delete poDS;
2009 : CPLError( CE_Failure, CPLE_AppDefined,
2010 : "The interleaving type of the file (%s) is not supported.",
2011 0 : pszInterleave );
2012 0 : return NULL;
2013 : }
2014 :
2015 102 : if (bIntOverflow)
2016 : {
2017 0 : delete poDS;
2018 : CPLError( CE_Failure, CPLE_AppDefined,
2019 0 : "Int overflow occured.");
2020 0 : return NULL;
2021 : }
2022 :
2023 : /* -------------------------------------------------------------------- */
2024 : /* Create band information objects. */
2025 : /* -------------------------------------------------------------------- */
2026 102 : poDS->nBands = nBands;
2027 274 : for( i = 0; i < poDS->nBands; i++ )
2028 : {
2029 : poDS->SetBand( i + 1,
2030 : new RawRasterBand(poDS, i + 1, poDS->fpImage,
2031 : nHeaderSize + nBandOffset * i,
2032 : nPixelOffset, nLineOffset, eType,
2033 172 : bNativeOrder, TRUE) );
2034 : }
2035 :
2036 : /* -------------------------------------------------------------------- */
2037 : /* Apply band names if we have them. */
2038 : /* Use wavelength for more descriptive information if possible */
2039 : /* -------------------------------------------------------------------- */
2040 102 : if( CSLFetchNameValue( poDS->papszHeader, "band_names" ) != NULL ||
2041 : CSLFetchNameValue( poDS->papszHeader, "wavelength" ) != NULL)
2042 : {
2043 : char **papszBandNames =
2044 : poDS->SplitList( CSLFetchNameValue( poDS->papszHeader,
2045 52 : "band_names" ) );
2046 : char **papszWL =
2047 : poDS->SplitList( CSLFetchNameValue( poDS->papszHeader,
2048 52 : "wavelength" ) );
2049 :
2050 52 : const char *pszWLUnits = NULL;
2051 52 : if (papszWL)
2052 : {
2053 : /* If WL information is present, process wavelength units */
2054 : pszWLUnits = CSLFetchNameValue( poDS->papszHeader,
2055 0 : "wavelength_units" );
2056 0 : if (pszWLUnits)
2057 : {
2058 : /* Don't show unknown or index units */
2059 0 : if (EQUAL(pszWLUnits,"Unknown") ||
2060 : EQUAL(pszWLUnits,"Index") )
2061 0 : pszWLUnits=0;
2062 : }
2063 : }
2064 :
2065 134 : for( i = 0; i < nBands; i++ )
2066 : {
2067 82 : CPLString osBandId, osBandName, osWavelength;
2068 :
2069 : /* First set up the wavelength names and units if available */
2070 82 : if (papszWL && CSLCount(papszWL) > i)
2071 : {
2072 0 : osWavelength = papszWL[i];
2073 0 : if (pszWLUnits)
2074 : {
2075 0 : osWavelength += " ";
2076 0 : osWavelength += pszWLUnits;
2077 : }
2078 : }
2079 :
2080 : /* Build the final name for this band */
2081 82 : if (papszBandNames && CSLCount(papszBandNames) > i)
2082 : {
2083 82 : osBandName = papszBandNames[i];
2084 82 : if (strlen(osWavelength) > 0)
2085 : {
2086 0 : osBandName += " (";
2087 0 : osBandName += osWavelength;
2088 0 : osBandName += ")";
2089 : }
2090 : }
2091 : else /* WL but no band names */
2092 0 : osBandName = osWavelength;
2093 :
2094 : /* Description is for internal GDAL usage */
2095 82 : poDS->GetRasterBand(i + 1)->SetDescription( osBandName );
2096 :
2097 : /* Metadata field named Band_1, etc. needed for ArcGIS integration */
2098 82 : osBandId = CPLSPrintf("Band_%i", i+1);
2099 82 : poDS->SetMetadataItem(osBandId, osBandName);
2100 : }
2101 52 : CSLDestroy( papszWL );
2102 52 : CSLDestroy( papszBandNames );
2103 : }
2104 : /* -------------------------------------------------------------------- */
2105 : /* Apply class names if we have them. */
2106 : /* -------------------------------------------------------------------- */
2107 102 : if( CSLFetchNameValue( poDS->papszHeader, "class_names" ) != NULL )
2108 : {
2109 : char **papszClassNames =
2110 : poDS->SplitList( CSLFetchNameValue( poDS->papszHeader,
2111 0 : "class_names" ) );
2112 :
2113 0 : poDS->GetRasterBand(1)->SetCategoryNames( papszClassNames );
2114 0 : CSLDestroy( papszClassNames );
2115 : }
2116 :
2117 : /* -------------------------------------------------------------------- */
2118 : /* Apply colormap if we have one. */
2119 : /* -------------------------------------------------------------------- */
2120 102 : if( CSLFetchNameValue( poDS->papszHeader, "class_lookup" ) != NULL )
2121 : {
2122 : char **papszClassColors =
2123 : poDS->SplitList( CSLFetchNameValue( poDS->papszHeader,
2124 0 : "class_lookup" ) );
2125 0 : int nColorValueCount = CSLCount(papszClassColors);
2126 0 : GDALColorTable oCT;
2127 :
2128 0 : for( i = 0; i*3+2 < nColorValueCount; i++ )
2129 : {
2130 : GDALColorEntry sEntry;
2131 :
2132 0 : sEntry.c1 = (short) atoi(papszClassColors[i*3+0]);
2133 0 : sEntry.c2 = (short) atoi(papszClassColors[i*3+1]);
2134 0 : sEntry.c3 = (short) atoi(papszClassColors[i*3+2]);
2135 0 : sEntry.c4 = 255;
2136 0 : oCT.SetColorEntry( i, &sEntry );
2137 : }
2138 :
2139 0 : CSLDestroy( papszClassColors );
2140 :
2141 0 : poDS->GetRasterBand(1)->SetColorTable( &oCT );
2142 0 : poDS->GetRasterBand(1)->SetColorInterpretation( GCI_PaletteIndex );
2143 : }
2144 :
2145 : /* -------------------------------------------------------------------- */
2146 : /* Set the nodata value if it is present */
2147 : /* -------------------------------------------------------------------- */
2148 102 : if( CSLFetchNameValue(poDS->papszHeader,"data_ignore_value" ) != NULL )
2149 : {
2150 0 : for( i = 0; i < poDS->nBands; i++ )
2151 : ((RawRasterBand*)poDS->GetRasterBand(i+1))->SetNoDataValue(atof(
2152 0 : CSLFetchNameValue(poDS->papszHeader,"data_ignore_value")));
2153 : }
2154 :
2155 : /* -------------------------------------------------------------------- */
2156 : /* Read the stats file if it is present */
2157 : /* -------------------------------------------------------------------- */
2158 102 : poDS->ProcessStatsFile();
2159 :
2160 : /* -------------------------------------------------------------------- */
2161 : /* Look for mapinfo */
2162 : /* -------------------------------------------------------------------- */
2163 102 : if( CSLFetchNameValue( poDS->papszHeader, "map_info" ) != NULL )
2164 : {
2165 : poDS->bFoundMapinfo =
2166 : poDS->ProcessMapinfo(
2167 52 : CSLFetchNameValue(poDS->papszHeader,"map_info") );
2168 : }
2169 :
2170 : /* -------------------------------------------------------------------- */
2171 : /* Look for RPC mapinfo */
2172 : /* -------------------------------------------------------------------- */
2173 102 : if( !poDS->bFoundMapinfo &&
2174 : CSLFetchNameValue( poDS->papszHeader, "rpc_info" ) != NULL )
2175 : {
2176 : poDS->ProcessRPCinfo(
2177 : CSLFetchNameValue(poDS->papszHeader,"rpc_info"),
2178 0 : poDS->nRasterXSize, poDS->nRasterYSize);
2179 : }
2180 :
2181 : /* -------------------------------------------------------------------- */
2182 : /* Initialize any PAM information. */
2183 : /* -------------------------------------------------------------------- */
2184 102 : poDS->SetDescription( poOpenInfo->pszFilename );
2185 102 : poDS->TryLoadXML();
2186 :
2187 : /* -------------------------------------------------------------------- */
2188 : /* Check for overviews. */
2189 : /* -------------------------------------------------------------------- */
2190 102 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
2191 :
2192 102 : return( poDS );
2193 : }
2194 :
2195 98 : int ENVIDataset::GetEnviType(GDALDataType eType)
2196 : {
2197 : int iENVIType;
2198 98 : switch( eType )
2199 : {
2200 : case GDT_Byte:
2201 44 : iENVIType = 1;
2202 44 : break;
2203 : case GDT_Int16:
2204 6 : iENVIType = 2;
2205 6 : break;
2206 : case GDT_Int32:
2207 6 : iENVIType = 3;
2208 6 : break;
2209 : case GDT_Float32:
2210 6 : iENVIType = 4;
2211 6 : break;
2212 : case GDT_Float64:
2213 6 : iENVIType = 5;
2214 6 : break;
2215 : case GDT_CFloat32:
2216 6 : iENVIType = 6;
2217 6 : break;
2218 : case GDT_CFloat64:
2219 6 : iENVIType = 9;
2220 6 : break;
2221 : case GDT_UInt16:
2222 6 : iENVIType = 12;
2223 6 : break;
2224 : case GDT_UInt32:
2225 6 : iENVIType = 13;
2226 6 : break;
2227 :
2228 : /* 14=Int64, 15=UInt64 */
2229 :
2230 : default:
2231 : CPLError( CE_Failure, CPLE_AppDefined,
2232 : "Attempt to create ENVI .hdr labelled dataset with an illegal\n"
2233 : "data type (%s).\n",
2234 6 : GDALGetDataTypeName(eType) );
2235 6 : return 1;
2236 : }
2237 92 : return iENVIType;
2238 : }
2239 :
2240 : /************************************************************************/
2241 : /* Create() */
2242 : /************************************************************************/
2243 :
2244 51 : GDALDataset *ENVIDataset::Create( const char * pszFilename,
2245 : int nXSize, int nYSize, int nBands,
2246 : GDALDataType eType,
2247 : char ** papszOptions )
2248 :
2249 : {
2250 : /* -------------------------------------------------------------------- */
2251 : /* Verify input options. */
2252 : /* -------------------------------------------------------------------- */
2253 51 : int iENVIType = GetEnviType(eType);
2254 51 : if (0 == iENVIType)
2255 0 : return 0;
2256 :
2257 : /* -------------------------------------------------------------------- */
2258 : /* Try to create the file. */
2259 : /* -------------------------------------------------------------------- */
2260 : VSILFILE *fp;
2261 :
2262 51 : fp = VSIFOpenL( pszFilename, "wb" );
2263 :
2264 51 : if( fp == NULL )
2265 : {
2266 : CPLError( CE_Failure, CPLE_OpenFailed,
2267 : "Attempt to create file `%s' failed.\n",
2268 2 : pszFilename );
2269 2 : return NULL;
2270 : }
2271 :
2272 : /* -------------------------------------------------------------------- */
2273 : /* Just write out a couple of bytes to establish the binary */
2274 : /* file, and then close it. */
2275 : /* -------------------------------------------------------------------- */
2276 49 : VSIFWriteL( (void *) "\0\0", 2, 1, fp );
2277 49 : VSIFCloseL( fp );
2278 :
2279 : /* -------------------------------------------------------------------- */
2280 : /* Create the .hdr filename. */
2281 : /* -------------------------------------------------------------------- */
2282 : const char *pszHDRFilename;
2283 : const char *pszSuffix;
2284 :
2285 49 : pszSuffix = CSLFetchNameValue( papszOptions, "SUFFIX" );
2286 49 : if ( pszSuffix && EQUALN( pszSuffix, "ADD", 3 ))
2287 0 : pszHDRFilename = CPLFormFilename( NULL, pszFilename, "hdr" );
2288 : else
2289 49 : pszHDRFilename = CPLResetExtension(pszFilename, "hdr" );
2290 :
2291 : /* -------------------------------------------------------------------- */
2292 : /* Open the file. */
2293 : /* -------------------------------------------------------------------- */
2294 49 : fp = VSIFOpenL( pszHDRFilename, "wt" );
2295 49 : if( fp == NULL )
2296 : {
2297 : CPLError( CE_Failure, CPLE_OpenFailed,
2298 : "Attempt to create file `%s' failed.\n",
2299 0 : pszHDRFilename );
2300 0 : return NULL;
2301 : }
2302 :
2303 : /* -------------------------------------------------------------------- */
2304 : /* Write out the header. */
2305 : /* -------------------------------------------------------------------- */
2306 : int iBigEndian;
2307 : const char *pszInterleaving;
2308 :
2309 : #ifdef CPL_LSB
2310 49 : iBigEndian = 0;
2311 : #else
2312 : iBigEndian = 1;
2313 : #endif
2314 :
2315 49 : VSIFPrintfL( fp, "ENVI\n" );
2316 : VSIFPrintfL( fp, "samples = %d\nlines = %d\nbands = %d\n",
2317 49 : nXSize, nYSize, nBands );
2318 49 : VSIFPrintfL( fp, "header offset = 0\nfile type = ENVI Standard\n" );
2319 49 : VSIFPrintfL( fp, "data type = %d\n", iENVIType );
2320 49 : pszInterleaving = CSLFetchNameValue( papszOptions, "INTERLEAVE" );
2321 49 : if ( pszInterleaving )
2322 : {
2323 0 : if ( EQUALN( pszInterleaving, "bip", 3 ) )
2324 0 : pszInterleaving = "bip"; // interleaved by pixel
2325 0 : else if ( EQUALN( pszInterleaving, "bil", 3 ) )
2326 0 : pszInterleaving = "bil"; // interleaved by line
2327 : else
2328 0 : pszInterleaving = "bsq"; // band sequental by default
2329 : }
2330 : else
2331 49 : pszInterleaving = "bsq";
2332 49 : VSIFPrintfL( fp, "interleave = %s\n", pszInterleaving);
2333 49 : VSIFPrintfL( fp, "byte order = %d\n", iBigEndian );
2334 :
2335 49 : VSIFCloseL( fp );
2336 :
2337 49 : return (GDALDataset *) GDALOpen( pszFilename, GA_Update );
2338 : }
2339 :
2340 : /************************************************************************/
2341 : /* GDALRegister_ENVI() */
2342 : /************************************************************************/
2343 :
2344 558 : void GDALRegister_ENVI()
2345 : {
2346 : GDALDriver *poDriver;
2347 :
2348 558 : if( GDALGetDriverByName( "ENVI" ) == NULL )
2349 : {
2350 537 : poDriver = new GDALDriver();
2351 :
2352 537 : poDriver->SetDescription( "ENVI" );
2353 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
2354 537 : "ENVI .hdr Labelled" );
2355 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
2356 537 : "frmt_various.html#ENVI" );
2357 537 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "" );
2358 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
2359 : "Byte Int16 UInt16 Int32 UInt32 "
2360 537 : "Float32 Float64 CFloat32 CFloat64" );
2361 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
2362 : "<CreationOptionList>"
2363 : " <Option name='SUFFIX' type='string-select'>"
2364 : " <Value>ADD</Value>"
2365 : " </Option>"
2366 : " <Option name='INTERLEAVE' type='string-select'>"
2367 : " <Value>BIP</Value>"
2368 : " <Value>BIL</Value>"
2369 : " <Value>BSQ</Value>"
2370 : " </Option>"
2371 537 : "</CreationOptionList>" );
2372 :
2373 537 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
2374 537 : poDriver->pfnOpen = ENVIDataset::Open;
2375 537 : poDriver->pfnCreate = ENVIDataset::Create;
2376 :
2377 537 : GetGDALDriverManager()->RegisterDriver( poDriver );
2378 : }
2379 558 : }
|
