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