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(GDALProgressFunc = nullptr, void * = nullptr) 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 : CSLConstList papszOptions);
72 : static GDALDataset *CreateCopy(const char *pszFilename,
73 : GDALDataset *poSrcDS, int bStrict,
74 : CSLConstList 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(CSLConstList 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(GDALProgressFunc, void *)
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(CSLConstList papszMetadata,
874 : const char *pszDomain)
875 : {
876 31 : if (pszDomain == nullptr || EQUAL(pszDomain, ""))
877 : {
878 31 : m_osCreator = CSLFetchNameValueDef(papszMetadata, "CREATOR", "");
879 31 : m_osCreated = CSLFetchNameValueDef(papszMetadata, "CREATED", "");
880 31 : SetHeaderDirty();
881 : }
882 31 : return GDALDataset::SetMetadata(papszMetadata, pszDomain);
883 : }
884 :
885 : /************************************************************************/
886 : /* SetMetadataItem() */
887 : /************************************************************************/
888 :
889 2 : CPLErr RRASTERDataset::SetMetadataItem(const char *pszName,
890 : const char *pszValue,
891 : const char *pszDomain)
892 : {
893 2 : if (pszDomain == nullptr || EQUAL(pszDomain, ""))
894 : {
895 2 : if (EQUAL(pszName, "CREATOR"))
896 : {
897 1 : m_osCreator = pszValue ? pszValue : "";
898 1 : SetHeaderDirty();
899 : }
900 2 : if (EQUAL(pszName, "CREATED"))
901 : {
902 1 : m_osCreated = pszValue ? pszValue : "";
903 1 : SetHeaderDirty();
904 : }
905 : }
906 2 : return GDALDataset::SetMetadataItem(pszName, pszValue, pszDomain);
907 : }
908 :
909 : /************************************************************************/
910 : /* Identify() */
911 : /************************************************************************/
912 :
913 61832 : int RRASTERDataset::Identify(GDALOpenInfo *poOpenInfo)
914 :
915 : {
916 66624 : if (poOpenInfo->nHeaderBytes < 40 || poOpenInfo->fpL == nullptr ||
917 4792 : !poOpenInfo->IsExtensionEqualToCI("grd"))
918 : {
919 61640 : return FALSE;
920 : }
921 :
922 : // We need to ingest enough bytes as the wkt= entry can be quite large
923 192 : if (poOpenInfo->nHeaderBytes <= 1024)
924 149 : poOpenInfo->TryToIngest(4096);
925 :
926 192 : if (strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
927 173 : "ncols") == nullptr ||
928 173 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
929 173 : "nrows") == nullptr ||
930 173 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
931 173 : "xmin") == nullptr ||
932 173 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
933 173 : "ymin") == nullptr ||
934 173 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
935 173 : "xmax") == nullptr ||
936 173 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
937 173 : "ymax") == nullptr ||
938 173 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
939 : "datatype") == nullptr)
940 : {
941 19 : return FALSE;
942 : }
943 :
944 173 : return TRUE;
945 : }
946 :
947 : /************************************************************************/
948 : /* ComputeSpacing() */
949 : /************************************************************************/
950 :
951 147 : bool RRASTERDataset::ComputeSpacings(const CPLString &osBandOrder, int nCols,
952 : int nRows, int l_nBands, GDALDataType eDT,
953 : int &nPixelOffset, int &nLineOffset,
954 : vsi_l_offset &nBandOffset)
955 : {
956 147 : nPixelOffset = 0;
957 147 : nLineOffset = 0;
958 147 : nBandOffset = 0;
959 147 : const int nPixelSize = GDALGetDataTypeSizeBytes(eDT);
960 147 : if (l_nBands == 1 || EQUAL(osBandOrder, "BIL"))
961 : {
962 137 : nPixelOffset = nPixelSize;
963 137 : if (l_nBands != 0 && nPixelSize != 0 &&
964 137 : nCols > INT_MAX / (l_nBands * nPixelSize))
965 : {
966 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many columns");
967 0 : return false;
968 : }
969 137 : nLineOffset = nPixelSize * nCols * l_nBands;
970 137 : nBandOffset = static_cast<vsi_l_offset>(nPixelSize) * nCols;
971 : }
972 10 : else if (EQUAL(osBandOrder, "BIP"))
973 : {
974 3 : if (l_nBands != 0 && nPixelSize != 0 &&
975 3 : nCols > INT_MAX / (l_nBands * nPixelSize))
976 : {
977 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many columns");
978 0 : return false;
979 : }
980 3 : nPixelOffset = nPixelSize * l_nBands;
981 3 : nLineOffset = nPixelSize * nCols * l_nBands;
982 3 : nBandOffset = nPixelSize;
983 : }
984 7 : else if (EQUAL(osBandOrder, "BSQ"))
985 : {
986 7 : if (nPixelSize != 0 && nCols > INT_MAX / nPixelSize)
987 : {
988 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many columns");
989 0 : return false;
990 : }
991 7 : nPixelOffset = nPixelSize;
992 7 : nLineOffset = nPixelSize * nCols;
993 7 : nBandOffset = static_cast<vsi_l_offset>(nLineOffset) * nRows;
994 : }
995 0 : else if (l_nBands > 1)
996 : {
997 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unknown bandorder");
998 0 : return false;
999 : }
1000 147 : return true;
1001 : }
1002 :
1003 : /************************************************************************/
1004 : /* CastToFloat() */
1005 : /************************************************************************/
1006 :
1007 0 : static float CastToFloat(double dfVal)
1008 : {
1009 0 : if (std::isinf(dfVal) || std::isnan(dfVal) ||
1010 0 : (dfVal >= -std::numeric_limits<float>::max() &&
1011 0 : dfVal <= std::numeric_limits<float>::max()))
1012 : {
1013 0 : return static_cast<float>(dfVal);
1014 : }
1015 0 : return std::numeric_limits<float>::quiet_NaN();
1016 : }
1017 :
1018 : /************************************************************************/
1019 : /* Open() */
1020 : /************************************************************************/
1021 :
1022 84 : GDALDataset *RRASTERDataset::Open(GDALOpenInfo *poOpenInfo)
1023 : {
1024 84 : if (!Identify(poOpenInfo))
1025 0 : return nullptr;
1026 :
1027 168 : auto poDS = std::make_unique<RRASTERDataset>();
1028 :
1029 84 : const char *pszLine = nullptr;
1030 84 : int nRows = 0;
1031 84 : int nCols = 0;
1032 84 : double dfXMin = 0.0;
1033 84 : double dfYMin = 0.0;
1034 84 : double dfXMax = 0.0;
1035 84 : double dfYMax = 0.0;
1036 84 : int l_nBands = 1;
1037 168 : CPLString osDataType;
1038 168 : CPLString osProjection;
1039 168 : CPLString osWkt;
1040 168 : CPLString osByteOrder;
1041 168 : CPLString osNoDataValue("NA");
1042 168 : CPLString osMinValue;
1043 168 : CPLString osMaxValue;
1044 168 : CPLString osLayerName;
1045 168 : CPLString osRatNames;
1046 168 : CPLString osRatTypes;
1047 168 : CPLString osRatValues;
1048 84 : bool bInLegend = false;
1049 84 : VSIRewindL(poOpenInfo->fpL);
1050 1592 : while ((pszLine = CPLReadLine2L(poOpenInfo->fpL, 1024 * 1024, nullptr)) !=
1051 : nullptr)
1052 : {
1053 1508 : if (pszLine[0] == '[')
1054 : {
1055 276 : bInLegend = EQUAL(pszLine, "[legend]");
1056 276 : continue;
1057 : }
1058 1232 : if (bInLegend)
1059 : {
1060 6 : poDS->m_osLegend += pszLine;
1061 6 : poDS->m_osLegend += "\n";
1062 : }
1063 1232 : char *pszKey = nullptr;
1064 1232 : const char *pszValue = CPLParseNameValue(pszLine, &pszKey);
1065 1232 : if (pszKey && pszValue)
1066 : {
1067 1232 : if (EQUAL(pszKey, "creator"))
1068 12 : poDS->m_osCreator = pszValue;
1069 1232 : if (EQUAL(pszKey, "created"))
1070 12 : poDS->m_osCreated = pszValue;
1071 1220 : else if (EQUAL(pszKey, "ncols"))
1072 84 : nCols = atoi(pszValue);
1073 1136 : else if (EQUAL(pszKey, "nrows"))
1074 84 : nRows = atoi(pszValue);
1075 1052 : else if (EQUAL(pszKey, "xmin"))
1076 84 : dfXMin = CPLAtof(pszValue);
1077 968 : else if (EQUAL(pszKey, "ymin"))
1078 84 : dfYMin = CPLAtof(pszValue);
1079 884 : else if (EQUAL(pszKey, "xmax"))
1080 84 : dfXMax = CPLAtof(pszValue);
1081 800 : else if (EQUAL(pszKey, "ymax"))
1082 84 : dfYMax = CPLAtof(pszValue);
1083 716 : else if (EQUAL(pszKey, "projection"))
1084 80 : osProjection = pszValue;
1085 636 : else if (EQUAL(pszKey, "wkt"))
1086 67 : osWkt = pszValue;
1087 569 : else if (EQUAL(pszKey, "nbands"))
1088 84 : l_nBands = atoi(pszValue);
1089 485 : else if (EQUAL(pszKey, "bandorder"))
1090 34 : poDS->m_osBandOrder = pszValue;
1091 451 : else if (EQUAL(pszKey, "datatype"))
1092 84 : osDataType = pszValue;
1093 367 : else if (EQUAL(pszKey, "byteorder"))
1094 84 : osByteOrder = pszValue;
1095 283 : else if (EQUAL(pszKey, "nodatavalue"))
1096 12 : osNoDataValue = pszValue;
1097 271 : else if (EQUAL(pszKey, "minvalue"))
1098 50 : osMinValue = pszValue;
1099 221 : else if (EQUAL(pszKey, "maxvalue"))
1100 50 : osMaxValue = pszValue;
1101 171 : else if (EQUAL(pszKey, "ratnames"))
1102 12 : osRatNames = pszValue;
1103 159 : else if (EQUAL(pszKey, "rattypes"))
1104 12 : osRatTypes = pszValue;
1105 147 : else if (EQUAL(pszKey, "ratvalues"))
1106 12 : osRatValues = pszValue;
1107 135 : else if (EQUAL(pszKey, "layername"))
1108 33 : osLayerName = pszValue;
1109 : }
1110 1232 : CPLFree(pszKey);
1111 : }
1112 84 : if (!GDALCheckDatasetDimensions(nCols, nRows))
1113 0 : return nullptr;
1114 84 : if (!GDALCheckBandCount(l_nBands, FALSE))
1115 0 : return nullptr;
1116 :
1117 84 : GDALDataType eDT = GDT_Unknown;
1118 84 : if (EQUAL(osDataType, "LOG1S"))
1119 0 : eDT = GDT_UInt8; // mapping TBC
1120 84 : else if (EQUAL(osDataType, "INT1S"))
1121 6 : eDT = GDT_Int8;
1122 78 : else if (EQUAL(osDataType, "INT2S"))
1123 5 : eDT = GDT_Int16;
1124 73 : else if (EQUAL(osDataType, "INT4S"))
1125 5 : eDT = GDT_Int32;
1126 : // else if( EQUAL(osDataType, "INT8S") )
1127 : // eDT = GDT_UInt64; // unhandled
1128 68 : else if (EQUAL(osDataType, "INT1U"))
1129 49 : eDT = GDT_UInt8;
1130 19 : else if (EQUAL(osDataType, "INT2U"))
1131 5 : eDT = GDT_UInt16;
1132 14 : else if (EQUAL(osDataType, "INT4U")) // Not documented
1133 5 : eDT = GDT_UInt32;
1134 9 : else if (EQUAL(osDataType, "FLT4S"))
1135 5 : eDT = GDT_Float32;
1136 4 : else if (EQUAL(osDataType, "FLT8S"))
1137 4 : eDT = GDT_Float64;
1138 : else
1139 : {
1140 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unhandled datatype=%s",
1141 : osDataType.c_str());
1142 0 : return nullptr;
1143 : }
1144 84 : if (l_nBands > 1 && poDS->m_osBandOrder.empty())
1145 : {
1146 0 : CPLError(CE_Failure, CPLE_AppDefined, "Missing 'bandorder'");
1147 0 : return nullptr;
1148 : }
1149 :
1150 84 : bool bNativeOrder = true;
1151 84 : if (EQUAL(osByteOrder, "little"))
1152 : {
1153 84 : bNativeOrder = CPL_TO_BOOL(CPL_IS_LSB);
1154 : }
1155 0 : else if (EQUAL(osByteOrder, "big"))
1156 : {
1157 0 : bNativeOrder = CPL_TO_BOOL(!CPL_IS_LSB);
1158 : }
1159 0 : else if (!EQUAL(osByteOrder, ""))
1160 : {
1161 0 : CPLError(CE_Warning, CPLE_AppDefined,
1162 : "Unhandled byteorder=%s. Assuming native order",
1163 : osByteOrder.c_str());
1164 : }
1165 :
1166 84 : int nPixelOffset = 0;
1167 84 : int nLineOffset = 0;
1168 84 : vsi_l_offset nBandOffset = 0;
1169 84 : if (!ComputeSpacings(poDS->m_osBandOrder, nCols, nRows, l_nBands, eDT,
1170 : nPixelOffset, nLineOffset, nBandOffset))
1171 : {
1172 0 : return nullptr;
1173 : }
1174 :
1175 168 : CPLString osDirname(CPLGetDirnameSafe(poOpenInfo->pszFilename));
1176 168 : CPLString osBasename(CPLGetBasenameSafe(poOpenInfo->pszFilename));
1177 168 : CPLString osGRDExtension(poOpenInfo->osExtension);
1178 168 : CPLString osGRIExtension((osGRDExtension[0] == 'g') ? "gri" : "GRI");
1179 84 : char **papszSiblings = poOpenInfo->GetSiblingFiles();
1180 84 : if (papszSiblings)
1181 : {
1182 84 : int iFile = CSLFindString(
1183 : papszSiblings,
1184 168 : CPLFormFilenameSafe(nullptr, osBasename, osGRIExtension).c_str());
1185 84 : if (iFile < 0)
1186 0 : return nullptr;
1187 84 : poDS->m_osGriFilename =
1188 168 : CPLFormFilenameSafe(osDirname, papszSiblings[iFile], nullptr);
1189 : }
1190 : else
1191 : {
1192 0 : poDS->m_osGriFilename =
1193 0 : CPLFormFilenameSafe(osDirname, osBasename, osGRIExtension);
1194 : }
1195 :
1196 : VSILFILE *fpImage =
1197 84 : VSIFOpenL(poDS->m_osGriFilename,
1198 84 : (poOpenInfo->eAccess == GA_Update) ? "rb+" : "rb");
1199 84 : if (fpImage == nullptr)
1200 : {
1201 0 : CPLError(CE_Failure, CPLE_FileIO, "Cannot open %s",
1202 0 : poDS->m_osGriFilename.c_str());
1203 0 : return nullptr;
1204 : }
1205 :
1206 84 : if (!RAWDatasetCheckMemoryUsage(nCols, nRows, l_nBands,
1207 : GDALGetDataTypeSizeBytes(eDT), nPixelOffset,
1208 : nLineOffset, 0, nBandOffset, fpImage))
1209 : {
1210 0 : VSIFCloseL(fpImage);
1211 0 : return nullptr;
1212 : }
1213 :
1214 84 : poDS->eAccess = poOpenInfo->eAccess;
1215 84 : poDS->nRasterXSize = nCols;
1216 84 : poDS->nRasterYSize = nRows;
1217 84 : poDS->m_bGeoTransformValid = true;
1218 84 : poDS->m_gt[0] = dfXMin;
1219 84 : poDS->m_gt[1] = (dfXMax - dfXMin) / nCols;
1220 84 : poDS->m_gt[2] = 0.0;
1221 84 : poDS->m_gt[3] = dfYMax;
1222 84 : poDS->m_gt[4] = 0.0;
1223 84 : poDS->m_gt[5] = -(dfYMax - dfYMin) / nRows;
1224 84 : poDS->m_fpImage = fpImage;
1225 84 : poDS->m_bNativeOrder = bNativeOrder;
1226 :
1227 84 : if (osWkt.empty())
1228 : {
1229 17 : if (!osProjection.empty())
1230 : {
1231 13 : poDS->m_oSRS.importFromProj4(osProjection.c_str());
1232 : }
1233 : }
1234 : else
1235 : {
1236 67 : poDS->m_oSRS.importFromWkt(osWkt.c_str());
1237 : }
1238 :
1239 84 : if (!poDS->m_osCreator.empty())
1240 12 : poDS->GDALDataset::SetMetadataItem("CREATOR", poDS->m_osCreator);
1241 :
1242 84 : if (!poDS->m_osCreated.empty())
1243 12 : poDS->GDALDataset::SetMetadataItem("CREATED", poDS->m_osCreated);
1244 :
1245 : // Instantiate RAT
1246 84 : if (!osRatNames.empty() && !osRatTypes.empty() && !osRatValues.empty())
1247 : {
1248 24 : CPLStringList aosRatNames(CSLTokenizeString2(osRatNames, ":", 0));
1249 24 : CPLStringList aosRatTypes(CSLTokenizeString2(osRatTypes, ":", 0));
1250 24 : CPLStringList aosRatValues(CSLTokenizeString2(osRatValues, ":", 0));
1251 12 : if (aosRatNames.size() >= 1 &&
1252 24 : aosRatNames.size() == aosRatTypes.size() &&
1253 12 : (aosRatValues.size() % aosRatNames.size()) == 0)
1254 : {
1255 12 : bool bIsCompatibleOfCT = false;
1256 12 : const int nValues = aosRatValues.size() / aosRatNames.size();
1257 20 : if ((aosRatNames.size() == 4 || aosRatNames.size() == 5) &&
1258 8 : EQUAL(aosRatNames[1], "red") &&
1259 8 : EQUAL(aosRatNames[2], "green") &&
1260 8 : EQUAL(aosRatNames[3], "blue") &&
1261 8 : (aosRatNames.size() == 4 || EQUAL(aosRatNames[4], "alpha")) &&
1262 8 : EQUAL(aosRatTypes[0], "integer") &&
1263 8 : EQUAL(aosRatTypes[1], "integer") &&
1264 8 : EQUAL(aosRatTypes[2], "integer") &&
1265 32 : EQUAL(aosRatTypes[3], "integer") &&
1266 8 : (aosRatTypes.size() == 4 || EQUAL(aosRatTypes[4], "integer")))
1267 : {
1268 8 : bIsCompatibleOfCT = true;
1269 8 : poDS->m_poCT.reset(new GDALColorTable());
1270 24 : for (int i = 0; i < nValues; i++)
1271 : {
1272 16 : const int nIndex = atoi(aosRatValues[i]);
1273 16 : if (nIndex >= 0 && nIndex < 65536)
1274 : {
1275 16 : const int nRed = atoi(aosRatValues[nValues + i]);
1276 16 : const int nGreen = atoi(aosRatValues[2 * nValues + i]);
1277 16 : const int nBlue = atoi(aosRatValues[3 * nValues + i]);
1278 : const int nAlpha =
1279 16 : aosRatTypes.size() == 4
1280 16 : ? 255
1281 8 : : atoi(aosRatValues[4 * nValues + i]);
1282 : const GDALColorEntry oEntry = {
1283 : static_cast<short>(nRed),
1284 : static_cast<short>(nGreen),
1285 : static_cast<short>(nBlue),
1286 16 : static_cast<short>(nAlpha)};
1287 :
1288 16 : poDS->m_poCT->SetColorEntry(nIndex, &oEntry);
1289 : }
1290 : else
1291 : {
1292 0 : bIsCompatibleOfCT = false;
1293 0 : poDS->m_poCT.reset();
1294 0 : break;
1295 : }
1296 : }
1297 : }
1298 :
1299 : // cppcheck-suppress knownConditionTrueFalse
1300 12 : if (!bIsCompatibleOfCT)
1301 : {
1302 4 : poDS->m_poRAT.reset(new GDALDefaultRasterAttributeTable());
1303 44 : for (int i = 0; i < aosRatNames.size(); i++)
1304 : {
1305 120 : poDS->m_poRAT->CreateColumn(
1306 40 : aosRatNames[i],
1307 40 : EQUAL(aosRatTypes[i], "integer") ? GFT_Integer
1308 12 : : EQUAL(aosRatTypes[i], "numeric") ? GFT_Real
1309 : : GFT_String,
1310 40 : EQUAL(aosRatNames[i], "red") ? GFU_Red
1311 68 : : EQUAL(aosRatNames[i], "green") ? GFU_Green
1312 60 : : EQUAL(aosRatNames[i], "blue") ? GFU_Blue
1313 52 : : EQUAL(aosRatNames[i], "alpha") ? GFU_Alpha
1314 44 : : EQUAL(aosRatNames[i], "name") ? GFU_Name
1315 20 : : EQUAL(aosRatNames[i], "pixelcount") ? GFU_PixelCount
1316 40 : : GFU_Generic);
1317 : }
1318 12 : for (int i = 0; i < nValues; i++)
1319 : {
1320 88 : for (int j = 0; j < aosRatTypes.size(); j++)
1321 : {
1322 80 : if (poDS->m_poRAT->GetTypeOfCol(j) == GFT_Integer)
1323 : {
1324 56 : poDS->m_poRAT->SetValue(
1325 56 : i, j, atoi(aosRatValues[j * nValues + i]));
1326 : }
1327 24 : else if (poDS->m_poRAT->GetTypeOfCol(j) == GFT_Real)
1328 : {
1329 16 : poDS->m_poRAT->SetValue(
1330 8 : i, j, CPLAtof(aosRatValues[j * nValues + i]));
1331 : }
1332 : else
1333 : {
1334 16 : poDS->m_poRAT->SetValue(
1335 16 : i, j, aosRatValues[j * nValues + i]);
1336 : }
1337 : }
1338 : }
1339 : }
1340 : }
1341 : }
1342 :
1343 168 : CPLStringList aosMinValues(CSLTokenizeString2(osMinValue, ":", 0));
1344 168 : CPLStringList aosMaxValues(CSLTokenizeString2(osMaxValue, ":", 0));
1345 :
1346 168 : CPLStringList aosLayerNames(CSLTokenizeString2(osLayerName, ":", 0));
1347 237 : for (int i = 1; i <= l_nBands; i++)
1348 : {
1349 : auto poBand = std::make_unique<RRASTERRasterBand>(
1350 153 : poDS.get(), i, fpImage, nBandOffset * (i - 1), nPixelOffset,
1351 153 : nLineOffset, eDT, bNativeOrder);
1352 153 : if (!poBand->IsValid())
1353 : {
1354 0 : return nullptr;
1355 : }
1356 153 : if (!osNoDataValue.empty() && !EQUAL(osNoDataValue, "NA"))
1357 : {
1358 10 : double dfNoDataValue = CPLAtof(osNoDataValue);
1359 10 : if (eDT == GDT_Float32)
1360 0 : dfNoDataValue = CastToFloat(dfNoDataValue);
1361 10 : poBand->m_bHasNoDataValue = true;
1362 10 : poBand->m_dfNoDataValue = dfNoDataValue;
1363 : }
1364 227 : if (i - 1 < static_cast<int>(aosMinValues.size()) &&
1365 74 : i - 1 < static_cast<int>(aosMaxValues.size()))
1366 : {
1367 148 : poBand->SetMinMax(CPLAtof(aosMinValues[i - 1]),
1368 74 : CPLAtof(aosMaxValues[i - 1]));
1369 : }
1370 153 : if (i - 1 < static_cast<int>(aosLayerNames.size()))
1371 : {
1372 156 : const CPLString osName(aosLayerNames[i - 1]);
1373 78 : poBand->GDALRasterBand::SetDescription(osName);
1374 78 : if (EQUAL(osName, "red"))
1375 15 : poBand->SetColorInterpretation(GCI_RedBand);
1376 63 : else if (EQUAL(osName, "green"))
1377 15 : poBand->SetColorInterpretation(GCI_GreenBand);
1378 48 : else if (EQUAL(osName, "blue"))
1379 15 : poBand->SetColorInterpretation(GCI_BlueBand);
1380 33 : else if (EQUAL(osName, "alpha"))
1381 15 : poBand->SetColorInterpretation(GCI_AlphaBand);
1382 : }
1383 153 : poBand->m_poRAT = poDS->m_poRAT;
1384 153 : poBand->m_poCT = poDS->m_poCT;
1385 153 : if (poBand->m_poCT)
1386 8 : poBand->SetColorInterpretation(GCI_PaletteIndex);
1387 153 : poDS->SetBand(i, std::move(poBand));
1388 : }
1389 :
1390 84 : return poDS.release();
1391 : }
1392 :
1393 : /************************************************************************/
1394 : /* Create() */
1395 : /************************************************************************/
1396 :
1397 80 : GDALDataset *RRASTERDataset::Create(const char *pszFilename, int nXSize,
1398 : int nYSize, int nBandsIn,
1399 : GDALDataType eType,
1400 : CSLConstList papszOptions)
1401 :
1402 : {
1403 : // Verify input options.
1404 80 : if (nBandsIn <= 0)
1405 : {
1406 1 : CPLError(CE_Failure, CPLE_NotSupported,
1407 : "RRASTER driver does not support %d bands.", nBandsIn);
1408 1 : return nullptr;
1409 : }
1410 :
1411 79 : if (eType != GDT_UInt8 && eType != GDT_Int8 && eType != GDT_UInt16 &&
1412 32 : eType != GDT_Int16 && eType != GDT_Int32 && eType != GDT_UInt32 &&
1413 20 : eType != GDT_Float32 && eType != GDT_Float64)
1414 : {
1415 16 : CPLError(CE_Failure, CPLE_NotSupported, "Unsupported data type (%s).",
1416 : GDALGetDataTypeName(eType));
1417 16 : return nullptr;
1418 : }
1419 :
1420 126 : CPLString osGRDExtension(CPLGetExtensionSafe(pszFilename));
1421 63 : if (!EQUAL(osGRDExtension, "grd"))
1422 : {
1423 0 : CPLError(CE_Failure, CPLE_NotSupported,
1424 : "RRASTER driver only supports grd extension");
1425 0 : return nullptr;
1426 : }
1427 :
1428 63 : int nPixelOffset = 0;
1429 63 : int nLineOffset = 0;
1430 63 : vsi_l_offset nBandOffset = 0;
1431 : CPLString osBandOrder(
1432 126 : CSLFetchNameValueDef(papszOptions, "INTERLEAVE", "BIL"));
1433 63 : if (!ComputeSpacings(osBandOrder, nXSize, nYSize, nBandsIn, eType,
1434 : nPixelOffset, nLineOffset, nBandOffset))
1435 : {
1436 0 : return nullptr;
1437 : }
1438 :
1439 63 : const std::string osGRIExtension((osGRDExtension[0] == 'g') ? "gri"
1440 189 : : "GRI");
1441 : const std::string osGriFilename(
1442 126 : CPLResetExtensionSafe(pszFilename, osGRIExtension.c_str()));
1443 :
1444 : // Try to create the file.
1445 63 : VSILFILE *fpImage = VSIFOpenL(osGriFilename.c_str(), "wb+");
1446 :
1447 63 : if (fpImage == nullptr)
1448 : {
1449 3 : CPLError(CE_Failure, CPLE_OpenFailed,
1450 : "Attempt to create file `%s' failed.", osGriFilename.c_str());
1451 3 : return nullptr;
1452 : }
1453 :
1454 60 : RRASTERDataset *poDS = new RRASTERDataset;
1455 60 : poDS->eAccess = GA_Update;
1456 60 : poDS->m_bHeaderDirty = true;
1457 60 : poDS->m_osGriFilename = osGriFilename;
1458 60 : poDS->nRasterXSize = nXSize;
1459 60 : poDS->nRasterYSize = nYSize;
1460 60 : poDS->m_fpImage = fpImage;
1461 60 : poDS->m_bNativeOrder = true;
1462 60 : poDS->m_osBandOrder = osBandOrder.toupper();
1463 60 : poDS->m_bInitRaster = CPLFetchBool(papszOptions, "@INIT_RASTER", true);
1464 :
1465 60 : const char *pszPixelType = CSLFetchNameValue(papszOptions, "PIXELTYPE");
1466 61 : const bool bByteSigned = (eType == GDT_UInt8 && pszPixelType &&
1467 1 : EQUAL(pszPixelType, "SIGNEDBYTE"));
1468 60 : if (bByteSigned)
1469 1 : poDS->m_bSignedByte = true;
1470 :
1471 169 : for (int i = 1; i <= nBandsIn; i++)
1472 : {
1473 : RRASTERRasterBand *poBand =
1474 109 : new RRASTERRasterBand(poDS, i, fpImage, nBandOffset * (i - 1),
1475 109 : nPixelOffset, nLineOffset, eType, true);
1476 109 : poDS->SetBand(i, poBand);
1477 : }
1478 :
1479 60 : return poDS;
1480 : }
1481 :
1482 : /************************************************************************/
1483 : /* CreateCopy() */
1484 : /************************************************************************/
1485 :
1486 47 : GDALDataset *RRASTERDataset::CreateCopy(const char *pszFilename,
1487 : GDALDataset *poSrcDS, int bStrict,
1488 : CSLConstList papszOptions,
1489 : GDALProgressFunc pfnProgress,
1490 : void *pProgressData)
1491 :
1492 : {
1493 : // Proceed with normal copying using the default createcopy operators.
1494 : GDALDriver *poDriver =
1495 47 : reinterpret_cast<GDALDriver *>(GDALGetDriverByName("RRASTER"));
1496 :
1497 47 : char **papszAdjustedOptions = CSLDuplicate(papszOptions);
1498 : papszAdjustedOptions =
1499 47 : CSLSetNameValue(papszAdjustedOptions, "@INIT_RASTER", "NO");
1500 47 : GDALDataset *poOutDS = poDriver->DefaultCreateCopy(
1501 : pszFilename, poSrcDS, bStrict, papszAdjustedOptions, pfnProgress,
1502 : pProgressData);
1503 47 : CSLDestroy(papszAdjustedOptions);
1504 :
1505 47 : if (poOutDS != nullptr)
1506 37 : poOutDS->FlushCache(false);
1507 :
1508 47 : return poOutDS;
1509 : }
1510 :
1511 : /************************************************************************/
1512 : /* GDALRegister_RRASTER() */
1513 : /************************************************************************/
1514 :
1515 2059 : void GDALRegister_RRASTER()
1516 :
1517 : {
1518 2059 : if (GDALGetDriverByName("RRASTER") != nullptr)
1519 283 : return;
1520 :
1521 1776 : GDALDriver *poDriver = new GDALDriver();
1522 :
1523 1776 : poDriver->SetDescription("RRASTER");
1524 1776 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1525 1776 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "grd");
1526 1776 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "R Raster");
1527 1776 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
1528 1776 : "drivers/raster/rraster.html");
1529 1776 : poDriver->SetMetadataItem(
1530 : GDAL_DMD_CREATIONDATATYPES,
1531 1776 : "Byte Int8 Int16 UInt16 Int32 UInt32 Float32 Float64");
1532 :
1533 1776 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1534 :
1535 1776 : poDriver->SetMetadataItem(
1536 : GDAL_DMD_CREATIONOPTIONLIST,
1537 : "<CreationOptionList>"
1538 : " <Option name='PIXELTYPE' type='string' description='(deprecated, "
1539 : "use Int8) "
1540 : "By setting this to SIGNEDBYTE, a new Byte file can be forced to be "
1541 : "written "
1542 : "as signed byte'/>"
1543 : " <Option name='INTERLEAVE' type='string-select' default='BIL'>"
1544 : " <Value>BIP</Value>"
1545 : " <Value>BIL</Value>"
1546 : " <Value>BSQ</Value>"
1547 : " </Option>"
1548 1776 : "</CreationOptionList>");
1549 :
1550 1776 : poDriver->pfnOpen = RRASTERDataset::Open;
1551 1776 : poDriver->pfnIdentify = RRASTERDataset::Identify;
1552 1776 : poDriver->pfnCreate = RRASTERDataset::Create;
1553 1776 : poDriver->pfnCreateCopy = RRASTERDataset::CreateCopy;
1554 :
1555 1776 : poDriver->SetMetadataItem(GDAL_DCAP_UPDATE, "YES");
1556 1776 : poDriver->SetMetadataItem(GDAL_DMD_UPDATE_ITEMS, "GeoTransform SRS NoData "
1557 : "RasterValues "
1558 1776 : "DatasetMetadata");
1559 :
1560 1776 : GetGDALDriverManager()->RegisterDriver(poDriver);
1561 : }
|
