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