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