Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: Implements R Raster Format.
5 : * Author: Even Rouault, <even dot rouault at spatialys dot com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2016, Even Rouault, <even dot rouault at spatialys dot com>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "cpl_port.h"
14 :
15 : #include "gdal_frmts.h"
16 : #include "gdal_rat.h"
17 : #include "gdal_priv.h"
18 :
19 : #include "rawdataset.h"
20 : #include "ogr_spatialref.h"
21 :
22 : #include <algorithm>
23 : #include <cmath>
24 : #include <limits>
25 : #include <memory>
26 :
27 : /************************************************************************/
28 : /* ==================================================================== */
29 : /* RRASTERDataset */
30 : /* ==================================================================== */
31 : /************************************************************************/
32 :
33 : class RRASTERDataset final : public RawDataset
34 : {
35 : bool m_bHeaderDirty = false;
36 : CPLString m_osGriFilename{};
37 : bool m_bGeoTransformValid = false;
38 : GDALGeoTransform m_gt{};
39 : VSILFILE *m_fpImage = nullptr;
40 : OGRSpatialReference m_oSRS{};
41 : std::shared_ptr<GDALRasterAttributeTable> m_poRAT{};
42 : std::shared_ptr<GDALColorTable> m_poCT{};
43 : bool m_bNativeOrder = true;
44 : CPLString m_osCreator{};
45 : CPLString m_osCreated{};
46 : CPLString m_osBandOrder{};
47 : CPLString m_osLegend{};
48 : bool m_bInitRaster = false;
49 : bool m_bSignedByte = false;
50 :
51 : static bool ComputeSpacings(const CPLString &osBandOrder, int nCols,
52 : int nRows, int l_nBands, GDALDataType eDT,
53 : int &nPixelOffset, int &nLineOffset,
54 : vsi_l_offset &nBandOffset);
55 : void RewriteHeader();
56 :
57 : CPL_DISALLOW_COPY_ASSIGN(RRASTERDataset)
58 :
59 : CPLErr Close() override;
60 :
61 : public:
62 : RRASTERDataset();
63 : ~RRASTERDataset() override;
64 :
65 : char **GetFileList() override;
66 :
67 : static GDALDataset *Open(GDALOpenInfo *);
68 : static int Identify(GDALOpenInfo *);
69 : static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
70 : int nBandsIn, GDALDataType eType,
71 : char **papszOptions);
72 : static GDALDataset *CreateCopy(const char *pszFilename,
73 : GDALDataset *poSrcDS, int bStrict,
74 : char **papszOptions,
75 : GDALProgressFunc pfnProgress,
76 : void *pProgressData);
77 :
78 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
79 : CPLErr SetGeoTransform(const GDALGeoTransform >) override;
80 : const OGRSpatialReference *GetSpatialRef() const override;
81 : CPLErr SetSpatialRef(const OGRSpatialReference *poSRS) override;
82 :
83 : CPLErr SetMetadata(char **papszMetadata,
84 : const char *pszDomain = "") override;
85 : CPLErr SetMetadataItem(const char *pszName, const char *pszValue,
86 : const char *pszDomain = "") override;
87 :
88 179 : void SetHeaderDirty()
89 : {
90 179 : m_bHeaderDirty = true;
91 179 : }
92 :
93 : void InitImageIfNeeded();
94 :
95 211 : inline bool IsSignedByte() const
96 : {
97 211 : return m_bSignedByte;
98 : }
99 : };
100 :
101 : /************************************************************************/
102 : /* ==================================================================== */
103 : /* RRASTERRasterBand */
104 : /* ==================================================================== */
105 : /************************************************************************/
106 :
107 : class RRASTERRasterBand final : public RawRasterBand
108 : {
109 : friend class RRASTERDataset;
110 :
111 : bool m_bHasNoDataValue = false;
112 : double m_dfNoDataValue = 0.0;
113 : double m_dfMin = std::numeric_limits<double>::infinity();
114 : double m_dfMax = -std::numeric_limits<double>::infinity();
115 : std::shared_ptr<GDALRasterAttributeTable> m_poRAT{};
116 : std::shared_ptr<GDALColorTable> m_poCT{};
117 :
118 : CPL_DISALLOW_COPY_ASSIGN(RRASTERRasterBand)
119 :
120 : public:
121 : RRASTERRasterBand(GDALDataset *poDS, int nBand, VSILFILE *fpRaw,
122 : vsi_l_offset nImgOffset, int nPixelOffset,
123 : int nLineOffset, GDALDataType eDataType,
124 : int bNativeOrder);
125 :
126 : void SetMinMax(double dfMin, double dfMax);
127 : double GetMinimum(int *pbSuccess = nullptr) override;
128 : double GetMaximum(int *pbSuccess = nullptr) override;
129 :
130 : double GetNoDataValue(int *pbSuccess = nullptr) override;
131 : CPLErr SetNoDataValue(double dfNoData) override;
132 :
133 : GDALColorTable *GetColorTable() override;
134 : CPLErr SetColorTable(GDALColorTable *poNewCT) override;
135 :
136 : GDALRasterAttributeTable *GetDefaultRAT() override;
137 : CPLErr SetDefaultRAT(const GDALRasterAttributeTable *poRAT) override;
138 :
139 : void SetDescription(const char *pszDesc) override;
140 :
141 : protected:
142 : CPLErr IWriteBlock(int, int, void *) override;
143 : CPLErr IRasterIO(GDALRWFlag, int, int, int, int, void *, int, int,
144 : GDALDataType, GSpacing nPixelSpace, GSpacing nLineSpace,
145 : GDALRasterIOExtraArg *psExtraArg) override;
146 : };
147 :
148 : /************************************************************************/
149 : /* RRASTERDataset() */
150 : /************************************************************************/
151 :
152 262 : RRASTERRasterBand::RRASTERRasterBand(GDALDataset *poDSIn, int nBandIn,
153 : VSILFILE *fpRawIn,
154 : vsi_l_offset nImgOffsetIn,
155 : int nPixelOffsetIn, int nLineOffsetIn,
156 : GDALDataType eDataTypeIn,
157 262 : int bNativeOrderIn)
158 : : RawRasterBand(poDSIn, nBandIn, fpRawIn, nImgOffsetIn, nPixelOffsetIn,
159 : nLineOffsetIn, eDataTypeIn, bNativeOrderIn,
160 262 : RawRasterBand::OwnFP::NO)
161 : {
162 262 : }
163 :
164 : /************************************************************************/
165 : /* SetMinMax() */
166 : /************************************************************************/
167 :
168 74 : void RRASTERRasterBand::SetMinMax(double dfMin, double dfMax)
169 : {
170 74 : m_dfMin = dfMin;
171 74 : m_dfMax = dfMax;
172 74 : }
173 :
174 : /************************************************************************/
175 : /* GetMinimum() */
176 : /************************************************************************/
177 :
178 29 : double RRASTERRasterBand::GetMinimum(int *pbSuccess)
179 : {
180 29 : if (m_dfMin <= m_dfMax)
181 : {
182 22 : if (pbSuccess)
183 22 : *pbSuccess = TRUE;
184 22 : return m_dfMin;
185 : }
186 7 : return RawRasterBand::GetMinimum(pbSuccess);
187 : }
188 :
189 : /************************************************************************/
190 : /* GetMaximum() */
191 : /************************************************************************/
192 :
193 28 : double RRASTERRasterBand::GetMaximum(int *pbSuccess)
194 : {
195 28 : if (m_dfMin <= m_dfMax)
196 : {
197 21 : if (pbSuccess)
198 21 : *pbSuccess = TRUE;
199 21 : return m_dfMax;
200 : }
201 7 : return RawRasterBand::GetMaximum(pbSuccess);
202 : }
203 :
204 : /************************************************************************/
205 : /* GetColorTable() */
206 : /************************************************************************/
207 :
208 88 : GDALColorTable *RRASTERRasterBand::GetColorTable()
209 : {
210 88 : return m_poCT.get();
211 : }
212 :
213 : /************************************************************************/
214 : /* SetColorTable() */
215 : /************************************************************************/
216 :
217 3 : CPLErr RRASTERRasterBand::SetColorTable(GDALColorTable *poNewCT)
218 : {
219 3 : RRASTERDataset *poGDS = static_cast<RRASTERDataset *>(poDS);
220 3 : if (poGDS->GetAccess() != GA_Update)
221 0 : return CE_Failure;
222 :
223 3 : if (poNewCT == nullptr)
224 1 : m_poCT.reset();
225 : else
226 2 : m_poCT.reset(poNewCT->Clone());
227 :
228 3 : poGDS->SetHeaderDirty();
229 :
230 3 : return CE_None;
231 : }
232 :
233 : /************************************************************************/
234 : /* GetDefaultRAT() */
235 : /************************************************************************/
236 :
237 112 : GDALRasterAttributeTable *RRASTERRasterBand::GetDefaultRAT()
238 : {
239 112 : return m_poRAT.get();
240 : }
241 :
242 : /************************************************************************/
243 : /* SetDefaultRAT() */
244 : /************************************************************************/
245 :
246 2 : CPLErr RRASTERRasterBand::SetDefaultRAT(const GDALRasterAttributeTable *poRAT)
247 : {
248 2 : RRASTERDataset *poGDS = static_cast<RRASTERDataset *>(poDS);
249 2 : if (poGDS->GetAccess() != GA_Update)
250 0 : return CE_Failure;
251 :
252 2 : if (poRAT == nullptr)
253 1 : m_poRAT.reset();
254 : else
255 1 : m_poRAT.reset(poRAT->Clone());
256 :
257 2 : poGDS->SetHeaderDirty();
258 :
259 2 : return CE_None;
260 : }
261 :
262 : /************************************************************************/
263 : /* SetDescription() */
264 : /************************************************************************/
265 :
266 20 : void RRASTERRasterBand::SetDescription(const char *pszDesc)
267 : {
268 20 : RRASTERDataset *poGDS = static_cast<RRASTERDataset *>(poDS);
269 20 : if (poGDS->GetAccess() != GA_Update)
270 0 : return;
271 :
272 20 : GDALRasterBand::SetDescription(pszDesc);
273 :
274 20 : poGDS->SetHeaderDirty();
275 : }
276 :
277 : /************************************************************************/
278 : /* GetNoDataValue() */
279 : /************************************************************************/
280 :
281 392 : double RRASTERRasterBand::GetNoDataValue(int *pbSuccess)
282 : {
283 392 : if (pbSuccess)
284 392 : *pbSuccess = m_bHasNoDataValue;
285 392 : return m_dfNoDataValue;
286 : }
287 :
288 : /************************************************************************/
289 : /* SetNoDataValue() */
290 : /************************************************************************/
291 :
292 3 : CPLErr RRASTERRasterBand::SetNoDataValue(double dfNoData)
293 : {
294 3 : RRASTERDataset *poGDS = static_cast<RRASTERDataset *>(poDS);
295 3 : if (poGDS->GetAccess() != GA_Update)
296 0 : return CE_Failure;
297 :
298 3 : m_bHasNoDataValue = true;
299 3 : m_dfNoDataValue = dfNoData;
300 3 : poGDS->SetHeaderDirty();
301 3 : return CE_None;
302 : }
303 :
304 : /************************************************************************/
305 : /* GetMinMax() */
306 : /************************************************************************/
307 :
308 : template <class T>
309 211 : static void GetMinMax(const T *buffer, int nBufXSize, int nBufYSize,
310 : GSpacing nPixelSpace, GSpacing nLineSpace,
311 : double dfNoDataValue, double &dfMin, double &dfMax)
312 : {
313 1091 : for (int iY = 0; iY < nBufYSize; iY++)
314 : {
315 42541 : for (int iX = 0; iX < nBufXSize; iX++)
316 : {
317 41661 : const double dfVal = buffer[iY * nLineSpace + iX * nPixelSpace];
318 41661 : if (dfVal != dfNoDataValue && !std::isnan(dfVal))
319 : {
320 41660 : dfMin = std::min(dfMin, dfVal);
321 41660 : dfMax = std::max(dfMax, dfVal);
322 : }
323 : }
324 : }
325 211 : }
326 :
327 211 : static void GetMinMax(const void *pBuffer, GDALDataType eDT, int nBufXSize,
328 : int nBufYSize, GSpacing nPixelSpace, GSpacing nLineSpace,
329 : double dfNoDataValue, double &dfMin, double &dfMax)
330 : {
331 211 : switch (eDT)
332 : {
333 197 : case GDT_UInt8:
334 197 : GetMinMax(static_cast<const GByte *>(pBuffer), nBufXSize, nBufYSize,
335 : nPixelSpace, nLineSpace, dfNoDataValue, dfMin, dfMax);
336 197 : break;
337 2 : case GDT_Int8:
338 2 : GetMinMax(static_cast<const GInt8 *>(pBuffer), nBufXSize, nBufYSize,
339 : nPixelSpace, nLineSpace, dfNoDataValue, dfMin, dfMax);
340 2 : break;
341 2 : case GDT_UInt16:
342 2 : GetMinMax(static_cast<const GUInt16 *>(pBuffer), nBufXSize,
343 : nBufYSize, nPixelSpace, nLineSpace, dfNoDataValue, dfMin,
344 : dfMax);
345 2 : break;
346 2 : case GDT_Int16:
347 2 : GetMinMax(static_cast<const GInt16 *>(pBuffer), nBufXSize,
348 : nBufYSize, nPixelSpace, nLineSpace, dfNoDataValue, dfMin,
349 : dfMax);
350 2 : break;
351 2 : case GDT_UInt32:
352 2 : GetMinMax(static_cast<const GUInt32 *>(pBuffer), nBufXSize,
353 : nBufYSize, nPixelSpace, nLineSpace, dfNoDataValue, dfMin,
354 : dfMax);
355 2 : break;
356 2 : case GDT_Int32:
357 2 : GetMinMax(static_cast<const GInt32 *>(pBuffer), nBufXSize,
358 : nBufYSize, nPixelSpace, nLineSpace, dfNoDataValue, dfMin,
359 : dfMax);
360 2 : break;
361 2 : case GDT_Float32:
362 2 : GetMinMax(static_cast<const float *>(pBuffer), nBufXSize, nBufYSize,
363 : nPixelSpace, nLineSpace, dfNoDataValue, dfMin, dfMax);
364 2 : break;
365 2 : case GDT_Float64:
366 2 : GetMinMax(static_cast<const double *>(pBuffer), nBufXSize,
367 : nBufYSize, nPixelSpace, nLineSpace, dfNoDataValue, dfMin,
368 : dfMax);
369 2 : break;
370 0 : default:
371 0 : CPLAssert(false);
372 : break;
373 : }
374 211 : }
375 :
376 : /************************************************************************/
377 : /* IWriteBlock() */
378 : /************************************************************************/
379 :
380 149 : CPLErr RRASTERRasterBand::IWriteBlock(int nBlockXOff, int nBlockYOff,
381 : void *pImage)
382 : {
383 149 : RRASTERDataset *poGDS = static_cast<RRASTERDataset *>(poDS);
384 149 : poGDS->InitImageIfNeeded();
385 :
386 149 : int bGotNoDataValue = false;
387 149 : double dfNoDataValue = GetNoDataValue(&bGotNoDataValue);
388 149 : if (!bGotNoDataValue)
389 148 : dfNoDataValue = std::numeric_limits<double>::quiet_NaN();
390 149 : GetMinMax(pImage, poGDS->IsSignedByte() ? GDT_Int8 : eDataType, nBlockXSize,
391 149 : nBlockYSize, 1, nBlockXSize, dfNoDataValue, m_dfMin, m_dfMax);
392 298 : return RawRasterBand::IWriteBlock(nBlockXOff, nBlockYOff, pImage);
393 : }
394 :
395 : /************************************************************************/
396 : /* IRasterIO() */
397 : /************************************************************************/
398 :
399 325 : CPLErr RRASTERRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
400 : int nXSize, int nYSize, void *pData,
401 : int nBufXSize, int nBufYSize,
402 : GDALDataType eBufType, GSpacing nPixelSpace,
403 : GSpacing nLineSpace,
404 : GDALRasterIOExtraArg *psExtraArg)
405 :
406 : {
407 325 : if (eRWFlag == GF_Write)
408 : {
409 62 : RRASTERDataset *poGDS = static_cast<RRASTERDataset *>(poDS);
410 62 : poGDS->InitImageIfNeeded();
411 :
412 62 : const int nDTSize = std::max(1, GDALGetDataTypeSizeBytes(eDataType));
413 62 : int bGotNoDataValue = false;
414 62 : double dfNoDataValue = GetNoDataValue(&bGotNoDataValue);
415 62 : if (!bGotNoDataValue)
416 61 : dfNoDataValue = std::numeric_limits<double>::quiet_NaN();
417 62 : GetMinMax(pData, poGDS->IsSignedByte() ? GDT_Int8 : eDataType,
418 62 : nBufXSize, nBufYSize, nPixelSpace / nDTSize,
419 62 : nLineSpace / nDTSize, dfNoDataValue, m_dfMin, m_dfMax);
420 : }
421 325 : return RawRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
422 : pData, nBufXSize, nBufYSize, eBufType,
423 325 : nPixelSpace, nLineSpace, psExtraArg);
424 : }
425 :
426 : /************************************************************************/
427 : /* RRASTERDataset() */
428 : /************************************************************************/
429 :
430 144 : RRASTERDataset::RRASTERDataset()
431 : {
432 144 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
433 144 : }
434 :
435 : /************************************************************************/
436 : /* ~RRASTERDataset() */
437 : /************************************************************************/
438 :
439 288 : RRASTERDataset::~RRASTERDataset()
440 :
441 : {
442 144 : RRASTERDataset::Close();
443 288 : }
444 :
445 : /************************************************************************/
446 : /* Close() */
447 : /************************************************************************/
448 :
449 281 : CPLErr RRASTERDataset::Close()
450 : {
451 281 : CPLErr eErr = CE_None;
452 281 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
453 : {
454 144 : if (m_fpImage)
455 : {
456 144 : InitImageIfNeeded();
457 144 : if (RRASTERDataset::FlushCache(true) != CE_None)
458 2 : eErr = CE_Failure;
459 144 : if (VSIFCloseL(m_fpImage) != CE_None)
460 0 : eErr = CE_Failure;
461 : }
462 144 : if (m_bHeaderDirty)
463 61 : RewriteHeader();
464 :
465 144 : if (GDALPamDataset::Close() != CE_None)
466 0 : eErr = CE_Failure;
467 : }
468 281 : return eErr;
469 : }
470 :
471 : /************************************************************************/
472 : /* InitImageIfNeeded() */
473 : /************************************************************************/
474 :
475 355 : void RRASTERDataset::InitImageIfNeeded()
476 : {
477 355 : CPLAssert(m_fpImage);
478 355 : if (!m_bInitRaster)
479 334 : return;
480 21 : m_bInitRaster = false;
481 21 : int bGotNoDataValue = false;
482 21 : double dfNoDataValue = GetRasterBand(1)->GetNoDataValue(&bGotNoDataValue);
483 21 : const GDALDataType eDT = GetRasterBand(1)->GetRasterDataType();
484 21 : const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
485 21 : if (dfNoDataValue == 0.0)
486 : {
487 20 : VSIFTruncateL(m_fpImage, static_cast<vsi_l_offset>(nRasterXSize) *
488 20 : nRasterYSize * nBands * nDTSize);
489 : }
490 : else
491 : {
492 : GByte abyNoDataValue[16];
493 1 : GDALCopyWords(&dfNoDataValue, GDT_Float64, 0, abyNoDataValue, eDT, 0,
494 : 1);
495 2 : for (GUIntBig i = 0;
496 2 : i < static_cast<GUIntBig>(nRasterXSize) * nRasterYSize * nBands;
497 : i++)
498 : {
499 1 : VSIFWriteL(abyNoDataValue, 1, nDTSize, m_fpImage);
500 : }
501 : }
502 : }
503 :
504 : /************************************************************************/
505 : /* RewriteHeader() */
506 : /************************************************************************/
507 :
508 61 : void RRASTERDataset::RewriteHeader()
509 : {
510 61 : VSILFILE *fp = VSIFOpenL(GetDescription(), "wb");
511 61 : if (!fp)
512 0 : return;
513 :
514 61 : VSIFPrintfL(fp, "[general]\n");
515 61 : if (!m_osCreator.empty())
516 4 : VSIFPrintfL(fp, "creator=%s\n", m_osCreator.c_str());
517 61 : if (!m_osCreated.empty())
518 4 : VSIFPrintfL(fp, "created=%s\n", m_osCreated.c_str());
519 :
520 61 : VSIFPrintfL(fp, "[data]\n");
521 61 : GDALDataType eDT = GetRasterBand(1)->GetRasterDataType();
522 119 : VSIFPrintfL(fp, "datatype=%s\n",
523 58 : (eDT == GDT_Int8 || m_bSignedByte) ? "INT1S"
524 : : (eDT == GDT_UInt8) ? "INT1U"
525 : : (eDT == GDT_UInt16) ? "INT2U"
526 : : (eDT == GDT_UInt32) ? "INT4U"
527 : : (eDT == GDT_Int16) ? "INT2S"
528 : : (eDT == GDT_Int32) ? "INT4S"
529 : : (eDT == GDT_Float32) ? "FLT4S"
530 : :
531 : /*(eDT == GDT_Float64) ?*/ "FLT8S");
532 :
533 61 : int bGotNoDataValue = false;
534 61 : double dfNoDataValue = GetRasterBand(1)->GetNoDataValue(&bGotNoDataValue);
535 61 : if (bGotNoDataValue)
536 3 : VSIFPrintfL(fp, "nodatavalue=%.17g\n", dfNoDataValue);
537 :
538 61 : VSIFPrintfL(fp, "byteorder=%s\n",
539 61 : (m_bNativeOrder ^ CPL_IS_LSB) ? "big" : "little");
540 61 : VSIFPrintfL(fp, "nbands=%d\n", nBands);
541 61 : if (nBands > 1)
542 20 : VSIFPrintfL(fp, "bandorder=%s\n", m_osBandOrder.c_str());
543 122 : CPLString osMinValue, osMaxValue;
544 123 : for (int i = 1; i <= nBands; i++)
545 : {
546 : RRASTERRasterBand *poBand =
547 83 : static_cast<RRASTERRasterBand *>(GetRasterBand(i));
548 83 : if (i > 1)
549 : {
550 22 : osMinValue += ":";
551 22 : osMaxValue += ":";
552 : }
553 83 : if (poBand->m_dfMin > poBand->m_dfMax)
554 : {
555 21 : osMinValue.clear();
556 21 : break;
557 : }
558 62 : osMinValue += CPLSPrintf("%.17g", poBand->m_dfMin);
559 62 : osMaxValue += CPLSPrintf("%.17g", poBand->m_dfMax);
560 : }
561 61 : if (!osMinValue.empty())
562 : {
563 40 : VSIFPrintfL(fp, "minvalue=%s\n", osMinValue.c_str());
564 40 : VSIFPrintfL(fp, "maxvalue=%s\n", osMaxValue.c_str());
565 : }
566 :
567 61 : GDALColorTable *poCT = GetRasterBand(1)->GetColorTable();
568 61 : GDALRasterAttributeTable *poRAT = GetRasterBand(1)->GetDefaultRAT();
569 61 : if (poCT == nullptr && poRAT == nullptr)
570 : {
571 58 : VSIFPrintfL(fp, "categorical=FALSE\n");
572 : }
573 : else
574 : {
575 3 : VSIFPrintfL(fp, "categorical=TRUE\n");
576 3 : if (poCT && poRAT)
577 : {
578 0 : CPLError(CE_Warning, CPLE_AppDefined,
579 : "Both color table and raster attribute table defined. "
580 : "Writing only the later");
581 : }
582 3 : if (poRAT)
583 : {
584 2 : CPLString osRatNames;
585 2 : CPLString osRatTypes;
586 11 : for (int i = 0; i < poRAT->GetColumnCount(); i++)
587 : {
588 10 : if (!osRatNames.empty())
589 : {
590 9 : osRatNames += ":";
591 9 : osRatTypes += ":";
592 : }
593 : osRatNames +=
594 10 : CPLString(poRAT->GetNameOfCol(i)).replaceAll(':', '.');
595 10 : GDALRATFieldType eColType = poRAT->GetTypeOfCol(i);
596 10 : if (eColType == GFT_Integer)
597 7 : osRatTypes += "integer";
598 3 : else if (eColType == GFT_Real)
599 1 : osRatTypes += "numeric";
600 : else
601 2 : osRatTypes += "character";
602 : }
603 1 : VSIFPrintfL(fp, "ratnames=%s\n", osRatNames.c_str());
604 1 : VSIFPrintfL(fp, "rattypes=%s\n", osRatTypes.c_str());
605 2 : CPLString osRatValues;
606 11 : for (int i = 0; i < poRAT->GetColumnCount(); i++)
607 : {
608 10 : GDALRATFieldType eColType = poRAT->GetTypeOfCol(i);
609 30 : for (int j = 0; j < poRAT->GetRowCount(); j++)
610 : {
611 20 : if (i != 0 || j != 0)
612 19 : osRatValues += ":";
613 20 : if (eColType == GFT_Integer)
614 : {
615 : osRatValues +=
616 14 : CPLSPrintf("%d", poRAT->GetValueAsInt(j, i));
617 : }
618 6 : else if (eColType == GFT_Real)
619 : {
620 : osRatValues +=
621 2 : CPLSPrintf("%.17g", poRAT->GetValueAsDouble(j, i));
622 : }
623 : else
624 : {
625 4 : const char *pszVal = poRAT->GetValueAsString(j, i);
626 4 : if (pszVal)
627 : {
628 : osRatValues +=
629 4 : CPLString(pszVal).replaceAll(':', '.');
630 : }
631 : }
632 : }
633 : }
634 1 : VSIFPrintfL(fp, "ratvalues=%s\n", osRatValues.c_str());
635 : }
636 : else
637 : {
638 2 : bool bNeedsAlpha = false;
639 4 : for (int i = 0; i < poCT->GetColorEntryCount(); i++)
640 : {
641 3 : if (poCT->GetColorEntry(i)->c4 != 255)
642 : {
643 1 : bNeedsAlpha = true;
644 1 : break;
645 : }
646 : }
647 2 : if (!bNeedsAlpha)
648 : {
649 1 : VSIFPrintfL(fp, "ratnames=%s\n", "ID:red:green:blue");
650 1 : VSIFPrintfL(fp, "rattypes=%s\n",
651 : "integer:integer:integer:integer");
652 : }
653 : else
654 : {
655 1 : VSIFPrintfL(fp, "ratnames=%s\n", "ID:red:green:blue:alpha");
656 1 : VSIFPrintfL(fp, "rattypes=%s\n",
657 : "integer:integer:integer:integer:integer");
658 : }
659 :
660 4 : CPLString osRatID;
661 4 : CPLString osRatR;
662 4 : CPLString osRatG;
663 4 : CPLString osRatB;
664 4 : CPLString osRatA;
665 6 : for (int i = 0; i < poCT->GetColorEntryCount(); i++)
666 : {
667 4 : const GDALColorEntry *psEntry = poCT->GetColorEntry(i);
668 4 : if (i > 0)
669 : {
670 2 : osRatID += ":";
671 2 : osRatR += ":";
672 2 : osRatG += ":";
673 2 : osRatB += ":";
674 2 : osRatA += ":";
675 : }
676 4 : osRatID += CPLSPrintf("%d", i);
677 4 : osRatR += CPLSPrintf("%d", psEntry->c1);
678 4 : osRatG += CPLSPrintf("%d", psEntry->c2);
679 4 : osRatB += CPLSPrintf("%d", psEntry->c3);
680 4 : osRatA += CPLSPrintf("%d", psEntry->c4);
681 : }
682 2 : if (!bNeedsAlpha)
683 : {
684 1 : VSIFPrintfL(fp, "ratvalues=%s:%s:%s:%s\n", osRatID.c_str(),
685 : osRatR.c_str(), osRatG.c_str(), osRatB.c_str());
686 : }
687 : else
688 : {
689 1 : VSIFPrintfL(fp, "ratvalues=%s:%s:%s:%s:%s\n", osRatID.c_str(),
690 : osRatR.c_str(), osRatG.c_str(), osRatB.c_str(),
691 : osRatA.c_str());
692 : }
693 : }
694 : }
695 :
696 61 : if (!m_osLegend.empty())
697 0 : VSIFPrintfL(fp, "[legend]\n%s", m_osLegend.c_str());
698 :
699 122 : CPLString osLayerName;
700 61 : bool bGotSignificantBandDesc = false;
701 171 : for (int i = 1; i <= nBands; i++)
702 : {
703 110 : GDALRasterBand *poBand = GetRasterBand(i);
704 110 : const char *pszDesc = poBand->GetDescription();
705 110 : if (EQUAL(pszDesc, ""))
706 : {
707 90 : GDALColorInterp eInterp = poBand->GetColorInterpretation();
708 90 : if (eInterp == GCI_RedBand)
709 : {
710 1 : bGotSignificantBandDesc = true;
711 1 : pszDesc = "red";
712 : }
713 89 : else if (eInterp == GCI_GreenBand)
714 : {
715 1 : bGotSignificantBandDesc = true;
716 1 : pszDesc = "green";
717 : }
718 88 : else if (eInterp == GCI_BlueBand)
719 : {
720 1 : bGotSignificantBandDesc = true;
721 1 : pszDesc = "blue";
722 : }
723 87 : else if (eInterp == GCI_AlphaBand)
724 : {
725 1 : bGotSignificantBandDesc = true;
726 1 : pszDesc = "alpha";
727 : }
728 : else
729 : {
730 86 : pszDesc = CPLSPrintf("Band%d", i);
731 : }
732 : }
733 : else
734 : {
735 20 : bGotSignificantBandDesc = true;
736 : }
737 110 : if (i > 1)
738 49 : osLayerName += ":";
739 110 : osLayerName += CPLString(pszDesc).replaceAll(':', '.');
740 : }
741 61 : if (bGotSignificantBandDesc)
742 : {
743 9 : VSIFPrintfL(fp, "[description]\n");
744 9 : VSIFPrintfL(fp, "layername=%s\n", osLayerName.c_str());
745 : }
746 :
747 : // Put the [georeference] section at end. Otherwise the wkt= entry
748 : // could cause GDAL < 3.5 to be unable to open the file due to the
749 : // Identify() function not using enough bytes
750 61 : VSIFPrintfL(fp, "[georeference]\n");
751 61 : VSIFPrintfL(fp, "nrows=%d\n", nRasterYSize);
752 61 : VSIFPrintfL(fp, "ncols=%d\n", nRasterXSize);
753 :
754 61 : VSIFPrintfL(fp, "xmin=%.17g\n", m_gt[0]);
755 61 : VSIFPrintfL(fp, "ymin=%.17g\n", m_gt[3] + nRasterYSize * m_gt[5]);
756 61 : VSIFPrintfL(fp, "xmax=%.17g\n", m_gt[0] + nRasterXSize * m_gt[1]);
757 61 : VSIFPrintfL(fp, "ymax=%.17g\n", m_gt[3]);
758 :
759 61 : if (!m_oSRS.IsEmpty())
760 : {
761 59 : char *pszProj4 = nullptr;
762 59 : m_oSRS.exportToProj4(&pszProj4);
763 59 : if (pszProj4)
764 : {
765 59 : VSIFPrintfL(fp, "projection=%s\n", pszProj4);
766 59 : VSIFree(pszProj4);
767 : }
768 59 : char *pszWKT = nullptr;
769 59 : const char *const apszOptions[] = {"FORMAT=WKT2_2019", nullptr};
770 59 : m_oSRS.exportToWkt(&pszWKT, apszOptions);
771 59 : if (pszWKT)
772 : {
773 59 : VSIFPrintfL(fp, "wkt=%s\n", pszWKT);
774 59 : VSIFree(pszWKT);
775 : }
776 : }
777 :
778 61 : VSIFCloseL(fp);
779 : }
780 :
781 : /************************************************************************/
782 : /* GetFileList() */
783 : /************************************************************************/
784 :
785 37 : char **RRASTERDataset::GetFileList()
786 :
787 : {
788 37 : char **papszFileList = RawDataset::GetFileList();
789 :
790 37 : papszFileList = CSLAddString(papszFileList, m_osGriFilename);
791 :
792 37 : return papszFileList;
793 : }
794 :
795 : /************************************************************************/
796 : /* GetGeoTransform() */
797 : /************************************************************************/
798 :
799 120 : CPLErr RRASTERDataset::GetGeoTransform(GDALGeoTransform >) const
800 : {
801 120 : if (m_bGeoTransformValid)
802 : {
803 120 : gt = m_gt;
804 120 : return CE_None;
805 : }
806 0 : return CE_Failure;
807 : }
808 :
809 : /************************************************************************/
810 : /* SetGeoTransform() */
811 : /************************************************************************/
812 :
813 59 : CPLErr RRASTERDataset::SetGeoTransform(const GDALGeoTransform >)
814 :
815 : {
816 59 : if (GetAccess() != GA_Update)
817 : {
818 0 : CPLError(CE_Failure, CPLE_NotSupported,
819 : "Cannot set geotransform on a read-only dataset");
820 0 : return CE_Failure;
821 : }
822 :
823 : // We only support non-rotated images with info in the .HDR file.
824 59 : if (gt[2] != 0.0 || gt[4] != 0.0)
825 : {
826 0 : CPLError(CE_Warning, CPLE_NotSupported,
827 : "Rotated / skewed images not supported");
828 0 : return GDALPamDataset::SetGeoTransform(gt);
829 : }
830 :
831 : // Record new geotransform.
832 59 : m_bGeoTransformValid = true;
833 59 : m_gt = gt;
834 59 : SetHeaderDirty();
835 :
836 59 : return CE_None;
837 : }
838 :
839 : /************************************************************************/
840 : /* GetSpatialRef() */
841 : /************************************************************************/
842 :
843 35 : const OGRSpatialReference *RRASTERDataset::GetSpatialRef() const
844 : {
845 35 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
846 : }
847 :
848 : /************************************************************************/
849 : /* SetSpatialRef() */
850 : /************************************************************************/
851 :
852 59 : CPLErr RRASTERDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
853 : {
854 59 : if (GetAccess() != GA_Update)
855 : {
856 0 : CPLError(CE_Failure, CPLE_NotSupported,
857 : "Cannot set projection on a read-only dataset");
858 0 : return CE_Failure;
859 : }
860 :
861 59 : m_oSRS.Clear();
862 59 : if (poSRS)
863 59 : m_oSRS = *poSRS;
864 59 : SetHeaderDirty();
865 :
866 59 : return CE_None;
867 : }
868 :
869 : /************************************************************************/
870 : /* SetMetadata() */
871 : /************************************************************************/
872 :
873 31 : CPLErr RRASTERDataset::SetMetadata(char **papszMetadata, const char *pszDomain)
874 : {
875 31 : if (pszDomain == nullptr || EQUAL(pszDomain, ""))
876 : {
877 31 : m_osCreator = CSLFetchNameValueDef(papszMetadata, "CREATOR", "");
878 31 : m_osCreated = CSLFetchNameValueDef(papszMetadata, "CREATED", "");
879 31 : SetHeaderDirty();
880 : }
881 31 : return GDALDataset::SetMetadata(papszMetadata, pszDomain);
882 : }
883 :
884 : /************************************************************************/
885 : /* SetMetadataItem() */
886 : /************************************************************************/
887 :
888 2 : CPLErr RRASTERDataset::SetMetadataItem(const char *pszName,
889 : const char *pszValue,
890 : const char *pszDomain)
891 : {
892 2 : if (pszDomain == nullptr || EQUAL(pszDomain, ""))
893 : {
894 2 : if (EQUAL(pszName, "CREATOR"))
895 : {
896 1 : m_osCreator = pszValue ? pszValue : "";
897 1 : SetHeaderDirty();
898 : }
899 2 : if (EQUAL(pszName, "CREATED"))
900 : {
901 1 : m_osCreated = pszValue ? pszValue : "";
902 1 : SetHeaderDirty();
903 : }
904 : }
905 2 : return GDALDataset::SetMetadataItem(pszName, pszValue, pszDomain);
906 : }
907 :
908 : /************************************************************************/
909 : /* Identify() */
910 : /************************************************************************/
911 :
912 59539 : int RRASTERDataset::Identify(GDALOpenInfo *poOpenInfo)
913 :
914 : {
915 63231 : if (poOpenInfo->nHeaderBytes < 40 || poOpenInfo->fpL == nullptr ||
916 3692 : !poOpenInfo->IsExtensionEqualToCI("grd"))
917 : {
918 59343 : return FALSE;
919 : }
920 :
921 : // We need to ingest enough bytes as the wkt= entry can be quite large
922 196 : if (poOpenInfo->nHeaderBytes <= 1024)
923 149 : poOpenInfo->TryToIngest(4096);
924 :
925 192 : if (strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
926 173 : "ncols") == nullptr ||
927 173 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
928 173 : "nrows") == nullptr ||
929 173 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
930 173 : "xmin") == nullptr ||
931 173 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
932 173 : "ymin") == nullptr ||
933 173 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
934 173 : "xmax") == nullptr ||
935 173 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
936 173 : "ymax") == nullptr ||
937 173 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
938 : "datatype") == nullptr)
939 : {
940 19 : return FALSE;
941 : }
942 :
943 173 : return TRUE;
944 : }
945 :
946 : /************************************************************************/
947 : /* ComputeSpacing() */
948 : /************************************************************************/
949 :
950 147 : bool RRASTERDataset::ComputeSpacings(const CPLString &osBandOrder, int nCols,
951 : int nRows, int l_nBands, GDALDataType eDT,
952 : int &nPixelOffset, int &nLineOffset,
953 : vsi_l_offset &nBandOffset)
954 : {
955 147 : nPixelOffset = 0;
956 147 : nLineOffset = 0;
957 147 : nBandOffset = 0;
958 147 : const int nPixelSize = GDALGetDataTypeSizeBytes(eDT);
959 147 : if (l_nBands == 1 || EQUAL(osBandOrder, "BIL"))
960 : {
961 137 : nPixelOffset = nPixelSize;
962 137 : if (l_nBands != 0 && nPixelSize != 0 &&
963 137 : nCols > INT_MAX / (l_nBands * nPixelSize))
964 : {
965 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many columns");
966 0 : return false;
967 : }
968 137 : nLineOffset = nPixelSize * nCols * l_nBands;
969 137 : nBandOffset = static_cast<vsi_l_offset>(nPixelSize) * nCols;
970 : }
971 10 : else if (EQUAL(osBandOrder, "BIP"))
972 : {
973 3 : if (l_nBands != 0 && nPixelSize != 0 &&
974 3 : nCols > INT_MAX / (l_nBands * nPixelSize))
975 : {
976 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many columns");
977 0 : return false;
978 : }
979 3 : nPixelOffset = nPixelSize * l_nBands;
980 3 : nLineOffset = nPixelSize * nCols * l_nBands;
981 3 : nBandOffset = nPixelSize;
982 : }
983 7 : else if (EQUAL(osBandOrder, "BSQ"))
984 : {
985 7 : if (nPixelSize != 0 && nCols > INT_MAX / nPixelSize)
986 : {
987 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many columns");
988 0 : return false;
989 : }
990 7 : nPixelOffset = nPixelSize;
991 7 : nLineOffset = nPixelSize * nCols;
992 7 : nBandOffset = static_cast<vsi_l_offset>(nLineOffset) * nRows;
993 : }
994 0 : else if (l_nBands > 1)
995 : {
996 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unknown bandorder");
997 0 : return false;
998 : }
999 147 : return true;
1000 : }
1001 :
1002 : /************************************************************************/
1003 : /* CastToFloat() */
1004 : /************************************************************************/
1005 :
1006 0 : static float CastToFloat(double dfVal)
1007 : {
1008 0 : if (std::isinf(dfVal) || std::isnan(dfVal) ||
1009 0 : (dfVal >= -std::numeric_limits<float>::max() &&
1010 0 : dfVal <= std::numeric_limits<float>::max()))
1011 : {
1012 0 : return static_cast<float>(dfVal);
1013 : }
1014 0 : return std::numeric_limits<float>::quiet_NaN();
1015 : }
1016 :
1017 : /************************************************************************/
1018 : /* Open() */
1019 : /************************************************************************/
1020 :
1021 84 : GDALDataset *RRASTERDataset::Open(GDALOpenInfo *poOpenInfo)
1022 : {
1023 84 : if (!Identify(poOpenInfo))
1024 0 : return nullptr;
1025 :
1026 168 : auto poDS = std::make_unique<RRASTERDataset>();
1027 :
1028 84 : const char *pszLine = nullptr;
1029 84 : int nRows = 0;
1030 84 : int nCols = 0;
1031 84 : double dfXMin = 0.0;
1032 84 : double dfYMin = 0.0;
1033 84 : double dfXMax = 0.0;
1034 84 : double dfYMax = 0.0;
1035 84 : int l_nBands = 1;
1036 168 : CPLString osDataType;
1037 168 : CPLString osProjection;
1038 168 : CPLString osWkt;
1039 168 : CPLString osByteOrder;
1040 168 : CPLString osNoDataValue("NA");
1041 168 : CPLString osMinValue;
1042 168 : CPLString osMaxValue;
1043 168 : CPLString osLayerName;
1044 168 : CPLString osRatNames;
1045 168 : CPLString osRatTypes;
1046 168 : CPLString osRatValues;
1047 84 : bool bInLegend = false;
1048 84 : VSIRewindL(poOpenInfo->fpL);
1049 1592 : while ((pszLine = CPLReadLine2L(poOpenInfo->fpL, 1024 * 1024, nullptr)) !=
1050 : nullptr)
1051 : {
1052 1508 : if (pszLine[0] == '[')
1053 : {
1054 276 : bInLegend = EQUAL(pszLine, "[legend]");
1055 276 : continue;
1056 : }
1057 1232 : if (bInLegend)
1058 : {
1059 6 : poDS->m_osLegend += pszLine;
1060 6 : poDS->m_osLegend += "\n";
1061 : }
1062 1232 : char *pszKey = nullptr;
1063 1232 : const char *pszValue = CPLParseNameValue(pszLine, &pszKey);
1064 1232 : if (pszKey && pszValue)
1065 : {
1066 1232 : if (EQUAL(pszKey, "creator"))
1067 12 : poDS->m_osCreator = pszValue;
1068 1232 : if (EQUAL(pszKey, "created"))
1069 12 : poDS->m_osCreated = pszValue;
1070 1220 : else if (EQUAL(pszKey, "ncols"))
1071 84 : nCols = atoi(pszValue);
1072 1136 : else if (EQUAL(pszKey, "nrows"))
1073 84 : nRows = atoi(pszValue);
1074 1052 : else if (EQUAL(pszKey, "xmin"))
1075 84 : dfXMin = CPLAtof(pszValue);
1076 968 : else if (EQUAL(pszKey, "ymin"))
1077 84 : dfYMin = CPLAtof(pszValue);
1078 884 : else if (EQUAL(pszKey, "xmax"))
1079 84 : dfXMax = CPLAtof(pszValue);
1080 800 : else if (EQUAL(pszKey, "ymax"))
1081 84 : dfYMax = CPLAtof(pszValue);
1082 716 : else if (EQUAL(pszKey, "projection"))
1083 80 : osProjection = pszValue;
1084 636 : else if (EQUAL(pszKey, "wkt"))
1085 67 : osWkt = pszValue;
1086 569 : else if (EQUAL(pszKey, "nbands"))
1087 84 : l_nBands = atoi(pszValue);
1088 485 : else if (EQUAL(pszKey, "bandorder"))
1089 34 : poDS->m_osBandOrder = pszValue;
1090 451 : else if (EQUAL(pszKey, "datatype"))
1091 84 : osDataType = pszValue;
1092 367 : else if (EQUAL(pszKey, "byteorder"))
1093 84 : osByteOrder = pszValue;
1094 283 : else if (EQUAL(pszKey, "nodatavalue"))
1095 12 : osNoDataValue = pszValue;
1096 271 : else if (EQUAL(pszKey, "minvalue"))
1097 50 : osMinValue = pszValue;
1098 221 : else if (EQUAL(pszKey, "maxvalue"))
1099 50 : osMaxValue = pszValue;
1100 171 : else if (EQUAL(pszKey, "ratnames"))
1101 12 : osRatNames = pszValue;
1102 159 : else if (EQUAL(pszKey, "rattypes"))
1103 12 : osRatTypes = pszValue;
1104 147 : else if (EQUAL(pszKey, "ratvalues"))
1105 12 : osRatValues = pszValue;
1106 135 : else if (EQUAL(pszKey, "layername"))
1107 33 : osLayerName = pszValue;
1108 : }
1109 1232 : CPLFree(pszKey);
1110 : }
1111 84 : if (!GDALCheckDatasetDimensions(nCols, nRows))
1112 0 : return nullptr;
1113 84 : if (!GDALCheckBandCount(l_nBands, FALSE))
1114 0 : return nullptr;
1115 :
1116 84 : GDALDataType eDT = GDT_Unknown;
1117 84 : if (EQUAL(osDataType, "LOG1S"))
1118 0 : eDT = GDT_UInt8; // mapping TBC
1119 84 : else if (EQUAL(osDataType, "INT1S"))
1120 6 : eDT = GDT_Int8;
1121 78 : else if (EQUAL(osDataType, "INT2S"))
1122 5 : eDT = GDT_Int16;
1123 73 : else if (EQUAL(osDataType, "INT4S"))
1124 5 : eDT = GDT_Int32;
1125 : // else if( EQUAL(osDataType, "INT8S") )
1126 : // eDT = GDT_UInt64; // unhandled
1127 68 : else if (EQUAL(osDataType, "INT1U"))
1128 49 : eDT = GDT_UInt8;
1129 19 : else if (EQUAL(osDataType, "INT2U"))
1130 5 : eDT = GDT_UInt16;
1131 14 : else if (EQUAL(osDataType, "INT4U")) // Not documented
1132 5 : eDT = GDT_UInt32;
1133 9 : else if (EQUAL(osDataType, "FLT4S"))
1134 5 : eDT = GDT_Float32;
1135 4 : else if (EQUAL(osDataType, "FLT8S"))
1136 4 : eDT = GDT_Float64;
1137 : else
1138 : {
1139 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unhandled datatype=%s",
1140 : osDataType.c_str());
1141 0 : return nullptr;
1142 : }
1143 84 : if (l_nBands > 1 && poDS->m_osBandOrder.empty())
1144 : {
1145 0 : CPLError(CE_Failure, CPLE_AppDefined, "Missing 'bandorder'");
1146 0 : return nullptr;
1147 : }
1148 :
1149 84 : bool bNativeOrder = true;
1150 84 : if (EQUAL(osByteOrder, "little"))
1151 : {
1152 84 : bNativeOrder = CPL_TO_BOOL(CPL_IS_LSB);
1153 : }
1154 0 : else if (EQUAL(osByteOrder, "big"))
1155 : {
1156 0 : bNativeOrder = CPL_TO_BOOL(!CPL_IS_LSB);
1157 : }
1158 0 : else if (!EQUAL(osByteOrder, ""))
1159 : {
1160 0 : CPLError(CE_Warning, CPLE_AppDefined,
1161 : "Unhandled byteorder=%s. Assuming native order",
1162 : osByteOrder.c_str());
1163 : }
1164 :
1165 84 : int nPixelOffset = 0;
1166 84 : int nLineOffset = 0;
1167 84 : vsi_l_offset nBandOffset = 0;
1168 84 : if (!ComputeSpacings(poDS->m_osBandOrder, nCols, nRows, l_nBands, eDT,
1169 : nPixelOffset, nLineOffset, nBandOffset))
1170 : {
1171 0 : return nullptr;
1172 : }
1173 :
1174 168 : CPLString osDirname(CPLGetDirnameSafe(poOpenInfo->pszFilename));
1175 168 : CPLString osBasename(CPLGetBasenameSafe(poOpenInfo->pszFilename));
1176 168 : CPLString osGRDExtension(poOpenInfo->osExtension);
1177 168 : CPLString osGRIExtension((osGRDExtension[0] == 'g') ? "gri" : "GRI");
1178 84 : char **papszSiblings = poOpenInfo->GetSiblingFiles();
1179 84 : if (papszSiblings)
1180 : {
1181 84 : int iFile = CSLFindString(
1182 : papszSiblings,
1183 168 : CPLFormFilenameSafe(nullptr, osBasename, osGRIExtension).c_str());
1184 84 : if (iFile < 0)
1185 0 : return nullptr;
1186 84 : poDS->m_osGriFilename =
1187 168 : CPLFormFilenameSafe(osDirname, papszSiblings[iFile], nullptr);
1188 : }
1189 : else
1190 : {
1191 0 : poDS->m_osGriFilename =
1192 0 : CPLFormFilenameSafe(osDirname, osBasename, osGRIExtension);
1193 : }
1194 :
1195 : VSILFILE *fpImage =
1196 84 : VSIFOpenL(poDS->m_osGriFilename,
1197 84 : (poOpenInfo->eAccess == GA_Update) ? "rb+" : "rb");
1198 84 : if (fpImage == nullptr)
1199 : {
1200 0 : CPLError(CE_Failure, CPLE_FileIO, "Cannot open %s",
1201 0 : poDS->m_osGriFilename.c_str());
1202 0 : return nullptr;
1203 : }
1204 :
1205 84 : if (!RAWDatasetCheckMemoryUsage(nCols, nRows, l_nBands,
1206 : GDALGetDataTypeSizeBytes(eDT), nPixelOffset,
1207 : nLineOffset, 0, nBandOffset, fpImage))
1208 : {
1209 0 : VSIFCloseL(fpImage);
1210 0 : return nullptr;
1211 : }
1212 :
1213 84 : poDS->eAccess = poOpenInfo->eAccess;
1214 84 : poDS->nRasterXSize = nCols;
1215 84 : poDS->nRasterYSize = nRows;
1216 84 : poDS->m_bGeoTransformValid = true;
1217 84 : poDS->m_gt[0] = dfXMin;
1218 84 : poDS->m_gt[1] = (dfXMax - dfXMin) / nCols;
1219 84 : poDS->m_gt[2] = 0.0;
1220 84 : poDS->m_gt[3] = dfYMax;
1221 84 : poDS->m_gt[4] = 0.0;
1222 84 : poDS->m_gt[5] = -(dfYMax - dfYMin) / nRows;
1223 84 : poDS->m_fpImage = fpImage;
1224 84 : poDS->m_bNativeOrder = bNativeOrder;
1225 :
1226 84 : if (osWkt.empty())
1227 : {
1228 17 : if (!osProjection.empty())
1229 : {
1230 13 : poDS->m_oSRS.importFromProj4(osProjection.c_str());
1231 : }
1232 : }
1233 : else
1234 : {
1235 67 : poDS->m_oSRS.importFromWkt(osWkt.c_str());
1236 : }
1237 :
1238 84 : if (!poDS->m_osCreator.empty())
1239 12 : poDS->GDALDataset::SetMetadataItem("CREATOR", poDS->m_osCreator);
1240 :
1241 84 : if (!poDS->m_osCreated.empty())
1242 12 : poDS->GDALDataset::SetMetadataItem("CREATED", poDS->m_osCreated);
1243 :
1244 : // Instantiate RAT
1245 84 : if (!osRatNames.empty() && !osRatTypes.empty() && !osRatValues.empty())
1246 : {
1247 24 : CPLStringList aosRatNames(CSLTokenizeString2(osRatNames, ":", 0));
1248 24 : CPLStringList aosRatTypes(CSLTokenizeString2(osRatTypes, ":", 0));
1249 24 : CPLStringList aosRatValues(CSLTokenizeString2(osRatValues, ":", 0));
1250 12 : if (aosRatNames.size() >= 1 &&
1251 24 : aosRatNames.size() == aosRatTypes.size() &&
1252 12 : (aosRatValues.size() % aosRatNames.size()) == 0)
1253 : {
1254 12 : bool bIsCompatibleOfCT = false;
1255 12 : const int nValues = aosRatValues.size() / aosRatNames.size();
1256 20 : if ((aosRatNames.size() == 4 || aosRatNames.size() == 5) &&
1257 8 : EQUAL(aosRatNames[1], "red") &&
1258 8 : EQUAL(aosRatNames[2], "green") &&
1259 8 : EQUAL(aosRatNames[3], "blue") &&
1260 8 : (aosRatNames.size() == 4 || EQUAL(aosRatNames[4], "alpha")) &&
1261 8 : EQUAL(aosRatTypes[0], "integer") &&
1262 8 : EQUAL(aosRatTypes[1], "integer") &&
1263 8 : EQUAL(aosRatTypes[2], "integer") &&
1264 32 : EQUAL(aosRatTypes[3], "integer") &&
1265 8 : (aosRatTypes.size() == 4 || EQUAL(aosRatTypes[4], "integer")))
1266 : {
1267 8 : bIsCompatibleOfCT = true;
1268 8 : poDS->m_poCT.reset(new GDALColorTable());
1269 24 : for (int i = 0; i < nValues; i++)
1270 : {
1271 16 : const int nIndex = atoi(aosRatValues[i]);
1272 16 : if (nIndex >= 0 && nIndex < 65536)
1273 : {
1274 16 : const int nRed = atoi(aosRatValues[nValues + i]);
1275 16 : const int nGreen = atoi(aosRatValues[2 * nValues + i]);
1276 16 : const int nBlue = atoi(aosRatValues[3 * nValues + i]);
1277 : const int nAlpha =
1278 16 : aosRatTypes.size() == 4
1279 16 : ? 255
1280 8 : : atoi(aosRatValues[4 * nValues + i]);
1281 : const GDALColorEntry oEntry = {
1282 : static_cast<short>(nRed),
1283 : static_cast<short>(nGreen),
1284 : static_cast<short>(nBlue),
1285 16 : static_cast<short>(nAlpha)};
1286 :
1287 16 : poDS->m_poCT->SetColorEntry(nIndex, &oEntry);
1288 : }
1289 : else
1290 : {
1291 0 : bIsCompatibleOfCT = false;
1292 0 : poDS->m_poCT.reset();
1293 0 : break;
1294 : }
1295 : }
1296 : }
1297 :
1298 : // cppcheck-suppress knownConditionTrueFalse
1299 12 : if (!bIsCompatibleOfCT)
1300 : {
1301 4 : poDS->m_poRAT.reset(new GDALDefaultRasterAttributeTable());
1302 44 : for (int i = 0; i < aosRatNames.size(); i++)
1303 : {
1304 120 : poDS->m_poRAT->CreateColumn(
1305 40 : aosRatNames[i],
1306 40 : EQUAL(aosRatTypes[i], "integer") ? GFT_Integer
1307 12 : : EQUAL(aosRatTypes[i], "numeric") ? GFT_Real
1308 : : GFT_String,
1309 40 : EQUAL(aosRatNames[i], "red") ? GFU_Red
1310 68 : : EQUAL(aosRatNames[i], "green") ? GFU_Green
1311 60 : : EQUAL(aosRatNames[i], "blue") ? GFU_Blue
1312 52 : : EQUAL(aosRatNames[i], "alpha") ? GFU_Alpha
1313 44 : : EQUAL(aosRatNames[i], "name") ? GFU_Name
1314 20 : : EQUAL(aosRatNames[i], "pixelcount") ? GFU_PixelCount
1315 40 : : GFU_Generic);
1316 : }
1317 12 : for (int i = 0; i < nValues; i++)
1318 : {
1319 88 : for (int j = 0; j < aosRatTypes.size(); j++)
1320 : {
1321 80 : if (poDS->m_poRAT->GetTypeOfCol(j) == GFT_Integer)
1322 : {
1323 56 : poDS->m_poRAT->SetValue(
1324 56 : i, j, atoi(aosRatValues[j * nValues + i]));
1325 : }
1326 24 : else if (poDS->m_poRAT->GetTypeOfCol(j) == GFT_Real)
1327 : {
1328 16 : poDS->m_poRAT->SetValue(
1329 8 : i, j, CPLAtof(aosRatValues[j * nValues + i]));
1330 : }
1331 : else
1332 : {
1333 16 : poDS->m_poRAT->SetValue(
1334 16 : i, j, aosRatValues[j * nValues + i]);
1335 : }
1336 : }
1337 : }
1338 : }
1339 : }
1340 : }
1341 :
1342 168 : CPLStringList aosMinValues(CSLTokenizeString2(osMinValue, ":", 0));
1343 168 : CPLStringList aosMaxValues(CSLTokenizeString2(osMaxValue, ":", 0));
1344 :
1345 168 : CPLStringList aosLayerNames(CSLTokenizeString2(osLayerName, ":", 0));
1346 237 : for (int i = 1; i <= l_nBands; i++)
1347 : {
1348 : auto poBand = std::make_unique<RRASTERRasterBand>(
1349 153 : poDS.get(), i, fpImage, nBandOffset * (i - 1), nPixelOffset,
1350 153 : nLineOffset, eDT, bNativeOrder);
1351 153 : if (!poBand->IsValid())
1352 : {
1353 0 : return nullptr;
1354 : }
1355 153 : if (!osNoDataValue.empty() && !EQUAL(osNoDataValue, "NA"))
1356 : {
1357 10 : double dfNoDataValue = CPLAtof(osNoDataValue);
1358 10 : if (eDT == GDT_Float32)
1359 0 : dfNoDataValue = CastToFloat(dfNoDataValue);
1360 10 : poBand->m_bHasNoDataValue = true;
1361 10 : poBand->m_dfNoDataValue = dfNoDataValue;
1362 : }
1363 227 : if (i - 1 < static_cast<int>(aosMinValues.size()) &&
1364 74 : i - 1 < static_cast<int>(aosMaxValues.size()))
1365 : {
1366 148 : poBand->SetMinMax(CPLAtof(aosMinValues[i - 1]),
1367 74 : CPLAtof(aosMaxValues[i - 1]));
1368 : }
1369 153 : if (i - 1 < static_cast<int>(aosLayerNames.size()))
1370 : {
1371 156 : const CPLString osName(aosLayerNames[i - 1]);
1372 78 : poBand->GDALRasterBand::SetDescription(osName);
1373 78 : if (EQUAL(osName, "red"))
1374 15 : poBand->SetColorInterpretation(GCI_RedBand);
1375 63 : else if (EQUAL(osName, "green"))
1376 15 : poBand->SetColorInterpretation(GCI_GreenBand);
1377 48 : else if (EQUAL(osName, "blue"))
1378 15 : poBand->SetColorInterpretation(GCI_BlueBand);
1379 33 : else if (EQUAL(osName, "alpha"))
1380 15 : poBand->SetColorInterpretation(GCI_AlphaBand);
1381 : }
1382 153 : poBand->m_poRAT = poDS->m_poRAT;
1383 153 : poBand->m_poCT = poDS->m_poCT;
1384 153 : if (poBand->m_poCT)
1385 8 : poBand->SetColorInterpretation(GCI_PaletteIndex);
1386 153 : poDS->SetBand(i, std::move(poBand));
1387 : }
1388 :
1389 84 : return poDS.release();
1390 : }
1391 :
1392 : /************************************************************************/
1393 : /* Create() */
1394 : /************************************************************************/
1395 :
1396 80 : GDALDataset *RRASTERDataset::Create(const char *pszFilename, int nXSize,
1397 : int nYSize, int nBandsIn,
1398 : GDALDataType eType, char **papszOptions)
1399 :
1400 : {
1401 : // Verify input options.
1402 80 : if (nBandsIn <= 0)
1403 : {
1404 1 : CPLError(CE_Failure, CPLE_NotSupported,
1405 : "RRASTER driver does not support %d bands.", nBandsIn);
1406 1 : return nullptr;
1407 : }
1408 :
1409 79 : if (eType != GDT_UInt8 && eType != GDT_Int8 && eType != GDT_UInt16 &&
1410 32 : eType != GDT_Int16 && eType != GDT_Int32 && eType != GDT_UInt32 &&
1411 20 : eType != GDT_Float32 && eType != GDT_Float64)
1412 : {
1413 16 : CPLError(CE_Failure, CPLE_NotSupported, "Unsupported data type (%s).",
1414 : GDALGetDataTypeName(eType));
1415 16 : return nullptr;
1416 : }
1417 :
1418 126 : CPLString osGRDExtension(CPLGetExtensionSafe(pszFilename));
1419 63 : if (!EQUAL(osGRDExtension, "grd"))
1420 : {
1421 0 : CPLError(CE_Failure, CPLE_NotSupported,
1422 : "RRASTER driver only supports grd extension");
1423 0 : return nullptr;
1424 : }
1425 :
1426 63 : int nPixelOffset = 0;
1427 63 : int nLineOffset = 0;
1428 63 : vsi_l_offset nBandOffset = 0;
1429 : CPLString osBandOrder(
1430 126 : CSLFetchNameValueDef(papszOptions, "INTERLEAVE", "BIL"));
1431 63 : if (!ComputeSpacings(osBandOrder, nXSize, nYSize, nBandsIn, eType,
1432 : nPixelOffset, nLineOffset, nBandOffset))
1433 : {
1434 0 : return nullptr;
1435 : }
1436 :
1437 63 : const std::string osGRIExtension((osGRDExtension[0] == 'g') ? "gri"
1438 189 : : "GRI");
1439 : const std::string osGriFilename(
1440 126 : CPLResetExtensionSafe(pszFilename, osGRIExtension.c_str()));
1441 :
1442 : // Try to create the file.
1443 63 : VSILFILE *fpImage = VSIFOpenL(osGriFilename.c_str(), "wb+");
1444 :
1445 63 : if (fpImage == nullptr)
1446 : {
1447 3 : CPLError(CE_Failure, CPLE_OpenFailed,
1448 : "Attempt to create file `%s' failed.", osGriFilename.c_str());
1449 3 : return nullptr;
1450 : }
1451 :
1452 60 : RRASTERDataset *poDS = new RRASTERDataset;
1453 60 : poDS->eAccess = GA_Update;
1454 60 : poDS->m_bHeaderDirty = true;
1455 60 : poDS->m_osGriFilename = osGriFilename;
1456 60 : poDS->nRasterXSize = nXSize;
1457 60 : poDS->nRasterYSize = nYSize;
1458 60 : poDS->m_fpImage = fpImage;
1459 60 : poDS->m_bNativeOrder = true;
1460 60 : poDS->m_osBandOrder = osBandOrder.toupper();
1461 60 : poDS->m_bInitRaster = CPLFetchBool(papszOptions, "@INIT_RASTER", true);
1462 :
1463 60 : const char *pszPixelType = CSLFetchNameValue(papszOptions, "PIXELTYPE");
1464 61 : const bool bByteSigned = (eType == GDT_UInt8 && pszPixelType &&
1465 1 : EQUAL(pszPixelType, "SIGNEDBYTE"));
1466 60 : if (bByteSigned)
1467 1 : poDS->m_bSignedByte = true;
1468 :
1469 169 : for (int i = 1; i <= nBandsIn; i++)
1470 : {
1471 : RRASTERRasterBand *poBand =
1472 109 : new RRASTERRasterBand(poDS, i, fpImage, nBandOffset * (i - 1),
1473 109 : nPixelOffset, nLineOffset, eType, true);
1474 109 : poDS->SetBand(i, poBand);
1475 : }
1476 :
1477 60 : return poDS;
1478 : }
1479 :
1480 : /************************************************************************/
1481 : /* CreateCopy() */
1482 : /************************************************************************/
1483 :
1484 47 : GDALDataset *RRASTERDataset::CreateCopy(const char *pszFilename,
1485 : GDALDataset *poSrcDS, int bStrict,
1486 : char **papszOptions,
1487 : GDALProgressFunc pfnProgress,
1488 : void *pProgressData)
1489 :
1490 : {
1491 : // Proceed with normal copying using the default createcopy operators.
1492 : GDALDriver *poDriver =
1493 47 : reinterpret_cast<GDALDriver *>(GDALGetDriverByName("RRASTER"));
1494 :
1495 47 : char **papszAdjustedOptions = CSLDuplicate(papszOptions);
1496 : papszAdjustedOptions =
1497 47 : CSLSetNameValue(papszAdjustedOptions, "@INIT_RASTER", "NO");
1498 47 : GDALDataset *poOutDS = poDriver->DefaultCreateCopy(
1499 : pszFilename, poSrcDS, bStrict, papszAdjustedOptions, pfnProgress,
1500 : pProgressData);
1501 47 : CSLDestroy(papszAdjustedOptions);
1502 :
1503 47 : if (poOutDS != nullptr)
1504 37 : poOutDS->FlushCache(false);
1505 :
1506 47 : return poOutDS;
1507 : }
1508 :
1509 : /************************************************************************/
1510 : /* GDALRegister_RRASTER() */
1511 : /************************************************************************/
1512 :
1513 2054 : void GDALRegister_RRASTER()
1514 :
1515 : {
1516 2054 : if (GDALGetDriverByName("RRASTER") != nullptr)
1517 283 : return;
1518 :
1519 1771 : GDALDriver *poDriver = new GDALDriver();
1520 :
1521 1771 : poDriver->SetDescription("RRASTER");
1522 1771 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1523 1771 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "grd");
1524 1771 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "R Raster");
1525 1771 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
1526 1771 : "drivers/raster/rraster.html");
1527 1771 : poDriver->SetMetadataItem(
1528 : GDAL_DMD_CREATIONDATATYPES,
1529 1771 : "Byte Int8 Int16 UInt16 Int32 UInt32 Float32 Float64");
1530 :
1531 1771 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1532 :
1533 1771 : poDriver->SetMetadataItem(
1534 : GDAL_DMD_CREATIONOPTIONLIST,
1535 : "<CreationOptionList>"
1536 : " <Option name='PIXELTYPE' type='string' description='(deprecated, "
1537 : "use Int8) "
1538 : "By setting this to SIGNEDBYTE, a new Byte file can be forced to be "
1539 : "written "
1540 : "as signed byte'/>"
1541 : " <Option name='INTERLEAVE' type='string-select' default='BIL'>"
1542 : " <Value>BIP</Value>"
1543 : " <Value>BIL</Value>"
1544 : " <Value>BSQ</Value>"
1545 : " </Option>"
1546 1771 : "</CreationOptionList>");
1547 :
1548 1771 : poDriver->pfnOpen = RRASTERDataset::Open;
1549 1771 : poDriver->pfnIdentify = RRASTERDataset::Identify;
1550 1771 : poDriver->pfnCreate = RRASTERDataset::Create;
1551 1771 : poDriver->pfnCreateCopy = RRASTERDataset::CreateCopy;
1552 :
1553 1771 : poDriver->SetMetadataItem(GDAL_DCAP_UPDATE, "YES");
1554 1771 : poDriver->SetMetadataItem(GDAL_DMD_UPDATE_ITEMS, "GeoTransform SRS NoData "
1555 : "RasterValues "
1556 1771 : "DatasetMetadata");
1557 :
1558 1771 : GetGDALDriverManager()->RegisterDriver(poDriver);
1559 : }
|
