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