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