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