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