Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: ESRI .hdr Driver
4 : * Purpose: Implementation of EHdrDataset
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 1999, Frank Warmerdam <warmerdam@pobox.com>
9 : * Copyright (c) 2007-2013, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_port.h"
15 : #include "ehdrdataset.h"
16 : #include "rawdataset.h"
17 :
18 : #include <algorithm>
19 : #include <cctype>
20 : #include <cerrno>
21 : #include <climits>
22 : #include <cmath>
23 : #include <cstddef>
24 : #include <cstdio>
25 : #include <cstdlib>
26 : #include <cstring>
27 : #if HAVE_FCNTL_H
28 : #include <fcntl.h>
29 : #endif
30 :
31 : #include <limits>
32 :
33 : #include "cpl_conv.h"
34 : #include "cpl_error.h"
35 : #include "cpl_progress.h"
36 : #include "cpl_string.h"
37 : #include "cpl_vsi.h"
38 : #include "gdal.h"
39 : #include "gdal_frmts.h"
40 : #include "gdal_pam.h"
41 : #include "gdal_priv.h"
42 : #include "ogr_core.h"
43 : #include "ogr_spatialref.h"
44 :
45 : constexpr int HAS_MIN_FLAG = 0x1;
46 : constexpr int HAS_MAX_FLAG = 0x2;
47 : constexpr int HAS_MEAN_FLAG = 0x4;
48 : constexpr int HAS_STDDEV_FLAG = 0x8;
49 : constexpr int HAS_ALL_FLAGS =
50 : HAS_MIN_FLAG | HAS_MAX_FLAG | HAS_MEAN_FLAG | HAS_STDDEV_FLAG;
51 :
52 : /************************************************************************/
53 : /* EHdrRasterBand() */
54 : /************************************************************************/
55 :
56 188 : EHdrRasterBand::EHdrRasterBand(GDALDataset *poDSIn, int nBandIn,
57 : VSILFILE *fpRawIn, vsi_l_offset nImgOffsetIn,
58 : int nPixelOffsetIn, int nLineOffsetIn,
59 : GDALDataType eDataTypeIn,
60 : RawRasterBand::ByteOrder eByteOrderIn,
61 188 : int nBitsIn)
62 : : RawRasterBand(poDSIn, nBandIn, fpRawIn, nImgOffsetIn, nPixelOffsetIn,
63 : nLineOffsetIn, eDataTypeIn, eByteOrderIn,
64 : RawRasterBand::OwnFP::NO),
65 : nBits(nBitsIn), nStartBit(0), nPixelOffsetBits(0), nLineOffsetBits(0),
66 : bNoDataSet(FALSE), dfNoData(0.0), dfMin(0.0), dfMax(0.0), dfMean(0.0),
67 188 : dfStdDev(0.0), minmaxmeanstddev(0)
68 : {
69 188 : m_bValid = RawRasterBand::IsValid();
70 :
71 188 : EHdrDataset *poEDS = reinterpret_cast<EHdrDataset *>(poDS);
72 :
73 188 : if (nBits < 8)
74 : {
75 0 : const int nSkipBytes = atoi(poEDS->GetKeyValue("SKIPBYTES"));
76 0 : if (nSkipBytes < 0 || nSkipBytes > std::numeric_limits<int>::max() / 8)
77 : {
78 0 : m_bValid = false;
79 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid SKIPBYTES: %d",
80 : nSkipBytes);
81 0 : nStartBit = 0;
82 : }
83 : else
84 : {
85 0 : nStartBit = static_cast<vsi_l_offset>(nSkipBytes) * 8;
86 : }
87 0 : if (nBand >= 2)
88 : {
89 : GIntBig nBandRowBytes =
90 0 : CPLAtoGIntBig(poEDS->GetKeyValue("BANDROWBYTES"));
91 0 : if (nBandRowBytes < 0)
92 : {
93 0 : m_bValid = false;
94 0 : CPLError(CE_Failure, CPLE_AppDefined,
95 : "Invalid BANDROWBYTES: " CPL_FRMT_GIB, nBandRowBytes);
96 0 : nBandRowBytes = 0;
97 : }
98 0 : vsi_l_offset nRowBytes = 0;
99 0 : if (nBandRowBytes == 0)
100 0 : nRowBytes =
101 0 : (static_cast<vsi_l_offset>(nBits) * poDS->GetRasterXSize() +
102 : 7) /
103 : 8;
104 : else
105 0 : nRowBytes = static_cast<vsi_l_offset>(nBandRowBytes);
106 :
107 0 : nStartBit += nRowBytes * (nBand - 1) * 8;
108 : }
109 :
110 0 : nPixelOffsetBits = nBits;
111 : GIntBig nTotalRowBytes =
112 0 : CPLAtoGIntBig(poEDS->GetKeyValue("TOTALROWBYTES"));
113 0 : if (nTotalRowBytes < 0 ||
114 0 : nTotalRowBytes > GINTBIG_MAX / 8 / poDS->GetRasterYSize())
115 : {
116 0 : m_bValid = false;
117 0 : CPLError(CE_Failure, CPLE_AppDefined,
118 : "Invalid TOTALROWBYTES: " CPL_FRMT_GIB, nTotalRowBytes);
119 0 : nTotalRowBytes = 0;
120 : }
121 0 : if (nTotalRowBytes > 0)
122 0 : nLineOffsetBits = static_cast<vsi_l_offset>(nTotalRowBytes * 8);
123 : else
124 0 : nLineOffsetBits = static_cast<vsi_l_offset>(nPixelOffsetBits) *
125 0 : poDS->GetRasterXSize();
126 :
127 0 : nBlockXSize = poDS->GetRasterXSize();
128 0 : nBlockYSize = 1;
129 :
130 0 : SetMetadataItem("NBITS", CPLString().Printf("%d", nBits),
131 : "IMAGE_STRUCTURE");
132 : }
133 188 : }
134 :
135 : /************************************************************************/
136 : /* IReadBlock() */
137 : /************************************************************************/
138 :
139 5351 : CPLErr EHdrRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage)
140 :
141 : {
142 5351 : if (nBits >= 8)
143 5351 : return RawRasterBand::IReadBlock(nBlockXOff, nBlockYOff, pImage);
144 :
145 : // Establish desired position.
146 0 : const vsi_l_offset nLineStart =
147 0 : (nStartBit + nLineOffsetBits * nBlockYOff) / 8;
148 0 : int iBitOffset =
149 0 : static_cast<int>((nStartBit + nLineOffsetBits * nBlockYOff) % 8);
150 0 : const vsi_l_offset nLineEnd =
151 0 : (nStartBit + nLineOffsetBits * nBlockYOff +
152 0 : static_cast<vsi_l_offset>(nPixelOffsetBits) * nBlockXSize - 1) /
153 : 8;
154 0 : const vsi_l_offset nLineBytesBig = nLineEnd - nLineStart + 1;
155 0 : if (nLineBytesBig >
156 0 : static_cast<vsi_l_offset>(std::numeric_limits<int>::max()))
157 0 : return CE_Failure;
158 0 : const unsigned int nLineBytes = static_cast<unsigned int>(nLineBytesBig);
159 :
160 : // Read data into buffer.
161 0 : GByte *pabyBuffer = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nLineBytes));
162 0 : if (pabyBuffer == nullptr)
163 0 : return CE_Failure;
164 :
165 0 : if (VSIFSeekL(GetFPL(), nLineStart, SEEK_SET) != 0 ||
166 0 : VSIFReadL(pabyBuffer, 1, nLineBytes, GetFPL()) != nLineBytes)
167 : {
168 0 : CPLError(CE_Failure, CPLE_FileIO,
169 : "Failed to read %u bytes at offset %lu.\n%s", nLineBytes,
170 0 : static_cast<unsigned long>(nLineStart), VSIStrerror(errno));
171 0 : CPLFree(pabyBuffer);
172 0 : return CE_Failure;
173 : }
174 :
175 : // Copy data, promoting to 8bit.
176 0 : for (int iX = 0, iPixel = 0; iX < nBlockXSize; iX++)
177 : {
178 0 : int nOutWord = 0;
179 :
180 0 : for (int iBit = 0; iBit < nBits; iBit++)
181 : {
182 0 : if (pabyBuffer[iBitOffset >> 3] & (0x80 >> (iBitOffset & 7)))
183 0 : nOutWord |= (1 << (nBits - 1 - iBit));
184 0 : iBitOffset++;
185 : }
186 :
187 0 : iBitOffset = iBitOffset + nPixelOffsetBits - nBits;
188 :
189 0 : reinterpret_cast<GByte *>(pImage)[iPixel++] =
190 : static_cast<GByte>(nOutWord);
191 : }
192 :
193 0 : CPLFree(pabyBuffer);
194 :
195 0 : return CE_None;
196 : }
197 :
198 : /************************************************************************/
199 : /* IWriteBlock() */
200 : /************************************************************************/
201 :
202 3429 : CPLErr EHdrRasterBand::IWriteBlock(int nBlockXOff, int nBlockYOff, void *pImage)
203 :
204 : {
205 3429 : if (nBits >= 8)
206 3429 : return RawRasterBand::IWriteBlock(nBlockXOff, nBlockYOff, pImage);
207 :
208 : // Establish desired position.
209 0 : const vsi_l_offset nLineStart =
210 0 : (nStartBit + nLineOffsetBits * nBlockYOff) / 8;
211 0 : int iBitOffset =
212 0 : static_cast<int>((nStartBit + nLineOffsetBits * nBlockYOff) % 8);
213 0 : const vsi_l_offset nLineEnd =
214 0 : (nStartBit + nLineOffsetBits * nBlockYOff +
215 0 : static_cast<vsi_l_offset>(nPixelOffsetBits) * nBlockXSize - 1) /
216 : 8;
217 0 : const vsi_l_offset nLineBytesBig = nLineEnd - nLineStart + 1;
218 0 : if (nLineBytesBig >
219 0 : static_cast<vsi_l_offset>(std::numeric_limits<int>::max()))
220 0 : return CE_Failure;
221 0 : const unsigned int nLineBytes = static_cast<unsigned int>(nLineBytesBig);
222 :
223 : // Read data into buffer.
224 0 : GByte *pabyBuffer = static_cast<GByte *>(VSI_CALLOC_VERBOSE(nLineBytes, 1));
225 0 : if (pabyBuffer == nullptr)
226 0 : return CE_Failure;
227 :
228 0 : if (VSIFSeekL(GetFPL(), nLineStart, SEEK_SET) != 0)
229 : {
230 0 : CPLError(CE_Failure, CPLE_FileIO,
231 : "Failed to read %u bytes at offset %lu.\n%s", nLineBytes,
232 0 : static_cast<unsigned long>(nLineStart), VSIStrerror(errno));
233 0 : CPLFree(pabyBuffer);
234 0 : return CE_Failure;
235 : }
236 :
237 0 : CPL_IGNORE_RET_VAL(VSIFReadL(pabyBuffer, nLineBytes, 1, GetFPL()));
238 :
239 : // Copy data, promoting to 8bit.
240 0 : for (int iX = 0, iPixel = 0; iX < nBlockXSize; iX++)
241 : {
242 0 : const int nOutWord = reinterpret_cast<GByte *>(pImage)[iPixel++];
243 :
244 0 : for (int iBit = 0; iBit < nBits; iBit++)
245 : {
246 0 : if (nOutWord & (1 << (nBits - 1 - iBit)))
247 0 : pabyBuffer[iBitOffset >> 3] |= (0x80 >> (iBitOffset & 7));
248 : else
249 0 : pabyBuffer[iBitOffset >> 3] &= ~((0x80 >> (iBitOffset & 7)));
250 :
251 0 : iBitOffset++;
252 : }
253 :
254 0 : iBitOffset = iBitOffset + nPixelOffsetBits - nBits;
255 : }
256 :
257 : // Write the data back out.
258 0 : if (VSIFSeekL(GetFPL(), nLineStart, SEEK_SET) != 0 ||
259 0 : VSIFWriteL(pabyBuffer, 1, nLineBytes, GetFPL()) != nLineBytes)
260 : {
261 0 : CPLError(CE_Failure, CPLE_FileIO,
262 : "Failed to write %u bytes at offset %lu.\n%s", nLineBytes,
263 0 : static_cast<unsigned long>(nLineStart), VSIStrerror(errno));
264 0 : return CE_Failure;
265 : }
266 :
267 0 : CPLFree(pabyBuffer);
268 :
269 0 : return CE_None;
270 : }
271 :
272 : /************************************************************************/
273 : /* IRasterIO() */
274 : /************************************************************************/
275 :
276 1744 : CPLErr EHdrRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
277 : int nXSize, int nYSize, void *pData,
278 : int nBufXSize, int nBufYSize,
279 : GDALDataType eBufType, GSpacing nPixelSpace,
280 : GSpacing nLineSpace,
281 : GDALRasterIOExtraArg *psExtraArg)
282 :
283 : {
284 : // Defer to RawRasterBand
285 1744 : if (nBits >= 8)
286 1744 : return RawRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
287 : pData, nBufXSize, nBufYSize, eBufType,
288 1744 : nPixelSpace, nLineSpace, psExtraArg);
289 :
290 : // Force use of IReadBlock() and IWriteBlock()
291 0 : return GDALRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
292 : pData, nBufXSize, nBufYSize, eBufType,
293 0 : nPixelSpace, nLineSpace, psExtraArg);
294 : }
295 :
296 : /************************************************************************/
297 : /* OSR_GDS() */
298 : /************************************************************************/
299 :
300 19 : static const char *OSR_GDS(char *pszResult, int nResultLen, char **papszNV,
301 : const char *pszField, const char *pszDefaultValue)
302 :
303 : {
304 19 : if (papszNV == nullptr || papszNV[0] == nullptr)
305 0 : return pszDefaultValue;
306 :
307 19 : int iLine = 0; // Used after for.
308 68 : for (; papszNV[iLine] != nullptr &&
309 50 : !EQUALN(papszNV[iLine], pszField, strlen(pszField));
310 : iLine++)
311 : {
312 : }
313 :
314 19 : if (papszNV[iLine] == nullptr)
315 18 : return pszDefaultValue;
316 :
317 1 : char **papszTokens = CSLTokenizeString(papszNV[iLine]);
318 :
319 1 : if (CSLCount(papszTokens) > 1)
320 1 : strncpy(pszResult, papszTokens[1], nResultLen - 1);
321 : else
322 0 : strncpy(pszResult, pszDefaultValue, nResultLen - 1);
323 1 : pszResult[nResultLen - 1] = '\0';
324 :
325 1 : CSLDestroy(papszTokens);
326 1 : return pszResult;
327 : }
328 :
329 : /************************************************************************/
330 : /* ==================================================================== */
331 : /* EHdrDataset */
332 : /* ==================================================================== */
333 : /************************************************************************/
334 :
335 : /************************************************************************/
336 : /* EHdrDataset() */
337 : /************************************************************************/
338 :
339 120 : EHdrDataset::EHdrDataset()
340 : : fpImage(nullptr), osHeaderExt("hdr"), bGotTransform(false),
341 120 : bHDRDirty(false), papszHDR(nullptr), bCLRDirty(false)
342 : {
343 120 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
344 120 : }
345 :
346 : /************************************************************************/
347 : /* ~EHdrDataset() */
348 : /************************************************************************/
349 :
350 240 : EHdrDataset::~EHdrDataset()
351 :
352 : {
353 120 : EHdrDataset::Close();
354 240 : }
355 :
356 : /************************************************************************/
357 : /* Close() */
358 : /************************************************************************/
359 :
360 235 : CPLErr EHdrDataset::Close()
361 : {
362 235 : CPLErr eErr = CE_None;
363 235 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
364 : {
365 120 : if (EHdrDataset::FlushCache(true) != CE_None)
366 0 : eErr = CE_Failure;
367 :
368 120 : if (nBands > 0 && GetAccess() == GA_Update)
369 : {
370 : int bNoDataSet;
371 : RawRasterBand *poBand =
372 53 : reinterpret_cast<RawRasterBand *>(GetRasterBand(1));
373 :
374 53 : const double dfNoData = poBand->GetNoDataValue(&bNoDataSet);
375 53 : if (bNoDataSet)
376 : {
377 2 : ResetKeyValue("NODATA", CPLString().Printf("%.8g", dfNoData));
378 : }
379 :
380 53 : if (bCLRDirty)
381 4 : RewriteCLR(poBand);
382 :
383 53 : if (bHDRDirty)
384 : {
385 40 : if (RewriteHDR() != CE_None)
386 2 : eErr = CE_Failure;
387 : }
388 : }
389 :
390 120 : if (fpImage)
391 : {
392 120 : if (VSIFCloseL(fpImage) != 0)
393 : {
394 0 : CPLError(CE_Failure, CPLE_FileIO, "I/O error");
395 0 : eErr = CE_Failure;
396 : }
397 : }
398 :
399 120 : CSLDestroy(papszHDR);
400 120 : if (GDALPamDataset::Close() != CE_None)
401 0 : eErr = CE_Failure;
402 : }
403 235 : return eErr;
404 : }
405 :
406 : /************************************************************************/
407 : /* GetKeyValue() */
408 : /************************************************************************/
409 :
410 0 : const char *EHdrDataset::GetKeyValue(const char *pszKey, const char *pszDefault)
411 :
412 : {
413 0 : for (int i = 0; papszHDR[i] != nullptr; i++)
414 : {
415 0 : if (EQUALN(pszKey, papszHDR[i], strlen(pszKey)) &&
416 0 : isspace(static_cast<unsigned char>(papszHDR[i][strlen(pszKey)])))
417 : {
418 0 : const char *pszValue = papszHDR[i] + strlen(pszKey);
419 0 : while (isspace(static_cast<unsigned char>(*pszValue)))
420 0 : pszValue++;
421 :
422 0 : return pszValue;
423 : }
424 : }
425 :
426 0 : return pszDefault;
427 : }
428 :
429 : /************************************************************************/
430 : /* ResetKeyValue() */
431 : /* */
432 : /* Replace or add the keyword with the indicated value in the */
433 : /* papszHDR list. */
434 : /************************************************************************/
435 :
436 158 : void EHdrDataset::ResetKeyValue(const char *pszKey, const char *pszValue)
437 :
438 : {
439 158 : if (strlen(pszValue) > 65)
440 : {
441 0 : CPLAssert(strlen(pszValue) <= 65);
442 0 : return;
443 : }
444 :
445 158 : char szNewLine[82] = {'\0'};
446 158 : snprintf(szNewLine, sizeof(szNewLine), "%-15s%s", pszKey, pszValue);
447 :
448 1782 : for (int i = CSLCount(papszHDR) - 1; i >= 0; i--)
449 : {
450 1624 : if (EQUALN(papszHDR[i], szNewLine, strlen(pszKey) + 1))
451 : {
452 0 : if (strcmp(papszHDR[i], szNewLine) != 0)
453 : {
454 0 : CPLFree(papszHDR[i]);
455 0 : papszHDR[i] = CPLStrdup(szNewLine);
456 0 : bHDRDirty = true;
457 : }
458 0 : return;
459 : }
460 : }
461 :
462 158 : bHDRDirty = true;
463 158 : papszHDR = CSLAddString(papszHDR, szNewLine);
464 : }
465 :
466 : /************************************************************************/
467 : /* RewriteCLR() */
468 : /************************************************************************/
469 :
470 4 : void EHdrDataset::RewriteCLR(GDALRasterBand *poBand) const
471 :
472 : {
473 4 : CPLString osCLRFilename = CPLResetExtensionSafe(GetDescription(), "clr");
474 4 : GDALColorTable *poTable = poBand->GetColorTable();
475 4 : GDALRasterAttributeTable *poRAT = poBand->GetDefaultRAT();
476 4 : if (poTable || poRAT)
477 : {
478 3 : VSILFILE *fp = VSIFOpenL(osCLRFilename, "wt");
479 3 : if (fp != nullptr)
480 : {
481 : // Write RAT in priority if both are defined
482 3 : if (poRAT)
483 : {
484 26 : for (int iEntry = 0; iEntry < poRAT->GetRowCount(); iEntry++)
485 : {
486 25 : CPLString oLine;
487 : oLine.Printf("%3d %3d %3d %3d\n",
488 25 : poRAT->GetValueAsInt(iEntry, 0),
489 25 : poRAT->GetValueAsInt(iEntry, 1),
490 25 : poRAT->GetValueAsInt(iEntry, 2),
491 25 : poRAT->GetValueAsInt(iEntry, 3));
492 50 : if (VSIFWriteL(reinterpret_cast<void *>(
493 25 : const_cast<char *>(oLine.c_str())),
494 25 : strlen(oLine), 1, fp) != 1)
495 : {
496 0 : CPLError(CE_Failure, CPLE_FileIO,
497 : "Error while write color table");
498 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
499 0 : return;
500 : }
501 : }
502 : }
503 : else
504 : {
505 8 : for (int iColor = 0; iColor < poTable->GetColorEntryCount();
506 : iColor++)
507 : {
508 : GDALColorEntry sEntry;
509 6 : poTable->GetColorEntryAsRGB(iColor, &sEntry);
510 :
511 : // I wish we had a way to mark transparency.
512 6 : CPLString oLine;
513 6 : oLine.Printf("%3d %3d %3d %3d\n", iColor, sEntry.c1,
514 6 : sEntry.c2, sEntry.c3);
515 12 : if (VSIFWriteL(reinterpret_cast<void *>(
516 6 : const_cast<char *>(oLine.c_str())),
517 6 : strlen(oLine), 1, fp) != 1)
518 : {
519 0 : CPLError(CE_Failure, CPLE_FileIO,
520 : "Error while write color table");
521 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
522 0 : return;
523 : }
524 : }
525 : }
526 3 : if (VSIFCloseL(fp) != 0)
527 : {
528 0 : CPLError(CE_Failure, CPLE_FileIO,
529 : "Error while write color table");
530 : }
531 : }
532 : else
533 : {
534 0 : CPLError(CE_Failure, CPLE_OpenFailed,
535 : "Unable to create color file %s.", osCLRFilename.c_str());
536 3 : }
537 : }
538 : else
539 : {
540 1 : VSIUnlink(osCLRFilename);
541 : }
542 : }
543 :
544 : /************************************************************************/
545 : /* SetSpatialRef() */
546 : /************************************************************************/
547 :
548 37 : CPLErr EHdrDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
549 :
550 : {
551 : // Reset coordinate system on the dataset.
552 37 : m_oSRS.Clear();
553 37 : if (poSRS == nullptr)
554 0 : return CE_None;
555 :
556 37 : m_oSRS = *poSRS;
557 : // Convert to ESRI WKT.
558 37 : char *pszESRI_SRS = nullptr;
559 37 : const char *const apszOptions[] = {"FORMAT=WKT1_ESRI", nullptr};
560 37 : m_oSRS.exportToWkt(&pszESRI_SRS, apszOptions);
561 :
562 37 : if (pszESRI_SRS)
563 : {
564 : // Write to .prj file.
565 : CPLString osPrjFilename =
566 37 : CPLResetExtensionSafe(GetDescription(), "prj");
567 37 : VSILFILE *fp = VSIFOpenL(osPrjFilename.c_str(), "wt");
568 37 : if (fp != nullptr)
569 : {
570 37 : size_t nCount = VSIFWriteL(pszESRI_SRS, strlen(pszESRI_SRS), 1, fp);
571 37 : nCount += VSIFWriteL("\n", 1, 1, fp);
572 37 : if (VSIFCloseL(fp) != 0 || nCount != 2)
573 : {
574 2 : CPLFree(pszESRI_SRS);
575 2 : return CE_Failure;
576 : }
577 : }
578 :
579 35 : CPLFree(pszESRI_SRS);
580 : }
581 :
582 35 : return CE_None;
583 : }
584 :
585 : /************************************************************************/
586 : /* GetGeoTransform() */
587 : /************************************************************************/
588 :
589 29 : CPLErr EHdrDataset::GetGeoTransform(GDALGeoTransform >) const
590 :
591 : {
592 29 : if (bGotTransform)
593 : {
594 29 : gt = m_gt;
595 29 : return CE_None;
596 : }
597 :
598 0 : return GDALPamDataset::GetGeoTransform(gt);
599 : }
600 :
601 : /************************************************************************/
602 : /* SetGeoTransform() */
603 : /************************************************************************/
604 :
605 39 : CPLErr EHdrDataset::SetGeoTransform(const GDALGeoTransform >)
606 :
607 : {
608 : // We only support non-rotated images with info in the .HDR file.
609 39 : if (gt[2] != 0.0 || gt[4] != 0.0)
610 : {
611 0 : return GDALPamDataset::SetGeoTransform(gt);
612 : }
613 :
614 : // Record new geotransform.
615 39 : bGotTransform = true;
616 39 : m_gt = gt;
617 :
618 : // Strip out all old geotransform keywords from HDR records.
619 385 : for (int i = CSLCount(papszHDR) - 1; i >= 0; i--)
620 : {
621 346 : if (STARTS_WITH_CI(papszHDR[i], "ul") ||
622 344 : STARTS_WITH_CI(papszHDR[i] + 1, "ll") ||
623 344 : STARTS_WITH_CI(papszHDR[i], "cell") ||
624 344 : STARTS_WITH_CI(papszHDR[i] + 1, "dim"))
625 : {
626 4 : papszHDR = CSLRemoveStrings(papszHDR, i, 1, nullptr);
627 : }
628 : }
629 :
630 : // Set the transformation information.
631 39 : CPLString oValue;
632 :
633 39 : oValue.Printf("%.15g", m_gt[0] + m_gt[1] * 0.5);
634 39 : ResetKeyValue("ULXMAP", oValue);
635 :
636 39 : oValue.Printf("%.15g", m_gt[3] + m_gt[5] * 0.5);
637 39 : ResetKeyValue("ULYMAP", oValue);
638 :
639 39 : oValue.Printf("%.15g", m_gt[1]);
640 39 : ResetKeyValue("XDIM", oValue);
641 :
642 39 : oValue.Printf("%.15g", fabs(m_gt[5]));
643 39 : ResetKeyValue("YDIM", oValue);
644 :
645 39 : return CE_None;
646 : }
647 :
648 : /************************************************************************/
649 : /* RewriteHDR() */
650 : /************************************************************************/
651 :
652 40 : CPLErr EHdrDataset::RewriteHDR()
653 :
654 : {
655 80 : const CPLString osPath = CPLGetPathSafe(GetDescription());
656 80 : const CPLString osName = CPLGetBasenameSafe(GetDescription());
657 : CPLString osHDRFilename =
658 80 : CPLFormCIFilenameSafe(osPath, osName, osHeaderExt);
659 :
660 : // Write .hdr file.
661 40 : VSILFILE *fp = VSIFOpenL(osHDRFilename, "wt");
662 :
663 40 : if (fp == nullptr)
664 : {
665 0 : CPLError(CE_Failure, CPLE_OpenFailed, "Failed to rewrite .hdr file %s.",
666 : osHDRFilename.c_str());
667 0 : return CE_Failure;
668 : }
669 :
670 541 : for (int i = 0; papszHDR[i] != nullptr; i++)
671 : {
672 503 : size_t nCount = VSIFWriteL(papszHDR[i], strlen(papszHDR[i]), 1, fp);
673 503 : nCount += VSIFWriteL("\n", 1, 1, fp);
674 503 : if (nCount != 2)
675 : {
676 2 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
677 2 : return CE_Failure;
678 : }
679 : }
680 :
681 38 : bHDRDirty = false;
682 :
683 38 : if (VSIFCloseL(fp) != 0)
684 0 : return CE_Failure;
685 :
686 38 : return CE_None;
687 : }
688 :
689 : /************************************************************************/
690 : /* RewriteSTX() */
691 : /************************************************************************/
692 :
693 4 : CPLErr EHdrDataset::RewriteSTX() const
694 : {
695 8 : const CPLString osPath = CPLGetPathSafe(GetDescription());
696 8 : const CPLString osName = CPLGetBasenameSafe(GetDescription());
697 : const CPLString osSTXFilename =
698 8 : CPLFormCIFilenameSafe(osPath, osName, "stx");
699 :
700 4 : VSILFILE *fp = VSIFOpenL(osSTXFilename, "wt");
701 4 : if (fp == nullptr)
702 : {
703 0 : CPLDebug("EHDR", "Failed to rewrite .stx file %s.",
704 : osSTXFilename.c_str());
705 0 : return CE_Failure;
706 : }
707 :
708 4 : bool bOK = true;
709 8 : for (int i = 0; bOK && i < nBands; ++i)
710 : {
711 4 : EHdrRasterBand *poBand =
712 4 : reinterpret_cast<EHdrRasterBand *>(papoBands[i]);
713 4 : bOK &= VSIFPrintfL(fp, "%d %.10f %.10f ", i + 1, poBand->dfMin,
714 4 : poBand->dfMax) >= 0;
715 4 : if (poBand->minmaxmeanstddev & HAS_MEAN_FLAG)
716 4 : bOK &= VSIFPrintfL(fp, "%.10f ", poBand->dfMean) >= 0;
717 : else
718 0 : bOK &= VSIFPrintfL(fp, "# ") >= 0;
719 :
720 4 : if (poBand->minmaxmeanstddev & HAS_STDDEV_FLAG)
721 4 : bOK &= VSIFPrintfL(fp, "%.10f\n", poBand->dfStdDev) >= 0;
722 : else
723 0 : bOK &= VSIFPrintfL(fp, "#\n") >= 0;
724 : }
725 :
726 4 : if (VSIFCloseL(fp) != 0)
727 0 : bOK = false;
728 :
729 4 : return bOK ? CE_None : CE_Failure;
730 : }
731 :
732 : /************************************************************************/
733 : /* ReadSTX() */
734 : /************************************************************************/
735 :
736 120 : CPLErr EHdrDataset::ReadSTX() const
737 : {
738 240 : const CPLString osPath = CPLGetPathSafe(GetDescription());
739 240 : const CPLString osName = CPLGetBasenameSafe(GetDescription());
740 : const CPLString osSTXFilename =
741 240 : CPLFormCIFilenameSafe(osPath, osName, "stx");
742 :
743 120 : VSILFILE *fp = VSIFOpenL(osSTXFilename, "rt");
744 120 : if (fp == nullptr)
745 117 : return CE_None;
746 :
747 3 : const char *pszLine = nullptr;
748 6 : while ((pszLine = CPLReadLineL(fp)) != nullptr)
749 : {
750 : char **papszTokens =
751 3 : CSLTokenizeStringComplex(pszLine, " \t", TRUE, FALSE);
752 3 : const int nTokens = CSLCount(papszTokens);
753 3 : if (nTokens >= 5)
754 : {
755 3 : const int i = atoi(papszTokens[0]);
756 3 : if (i > 0 && i <= nBands)
757 : {
758 3 : EHdrRasterBand *poBand =
759 3 : reinterpret_cast<EHdrRasterBand *>(papoBands[i - 1]);
760 3 : poBand->dfMin = CPLAtof(papszTokens[1]);
761 3 : poBand->dfMax = CPLAtof(papszTokens[2]);
762 :
763 3 : int bNoDataSet = FALSE;
764 3 : const double dfNoData = poBand->GetNoDataValue(&bNoDataSet);
765 3 : if (bNoDataSet && dfNoData == poBand->dfMin)
766 : {
767 : // Triggered by
768 : // /vsicurl/http://eros.usgs.gov/archive/nslrsda/GeoTowns/HongKong/srtm/n22e113.zip/n22e113.bil
769 0 : CPLDebug(
770 : "EHDr",
771 : "Ignoring .stx file where min == nodata. "
772 : "The nodata value should not be taken into account "
773 : "in minimum value computation.");
774 0 : CSLDestroy(papszTokens);
775 0 : papszTokens = nullptr;
776 0 : break;
777 : }
778 :
779 3 : poBand->minmaxmeanstddev = HAS_MIN_FLAG | HAS_MAX_FLAG;
780 : // Reads optional mean and stddev.
781 3 : if (!EQUAL(papszTokens[3], "#"))
782 : {
783 3 : poBand->dfMean = CPLAtof(papszTokens[3]);
784 3 : poBand->minmaxmeanstddev |= HAS_MEAN_FLAG;
785 : }
786 3 : if (!EQUAL(papszTokens[4], "#"))
787 : {
788 3 : poBand->dfStdDev = CPLAtof(papszTokens[4]);
789 3 : poBand->minmaxmeanstddev |= HAS_STDDEV_FLAG;
790 : }
791 :
792 3 : if (nTokens >= 6 && !EQUAL(papszTokens[5], "#"))
793 0 : poBand->SetMetadataItem("STRETCHMIN", papszTokens[5],
794 0 : "RENDERING_HINTS");
795 :
796 3 : if (nTokens >= 7 && !EQUAL(papszTokens[6], "#"))
797 0 : poBand->SetMetadataItem("STRETCHMAX", papszTokens[6],
798 0 : "RENDERING_HINTS");
799 : }
800 : }
801 :
802 3 : CSLDestroy(papszTokens);
803 : }
804 :
805 3 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
806 :
807 3 : return CE_None;
808 : }
809 :
810 : /************************************************************************/
811 : /* GetImageRepFilename() */
812 : /************************************************************************/
813 :
814 : // Check for IMAGE.REP (Spatiocarte Defense 1.0) or name_of_image.rep
815 : // if it is a GIS-GeoSPOT image.
816 : // For the specification of SPDF (in French), see
817 : // http://eden.ign.fr/download/pub/doc/emabgi/spdf10.pdf/download
818 :
819 103 : CPLString EHdrDataset::GetImageRepFilename(const char *pszFilename)
820 : {
821 :
822 206 : const CPLString osPath = CPLGetPathSafe(pszFilename);
823 206 : const CPLString osName = CPLGetBasenameSafe(pszFilename);
824 206 : CPLString osREPFilename = CPLFormCIFilenameSafe(osPath, osName, "rep");
825 :
826 : VSIStatBufL sStatBuf;
827 103 : if (VSIStatExL(osREPFilename.c_str(), &sStatBuf, VSI_STAT_EXISTS_FLAG) == 0)
828 0 : return osREPFilename;
829 :
830 206 : if (EQUAL(CPLGetFilename(pszFilename), "imspatio.bil") ||
831 103 : EQUAL(CPLGetFilename(pszFilename), "haspatio.bil"))
832 : {
833 : CPLString osImageRepFilename(
834 0 : CPLFormCIFilenameSafe(osPath, "image", "rep"));
835 0 : if (VSIStatExL(osImageRepFilename, &sStatBuf, VSI_STAT_EXISTS_FLAG) ==
836 : 0)
837 0 : return osImageRepFilename;
838 :
839 : // Try in the upper directories if not found in the BIL image directory.
840 0 : CPLString dirName(CPLGetDirnameSafe(osPath));
841 0 : if (CPLIsFilenameRelative(osPath.c_str()))
842 : {
843 0 : char *cwd = CPLGetCurrentDir();
844 0 : if (cwd)
845 : {
846 0 : dirName = CPLFormFilenameSafe(cwd, dirName.c_str(), nullptr);
847 0 : CPLFree(cwd);
848 : }
849 : }
850 0 : while (dirName[0] != 0 && EQUAL(dirName, ".") == FALSE &&
851 0 : EQUAL(dirName, "/") == FALSE)
852 : {
853 : osImageRepFilename =
854 0 : CPLFormCIFilenameSafe(dirName.c_str(), "image", "rep");
855 0 : if (VSIStatExL(osImageRepFilename.c_str(), &sStatBuf,
856 0 : VSI_STAT_EXISTS_FLAG) == 0)
857 0 : return osImageRepFilename;
858 :
859 : // Don't try to recurse above the 'image' subdirectory.
860 0 : if (EQUAL(dirName, "image"))
861 : {
862 0 : break;
863 : }
864 0 : dirName = CPLString(CPLGetDirnameSafe(dirName));
865 : }
866 : }
867 103 : return CPLString();
868 : }
869 :
870 : /************************************************************************/
871 : /* GetFileList() */
872 : /************************************************************************/
873 :
874 22 : char **EHdrDataset::GetFileList()
875 :
876 : {
877 44 : const CPLString osPath = CPLGetPathSafe(GetDescription());
878 44 : const CPLString osName = CPLGetBasenameSafe(GetDescription());
879 :
880 : // Main data file, etc.
881 22 : char **papszFileList = GDALPamDataset::GetFileList();
882 :
883 : // Header file.
884 44 : CPLString osFilename = CPLFormCIFilenameSafe(osPath, osName, osHeaderExt);
885 22 : papszFileList = CSLAddString(papszFileList, osFilename);
886 :
887 : // Statistics file
888 22 : osFilename = CPLFormCIFilenameSafe(osPath, osName, "stx");
889 : VSIStatBufL sStatBuf;
890 22 : if (VSIStatExL(osFilename, &sStatBuf, VSI_STAT_EXISTS_FLAG) == 0)
891 2 : papszFileList = CSLAddString(papszFileList, osFilename);
892 :
893 : // color table file.
894 22 : osFilename = CPLFormCIFilenameSafe(osPath, osName, "clr");
895 22 : if (VSIStatExL(osFilename, &sStatBuf, VSI_STAT_EXISTS_FLAG) == 0)
896 2 : papszFileList = CSLAddString(papszFileList, osFilename);
897 :
898 : // projections file.
899 22 : osFilename = CPLFormCIFilenameSafe(osPath, osName, "prj");
900 22 : if (VSIStatExL(osFilename, &sStatBuf, VSI_STAT_EXISTS_FLAG) == 0)
901 8 : papszFileList = CSLAddString(papszFileList, osFilename);
902 :
903 22 : const CPLString imageRepFilename = GetImageRepFilename(GetDescription());
904 22 : if (!imageRepFilename.empty())
905 0 : papszFileList = CSLAddString(papszFileList, imageRepFilename.c_str());
906 :
907 44 : return papszFileList;
908 : }
909 :
910 : /************************************************************************/
911 : /* Open() */
912 : /************************************************************************/
913 :
914 30763 : GDALDataset *EHdrDataset::Open(GDALOpenInfo *poOpenInfo)
915 :
916 : {
917 30763 : return Open(poOpenInfo, true);
918 : }
919 :
920 30816 : GDALDataset *EHdrDataset::Open(GDALOpenInfo *poOpenInfo, bool bFileSizeCheck)
921 :
922 : {
923 : // Assume the caller is pointing to the binary (i.e. .bil) file.
924 30816 : if (poOpenInfo->nHeaderBytes < 2 || poOpenInfo->fpL == nullptr)
925 29469 : return nullptr;
926 :
927 : // Tear apart the filename to form a .HDR filename.
928 2694 : const CPLString osPath = CPLGetPathSafe(poOpenInfo->pszFilename);
929 2694 : const CPLString osName = CPLGetBasenameSafe(poOpenInfo->pszFilename);
930 :
931 1347 : const char *pszHeaderExt = "hdr";
932 1347 : if (poOpenInfo->IsExtensionEqualToCI("SRC") && osName.size() == 7 &&
933 0 : (osName[0] == 'e' || osName[0] == 'E' || osName[0] == 'w' ||
934 1347 : osName[0] == 'W') &&
935 0 : (osName[4] == 'n' || osName[4] == 'N' || osName[4] == 's' ||
936 0 : osName[4] == 'S'))
937 : {
938 : // It is a GTOPO30 or SRTM30 source file, whose header extension is .sch
939 : // see http://dds.cr.usgs.gov/srtm/version1/SRTM30/GTOPO30_Documentation
940 0 : pszHeaderExt = "sch";
941 : }
942 :
943 1347 : char **papszSiblingFiles = poOpenInfo->GetSiblingFiles();
944 2694 : CPLString osHDRFilename;
945 1347 : if (papszSiblingFiles)
946 : {
947 1342 : const int iFile = CSLFindString(
948 : papszSiblingFiles,
949 2684 : CPLFormFilenameSafe(nullptr, osName, pszHeaderExt).c_str());
950 1342 : if (iFile < 0) // Return if there is no corresponding .hdr file.
951 1200 : return nullptr;
952 :
953 : osHDRFilename =
954 142 : CPLFormFilenameSafe(osPath, papszSiblingFiles[iFile], nullptr);
955 : }
956 : else
957 : {
958 5 : osHDRFilename = CPLFormCIFilenameSafe(osPath, osName, pszHeaderExt);
959 : }
960 :
961 147 : const bool bSelectedHDR = EQUAL(osHDRFilename, poOpenInfo->pszFilename);
962 :
963 : // Do we have a .hdr file?
964 147 : VSILFILE *fp = VSIFOpenL(osHDRFilename, "r");
965 147 : if (fp == nullptr)
966 : {
967 5 : return nullptr;
968 : }
969 :
970 : // Is this file an ESRI header file? Read a few lines of text
971 : // searching for something starting with nrows or ncols.
972 142 : int nRows = -1;
973 142 : int nCols = -1;
974 142 : int l_nBands = 1;
975 142 : int nSkipBytes = 0;
976 142 : double dfULXMap = 0.5;
977 142 : double dfULYMap = 0.5;
978 142 : double dfYLLCorner = -123.456;
979 142 : int bCenter = TRUE;
980 142 : double dfXDim = 1.0;
981 142 : double dfYDim = 1.0;
982 142 : double dfNoData = 0.0;
983 142 : int nLineCount = 0;
984 142 : int bNoDataSet = FALSE;
985 142 : GDALDataType eDataType = GDT_Byte;
986 142 : int nBits = -1;
987 142 : char chByteOrder = 'M';
988 142 : char chPixelType = 'N'; // Not defined.
989 142 : char szLayout[10] = "BIL";
990 142 : char **papszHDR = nullptr;
991 142 : int bHasInternalProjection = FALSE;
992 142 : int bHasMin = FALSE;
993 142 : int bHasMax = FALSE;
994 142 : double dfMin = 0;
995 142 : double dfMax = 0;
996 :
997 142 : const char *pszLine = nullptr;
998 1643 : while ((pszLine = CPLReadLineL(fp)) != nullptr)
999 : {
1000 1502 : nLineCount++;
1001 :
1002 1502 : if (nLineCount > 50 || strlen(pszLine) > 1000)
1003 : break;
1004 :
1005 1501 : papszHDR = CSLAddString(papszHDR, pszLine);
1006 :
1007 : char **papszTokens =
1008 1501 : CSLTokenizeStringComplex(pszLine, " \t", TRUE, FALSE);
1009 1501 : if (CSLCount(papszTokens) < 2)
1010 : {
1011 45 : CSLDestroy(papszTokens);
1012 45 : continue;
1013 : }
1014 :
1015 1456 : if (EQUAL(papszTokens[0], "ncols"))
1016 : {
1017 120 : nCols = atoi(papszTokens[1]);
1018 : }
1019 1336 : else if (EQUAL(papszTokens[0], "nrows"))
1020 : {
1021 122 : nRows = atoi(papszTokens[1]);
1022 : }
1023 1214 : else if (EQUAL(papszTokens[0], "skipbytes"))
1024 : {
1025 0 : nSkipBytes = atoi(papszTokens[1]);
1026 : }
1027 1214 : else if (EQUAL(papszTokens[0], "ulxmap") ||
1028 1170 : EQUAL(papszTokens[0], "xllcorner") ||
1029 1166 : EQUAL(papszTokens[0], "xllcenter"))
1030 : {
1031 48 : dfULXMap = CPLAtofM(papszTokens[1]);
1032 48 : if (EQUAL(papszTokens[0], "xllcorner"))
1033 4 : bCenter = FALSE;
1034 : }
1035 1166 : else if (EQUAL(papszTokens[0], "ulymap"))
1036 : {
1037 44 : dfULYMap = CPLAtofM(papszTokens[1]);
1038 : }
1039 1122 : else if (EQUAL(papszTokens[0], "yllcorner") ||
1040 1118 : EQUAL(papszTokens[0], "yllcenter"))
1041 : {
1042 4 : dfYLLCorner = CPLAtofM(papszTokens[1]);
1043 4 : if (EQUAL(papszTokens[0], "yllcorner"))
1044 4 : bCenter = FALSE;
1045 : }
1046 1118 : else if (EQUAL(papszTokens[0], "xdim"))
1047 : {
1048 44 : dfXDim = CPLAtofM(papszTokens[1]);
1049 : }
1050 1074 : else if (EQUAL(papszTokens[0], "ydim"))
1051 : {
1052 44 : dfYDim = CPLAtofM(papszTokens[1]);
1053 : }
1054 1030 : else if (EQUAL(papszTokens[0], "cellsize"))
1055 : {
1056 4 : dfXDim = CPLAtofM(papszTokens[1]);
1057 4 : dfYDim = dfXDim;
1058 : }
1059 1026 : else if (EQUAL(papszTokens[0], "nbands"))
1060 : {
1061 114 : l_nBands = atoi(papszTokens[1]);
1062 : }
1063 912 : else if (EQUAL(papszTokens[0], "layout"))
1064 : {
1065 120 : snprintf(szLayout, sizeof(szLayout), "%s", papszTokens[1]);
1066 : }
1067 792 : else if (EQUAL(papszTokens[0], "NODATA_value") ||
1068 792 : EQUAL(papszTokens[0], "NODATA"))
1069 : {
1070 4 : dfNoData = CPLAtofM(papszTokens[1]);
1071 4 : bNoDataSet = TRUE;
1072 : }
1073 788 : else if (EQUAL(papszTokens[0], "NBITS"))
1074 : {
1075 114 : nBits = atoi(papszTokens[1]);
1076 : }
1077 674 : else if (EQUAL(papszTokens[0], "PIXELTYPE"))
1078 : {
1079 110 : chPixelType = static_cast<char>(
1080 110 : toupper(static_cast<unsigned char>(papszTokens[1][0])));
1081 : }
1082 564 : else if (EQUAL(papszTokens[0], "byteorder"))
1083 : {
1084 126 : chByteOrder = static_cast<char>(
1085 126 : toupper(static_cast<unsigned char>(papszTokens[1][0])));
1086 : }
1087 :
1088 : // http://www.worldclim.org/futdown.htm have the projection extensions
1089 438 : else if (EQUAL(papszTokens[0], "Projection"))
1090 : {
1091 3 : bHasInternalProjection = TRUE;
1092 : }
1093 435 : else if (EQUAL(papszTokens[0], "MinValue") ||
1094 434 : EQUAL(papszTokens[0], "MIN_VALUE"))
1095 : {
1096 1 : dfMin = CPLAtofM(papszTokens[1]);
1097 1 : bHasMin = TRUE;
1098 : }
1099 434 : else if (EQUAL(papszTokens[0], "MaxValue") ||
1100 433 : EQUAL(papszTokens[0], "MAX_VALUE"))
1101 : {
1102 1 : dfMax = CPLAtofM(papszTokens[1]);
1103 1 : bHasMax = TRUE;
1104 : }
1105 :
1106 1456 : CSLDestroy(papszTokens);
1107 : }
1108 :
1109 142 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
1110 :
1111 : // Did we get the required keywords? If not, return with this never having
1112 : // been considered to be a match. This isn't an error!
1113 142 : if (nRows == -1 || nCols == -1)
1114 : {
1115 22 : CSLDestroy(papszHDR);
1116 22 : return nullptr;
1117 : }
1118 :
1119 240 : if (!GDALCheckDatasetDimensions(nCols, nRows) ||
1120 120 : !GDALCheckBandCount(l_nBands, FALSE))
1121 : {
1122 0 : CSLDestroy(papszHDR);
1123 0 : return nullptr;
1124 : }
1125 :
1126 : // Has the caller selected the .hdr file to open?
1127 120 : if (bSelectedHDR)
1128 : {
1129 0 : CPLError(CE_Failure, CPLE_AppDefined,
1130 : "The selected file is an ESRI BIL header file, but to "
1131 : "open ESRI BIL datasets, the data file should be selected "
1132 : "instead of the .hdr file. Please try again selecting "
1133 : "the data file (often with the extension .bil) corresponding "
1134 : "to the header file: %s",
1135 : poOpenInfo->pszFilename);
1136 0 : CSLDestroy(papszHDR);
1137 0 : return nullptr;
1138 : }
1139 :
1140 : // If we aren't sure of the file type, check the data file size. If it is 4
1141 : // bytes or more per pixel then we assume it is floating point data.
1142 120 : if (nBits == -1 && chPixelType == 'N')
1143 : {
1144 : VSIStatBufL sStatBuf;
1145 6 : if (VSIStatL(poOpenInfo->pszFilename, &sStatBuf) == 0)
1146 : {
1147 6 : const vsi_l_offset nBytes =
1148 6 : sStatBuf.st_size / nCols / nRows / l_nBands;
1149 : // Exit now if nBytes value does not make sense to avoid later overflow in below multiplication
1150 6 : if (nBytes > 8)
1151 : {
1152 0 : CPLError(CE_Failure, CPLE_AppDefined,
1153 : "EHdr driver cannot infer NBITS value");
1154 0 : CSLDestroy(papszHDR);
1155 0 : return nullptr;
1156 : }
1157 6 : if (nBytes > 0 && nBytes != 3)
1158 2 : nBits = static_cast<int>(nBytes * 8);
1159 :
1160 6 : if (nBytes == 4)
1161 2 : chPixelType = 'F';
1162 : }
1163 : }
1164 :
1165 : // If the extension is FLT it is likely a floating point file.
1166 120 : if (chPixelType == 'N')
1167 : {
1168 8 : if (poOpenInfo->IsExtensionEqualToCI("FLT"))
1169 2 : chPixelType = 'F';
1170 : }
1171 :
1172 : // If we have a negative nodata value, assume that the
1173 : // pixel type is signed. This is necessary for datasets from
1174 : // http://www.worldclim.org/futdown.htm
1175 :
1176 120 : if (bNoDataSet && dfNoData < 0 && chPixelType == 'N')
1177 : {
1178 2 : chPixelType = 'S';
1179 : }
1180 :
1181 240 : auto poDS = std::make_unique<EHdrDataset>();
1182 :
1183 120 : poDS->osHeaderExt = pszHeaderExt;
1184 :
1185 120 : poDS->nRasterXSize = nCols;
1186 120 : poDS->nRasterYSize = nRows;
1187 120 : poDS->papszHDR = papszHDR;
1188 120 : std::swap(poDS->fpImage, poOpenInfo->fpL);
1189 120 : poDS->eAccess = poOpenInfo->eAccess;
1190 :
1191 : // Figure out the data type.
1192 120 : if (nBits == 16)
1193 : {
1194 21 : if (chPixelType == 'S')
1195 13 : eDataType = GDT_Int16;
1196 : else
1197 8 : eDataType = GDT_UInt16; // Default
1198 : }
1199 99 : else if (nBits == 32)
1200 : {
1201 34 : if (chPixelType == 'S')
1202 8 : eDataType = GDT_Int32;
1203 26 : else if (chPixelType == 'F')
1204 21 : eDataType = GDT_Float32;
1205 : else
1206 5 : eDataType = GDT_UInt32; // Default
1207 : }
1208 65 : else if (nBits >= 1 && nBits <= 8)
1209 : {
1210 61 : if (chPixelType == 'S')
1211 6 : eDataType = GDT_Int8;
1212 : else
1213 55 : eDataType = GDT_Byte;
1214 61 : nBits = 8;
1215 : }
1216 4 : else if (nBits == -1)
1217 : {
1218 4 : if (chPixelType == 'F')
1219 : {
1220 0 : eDataType = GDT_Float32;
1221 0 : nBits = 32;
1222 : }
1223 : else
1224 : {
1225 4 : eDataType = GDT_Byte;
1226 4 : nBits = 8;
1227 : }
1228 : }
1229 : else
1230 : {
1231 0 : CPLError(CE_Failure, CPLE_NotSupported,
1232 : "EHdr driver does not support %d NBITS value.", nBits);
1233 0 : return nullptr;
1234 : }
1235 :
1236 : // Compute the line offset.
1237 120 : const int nItemSize = GDALGetDataTypeSizeBytes(eDataType);
1238 120 : CPLAssert(nItemSize != 0);
1239 120 : CPLAssert(l_nBands != 0);
1240 :
1241 120 : int nPixelOffset = 0;
1242 120 : int nLineOffset = 0;
1243 120 : vsi_l_offset nBandOffset = 0;
1244 :
1245 120 : if (EQUAL(szLayout, "BIP"))
1246 : {
1247 0 : if (nCols > std::numeric_limits<int>::max() / (nItemSize * l_nBands))
1248 : {
1249 0 : CPLError(CE_Failure, CPLE_AppDefined, "Int overflow occurred.");
1250 0 : return nullptr;
1251 : }
1252 0 : nPixelOffset = nItemSize * l_nBands;
1253 0 : nLineOffset = nPixelOffset * nCols;
1254 0 : nBandOffset = static_cast<vsi_l_offset>(nItemSize);
1255 : }
1256 120 : else if (EQUAL(szLayout, "BSQ"))
1257 : {
1258 0 : if (nCols > std::numeric_limits<int>::max() / nItemSize)
1259 : {
1260 0 : CPLError(CE_Failure, CPLE_AppDefined, "Int overflow occurred.");
1261 0 : return nullptr;
1262 : }
1263 0 : nPixelOffset = nItemSize;
1264 0 : nLineOffset = nPixelOffset * nCols;
1265 0 : nBandOffset = static_cast<vsi_l_offset>(nLineOffset) * nRows;
1266 : }
1267 : else
1268 : {
1269 : // Assume BIL.
1270 120 : if (nCols > std::numeric_limits<int>::max() / (nItemSize * l_nBands))
1271 : {
1272 0 : CPLError(CE_Failure, CPLE_AppDefined, "Int overflow occurred.");
1273 0 : return nullptr;
1274 : }
1275 120 : nPixelOffset = nItemSize;
1276 120 : nLineOffset = nItemSize * l_nBands * nCols;
1277 120 : nBandOffset = static_cast<vsi_l_offset>(nItemSize) * nCols;
1278 : }
1279 :
1280 192 : if (nBits >= 8 && bFileSizeCheck &&
1281 72 : !RAWDatasetCheckMemoryUsage(
1282 72 : poDS->nRasterXSize, poDS->nRasterYSize, l_nBands, nItemSize,
1283 72 : nPixelOffset, nLineOffset, nSkipBytes, nBandOffset, poDS->fpImage))
1284 : {
1285 0 : return nullptr;
1286 : }
1287 :
1288 120 : poDS->SetDescription(poOpenInfo->pszFilename);
1289 120 : poDS->PamInitialize();
1290 :
1291 : // Create band information objects.
1292 308 : for (int i = 0; i < l_nBands; i++)
1293 : {
1294 : auto poBand = std::make_unique<EHdrRasterBand>(
1295 188 : poDS.get(), i + 1, poDS->fpImage, nSkipBytes + nBandOffset * i,
1296 : nPixelOffset, nLineOffset, eDataType,
1297 : chByteOrder == 'I' || chByteOrder == 'L'
1298 188 : ? RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN
1299 : : RawRasterBand::ByteOrder::ORDER_BIG_ENDIAN,
1300 188 : nBits);
1301 188 : if (!poBand->IsValid())
1302 0 : return nullptr;
1303 :
1304 188 : poBand->bNoDataSet = bNoDataSet;
1305 188 : poBand->dfNoData = dfNoData;
1306 :
1307 188 : if (bHasMin && bHasMax)
1308 : {
1309 1 : poBand->dfMin = dfMin;
1310 1 : poBand->dfMax = dfMax;
1311 1 : poBand->minmaxmeanstddev = HAS_MIN_FLAG | HAS_MAX_FLAG;
1312 : }
1313 :
1314 188 : poDS->SetBand(i + 1, std::move(poBand));
1315 : }
1316 :
1317 : // If we didn't get bounds in the .hdr, look for a worldfile.
1318 120 : if (dfYLLCorner != -123.456)
1319 : {
1320 4 : if (bCenter)
1321 0 : dfULYMap = dfYLLCorner + (nRows - 1) * dfYDim;
1322 : else
1323 4 : dfULYMap = dfYLLCorner + nRows * dfYDim;
1324 : }
1325 :
1326 120 : if (dfULXMap != 0.5 || dfULYMap != 0.5 || dfXDim != 1.0 || dfYDim != 1.0)
1327 : {
1328 48 : poDS->bGotTransform = true;
1329 :
1330 48 : if (bCenter)
1331 : {
1332 44 : poDS->m_gt[0] = dfULXMap - dfXDim * 0.5;
1333 44 : poDS->m_gt[1] = dfXDim;
1334 44 : poDS->m_gt[2] = 0.0;
1335 44 : poDS->m_gt[3] = dfULYMap + dfYDim * 0.5;
1336 44 : poDS->m_gt[4] = 0.0;
1337 44 : poDS->m_gt[5] = -dfYDim;
1338 : }
1339 : else
1340 : {
1341 4 : poDS->m_gt[0] = dfULXMap;
1342 4 : poDS->m_gt[1] = dfXDim;
1343 4 : poDS->m_gt[2] = 0.0;
1344 4 : poDS->m_gt[3] = dfULYMap;
1345 4 : poDS->m_gt[4] = 0.0;
1346 4 : poDS->m_gt[5] = -dfYDim;
1347 : }
1348 : }
1349 :
1350 120 : if (!poDS->bGotTransform)
1351 144 : poDS->bGotTransform = CPL_TO_BOOL(GDALReadWorldFile(
1352 144 : poOpenInfo->pszFilename, nullptr, poDS->m_gt.data()));
1353 :
1354 120 : if (!poDS->bGotTransform)
1355 144 : poDS->bGotTransform = CPL_TO_BOOL(GDALReadWorldFile(
1356 144 : poOpenInfo->pszFilename, "wld", poDS->m_gt.data()));
1357 :
1358 : // Check for a .prj file.
1359 240 : std::string osPrjFilename = CPLFormCIFilenameSafe(osPath, osName, "prj");
1360 :
1361 120 : fp = VSIFOpenL(osPrjFilename.c_str(), "r");
1362 :
1363 : // .hdr files from http://www.worldclim.org/futdown.htm have the projection
1364 : // info in the .hdr file itself.
1365 120 : if (fp == nullptr && bHasInternalProjection)
1366 : {
1367 1 : osPrjFilename = std::move(osHDRFilename);
1368 1 : fp = VSIFOpenL(osPrjFilename.c_str(), "r");
1369 : }
1370 :
1371 120 : if (fp != nullptr)
1372 : {
1373 39 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
1374 39 : fp = nullptr;
1375 :
1376 39 : char **papszLines = CSLLoad(osPrjFilename.c_str());
1377 :
1378 39 : if (poDS->m_oSRS.importFromESRI(papszLines) == OGRERR_NONE)
1379 : {
1380 : // If geographic values are in seconds, we must transform.
1381 : // Is there a code for minutes too?
1382 37 : char szResult[80] = {'\0'};
1383 56 : if (poDS->m_oSRS.IsGeographic() &&
1384 19 : EQUAL(OSR_GDS(szResult, sizeof(szResult), papszLines, "Units",
1385 : ""),
1386 : "DS"))
1387 : {
1388 0 : poDS->m_gt[0] /= 3600.0;
1389 0 : poDS->m_gt[1] /= 3600.0;
1390 0 : poDS->m_gt[2] /= 3600.0;
1391 0 : poDS->m_gt[3] /= 3600.0;
1392 0 : poDS->m_gt[4] /= 3600.0;
1393 0 : poDS->m_gt[5] /= 3600.0;
1394 : }
1395 : }
1396 : else
1397 : {
1398 2 : poDS->m_oSRS.Clear();
1399 : }
1400 :
1401 39 : CSLDestroy(papszLines);
1402 : }
1403 : else
1404 : {
1405 : // Check for IMAGE.REP (Spatiocarte Defense 1.0) or name_of_image.rep
1406 : // if it is a GIS-GeoSPOT image
1407 : // For the specification of SPDF (in French), see
1408 : // http://eden.ign.fr/download/pub/doc/emabgi/spdf10.pdf/download
1409 : const CPLString szImageRepFilename =
1410 162 : GetImageRepFilename(poOpenInfo->pszFilename);
1411 81 : if (!szImageRepFilename.empty())
1412 : {
1413 0 : fp = VSIFOpenL(szImageRepFilename.c_str(), "r");
1414 : }
1415 81 : if (fp != nullptr)
1416 : {
1417 0 : bool bUTM = false;
1418 0 : bool bWGS84 = false;
1419 0 : int bNorth = FALSE;
1420 0 : bool bSouth = false;
1421 0 : int utmZone = 0;
1422 :
1423 0 : while ((pszLine = CPLReadLineL(fp)) != nullptr)
1424 : {
1425 0 : if (STARTS_WITH(pszLine, "PROJ_ID") && strstr(pszLine, "UTM"))
1426 : {
1427 0 : bUTM = true;
1428 : }
1429 0 : else if (STARTS_WITH(pszLine, "PROJ_ZONE"))
1430 : {
1431 0 : const char *c = strchr(pszLine, '"');
1432 0 : if (c)
1433 : {
1434 0 : c++;
1435 0 : if (*c >= '0' && *c <= '9')
1436 : {
1437 0 : utmZone = atoi(c);
1438 0 : if (utmZone >= 1 && utmZone <= 60)
1439 : {
1440 0 : if (strstr(pszLine, "Nord") ||
1441 0 : strstr(pszLine, "NORD"))
1442 : {
1443 0 : bNorth = TRUE;
1444 : }
1445 0 : else if (strstr(pszLine, "Sud") ||
1446 0 : strstr(pszLine, "SUD"))
1447 : {
1448 0 : bSouth = true;
1449 : }
1450 : }
1451 : }
1452 : }
1453 : }
1454 0 : else if (STARTS_WITH(pszLine, "PROJ_CODE") &&
1455 0 : strstr(pszLine, "FR-MINDEF"))
1456 : {
1457 0 : const char *c = strchr(pszLine, 'A');
1458 0 : if (c)
1459 : {
1460 0 : c++;
1461 0 : if (*c >= '0' && *c <= '9')
1462 : {
1463 0 : utmZone = atoi(c);
1464 0 : if (utmZone >= 1 && utmZone <= 60)
1465 : {
1466 0 : if (c[1] == 'N' ||
1467 0 : (c[1] != '\0' && c[2] == 'N'))
1468 : {
1469 0 : bNorth = TRUE;
1470 : }
1471 0 : else if (c[1] == 'S' ||
1472 0 : (c[1] != '\0' && c[2] == 'S'))
1473 : {
1474 0 : bSouth = true;
1475 : }
1476 : }
1477 : }
1478 0 : }
1479 : }
1480 0 : else if (STARTS_WITH(pszLine, "HORIZ_DATUM") &&
1481 0 : (strstr(pszLine, "WGS 84") ||
1482 0 : strstr(pszLine, "WGS84")))
1483 : {
1484 0 : bWGS84 = true;
1485 : }
1486 0 : else if (STARTS_WITH(pszLine, "MAP_NUMBER"))
1487 : {
1488 0 : const char *c = strchr(pszLine, '"');
1489 0 : if (c)
1490 : {
1491 0 : char *pszMapNumber = CPLStrdup(c + 1);
1492 0 : char *c2 = strchr(pszMapNumber, '"');
1493 0 : if (c2)
1494 0 : *c2 = 0;
1495 0 : poDS->SetMetadataItem("SPDF_MAP_NUMBER", pszMapNumber);
1496 0 : CPLFree(pszMapNumber);
1497 : }
1498 : }
1499 0 : else if (STARTS_WITH(pszLine, "PRODUCTION_DATE"))
1500 : {
1501 0 : const char *c = pszLine + strlen("PRODUCTION_DATE");
1502 0 : while (*c == ' ')
1503 0 : c++;
1504 0 : if (*c)
1505 : {
1506 0 : poDS->SetMetadataItem("SPDF_PRODUCTION_DATE", c);
1507 : }
1508 : }
1509 : }
1510 :
1511 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
1512 :
1513 0 : if (utmZone >= 1 && utmZone <= 60 && bUTM && bWGS84 &&
1514 0 : (bNorth || bSouth))
1515 : {
1516 0 : char projCSStr[64] = {'\0'};
1517 0 : snprintf(projCSStr, sizeof(projCSStr), "WGS 84 / UTM zone %d%c",
1518 : utmZone, (bNorth) ? 'N' : 'S');
1519 :
1520 0 : poDS->m_oSRS.SetProjCS(projCSStr);
1521 0 : poDS->m_oSRS.SetWellKnownGeogCS("WGS84");
1522 0 : poDS->m_oSRS.SetUTM(utmZone, bNorth);
1523 0 : poDS->m_oSRS.SetAuthority("PROJCS", "EPSG",
1524 0 : (bNorth ? 32600 : 32700) + utmZone);
1525 0 : poDS->m_oSRS.AutoIdentifyEPSG();
1526 : }
1527 : else
1528 : {
1529 0 : CPLError(CE_Warning, CPLE_NotSupported,
1530 : "Cannot retrieve projection from IMAGE.REP");
1531 : }
1532 : }
1533 : }
1534 :
1535 : // Check for a color table.
1536 : const std::string osCLRFilename =
1537 240 : CPLFormCIFilenameSafe(osPath, osName, "clr");
1538 :
1539 : // Only read the .clr for byte, int16 or uint16 bands.
1540 120 : if (nItemSize <= 2)
1541 86 : fp = VSIFOpenL(osCLRFilename.c_str(), "r");
1542 : else
1543 34 : fp = nullptr;
1544 :
1545 120 : if (fp != nullptr)
1546 : {
1547 : std::shared_ptr<GDALRasterAttributeTable> poRat(
1548 7 : new GDALDefaultRasterAttributeTable());
1549 7 : poRat->CreateColumn("Value", GFT_Integer, GFU_Generic);
1550 7 : poRat->CreateColumn("Red", GFT_Integer, GFU_Red);
1551 7 : poRat->CreateColumn("Green", GFT_Integer, GFU_Green);
1552 7 : poRat->CreateColumn("Blue", GFT_Integer, GFU_Blue);
1553 :
1554 7 : poDS->m_poColorTable.reset(new GDALColorTable());
1555 :
1556 7 : bool bHasFoundNonCTValues = false;
1557 7 : int nRatRow = 0;
1558 :
1559 : while (true)
1560 : {
1561 94 : pszLine = CPLReadLineL(fp);
1562 94 : if (!pszLine)
1563 7 : break;
1564 :
1565 87 : if (*pszLine == '#' || *pszLine == '!')
1566 0 : continue;
1567 :
1568 : char **papszValues =
1569 87 : CSLTokenizeString2(pszLine, "\t ", CSLT_HONOURSTRINGS);
1570 :
1571 87 : if (CSLCount(papszValues) >= 4)
1572 : {
1573 87 : const int nIndex = atoi(papszValues[0]);
1574 87 : poRat->SetValue(nRatRow, 0, nIndex);
1575 87 : poRat->SetValue(nRatRow, 1, atoi(papszValues[1]));
1576 87 : poRat->SetValue(nRatRow, 2, atoi(papszValues[2]));
1577 87 : poRat->SetValue(nRatRow, 3, atoi(papszValues[3]));
1578 87 : nRatRow++;
1579 :
1580 87 : if (nIndex >= 0 && nIndex < 65536)
1581 : {
1582 72 : const GDALColorEntry oEntry = {
1583 72 : static_cast<short>(atoi(papszValues[1])), // Red
1584 72 : static_cast<short>(atoi(papszValues[2])), // Green
1585 72 : static_cast<short>(atoi(papszValues[3])), // Blue
1586 72 : 255};
1587 :
1588 72 : poDS->m_poColorTable->SetColorEntry(nIndex, &oEntry);
1589 : }
1590 : else
1591 : {
1592 : // Negative values are valid. At least we can find use of
1593 : // them here:
1594 : // http://www.ngdc.noaa.gov/mgg/topo/elev/esri/clr/
1595 : // But, there's no way of representing them with GDAL color
1596 : // table model.
1597 15 : if (!bHasFoundNonCTValues)
1598 3 : CPLDebug("EHdr", "Ignoring color index : %d", nIndex);
1599 15 : bHasFoundNonCTValues = true;
1600 : }
1601 : }
1602 :
1603 87 : CSLDestroy(papszValues);
1604 87 : }
1605 :
1606 7 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
1607 :
1608 7 : if (bHasFoundNonCTValues)
1609 : {
1610 3 : poDS->m_poRAT.swap(poRat);
1611 : }
1612 :
1613 14 : for (int i = 1; i <= poDS->nBands; i++)
1614 : {
1615 : EHdrRasterBand *poBand =
1616 7 : cpl::down_cast<EHdrRasterBand *>(poDS->GetRasterBand(i));
1617 7 : poBand->m_poColorTable = poDS->m_poColorTable;
1618 7 : poBand->m_poRAT = poDS->m_poRAT;
1619 7 : poBand->SetColorInterpretation(GCI_PaletteIndex);
1620 : }
1621 :
1622 7 : poDS->bCLRDirty = false;
1623 : }
1624 :
1625 : // Read statistics (.STX).
1626 120 : poDS->ReadSTX();
1627 :
1628 : // Initialize any PAM information.
1629 120 : poDS->TryLoadXML();
1630 :
1631 : // Check for overviews.
1632 120 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
1633 :
1634 120 : return poDS.release();
1635 : }
1636 :
1637 : /************************************************************************/
1638 : /* Create() */
1639 : /************************************************************************/
1640 :
1641 79 : GDALDataset *EHdrDataset::Create(const char *pszFilename, int nXSize,
1642 : int nYSize, int nBandsIn, GDALDataType eType,
1643 : char **papszParamList)
1644 :
1645 : {
1646 : // Verify input options.
1647 79 : if (nBandsIn <= 0)
1648 : {
1649 1 : CPLError(CE_Failure, CPLE_NotSupported,
1650 : "EHdr driver does not support %d bands.", nBandsIn);
1651 1 : return nullptr;
1652 : }
1653 :
1654 78 : if (eType != GDT_Byte && eType != GDT_Int8 && eType != GDT_Float32 &&
1655 30 : eType != GDT_UInt16 && eType != GDT_Int16 && eType != GDT_Int32 &&
1656 : eType != GDT_UInt32)
1657 : {
1658 19 : CPLError(CE_Failure, CPLE_AppDefined,
1659 : "Attempt to create ESRI .hdr labelled dataset with an illegal"
1660 : "data type (%s).",
1661 : GDALGetDataTypeName(eType));
1662 :
1663 19 : return nullptr;
1664 : }
1665 :
1666 : // Try to create the file.
1667 59 : VSILFILE *fp = VSIFOpenL(pszFilename, "wb");
1668 :
1669 59 : if (fp == nullptr)
1670 : {
1671 3 : CPLError(CE_Failure, CPLE_OpenFailed,
1672 : "Attempt to create file `%s' failed.", pszFilename);
1673 3 : return nullptr;
1674 : }
1675 :
1676 : // Just write out a couple of bytes to establish the binary
1677 : // file, and then close it.
1678 56 : bool bOK = VSIFWriteL(reinterpret_cast<void *>(const_cast<char *>("\0\0")),
1679 56 : 2, 1, fp) == 1;
1680 56 : if (VSIFCloseL(fp) != 0)
1681 0 : bOK = false;
1682 56 : fp = nullptr;
1683 56 : if (!bOK)
1684 1 : return nullptr;
1685 :
1686 : // Create the hdr filename.
1687 110 : const std::string osHdrFilename = CPLResetExtensionSafe(pszFilename, "hdr");
1688 :
1689 : // Open the file.
1690 55 : fp = VSIFOpenL(osHdrFilename.c_str(), "wt");
1691 55 : if (fp == nullptr)
1692 : {
1693 0 : CPLError(CE_Failure, CPLE_OpenFailed,
1694 : "Attempt to create file `%s' failed.", osHdrFilename.c_str());
1695 0 : return nullptr;
1696 : }
1697 :
1698 : // Decide how many bits the file should have.
1699 55 : const char *pszNBITS = CSLFetchNameValue(papszParamList, "NBITS");
1700 : const int nBits =
1701 55 : pszNBITS ? atoi(pszNBITS) : GDALGetDataTypeSizeBits(eType);
1702 55 : if (nBits <= 0 || nXSize > (static_cast<int64_t>(INT_MAX) * 8 - 7) / nBits)
1703 : {
1704 0 : CPLError(CE_Failure, CPLE_AppDefined,
1705 : "Invalid NBITS or too large image width");
1706 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
1707 0 : return nullptr;
1708 : }
1709 :
1710 55 : const int nRowBytes = (nBits * nXSize + 7) / 8;
1711 :
1712 55 : if (nRowBytes > INT_MAX / nBandsIn)
1713 : {
1714 0 : CPLError(CE_Failure, CPLE_AppDefined,
1715 : "Too large band count and/or image width");
1716 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
1717 0 : return nullptr;
1718 : }
1719 :
1720 : // Check for signed byte.
1721 55 : const char *pszPixelType = CSLFetchNameValue(papszParamList, "PIXELTYPE");
1722 55 : if (pszPixelType == nullptr)
1723 54 : pszPixelType = "";
1724 :
1725 : // Write out the raw definition for the dataset as a whole.
1726 55 : bOK &= VSIFPrintfL(fp, "BYTEORDER I\n") >= 0;
1727 55 : bOK &= VSIFPrintfL(fp, "LAYOUT BIL\n") >= 0;
1728 55 : bOK &= VSIFPrintfL(fp, "NROWS %d\n", nYSize) >= 0;
1729 55 : bOK &= VSIFPrintfL(fp, "NCOLS %d\n", nXSize) >= 0;
1730 55 : bOK &= VSIFPrintfL(fp, "NBANDS %d\n", nBandsIn) >= 0;
1731 55 : bOK &= VSIFPrintfL(fp, "NBITS %d\n", nBits) >= 0;
1732 55 : bOK &= VSIFPrintfL(fp, "BANDROWBYTES %d\n", nRowBytes) >= 0;
1733 55 : bOK &= VSIFPrintfL(fp, "TOTALROWBYTES %d\n", nRowBytes * nBandsIn) >= 0;
1734 :
1735 55 : if (eType == GDT_Float32)
1736 5 : bOK &= VSIFPrintfL(fp, "PIXELTYPE FLOAT\n") >= 0;
1737 50 : else if (eType == GDT_Int8 || eType == GDT_Int16 || eType == GDT_Int32)
1738 10 : bOK &= VSIFPrintfL(fp, "PIXELTYPE SIGNEDINT\n") >= 0;
1739 40 : else if (eType == GDT_Byte && EQUAL(pszPixelType, "SIGNEDBYTE"))
1740 1 : bOK &= VSIFPrintfL(fp, "PIXELTYPE SIGNEDINT\n") >= 0;
1741 : else
1742 39 : bOK &= VSIFPrintfL(fp, "PIXELTYPE UNSIGNEDINT\n") >= 0;
1743 :
1744 55 : if (VSIFCloseL(fp) != 0)
1745 0 : bOK = false;
1746 :
1747 55 : if (!bOK)
1748 0 : return nullptr;
1749 :
1750 110 : GDALOpenInfo oOpenInfo(pszFilename, GA_Update);
1751 55 : return Open(&oOpenInfo, false);
1752 : }
1753 :
1754 : /************************************************************************/
1755 : /* CreateCopy() */
1756 : /************************************************************************/
1757 :
1758 40 : GDALDataset *EHdrDataset::CreateCopy(const char *pszFilename,
1759 : GDALDataset *poSrcDS, int bStrict,
1760 : char **papszOptions,
1761 : GDALProgressFunc pfnProgress,
1762 : void *pProgressData)
1763 :
1764 : {
1765 40 : const int nBands = poSrcDS->GetRasterCount();
1766 40 : if (nBands == 0)
1767 : {
1768 1 : CPLError(CE_Failure, CPLE_NotSupported,
1769 : "EHdr driver does not support source dataset without any "
1770 : "bands.");
1771 1 : return nullptr;
1772 : }
1773 :
1774 39 : char **papszAdjustedOptions = CSLDuplicate(papszOptions);
1775 :
1776 : // Ensure we pass on NBITS and PIXELTYPE structure information.
1777 39 : auto poSrcBand = poSrcDS->GetRasterBand(1);
1778 39 : if (poSrcBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE") != nullptr &&
1779 0 : CSLFetchNameValue(papszOptions, "NBITS") == nullptr)
1780 : {
1781 0 : papszAdjustedOptions = CSLSetNameValue(
1782 : papszAdjustedOptions, "NBITS",
1783 0 : poSrcBand->GetMetadataItem("NBITS", "IMAGE_STRUCTURE"));
1784 : }
1785 :
1786 64 : if (poSrcBand->GetRasterDataType() == GDT_Byte &&
1787 25 : CSLFetchNameValue(papszOptions, "PIXELTYPE") == nullptr)
1788 : {
1789 25 : poSrcBand->EnablePixelTypeSignedByteWarning(false);
1790 : const char *pszPixelType =
1791 25 : poSrcBand->GetMetadataItem("PIXELTYPE", "IMAGE_STRUCTURE");
1792 25 : poSrcBand->EnablePixelTypeSignedByteWarning(true);
1793 25 : if (pszPixelType != nullptr)
1794 : {
1795 1 : papszAdjustedOptions = CSLSetNameValue(papszAdjustedOptions,
1796 : "PIXELTYPE", pszPixelType);
1797 : }
1798 : }
1799 :
1800 : // Proceed with normal copying using the default createcopy operators.
1801 : GDALDriver *poDriver =
1802 39 : reinterpret_cast<GDALDriver *>(GDALGetDriverByName("EHdr"));
1803 :
1804 39 : GDALDataset *poOutDS = poDriver->DefaultCreateCopy(
1805 : pszFilename, poSrcDS, bStrict, papszAdjustedOptions, pfnProgress,
1806 : pProgressData);
1807 39 : CSLDestroy(papszAdjustedOptions);
1808 :
1809 39 : if (poOutDS != nullptr)
1810 21 : poOutDS->FlushCache(false);
1811 :
1812 39 : return poOutDS;
1813 : }
1814 :
1815 : /************************************************************************/
1816 : /* GetNoDataValue() */
1817 : /************************************************************************/
1818 :
1819 102 : double EHdrRasterBand::GetNoDataValue(int *pbSuccess)
1820 : {
1821 102 : if (pbSuccess)
1822 102 : *pbSuccess = bNoDataSet;
1823 :
1824 102 : if (bNoDataSet)
1825 1 : return dfNoData;
1826 :
1827 101 : return RawRasterBand::GetNoDataValue(pbSuccess);
1828 : }
1829 :
1830 : /************************************************************************/
1831 : /* GetMinimum() */
1832 : /************************************************************************/
1833 :
1834 3 : double EHdrRasterBand::GetMinimum(int *pbSuccess)
1835 : {
1836 3 : if (pbSuccess != nullptr)
1837 3 : *pbSuccess = (minmaxmeanstddev & HAS_MIN_FLAG) != 0;
1838 :
1839 3 : if (minmaxmeanstddev & HAS_MIN_FLAG)
1840 2 : return dfMin;
1841 :
1842 1 : return RawRasterBand::GetMinimum(pbSuccess);
1843 : }
1844 :
1845 : /************************************************************************/
1846 : /* GetMaximum() */
1847 : /************************************************************************/
1848 :
1849 2 : double EHdrRasterBand::GetMaximum(int *pbSuccess)
1850 : {
1851 2 : if (pbSuccess != nullptr)
1852 2 : *pbSuccess = (minmaxmeanstddev & HAS_MAX_FLAG) != 0;
1853 :
1854 2 : if (minmaxmeanstddev & HAS_MAX_FLAG)
1855 1 : return dfMax;
1856 :
1857 1 : return RawRasterBand::GetMaximum(pbSuccess);
1858 : }
1859 :
1860 : /************************************************************************/
1861 : /* GetStatistics() */
1862 : /************************************************************************/
1863 :
1864 6 : CPLErr EHdrRasterBand::GetStatistics(int bApproxOK, int bForce, double *pdfMin,
1865 : double *pdfMax, double *pdfMean,
1866 : double *pdfStdDev)
1867 : {
1868 6 : if (!(GetMetadataItem("STATISTICS_APPROXIMATE") && !bApproxOK))
1869 : {
1870 5 : if ((minmaxmeanstddev & HAS_ALL_FLAGS) == HAS_ALL_FLAGS)
1871 : {
1872 1 : if (pdfMin)
1873 1 : *pdfMin = dfMin;
1874 1 : if (pdfMax)
1875 1 : *pdfMax = dfMax;
1876 1 : if (pdfMean)
1877 1 : *pdfMean = dfMean;
1878 1 : if (pdfStdDev)
1879 1 : *pdfStdDev = dfStdDev;
1880 1 : return CE_None;
1881 : }
1882 : }
1883 :
1884 5 : const CPLErr eErr = RawRasterBand::GetStatistics(
1885 : bApproxOK, bForce, &dfMin, &dfMax, &dfMean, &dfStdDev);
1886 5 : if (eErr != CE_None)
1887 2 : return eErr;
1888 :
1889 3 : EHdrDataset *poEDS = reinterpret_cast<EHdrDataset *>(poDS);
1890 :
1891 3 : minmaxmeanstddev = HAS_ALL_FLAGS;
1892 :
1893 3 : if (!bApproxOK && poEDS->RewriteSTX() != CE_None)
1894 0 : RawRasterBand::SetStatistics(dfMin, dfMax, dfMean, dfStdDev);
1895 :
1896 3 : if (pdfMin)
1897 3 : *pdfMin = dfMin;
1898 3 : if (pdfMax)
1899 3 : *pdfMax = dfMax;
1900 3 : if (pdfMean)
1901 3 : *pdfMean = dfMean;
1902 3 : if (pdfStdDev)
1903 3 : *pdfStdDev = dfStdDev;
1904 :
1905 3 : return CE_None;
1906 : }
1907 :
1908 : /************************************************************************/
1909 : /* SetStatistics() */
1910 : /************************************************************************/
1911 :
1912 3 : CPLErr EHdrRasterBand::SetStatistics(double dfMinIn, double dfMaxIn,
1913 : double dfMeanIn, double dfStdDevIn)
1914 : {
1915 : // Avoid churn if nothing is changing.
1916 3 : if (dfMin == dfMinIn && dfMax == dfMaxIn && dfMean == dfMeanIn &&
1917 1 : dfStdDev == dfStdDevIn)
1918 1 : return CE_None;
1919 :
1920 2 : dfMin = dfMinIn;
1921 2 : dfMax = dfMaxIn;
1922 2 : dfMean = dfMeanIn;
1923 2 : dfStdDev = dfStdDevIn;
1924 :
1925 : // marks stats valid
1926 2 : minmaxmeanstddev = HAS_ALL_FLAGS;
1927 :
1928 2 : EHdrDataset *poEDS = reinterpret_cast<EHdrDataset *>(poDS);
1929 :
1930 2 : if (GetMetadataItem("STATISTICS_APPROXIMATE") == nullptr)
1931 : {
1932 2 : if (GetMetadataItem("STATISTICS_MINIMUM"))
1933 : {
1934 0 : SetMetadataItem("STATISTICS_MINIMUM", nullptr);
1935 0 : SetMetadataItem("STATISTICS_MAXIMUM", nullptr);
1936 0 : SetMetadataItem("STATISTICS_MEAN", nullptr);
1937 0 : SetMetadataItem("STATISTICS_STDDEV", nullptr);
1938 : }
1939 2 : return poEDS->RewriteSTX();
1940 : }
1941 :
1942 0 : return RawRasterBand::SetStatistics(dfMinIn, dfMaxIn, dfMeanIn, dfStdDevIn);
1943 : }
1944 :
1945 : /************************************************************************/
1946 : /* GetColorTable() */
1947 : /************************************************************************/
1948 :
1949 13 : GDALColorTable *EHdrRasterBand::GetColorTable()
1950 : {
1951 13 : return m_poColorTable.get();
1952 : }
1953 :
1954 : /************************************************************************/
1955 : /* SetColorTable() */
1956 : /************************************************************************/
1957 :
1958 6 : CPLErr EHdrRasterBand::SetColorTable(GDALColorTable *poNewCT)
1959 : {
1960 6 : if (poNewCT == nullptr)
1961 2 : m_poColorTable.reset();
1962 : else
1963 4 : m_poColorTable.reset(poNewCT->Clone());
1964 :
1965 6 : reinterpret_cast<EHdrDataset *>(poDS)->bCLRDirty = true;
1966 :
1967 6 : return CE_None;
1968 : }
1969 :
1970 : /************************************************************************/
1971 : /* GetDefaultRAT() */
1972 : /************************************************************************/
1973 :
1974 13 : GDALRasterAttributeTable *EHdrRasterBand::GetDefaultRAT()
1975 : {
1976 13 : return m_poRAT.get();
1977 : }
1978 :
1979 : /************************************************************************/
1980 : /* SetDefaultRAT() */
1981 : /************************************************************************/
1982 :
1983 3 : CPLErr EHdrRasterBand::SetDefaultRAT(const GDALRasterAttributeTable *poRAT)
1984 : {
1985 3 : if (poRAT)
1986 : {
1987 3 : if (!(poRAT->GetColumnCount() == 4 &&
1988 1 : poRAT->GetTypeOfCol(0) == GFT_Integer &&
1989 1 : poRAT->GetTypeOfCol(1) == GFT_Integer &&
1990 1 : poRAT->GetTypeOfCol(2) == GFT_Integer &&
1991 1 : poRAT->GetTypeOfCol(3) == GFT_Integer &&
1992 1 : poRAT->GetUsageOfCol(0) == GFU_Generic &&
1993 1 : poRAT->GetUsageOfCol(1) == GFU_Red &&
1994 1 : poRAT->GetUsageOfCol(2) == GFU_Green &&
1995 1 : poRAT->GetUsageOfCol(3) == GFU_Blue))
1996 : {
1997 1 : CPLError(CE_Warning, CPLE_NotSupported,
1998 : "Unsupported type of RAT: "
1999 : "only value,R,G,B ones are supported");
2000 1 : return CE_Failure;
2001 : }
2002 : }
2003 :
2004 2 : if (poRAT == nullptr)
2005 1 : m_poRAT.reset();
2006 : else
2007 1 : m_poRAT.reset(poRAT->Clone());
2008 :
2009 2 : reinterpret_cast<EHdrDataset *>(poDS)->bCLRDirty = true;
2010 :
2011 2 : return CE_None;
2012 : }
2013 :
2014 : /************************************************************************/
2015 : /* GDALRegister_EHdr() */
2016 : /************************************************************************/
2017 :
2018 1935 : void GDALRegister_EHdr()
2019 :
2020 : {
2021 1935 : if (GDALGetDriverByName("EHdr") != nullptr)
2022 282 : return;
2023 :
2024 1653 : GDALDriver *poDriver = new GDALDriver();
2025 :
2026 1653 : poDriver->SetDescription("EHdr");
2027 1653 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
2028 1653 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "ESRI .hdr Labelled");
2029 1653 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/ehdr.html");
2030 1653 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "bil");
2031 1653 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
2032 1653 : "Byte Int8 Int16 UInt16 Int32 UInt32 Float32");
2033 :
2034 1653 : poDriver->SetMetadataItem(
2035 : GDAL_DMD_CREATIONOPTIONLIST,
2036 : "<CreationOptionList>"
2037 : " <Option name='NBITS' type='int' description='Special pixel bits "
2038 : "(1-7)'/>"
2039 : " <Option name='PIXELTYPE' type='string' description='By setting "
2040 : "this to SIGNEDBYTE, a new Byte file can be forced to be written as "
2041 : "signed byte'/>"
2042 1653 : "</CreationOptionList>");
2043 :
2044 1653 : poDriver->SetMetadataItem(GDAL_DCAP_UPDATE, "YES");
2045 1653 : poDriver->SetMetadataItem(GDAL_DMD_UPDATE_ITEMS, "GeoTransform SRS NoData "
2046 1653 : "RasterValues");
2047 :
2048 1653 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
2049 1653 : poDriver->pfnOpen = EHdrDataset::Open;
2050 1653 : poDriver->pfnCreate = EHdrDataset::Create;
2051 1653 : poDriver->pfnCreateCopy = EHdrDataset::CreateCopy;
2052 :
2053 1653 : GetGDALDriverManager()->RegisterDriver(poDriver);
2054 : }
|
