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