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