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