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_Byte:
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_Byte) ? "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 : #if CPL_IS_LSB
539 61 : VSIFPrintfL(fp, "byteorder=%s\n", m_bNativeOrder ? "little" : "big");
540 : #else
541 : VSIFPrintfL(fp, "byteorder=%s\n", !m_bNativeOrder ? "little" : "big");
542 : #endif
543 61 : VSIFPrintfL(fp, "nbands=%d\n", nBands);
544 61 : if (nBands > 1)
545 20 : VSIFPrintfL(fp, "bandorder=%s\n", m_osBandOrder.c_str());
546 122 : CPLString osMinValue, osMaxValue;
547 123 : for (int i = 1; i <= nBands; i++)
548 : {
549 : RRASTERRasterBand *poBand =
550 83 : static_cast<RRASTERRasterBand *>(GetRasterBand(i));
551 83 : if (i > 1)
552 : {
553 22 : osMinValue += ":";
554 22 : osMaxValue += ":";
555 : }
556 83 : if (poBand->m_dfMin > poBand->m_dfMax)
557 : {
558 21 : osMinValue.clear();
559 21 : break;
560 : }
561 62 : osMinValue += CPLSPrintf("%.17g", poBand->m_dfMin);
562 62 : osMaxValue += CPLSPrintf("%.17g", poBand->m_dfMax);
563 : }
564 61 : if (!osMinValue.empty())
565 : {
566 40 : VSIFPrintfL(fp, "minvalue=%s\n", osMinValue.c_str());
567 40 : VSIFPrintfL(fp, "maxvalue=%s\n", osMaxValue.c_str());
568 : }
569 :
570 61 : GDALColorTable *poCT = GetRasterBand(1)->GetColorTable();
571 61 : GDALRasterAttributeTable *poRAT = GetRasterBand(1)->GetDefaultRAT();
572 61 : if (poCT == nullptr && poRAT == nullptr)
573 : {
574 58 : VSIFPrintfL(fp, "categorical=FALSE\n");
575 : }
576 : else
577 : {
578 3 : VSIFPrintfL(fp, "categorical=TRUE\n");
579 3 : if (poCT && poRAT)
580 : {
581 0 : CPLError(CE_Warning, CPLE_AppDefined,
582 : "Both color table and raster attribute table defined. "
583 : "Writing only the later");
584 : }
585 3 : if (poRAT)
586 : {
587 2 : CPLString osRatNames;
588 2 : CPLString osRatTypes;
589 11 : for (int i = 0; i < poRAT->GetColumnCount(); i++)
590 : {
591 10 : if (!osRatNames.empty())
592 : {
593 9 : osRatNames += ":";
594 9 : osRatTypes += ":";
595 : }
596 : osRatNames +=
597 10 : CPLString(poRAT->GetNameOfCol(i)).replaceAll(':', '.');
598 10 : GDALRATFieldType eColType = poRAT->GetTypeOfCol(i);
599 10 : if (eColType == GFT_Integer)
600 7 : osRatTypes += "integer";
601 3 : else if (eColType == GFT_Real)
602 1 : osRatTypes += "numeric";
603 : else
604 2 : osRatTypes += "character";
605 : }
606 1 : VSIFPrintfL(fp, "ratnames=%s\n", osRatNames.c_str());
607 1 : VSIFPrintfL(fp, "rattypes=%s\n", osRatTypes.c_str());
608 2 : CPLString osRatValues;
609 11 : for (int i = 0; i < poRAT->GetColumnCount(); i++)
610 : {
611 10 : GDALRATFieldType eColType = poRAT->GetTypeOfCol(i);
612 30 : for (int j = 0; j < poRAT->GetRowCount(); j++)
613 : {
614 20 : if (i != 0 || j != 0)
615 19 : osRatValues += ":";
616 20 : if (eColType == GFT_Integer)
617 : {
618 : osRatValues +=
619 14 : CPLSPrintf("%d", poRAT->GetValueAsInt(j, i));
620 : }
621 6 : else if (eColType == GFT_Real)
622 : {
623 : osRatValues +=
624 2 : CPLSPrintf("%.17g", poRAT->GetValueAsDouble(j, i));
625 : }
626 : else
627 : {
628 4 : const char *pszVal = poRAT->GetValueAsString(j, i);
629 4 : if (pszVal)
630 : {
631 : osRatValues +=
632 4 : CPLString(pszVal).replaceAll(':', '.');
633 : }
634 : }
635 : }
636 : }
637 1 : VSIFPrintfL(fp, "ratvalues=%s\n", osRatValues.c_str());
638 : }
639 : else
640 : {
641 2 : bool bNeedsAlpha = false;
642 4 : for (int i = 0; i < poCT->GetColorEntryCount(); i++)
643 : {
644 3 : if (poCT->GetColorEntry(i)->c4 != 255)
645 : {
646 1 : bNeedsAlpha = true;
647 1 : break;
648 : }
649 : }
650 2 : if (!bNeedsAlpha)
651 : {
652 1 : VSIFPrintfL(fp, "ratnames=%s\n", "ID:red:green:blue");
653 1 : VSIFPrintfL(fp, "rattypes=%s\n",
654 : "integer:integer:integer:integer");
655 : }
656 : else
657 : {
658 1 : VSIFPrintfL(fp, "ratnames=%s\n", "ID:red:green:blue:alpha");
659 1 : VSIFPrintfL(fp, "rattypes=%s\n",
660 : "integer:integer:integer:integer:integer");
661 : }
662 :
663 4 : CPLString osRatID;
664 4 : CPLString osRatR;
665 4 : CPLString osRatG;
666 4 : CPLString osRatB;
667 4 : CPLString osRatA;
668 6 : for (int i = 0; i < poCT->GetColorEntryCount(); i++)
669 : {
670 4 : const GDALColorEntry *psEntry = poCT->GetColorEntry(i);
671 4 : if (i > 0)
672 : {
673 2 : osRatID += ":";
674 2 : osRatR += ":";
675 2 : osRatG += ":";
676 2 : osRatB += ":";
677 2 : osRatA += ":";
678 : }
679 4 : osRatID += CPLSPrintf("%d", i);
680 4 : osRatR += CPLSPrintf("%d", psEntry->c1);
681 4 : osRatG += CPLSPrintf("%d", psEntry->c2);
682 4 : osRatB += CPLSPrintf("%d", psEntry->c3);
683 4 : osRatA += CPLSPrintf("%d", psEntry->c4);
684 : }
685 2 : if (!bNeedsAlpha)
686 : {
687 1 : VSIFPrintfL(fp, "ratvalues=%s:%s:%s:%s\n", osRatID.c_str(),
688 : osRatR.c_str(), osRatG.c_str(), osRatB.c_str());
689 : }
690 : else
691 : {
692 1 : VSIFPrintfL(fp, "ratvalues=%s:%s:%s:%s:%s\n", osRatID.c_str(),
693 : osRatR.c_str(), osRatG.c_str(), osRatB.c_str(),
694 : osRatA.c_str());
695 : }
696 : }
697 : }
698 :
699 61 : if (!m_osLegend.empty())
700 0 : VSIFPrintfL(fp, "[legend]\n%s", m_osLegend.c_str());
701 :
702 122 : CPLString osLayerName;
703 61 : bool bGotSignificantBandDesc = false;
704 171 : for (int i = 1; i <= nBands; i++)
705 : {
706 110 : GDALRasterBand *poBand = GetRasterBand(i);
707 110 : const char *pszDesc = poBand->GetDescription();
708 110 : if (EQUAL(pszDesc, ""))
709 : {
710 90 : GDALColorInterp eInterp = poBand->GetColorInterpretation();
711 90 : if (eInterp == GCI_RedBand)
712 : {
713 1 : bGotSignificantBandDesc = true;
714 1 : pszDesc = "red";
715 : }
716 89 : else if (eInterp == GCI_GreenBand)
717 : {
718 1 : bGotSignificantBandDesc = true;
719 1 : pszDesc = "green";
720 : }
721 88 : else if (eInterp == GCI_BlueBand)
722 : {
723 1 : bGotSignificantBandDesc = true;
724 1 : pszDesc = "blue";
725 : }
726 87 : else if (eInterp == GCI_AlphaBand)
727 : {
728 1 : bGotSignificantBandDesc = true;
729 1 : pszDesc = "alpha";
730 : }
731 : else
732 : {
733 86 : pszDesc = CPLSPrintf("Band%d", i);
734 : }
735 : }
736 : else
737 : {
738 20 : bGotSignificantBandDesc = true;
739 : }
740 110 : if (i > 1)
741 49 : osLayerName += ":";
742 110 : osLayerName += CPLString(pszDesc).replaceAll(':', '.');
743 : }
744 61 : if (bGotSignificantBandDesc)
745 : {
746 9 : VSIFPrintfL(fp, "[description]\n");
747 9 : VSIFPrintfL(fp, "layername=%s\n", osLayerName.c_str());
748 : }
749 :
750 : // Put the [georeference] section at end. Otherwise the wkt= entry
751 : // could cause GDAL < 3.5 to be unable to open the file due to the
752 : // Identify() function not using enough bytes
753 61 : VSIFPrintfL(fp, "[georeference]\n");
754 61 : VSIFPrintfL(fp, "nrows=%d\n", nRasterYSize);
755 61 : VSIFPrintfL(fp, "ncols=%d\n", nRasterXSize);
756 :
757 61 : VSIFPrintfL(fp, "xmin=%.17g\n", m_gt[0]);
758 61 : VSIFPrintfL(fp, "ymin=%.17g\n", m_gt[3] + nRasterYSize * m_gt[5]);
759 61 : VSIFPrintfL(fp, "xmax=%.17g\n", m_gt[0] + nRasterXSize * m_gt[1]);
760 61 : VSIFPrintfL(fp, "ymax=%.17g\n", m_gt[3]);
761 :
762 61 : if (!m_oSRS.IsEmpty())
763 : {
764 59 : char *pszProj4 = nullptr;
765 59 : m_oSRS.exportToProj4(&pszProj4);
766 59 : if (pszProj4)
767 : {
768 59 : VSIFPrintfL(fp, "projection=%s\n", pszProj4);
769 59 : VSIFree(pszProj4);
770 : }
771 59 : char *pszWKT = nullptr;
772 59 : const char *const apszOptions[] = {"FORMAT=WKT2_2019", nullptr};
773 59 : m_oSRS.exportToWkt(&pszWKT, apszOptions);
774 59 : if (pszWKT)
775 : {
776 59 : VSIFPrintfL(fp, "wkt=%s\n", pszWKT);
777 59 : VSIFree(pszWKT);
778 : }
779 : }
780 :
781 61 : VSIFCloseL(fp);
782 : }
783 :
784 : /************************************************************************/
785 : /* GetFileList() */
786 : /************************************************************************/
787 :
788 37 : char **RRASTERDataset::GetFileList()
789 :
790 : {
791 37 : char **papszFileList = RawDataset::GetFileList();
792 :
793 37 : papszFileList = CSLAddString(papszFileList, m_osGriFilename);
794 :
795 37 : return papszFileList;
796 : }
797 :
798 : /************************************************************************/
799 : /* GetGeoTransform() */
800 : /************************************************************************/
801 :
802 120 : CPLErr RRASTERDataset::GetGeoTransform(GDALGeoTransform >) const
803 : {
804 120 : if (m_bGeoTransformValid)
805 : {
806 120 : gt = m_gt;
807 120 : return CE_None;
808 : }
809 0 : return CE_Failure;
810 : }
811 :
812 : /************************************************************************/
813 : /* SetGeoTransform() */
814 : /************************************************************************/
815 :
816 59 : CPLErr RRASTERDataset::SetGeoTransform(const GDALGeoTransform >)
817 :
818 : {
819 59 : if (GetAccess() != GA_Update)
820 : {
821 0 : CPLError(CE_Failure, CPLE_NotSupported,
822 : "Cannot set geotransform on a read-only dataset");
823 0 : return CE_Failure;
824 : }
825 :
826 : // We only support non-rotated images with info in the .HDR file.
827 59 : if (gt[2] != 0.0 || gt[4] != 0.0)
828 : {
829 0 : CPLError(CE_Warning, CPLE_NotSupported,
830 : "Rotated / skewed images not supported");
831 0 : return GDALPamDataset::SetGeoTransform(gt);
832 : }
833 :
834 : // Record new geotransform.
835 59 : m_bGeoTransformValid = true;
836 59 : m_gt = gt;
837 59 : SetHeaderDirty();
838 :
839 59 : return CE_None;
840 : }
841 :
842 : /************************************************************************/
843 : /* GetSpatialRef() */
844 : /************************************************************************/
845 :
846 35 : const OGRSpatialReference *RRASTERDataset::GetSpatialRef() const
847 : {
848 35 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
849 : }
850 :
851 : /************************************************************************/
852 : /* SetSpatialRef() */
853 : /************************************************************************/
854 :
855 59 : CPLErr RRASTERDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
856 : {
857 59 : if (GetAccess() != GA_Update)
858 : {
859 0 : CPLError(CE_Failure, CPLE_NotSupported,
860 : "Cannot set projection on a read-only dataset");
861 0 : return CE_Failure;
862 : }
863 :
864 59 : m_oSRS.Clear();
865 59 : if (poSRS)
866 59 : m_oSRS = *poSRS;
867 59 : SetHeaderDirty();
868 :
869 59 : return CE_None;
870 : }
871 :
872 : /************************************************************************/
873 : /* SetMetadata() */
874 : /************************************************************************/
875 :
876 31 : CPLErr RRASTERDataset::SetMetadata(char **papszMetadata, const char *pszDomain)
877 : {
878 31 : if (pszDomain == nullptr || EQUAL(pszDomain, ""))
879 : {
880 31 : m_osCreator = CSLFetchNameValueDef(papszMetadata, "CREATOR", "");
881 31 : m_osCreated = CSLFetchNameValueDef(papszMetadata, "CREATED", "");
882 31 : SetHeaderDirty();
883 : }
884 31 : return GDALDataset::SetMetadata(papszMetadata, pszDomain);
885 : }
886 :
887 : /************************************************************************/
888 : /* SetMetadataItem() */
889 : /************************************************************************/
890 :
891 2 : CPLErr RRASTERDataset::SetMetadataItem(const char *pszName,
892 : const char *pszValue,
893 : const char *pszDomain)
894 : {
895 2 : if (pszDomain == nullptr || EQUAL(pszDomain, ""))
896 : {
897 2 : if (EQUAL(pszName, "CREATOR"))
898 : {
899 1 : m_osCreator = pszValue ? pszValue : "";
900 1 : SetHeaderDirty();
901 : }
902 2 : if (EQUAL(pszName, "CREATED"))
903 : {
904 1 : m_osCreated = pszValue ? pszValue : "";
905 1 : SetHeaderDirty();
906 : }
907 : }
908 2 : return GDALDataset::SetMetadataItem(pszName, pszValue, pszDomain);
909 : }
910 :
911 : /************************************************************************/
912 : /* Identify() */
913 : /************************************************************************/
914 :
915 57191 : int RRASTERDataset::Identify(GDALOpenInfo *poOpenInfo)
916 :
917 : {
918 60747 : if (poOpenInfo->nHeaderBytes < 40 || poOpenInfo->fpL == nullptr ||
919 3556 : !poOpenInfo->IsExtensionEqualToCI("grd"))
920 : {
921 57000 : return FALSE;
922 : }
923 :
924 : // We need to ingest enough bytes as the wkt= entry can be quite large
925 191 : if (poOpenInfo->nHeaderBytes <= 1024)
926 148 : poOpenInfo->TryToIngest(4096);
927 :
928 191 : if (strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
929 173 : "ncols") == nullptr ||
930 173 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
931 173 : "nrows") == nullptr ||
932 173 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
933 173 : "xmin") == nullptr ||
934 173 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
935 173 : "ymin") == nullptr ||
936 173 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
937 173 : "xmax") == nullptr ||
938 173 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
939 173 : "ymax") == nullptr ||
940 173 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
941 : "datatype") == nullptr)
942 : {
943 18 : return FALSE;
944 : }
945 :
946 173 : return TRUE;
947 : }
948 :
949 : /************************************************************************/
950 : /* ComputeSpacing() */
951 : /************************************************************************/
952 :
953 147 : bool RRASTERDataset::ComputeSpacings(const CPLString &osBandOrder, int nCols,
954 : int nRows, int l_nBands, GDALDataType eDT,
955 : int &nPixelOffset, int &nLineOffset,
956 : vsi_l_offset &nBandOffset)
957 : {
958 147 : nPixelOffset = 0;
959 147 : nLineOffset = 0;
960 147 : nBandOffset = 0;
961 147 : const int nPixelSize = GDALGetDataTypeSizeBytes(eDT);
962 147 : if (l_nBands == 1 || EQUAL(osBandOrder, "BIL"))
963 : {
964 137 : nPixelOffset = nPixelSize;
965 137 : if (l_nBands != 0 && nPixelSize != 0 &&
966 137 : nCols > INT_MAX / (l_nBands * nPixelSize))
967 : {
968 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many columns");
969 0 : return false;
970 : }
971 137 : nLineOffset = nPixelSize * nCols * l_nBands;
972 137 : nBandOffset = static_cast<vsi_l_offset>(nPixelSize) * nCols;
973 : }
974 10 : else if (EQUAL(osBandOrder, "BIP"))
975 : {
976 3 : if (l_nBands != 0 && nPixelSize != 0 &&
977 3 : nCols > INT_MAX / (l_nBands * nPixelSize))
978 : {
979 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many columns");
980 0 : return false;
981 : }
982 3 : nPixelOffset = nPixelSize * l_nBands;
983 3 : nLineOffset = nPixelSize * nCols * l_nBands;
984 3 : nBandOffset = nPixelSize;
985 : }
986 7 : else if (EQUAL(osBandOrder, "BSQ"))
987 : {
988 7 : if (nPixelSize != 0 && nCols > INT_MAX / nPixelSize)
989 : {
990 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many columns");
991 0 : return false;
992 : }
993 7 : nPixelOffset = nPixelSize;
994 7 : nLineOffset = nPixelSize * nCols;
995 7 : nBandOffset = static_cast<vsi_l_offset>(nLineOffset) * nRows;
996 : }
997 0 : else if (l_nBands > 1)
998 : {
999 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unknown bandorder");
1000 0 : return false;
1001 : }
1002 147 : return true;
1003 : }
1004 :
1005 : /************************************************************************/
1006 : /* CastToFloat() */
1007 : /************************************************************************/
1008 :
1009 0 : static float CastToFloat(double dfVal)
1010 : {
1011 0 : if (std::isinf(dfVal) || std::isnan(dfVal) ||
1012 0 : (dfVal >= -std::numeric_limits<float>::max() &&
1013 0 : dfVal <= std::numeric_limits<float>::max()))
1014 : {
1015 0 : return static_cast<float>(dfVal);
1016 : }
1017 0 : return std::numeric_limits<float>::quiet_NaN();
1018 : }
1019 :
1020 : /************************************************************************/
1021 : /* Open() */
1022 : /************************************************************************/
1023 :
1024 84 : GDALDataset *RRASTERDataset::Open(GDALOpenInfo *poOpenInfo)
1025 : {
1026 84 : if (!Identify(poOpenInfo))
1027 0 : return nullptr;
1028 :
1029 168 : auto poDS = std::make_unique<RRASTERDataset>();
1030 :
1031 84 : const char *pszLine = nullptr;
1032 84 : int nRows = 0;
1033 84 : int nCols = 0;
1034 84 : double dfXMin = 0.0;
1035 84 : double dfYMin = 0.0;
1036 84 : double dfXMax = 0.0;
1037 84 : double dfYMax = 0.0;
1038 84 : int l_nBands = 1;
1039 168 : CPLString osDataType;
1040 168 : CPLString osProjection;
1041 168 : CPLString osWkt;
1042 168 : CPLString osByteOrder;
1043 168 : CPLString osNoDataValue("NA");
1044 168 : CPLString osMinValue;
1045 168 : CPLString osMaxValue;
1046 168 : CPLString osLayerName;
1047 168 : CPLString osRatNames;
1048 168 : CPLString osRatTypes;
1049 168 : CPLString osRatValues;
1050 84 : bool bInLegend = false;
1051 84 : VSIRewindL(poOpenInfo->fpL);
1052 1592 : while ((pszLine = CPLReadLine2L(poOpenInfo->fpL, 1024 * 1024, nullptr)) !=
1053 : nullptr)
1054 : {
1055 1508 : if (pszLine[0] == '[')
1056 : {
1057 276 : bInLegend = EQUAL(pszLine, "[legend]");
1058 276 : continue;
1059 : }
1060 1232 : if (bInLegend)
1061 : {
1062 6 : poDS->m_osLegend += pszLine;
1063 6 : poDS->m_osLegend += "\n";
1064 : }
1065 1232 : char *pszKey = nullptr;
1066 1232 : const char *pszValue = CPLParseNameValue(pszLine, &pszKey);
1067 1232 : if (pszKey && pszValue)
1068 : {
1069 1232 : if (EQUAL(pszKey, "creator"))
1070 12 : poDS->m_osCreator = pszValue;
1071 1232 : if (EQUAL(pszKey, "created"))
1072 12 : poDS->m_osCreated = pszValue;
1073 1220 : else if (EQUAL(pszKey, "ncols"))
1074 84 : nCols = atoi(pszValue);
1075 1136 : else if (EQUAL(pszKey, "nrows"))
1076 84 : nRows = atoi(pszValue);
1077 1052 : else if (EQUAL(pszKey, "xmin"))
1078 84 : dfXMin = CPLAtof(pszValue);
1079 968 : else if (EQUAL(pszKey, "ymin"))
1080 84 : dfYMin = CPLAtof(pszValue);
1081 884 : else if (EQUAL(pszKey, "xmax"))
1082 84 : dfXMax = CPLAtof(pszValue);
1083 800 : else if (EQUAL(pszKey, "ymax"))
1084 84 : dfYMax = CPLAtof(pszValue);
1085 716 : else if (EQUAL(pszKey, "projection"))
1086 80 : osProjection = pszValue;
1087 636 : else if (EQUAL(pszKey, "wkt"))
1088 67 : osWkt = pszValue;
1089 569 : else if (EQUAL(pszKey, "nbands"))
1090 84 : l_nBands = atoi(pszValue);
1091 485 : else if (EQUAL(pszKey, "bandorder"))
1092 34 : poDS->m_osBandOrder = pszValue;
1093 451 : else if (EQUAL(pszKey, "datatype"))
1094 84 : osDataType = pszValue;
1095 367 : else if (EQUAL(pszKey, "byteorder"))
1096 84 : osByteOrder = pszValue;
1097 283 : else if (EQUAL(pszKey, "nodatavalue"))
1098 12 : osNoDataValue = pszValue;
1099 271 : else if (EQUAL(pszKey, "minvalue"))
1100 50 : osMinValue = pszValue;
1101 221 : else if (EQUAL(pszKey, "maxvalue"))
1102 50 : osMaxValue = pszValue;
1103 171 : else if (EQUAL(pszKey, "ratnames"))
1104 12 : osRatNames = pszValue;
1105 159 : else if (EQUAL(pszKey, "rattypes"))
1106 12 : osRatTypes = pszValue;
1107 147 : else if (EQUAL(pszKey, "ratvalues"))
1108 12 : osRatValues = pszValue;
1109 135 : else if (EQUAL(pszKey, "layername"))
1110 33 : osLayerName = pszValue;
1111 : }
1112 1232 : CPLFree(pszKey);
1113 : }
1114 84 : if (!GDALCheckDatasetDimensions(nCols, nRows))
1115 0 : return nullptr;
1116 84 : if (!GDALCheckBandCount(l_nBands, FALSE))
1117 0 : return nullptr;
1118 :
1119 84 : GDALDataType eDT = GDT_Unknown;
1120 84 : if (EQUAL(osDataType, "LOG1S"))
1121 0 : eDT = GDT_Byte; // mapping TBC
1122 84 : else if (EQUAL(osDataType, "INT1S"))
1123 6 : eDT = GDT_Int8;
1124 78 : else if (EQUAL(osDataType, "INT2S"))
1125 5 : eDT = GDT_Int16;
1126 73 : else if (EQUAL(osDataType, "INT4S"))
1127 5 : eDT = GDT_Int32;
1128 : // else if( EQUAL(osDataType, "INT8S") )
1129 : // eDT = GDT_UInt64; // unhandled
1130 68 : else if (EQUAL(osDataType, "INT1U"))
1131 49 : eDT = GDT_Byte;
1132 19 : else if (EQUAL(osDataType, "INT2U"))
1133 5 : eDT = GDT_UInt16;
1134 14 : else if (EQUAL(osDataType, "INT4U")) // Not documented
1135 5 : eDT = GDT_UInt32;
1136 9 : else if (EQUAL(osDataType, "FLT4S"))
1137 5 : eDT = GDT_Float32;
1138 4 : else if (EQUAL(osDataType, "FLT8S"))
1139 4 : eDT = GDT_Float64;
1140 : else
1141 : {
1142 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unhandled datatype=%s",
1143 : osDataType.c_str());
1144 0 : return nullptr;
1145 : }
1146 84 : if (l_nBands > 1 && poDS->m_osBandOrder.empty())
1147 : {
1148 0 : CPLError(CE_Failure, CPLE_AppDefined, "Missing 'bandorder'");
1149 0 : return nullptr;
1150 : }
1151 :
1152 84 : bool bNativeOrder = true;
1153 84 : if (EQUAL(osByteOrder, "little"))
1154 : {
1155 84 : bNativeOrder = CPL_TO_BOOL(CPL_IS_LSB);
1156 : }
1157 0 : else if (EQUAL(osByteOrder, "big"))
1158 : {
1159 0 : bNativeOrder = CPL_TO_BOOL(!CPL_IS_LSB);
1160 : }
1161 0 : else if (!EQUAL(osByteOrder, ""))
1162 : {
1163 0 : CPLError(CE_Warning, CPLE_AppDefined,
1164 : "Unhandled byteorder=%s. Assuming native order",
1165 : osByteOrder.c_str());
1166 : }
1167 :
1168 84 : int nPixelOffset = 0;
1169 84 : int nLineOffset = 0;
1170 84 : vsi_l_offset nBandOffset = 0;
1171 84 : if (!ComputeSpacings(poDS->m_osBandOrder, nCols, nRows, l_nBands, eDT,
1172 : nPixelOffset, nLineOffset, nBandOffset))
1173 : {
1174 0 : return nullptr;
1175 : }
1176 :
1177 168 : CPLString osDirname(CPLGetDirnameSafe(poOpenInfo->pszFilename));
1178 168 : CPLString osBasename(CPLGetBasenameSafe(poOpenInfo->pszFilename));
1179 168 : CPLString osGRDExtension(poOpenInfo->osExtension);
1180 168 : CPLString osGRIExtension((osGRDExtension[0] == 'g') ? "gri" : "GRI");
1181 84 : char **papszSiblings = poOpenInfo->GetSiblingFiles();
1182 84 : if (papszSiblings)
1183 : {
1184 84 : int iFile = CSLFindString(
1185 : papszSiblings,
1186 168 : CPLFormFilenameSafe(nullptr, osBasename, osGRIExtension).c_str());
1187 84 : if (iFile < 0)
1188 0 : return nullptr;
1189 84 : poDS->m_osGriFilename =
1190 168 : CPLFormFilenameSafe(osDirname, papszSiblings[iFile], nullptr);
1191 : }
1192 : else
1193 : {
1194 0 : poDS->m_osGriFilename =
1195 0 : CPLFormFilenameSafe(osDirname, osBasename, osGRIExtension);
1196 : }
1197 :
1198 : VSILFILE *fpImage =
1199 84 : VSIFOpenL(poDS->m_osGriFilename,
1200 84 : (poOpenInfo->eAccess == GA_Update) ? "rb+" : "rb");
1201 84 : if (fpImage == nullptr)
1202 : {
1203 0 : CPLError(CE_Failure, CPLE_FileIO, "Cannot open %s",
1204 0 : poDS->m_osGriFilename.c_str());
1205 0 : return nullptr;
1206 : }
1207 :
1208 84 : if (!RAWDatasetCheckMemoryUsage(nCols, nRows, l_nBands,
1209 : GDALGetDataTypeSizeBytes(eDT), nPixelOffset,
1210 : nLineOffset, 0, nBandOffset, fpImage))
1211 : {
1212 0 : VSIFCloseL(fpImage);
1213 0 : return nullptr;
1214 : }
1215 :
1216 84 : poDS->eAccess = poOpenInfo->eAccess;
1217 84 : poDS->nRasterXSize = nCols;
1218 84 : poDS->nRasterYSize = nRows;
1219 84 : poDS->m_bGeoTransformValid = true;
1220 84 : poDS->m_gt[0] = dfXMin;
1221 84 : poDS->m_gt[1] = (dfXMax - dfXMin) / nCols;
1222 84 : poDS->m_gt[2] = 0.0;
1223 84 : poDS->m_gt[3] = dfYMax;
1224 84 : poDS->m_gt[4] = 0.0;
1225 84 : poDS->m_gt[5] = -(dfYMax - dfYMin) / nRows;
1226 84 : poDS->m_fpImage = fpImage;
1227 84 : poDS->m_bNativeOrder = bNativeOrder;
1228 :
1229 84 : if (osWkt.empty())
1230 : {
1231 17 : if (!osProjection.empty())
1232 : {
1233 13 : poDS->m_oSRS.importFromProj4(osProjection.c_str());
1234 : }
1235 : }
1236 : else
1237 : {
1238 67 : poDS->m_oSRS.importFromWkt(osWkt.c_str());
1239 : }
1240 :
1241 84 : if (!poDS->m_osCreator.empty())
1242 12 : poDS->GDALDataset::SetMetadataItem("CREATOR", poDS->m_osCreator);
1243 :
1244 84 : if (!poDS->m_osCreated.empty())
1245 12 : poDS->GDALDataset::SetMetadataItem("CREATED", poDS->m_osCreated);
1246 :
1247 : // Instantiate RAT
1248 84 : if (!osRatNames.empty() && !osRatTypes.empty() && !osRatValues.empty())
1249 : {
1250 24 : CPLStringList aosRatNames(CSLTokenizeString2(osRatNames, ":", 0));
1251 24 : CPLStringList aosRatTypes(CSLTokenizeString2(osRatTypes, ":", 0));
1252 24 : CPLStringList aosRatValues(CSLTokenizeString2(osRatValues, ":", 0));
1253 12 : if (aosRatNames.size() >= 1 &&
1254 24 : aosRatNames.size() == aosRatTypes.size() &&
1255 12 : (aosRatValues.size() % aosRatNames.size()) == 0)
1256 : {
1257 12 : bool bIsCompatibleOfCT = false;
1258 12 : const int nValues = aosRatValues.size() / aosRatNames.size();
1259 20 : if ((aosRatNames.size() == 4 || aosRatNames.size() == 5) &&
1260 8 : EQUAL(aosRatNames[1], "red") &&
1261 8 : EQUAL(aosRatNames[2], "green") &&
1262 8 : EQUAL(aosRatNames[3], "blue") &&
1263 8 : (aosRatNames.size() == 4 || EQUAL(aosRatNames[4], "alpha")) &&
1264 8 : EQUAL(aosRatTypes[0], "integer") &&
1265 8 : EQUAL(aosRatTypes[1], "integer") &&
1266 8 : EQUAL(aosRatTypes[2], "integer") &&
1267 32 : EQUAL(aosRatTypes[3], "integer") &&
1268 8 : (aosRatTypes.size() == 4 || EQUAL(aosRatTypes[4], "integer")))
1269 : {
1270 8 : bIsCompatibleOfCT = true;
1271 8 : poDS->m_poCT.reset(new GDALColorTable());
1272 24 : for (int i = 0; i < nValues; i++)
1273 : {
1274 16 : const int nIndex = atoi(aosRatValues[i]);
1275 16 : if (nIndex >= 0 && nIndex < 65536)
1276 : {
1277 16 : const int nRed = atoi(aosRatValues[nValues + i]);
1278 16 : const int nGreen = atoi(aosRatValues[2 * nValues + i]);
1279 16 : const int nBlue = atoi(aosRatValues[3 * nValues + i]);
1280 : const int nAlpha =
1281 16 : aosRatTypes.size() == 4
1282 16 : ? 255
1283 8 : : atoi(aosRatValues[4 * nValues + i]);
1284 : const GDALColorEntry oEntry = {
1285 : static_cast<short>(nRed),
1286 : static_cast<short>(nGreen),
1287 : static_cast<short>(nBlue),
1288 16 : static_cast<short>(nAlpha)};
1289 :
1290 16 : poDS->m_poCT->SetColorEntry(nIndex, &oEntry);
1291 : }
1292 : else
1293 : {
1294 0 : bIsCompatibleOfCT = false;
1295 0 : poDS->m_poCT.reset();
1296 0 : break;
1297 : }
1298 : }
1299 : }
1300 :
1301 : // cppcheck-suppress knownConditionTrueFalse
1302 12 : if (!bIsCompatibleOfCT)
1303 : {
1304 4 : poDS->m_poRAT.reset(new GDALDefaultRasterAttributeTable());
1305 44 : for (int i = 0; i < aosRatNames.size(); i++)
1306 : {
1307 120 : poDS->m_poRAT->CreateColumn(
1308 40 : aosRatNames[i],
1309 40 : EQUAL(aosRatTypes[i], "integer") ? GFT_Integer
1310 12 : : EQUAL(aosRatTypes[i], "numeric") ? GFT_Real
1311 : : GFT_String,
1312 40 : EQUAL(aosRatNames[i], "red") ? GFU_Red
1313 68 : : EQUAL(aosRatNames[i], "green") ? GFU_Green
1314 60 : : EQUAL(aosRatNames[i], "blue") ? GFU_Blue
1315 52 : : EQUAL(aosRatNames[i], "alpha") ? GFU_Alpha
1316 44 : : EQUAL(aosRatNames[i], "name") ? GFU_Name
1317 20 : : EQUAL(aosRatNames[i], "pixelcount") ? GFU_PixelCount
1318 40 : : GFU_Generic);
1319 : }
1320 12 : for (int i = 0; i < nValues; i++)
1321 : {
1322 88 : for (int j = 0; j < aosRatTypes.size(); j++)
1323 : {
1324 80 : if (poDS->m_poRAT->GetTypeOfCol(j) == GFT_Integer)
1325 : {
1326 56 : poDS->m_poRAT->SetValue(
1327 56 : i, j, atoi(aosRatValues[j * nValues + i]));
1328 : }
1329 24 : else if (poDS->m_poRAT->GetTypeOfCol(j) == GFT_Real)
1330 : {
1331 16 : poDS->m_poRAT->SetValue(
1332 8 : i, j, CPLAtof(aosRatValues[j * nValues + i]));
1333 : }
1334 : else
1335 : {
1336 16 : poDS->m_poRAT->SetValue(
1337 16 : i, j, aosRatValues[j * nValues + i]);
1338 : }
1339 : }
1340 : }
1341 : }
1342 : }
1343 : }
1344 :
1345 168 : CPLStringList aosMinValues(CSLTokenizeString2(osMinValue, ":", 0));
1346 168 : CPLStringList aosMaxValues(CSLTokenizeString2(osMaxValue, ":", 0));
1347 :
1348 168 : CPLStringList aosLayerNames(CSLTokenizeString2(osLayerName, ":", 0));
1349 237 : for (int i = 1; i <= l_nBands; i++)
1350 : {
1351 : auto poBand = std::make_unique<RRASTERRasterBand>(
1352 153 : poDS.get(), i, fpImage, nBandOffset * (i - 1), nPixelOffset,
1353 153 : nLineOffset, eDT, bNativeOrder);
1354 153 : if (!poBand->IsValid())
1355 : {
1356 0 : return nullptr;
1357 : }
1358 153 : if (!osNoDataValue.empty() && !EQUAL(osNoDataValue, "NA"))
1359 : {
1360 10 : double dfNoDataValue = CPLAtof(osNoDataValue);
1361 10 : if (eDT == GDT_Float32)
1362 0 : dfNoDataValue = CastToFloat(dfNoDataValue);
1363 10 : poBand->m_bHasNoDataValue = true;
1364 10 : poBand->m_dfNoDataValue = dfNoDataValue;
1365 : }
1366 227 : if (i - 1 < static_cast<int>(aosMinValues.size()) &&
1367 74 : i - 1 < static_cast<int>(aosMaxValues.size()))
1368 : {
1369 148 : poBand->SetMinMax(CPLAtof(aosMinValues[i - 1]),
1370 74 : CPLAtof(aosMaxValues[i - 1]));
1371 : }
1372 153 : if (i - 1 < static_cast<int>(aosLayerNames.size()))
1373 : {
1374 156 : const CPLString osName(aosLayerNames[i - 1]);
1375 78 : poBand->GDALRasterBand::SetDescription(osName);
1376 78 : if (EQUAL(osName, "red"))
1377 15 : poBand->SetColorInterpretation(GCI_RedBand);
1378 63 : else if (EQUAL(osName, "green"))
1379 15 : poBand->SetColorInterpretation(GCI_GreenBand);
1380 48 : else if (EQUAL(osName, "blue"))
1381 15 : poBand->SetColorInterpretation(GCI_BlueBand);
1382 33 : else if (EQUAL(osName, "alpha"))
1383 15 : poBand->SetColorInterpretation(GCI_AlphaBand);
1384 : }
1385 153 : poBand->m_poRAT = poDS->m_poRAT;
1386 153 : poBand->m_poCT = poDS->m_poCT;
1387 153 : if (poBand->m_poCT)
1388 8 : poBand->SetColorInterpretation(GCI_PaletteIndex);
1389 153 : poDS->SetBand(i, std::move(poBand));
1390 : }
1391 :
1392 84 : return poDS.release();
1393 : }
1394 :
1395 : /************************************************************************/
1396 : /* Create() */
1397 : /************************************************************************/
1398 :
1399 80 : GDALDataset *RRASTERDataset::Create(const char *pszFilename, int nXSize,
1400 : int nYSize, int nBandsIn,
1401 : GDALDataType eType, char **papszOptions)
1402 :
1403 : {
1404 : // Verify input options.
1405 80 : if (nBandsIn <= 0)
1406 : {
1407 1 : CPLError(CE_Failure, CPLE_NotSupported,
1408 : "RRASTER driver does not support %d bands.", nBandsIn);
1409 1 : return nullptr;
1410 : }
1411 :
1412 79 : if (eType != GDT_Byte && eType != GDT_Int8 && eType != GDT_UInt16 &&
1413 32 : eType != GDT_Int16 && eType != GDT_Int32 && eType != GDT_UInt32 &&
1414 20 : eType != GDT_Float32 && eType != GDT_Float64)
1415 : {
1416 16 : CPLError(CE_Failure, CPLE_NotSupported, "Unsupported data type (%s).",
1417 : GDALGetDataTypeName(eType));
1418 16 : return nullptr;
1419 : }
1420 :
1421 126 : CPLString osGRDExtension(CPLGetExtensionSafe(pszFilename));
1422 63 : if (!EQUAL(osGRDExtension, "grd"))
1423 : {
1424 0 : CPLError(CE_Failure, CPLE_NotSupported,
1425 : "RRASTER driver only supports grd extension");
1426 0 : return nullptr;
1427 : }
1428 :
1429 63 : int nPixelOffset = 0;
1430 63 : int nLineOffset = 0;
1431 63 : vsi_l_offset nBandOffset = 0;
1432 : CPLString osBandOrder(
1433 126 : CSLFetchNameValueDef(papszOptions, "INTERLEAVE", "BIL"));
1434 63 : if (!ComputeSpacings(osBandOrder, nXSize, nYSize, nBandsIn, eType,
1435 : nPixelOffset, nLineOffset, nBandOffset))
1436 : {
1437 0 : return nullptr;
1438 : }
1439 :
1440 63 : const std::string osGRIExtension((osGRDExtension[0] == 'g') ? "gri"
1441 189 : : "GRI");
1442 : const std::string osGriFilename(
1443 126 : CPLResetExtensionSafe(pszFilename, osGRIExtension.c_str()));
1444 :
1445 : // Try to create the file.
1446 63 : VSILFILE *fpImage = VSIFOpenL(osGriFilename.c_str(), "wb+");
1447 :
1448 63 : if (fpImage == nullptr)
1449 : {
1450 3 : CPLError(CE_Failure, CPLE_OpenFailed,
1451 : "Attempt to create file `%s' failed.", osGriFilename.c_str());
1452 3 : return nullptr;
1453 : }
1454 :
1455 60 : RRASTERDataset *poDS = new RRASTERDataset;
1456 60 : poDS->eAccess = GA_Update;
1457 60 : poDS->m_bHeaderDirty = true;
1458 60 : poDS->m_osGriFilename = osGriFilename;
1459 60 : poDS->nRasterXSize = nXSize;
1460 60 : poDS->nRasterYSize = nYSize;
1461 60 : poDS->m_fpImage = fpImage;
1462 60 : poDS->m_bNativeOrder = true;
1463 60 : poDS->m_osBandOrder = osBandOrder.toupper();
1464 60 : poDS->m_bInitRaster = CPLFetchBool(papszOptions, "@INIT_RASTER", true);
1465 :
1466 60 : const char *pszPixelType = CSLFetchNameValue(papszOptions, "PIXELTYPE");
1467 61 : const bool bByteSigned = (eType == GDT_Byte && pszPixelType &&
1468 1 : EQUAL(pszPixelType, "SIGNEDBYTE"));
1469 60 : if (bByteSigned)
1470 1 : poDS->m_bSignedByte = true;
1471 :
1472 169 : for (int i = 1; i <= nBandsIn; i++)
1473 : {
1474 : RRASTERRasterBand *poBand =
1475 109 : new RRASTERRasterBand(poDS, i, fpImage, nBandOffset * (i - 1),
1476 109 : nPixelOffset, nLineOffset, eType, true);
1477 109 : poDS->SetBand(i, poBand);
1478 : }
1479 :
1480 60 : return poDS;
1481 : }
1482 :
1483 : /************************************************************************/
1484 : /* CreateCopy() */
1485 : /************************************************************************/
1486 :
1487 47 : GDALDataset *RRASTERDataset::CreateCopy(const char *pszFilename,
1488 : GDALDataset *poSrcDS, int bStrict,
1489 : char **papszOptions,
1490 : GDALProgressFunc pfnProgress,
1491 : void *pProgressData)
1492 :
1493 : {
1494 : // Proceed with normal copying using the default createcopy operators.
1495 : GDALDriver *poDriver =
1496 47 : reinterpret_cast<GDALDriver *>(GDALGetDriverByName("RRASTER"));
1497 :
1498 47 : char **papszAdjustedOptions = CSLDuplicate(papszOptions);
1499 : papszAdjustedOptions =
1500 47 : CSLSetNameValue(papszAdjustedOptions, "@INIT_RASTER", "NO");
1501 47 : GDALDataset *poOutDS = poDriver->DefaultCreateCopy(
1502 : pszFilename, poSrcDS, bStrict, papszAdjustedOptions, pfnProgress,
1503 : pProgressData);
1504 47 : CSLDestroy(papszAdjustedOptions);
1505 :
1506 47 : if (poOutDS != nullptr)
1507 37 : poOutDS->FlushCache(false);
1508 :
1509 47 : return poOutDS;
1510 : }
1511 :
1512 : /************************************************************************/
1513 : /* GDALRegister_RRASTER() */
1514 : /************************************************************************/
1515 :
1516 1911 : void GDALRegister_RRASTER()
1517 :
1518 : {
1519 1911 : if (GDALGetDriverByName("RRASTER") != nullptr)
1520 282 : return;
1521 :
1522 1629 : GDALDriver *poDriver = new GDALDriver();
1523 :
1524 1629 : poDriver->SetDescription("RRASTER");
1525 1629 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1526 1629 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "grd");
1527 1629 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "R Raster");
1528 1629 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
1529 1629 : "drivers/raster/rraster.html");
1530 1629 : poDriver->SetMetadataItem(
1531 : GDAL_DMD_CREATIONDATATYPES,
1532 1629 : "Byte Int8 Int16 UInt16 Int32 UInt32 Float32 Float64");
1533 :
1534 1629 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1535 :
1536 1629 : poDriver->SetMetadataItem(
1537 : GDAL_DMD_CREATIONOPTIONLIST,
1538 : "<CreationOptionList>"
1539 : " <Option name='PIXELTYPE' type='string' description='(deprecated, "
1540 : "use Int8) "
1541 : "By setting this to SIGNEDBYTE, a new Byte file can be forced to be "
1542 : "written "
1543 : "as signed byte'/>"
1544 : " <Option name='INTERLEAVE' type='string-select' default='BIL'>"
1545 : " <Value>BIP</Value>"
1546 : " <Value>BIL</Value>"
1547 : " <Value>BSQ</Value>"
1548 : " </Option>"
1549 1629 : "</CreationOptionList>");
1550 :
1551 1629 : poDriver->pfnOpen = RRASTERDataset::Open;
1552 1629 : poDriver->pfnIdentify = RRASTERDataset::Identify;
1553 1629 : poDriver->pfnCreate = RRASTERDataset::Create;
1554 1629 : poDriver->pfnCreateCopy = RRASTERDataset::CreateCopy;
1555 :
1556 1629 : poDriver->SetMetadataItem(GDAL_DCAP_UPDATE, "YES");
1557 1629 : poDriver->SetMetadataItem(GDAL_DMD_UPDATE_ITEMS, "GeoTransform SRS NoData "
1558 : "RasterValues "
1559 1629 : "DatasetMetadata");
1560 :
1561 1629 : GetGDALDriverManager()->RegisterDriver(poDriver);
1562 : }
|
