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 : double m_adfGeoTransform[6]{0, 1, 0, 0, 0, -1};
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(double *) override;
79 : CPLErr SetGeoTransform(double *padfGeoTransform) 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 173 : void SetHeaderDirty()
89 : {
90 173 : m_bHeaderDirty = true;
91 173 : }
92 :
93 : void InitImageIfNeeded();
94 :
95 80 : inline bool IsSignedByte() const
96 : {
97 80 : 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 260 : RRASTERRasterBand::RRASTERRasterBand(GDALDataset *poDSIn, int nBandIn,
153 : VSILFILE *fpRawIn,
154 : vsi_l_offset nImgOffsetIn,
155 : int nPixelOffsetIn, int nLineOffsetIn,
156 : GDALDataType eDataTypeIn,
157 260 : int bNativeOrderIn)
158 : : RawRasterBand(poDSIn, nBandIn, fpRawIn, nImgOffsetIn, nPixelOffsetIn,
159 : nLineOffsetIn, eDataTypeIn, bNativeOrderIn,
160 260 : RawRasterBand::OwnFP::NO)
161 : {
162 260 : }
163 :
164 : /************************************************************************/
165 : /* SetMinMax() */
166 : /************************************************************************/
167 :
168 73 : void RRASTERRasterBand::SetMinMax(double dfMin, double dfMax)
169 : {
170 73 : m_dfMin = dfMin;
171 73 : m_dfMax = dfMax;
172 73 : }
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 86 : GDALColorTable *RRASTERRasterBand::GetColorTable()
209 : {
210 86 : 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 110 : GDALRasterAttributeTable *RRASTERRasterBand::GetDefaultRAT()
238 : {
239 110 : 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 259 : double RRASTERRasterBand::GetNoDataValue(int *pbSuccess)
282 : {
283 259 : if (pbSuccess)
284 259 : *pbSuccess = m_bHasNoDataValue;
285 259 : return m_dfNoDataValue;
286 : }
287 :
288 : /************************************************************************/
289 : /* SetNoDataValue() */
290 : /************************************************************************/
291 :
292 2 : CPLErr RRASTERRasterBand::SetNoDataValue(double dfNoData)
293 : {
294 2 : RRASTERDataset *poGDS = static_cast<RRASTERDataset *>(poDS);
295 2 : if (poGDS->GetAccess() != GA_Update)
296 0 : return CE_Failure;
297 :
298 2 : m_bHasNoDataValue = true;
299 2 : m_dfNoDataValue = dfNoData;
300 2 : poGDS->SetHeaderDirty();
301 2 : return CE_None;
302 : }
303 :
304 : /************************************************************************/
305 : /* GetMinMax() */
306 : /************************************************************************/
307 :
308 : template <class T>
309 80 : static void GetMinMax(const T *buffer, int nBufXSize, int nBufYSize,
310 : GSpacing nPixelSpace, GSpacing nLineSpace,
311 : double dfNoDataValue, double &dfMin, double &dfMax)
312 : {
313 702 : for (int iY = 0; iY < nBufYSize; iY++)
314 : {
315 9386 : for (int iX = 0; iX < nBufXSize; iX++)
316 : {
317 8764 : const double dfVal = buffer[iY * nLineSpace + iX * nPixelSpace];
318 8764 : if (dfVal != dfNoDataValue && !std::isnan(dfVal))
319 : {
320 8764 : dfMin = std::min(dfMin, dfVal);
321 8764 : dfMax = std::max(dfMax, dfVal);
322 : }
323 : }
324 : }
325 80 : }
326 :
327 80 : 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 80 : switch (eDT)
332 : {
333 66 : case GDT_Byte:
334 66 : GetMinMax(static_cast<const GByte *>(pBuffer), nBufXSize, nBufYSize,
335 : nPixelSpace, nLineSpace, dfNoDataValue, dfMin, dfMax);
336 66 : 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 80 : }
375 :
376 : /************************************************************************/
377 : /* IWriteBlock() */
378 : /************************************************************************/
379 :
380 20 : CPLErr RRASTERRasterBand::IWriteBlock(int nBlockXOff, int nBlockYOff,
381 : void *pImage)
382 : {
383 20 : RRASTERDataset *poGDS = static_cast<RRASTERDataset *>(poDS);
384 20 : poGDS->InitImageIfNeeded();
385 :
386 20 : int bGotNoDataValue = false;
387 20 : double dfNoDataValue = GetNoDataValue(&bGotNoDataValue);
388 20 : if (!bGotNoDataValue)
389 20 : dfNoDataValue = std::numeric_limits<double>::quiet_NaN();
390 20 : GetMinMax(pImage, poGDS->IsSignedByte() ? GDT_Int8 : eDataType, nBlockXSize,
391 20 : nBlockYSize, 1, nBlockXSize, dfNoDataValue, m_dfMin, m_dfMax);
392 40 : return RawRasterBand::IWriteBlock(nBlockXOff, nBlockYOff, pImage);
393 : }
394 :
395 : /************************************************************************/
396 : /* IRasterIO() */
397 : /************************************************************************/
398 :
399 323 : 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 323 : if (eRWFlag == GF_Write)
408 : {
409 60 : RRASTERDataset *poGDS = static_cast<RRASTERDataset *>(poDS);
410 60 : poGDS->InitImageIfNeeded();
411 :
412 60 : const int nDTSize = std::max(1, GDALGetDataTypeSizeBytes(eDataType));
413 60 : int bGotNoDataValue = false;
414 60 : double dfNoDataValue = GetNoDataValue(&bGotNoDataValue);
415 60 : if (!bGotNoDataValue)
416 59 : dfNoDataValue = std::numeric_limits<double>::quiet_NaN();
417 60 : GetMinMax(pData, poGDS->IsSignedByte() ? GDT_Int8 : eDataType,
418 60 : nBufXSize, nBufYSize, nPixelSpace / nDTSize,
419 60 : nLineSpace / nDTSize, dfNoDataValue, m_dfMin, m_dfMax);
420 : }
421 323 : return RawRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
422 : pData, nBufXSize, nBufYSize, eBufType,
423 323 : nPixelSpace, nLineSpace, psExtraArg);
424 : }
425 :
426 : /************************************************************************/
427 : /* RRASTERDataset() */
428 : /************************************************************************/
429 :
430 142 : RRASTERDataset::RRASTERDataset()
431 : {
432 142 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
433 142 : }
434 :
435 : /************************************************************************/
436 : /* ~RRASTERDataset() */
437 : /************************************************************************/
438 :
439 284 : RRASTERDataset::~RRASTERDataset()
440 :
441 : {
442 142 : RRASTERDataset::Close();
443 284 : }
444 :
445 : /************************************************************************/
446 : /* Close() */
447 : /************************************************************************/
448 :
449 277 : CPLErr RRASTERDataset::Close()
450 : {
451 277 : CPLErr eErr = CE_None;
452 277 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
453 : {
454 142 : if (m_fpImage)
455 : {
456 142 : InitImageIfNeeded();
457 142 : if (RRASTERDataset::FlushCache(true) != CE_None)
458 2 : eErr = CE_Failure;
459 142 : if (VSIFCloseL(m_fpImage) != CE_None)
460 0 : eErr = CE_Failure;
461 : }
462 142 : if (m_bHeaderDirty)
463 59 : RewriteHeader();
464 :
465 142 : if (GDALPamDataset::Close() != CE_None)
466 0 : eErr = CE_Failure;
467 : }
468 277 : return eErr;
469 : }
470 :
471 : /************************************************************************/
472 : /* InitImageIfNeeded() */
473 : /************************************************************************/
474 :
475 222 : void RRASTERDataset::InitImageIfNeeded()
476 : {
477 222 : CPLAssert(m_fpImage);
478 222 : if (!m_bInitRaster)
479 201 : 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 59 : void RRASTERDataset::RewriteHeader()
509 : {
510 59 : VSILFILE *fp = VSIFOpenL(GetDescription(), "wb");
511 59 : if (!fp)
512 0 : return;
513 :
514 59 : VSIFPrintfL(fp, "[general]\n");
515 59 : if (!m_osCreator.empty())
516 4 : VSIFPrintfL(fp, "creator=%s\n", m_osCreator.c_str());
517 59 : if (!m_osCreated.empty())
518 4 : VSIFPrintfL(fp, "created=%s\n", m_osCreated.c_str());
519 :
520 59 : VSIFPrintfL(fp, "[data]\n");
521 59 : GDALDataType eDT = GetRasterBand(1)->GetRasterDataType();
522 115 : VSIFPrintfL(fp, "datatype=%s\n",
523 56 : (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 59 : int bGotNoDataValue = false;
534 59 : double dfNoDataValue = GetRasterBand(1)->GetNoDataValue(&bGotNoDataValue);
535 59 : if (bGotNoDataValue)
536 2 : VSIFPrintfL(fp, "nodatavalue=%.17g\n", dfNoDataValue);
537 :
538 : #if CPL_IS_LSB
539 59 : VSIFPrintfL(fp, "byteorder=%s\n", m_bNativeOrder ? "little" : "big");
540 : #else
541 : VSIFPrintfL(fp, "byteorder=%s\n", !m_bNativeOrder ? "little" : "big");
542 : #endif
543 59 : VSIFPrintfL(fp, "nbands=%d\n", nBands);
544 59 : if (nBands > 1)
545 20 : VSIFPrintfL(fp, "bandorder=%s\n", m_osBandOrder.c_str());
546 118 : CPLString osMinValue, osMaxValue;
547 119 : for (int i = 1; i <= nBands; i++)
548 : {
549 : RRASTERRasterBand *poBand =
550 81 : static_cast<RRASTERRasterBand *>(GetRasterBand(i));
551 81 : if (i > 1)
552 : {
553 22 : osMinValue += ":";
554 22 : osMaxValue += ":";
555 : }
556 81 : if (poBand->m_dfMin > poBand->m_dfMax)
557 : {
558 21 : osMinValue.clear();
559 21 : break;
560 : }
561 60 : osMinValue += CPLSPrintf("%.17g", poBand->m_dfMin);
562 60 : osMaxValue += CPLSPrintf("%.17g", poBand->m_dfMax);
563 : }
564 59 : if (!osMinValue.empty())
565 : {
566 38 : VSIFPrintfL(fp, "minvalue=%s\n", osMinValue.c_str());
567 38 : VSIFPrintfL(fp, "maxvalue=%s\n", osMaxValue.c_str());
568 : }
569 :
570 59 : GDALColorTable *poCT = GetRasterBand(1)->GetColorTable();
571 59 : GDALRasterAttributeTable *poRAT = GetRasterBand(1)->GetDefaultRAT();
572 59 : if (poCT == nullptr && poRAT == nullptr)
573 : {
574 56 : 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 59 : if (!m_osLegend.empty())
700 0 : VSIFPrintfL(fp, "[legend]\n%s", m_osLegend.c_str());
701 :
702 118 : CPLString osLayerName;
703 59 : bool bGotSignificantBandDesc = false;
704 167 : for (int i = 1; i <= nBands; i++)
705 : {
706 108 : GDALRasterBand *poBand = GetRasterBand(i);
707 108 : const char *pszDesc = poBand->GetDescription();
708 108 : if (EQUAL(pszDesc, ""))
709 : {
710 88 : GDALColorInterp eInterp = poBand->GetColorInterpretation();
711 88 : if (eInterp == GCI_RedBand)
712 : {
713 1 : bGotSignificantBandDesc = true;
714 1 : pszDesc = "red";
715 : }
716 87 : else if (eInterp == GCI_GreenBand)
717 : {
718 1 : bGotSignificantBandDesc = true;
719 1 : pszDesc = "green";
720 : }
721 86 : else if (eInterp == GCI_BlueBand)
722 : {
723 1 : bGotSignificantBandDesc = true;
724 1 : pszDesc = "blue";
725 : }
726 85 : else if (eInterp == GCI_AlphaBand)
727 : {
728 1 : bGotSignificantBandDesc = true;
729 1 : pszDesc = "alpha";
730 : }
731 : else
732 : {
733 84 : pszDesc = CPLSPrintf("Band%d", i);
734 : }
735 : }
736 : else
737 : {
738 20 : bGotSignificantBandDesc = true;
739 : }
740 108 : if (i > 1)
741 49 : osLayerName += ":";
742 108 : osLayerName += CPLString(pszDesc).replaceAll(':', '.');
743 : }
744 59 : 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 59 : VSIFPrintfL(fp, "[georeference]\n");
754 59 : VSIFPrintfL(fp, "nrows=%d\n", nRasterYSize);
755 59 : VSIFPrintfL(fp, "ncols=%d\n", nRasterXSize);
756 :
757 59 : VSIFPrintfL(fp, "xmin=%.17g\n", m_adfGeoTransform[0]);
758 59 : VSIFPrintfL(fp, "ymin=%.17g\n",
759 59 : m_adfGeoTransform[3] + nRasterYSize * m_adfGeoTransform[5]);
760 59 : VSIFPrintfL(fp, "xmax=%.17g\n",
761 59 : m_adfGeoTransform[0] + nRasterXSize * m_adfGeoTransform[1]);
762 59 : VSIFPrintfL(fp, "ymax=%.17g\n", m_adfGeoTransform[3]);
763 :
764 59 : if (!m_oSRS.IsEmpty())
765 : {
766 57 : char *pszProj4 = nullptr;
767 57 : m_oSRS.exportToProj4(&pszProj4);
768 57 : if (pszProj4)
769 : {
770 57 : VSIFPrintfL(fp, "projection=%s\n", pszProj4);
771 57 : VSIFree(pszProj4);
772 : }
773 57 : char *pszWKT = nullptr;
774 57 : const char *const apszOptions[] = {"FORMAT=WKT2_2019", nullptr};
775 57 : m_oSRS.exportToWkt(&pszWKT, apszOptions);
776 57 : if (pszWKT)
777 : {
778 57 : VSIFPrintfL(fp, "wkt=%s\n", pszWKT);
779 57 : VSIFree(pszWKT);
780 : }
781 : }
782 :
783 59 : VSIFCloseL(fp);
784 : }
785 :
786 : /************************************************************************/
787 : /* GetFileList() */
788 : /************************************************************************/
789 :
790 37 : char **RRASTERDataset::GetFileList()
791 :
792 : {
793 37 : char **papszFileList = RawDataset::GetFileList();
794 :
795 37 : papszFileList = CSLAddString(papszFileList, m_osGriFilename);
796 :
797 37 : return papszFileList;
798 : }
799 :
800 : /************************************************************************/
801 : /* GetGeoTransform() */
802 : /************************************************************************/
803 :
804 120 : CPLErr RRASTERDataset::GetGeoTransform(double *padfGeoTransform)
805 : {
806 120 : if (m_bGeoTransformValid)
807 : {
808 120 : memcpy(padfGeoTransform, m_adfGeoTransform, 6 * sizeof(double));
809 120 : return CE_None;
810 : }
811 0 : return CE_Failure;
812 : }
813 :
814 : /************************************************************************/
815 : /* SetGeoTransform() */
816 : /************************************************************************/
817 :
818 57 : CPLErr RRASTERDataset::SetGeoTransform(double *padfGeoTransform)
819 :
820 : {
821 57 : if (GetAccess() != GA_Update)
822 : {
823 0 : CPLError(CE_Failure, CPLE_NotSupported,
824 : "Cannot set geotransform on a read-only dataset");
825 0 : return CE_Failure;
826 : }
827 :
828 : // We only support non-rotated images with info in the .HDR file.
829 57 : if (padfGeoTransform[2] != 0.0 || padfGeoTransform[4] != 0.0)
830 : {
831 0 : CPLError(CE_Warning, CPLE_NotSupported,
832 : "Rotated / skewed images not supported");
833 0 : return GDALPamDataset::SetGeoTransform(padfGeoTransform);
834 : }
835 :
836 : // Record new geotransform.
837 57 : m_bGeoTransformValid = true;
838 57 : memcpy(m_adfGeoTransform, padfGeoTransform, sizeof(double) * 6);
839 57 : SetHeaderDirty();
840 :
841 57 : return CE_None;
842 : }
843 :
844 : /************************************************************************/
845 : /* GetSpatialRef() */
846 : /************************************************************************/
847 :
848 35 : const OGRSpatialReference *RRASTERDataset::GetSpatialRef() const
849 : {
850 35 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
851 : }
852 :
853 : /************************************************************************/
854 : /* SetSpatialRef() */
855 : /************************************************************************/
856 :
857 57 : CPLErr RRASTERDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
858 : {
859 57 : if (GetAccess() != GA_Update)
860 : {
861 0 : CPLError(CE_Failure, CPLE_NotSupported,
862 : "Cannot set projection on a read-only dataset");
863 0 : return CE_Failure;
864 : }
865 :
866 57 : m_oSRS.Clear();
867 57 : if (poSRS)
868 57 : m_oSRS = *poSRS;
869 57 : SetHeaderDirty();
870 :
871 57 : return CE_None;
872 : }
873 :
874 : /************************************************************************/
875 : /* SetMetadata() */
876 : /************************************************************************/
877 :
878 30 : CPLErr RRASTERDataset::SetMetadata(char **papszMetadata, const char *pszDomain)
879 : {
880 30 : if (pszDomain == nullptr || EQUAL(pszDomain, ""))
881 : {
882 30 : m_osCreator = CSLFetchNameValueDef(papszMetadata, "CREATOR", "");
883 30 : m_osCreated = CSLFetchNameValueDef(papszMetadata, "CREATED", "");
884 30 : SetHeaderDirty();
885 : }
886 30 : return GDALDataset::SetMetadata(papszMetadata, pszDomain);
887 : }
888 :
889 : /************************************************************************/
890 : /* SetMetadataItem() */
891 : /************************************************************************/
892 :
893 2 : CPLErr RRASTERDataset::SetMetadataItem(const char *pszName,
894 : const char *pszValue,
895 : const char *pszDomain)
896 : {
897 2 : if (pszDomain == nullptr || EQUAL(pszDomain, ""))
898 : {
899 2 : if (EQUAL(pszName, "CREATOR"))
900 : {
901 1 : m_osCreator = pszValue ? pszValue : "";
902 1 : SetHeaderDirty();
903 : }
904 2 : if (EQUAL(pszName, "CREATED"))
905 : {
906 1 : m_osCreated = pszValue ? pszValue : "";
907 1 : SetHeaderDirty();
908 : }
909 : }
910 2 : return GDALDataset::SetMetadataItem(pszName, pszValue, pszDomain);
911 : }
912 :
913 : /************************************************************************/
914 : /* Identify() */
915 : /************************************************************************/
916 :
917 51481 : int RRASTERDataset::Identify(GDALOpenInfo *poOpenInfo)
918 :
919 : {
920 54860 : if (poOpenInfo->nHeaderBytes < 40 || poOpenInfo->fpL == nullptr ||
921 3379 : !EQUAL(CPLGetExtension(poOpenInfo->pszFilename), "grd"))
922 : {
923 51285 : return FALSE;
924 : }
925 :
926 : // We need to ingest enough bytes as the wkt= entry can be quite large
927 196 : if (poOpenInfo->nHeaderBytes <= 1024)
928 147 : poOpenInfo->TryToIngest(4096);
929 :
930 190 : if (strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
931 171 : "ncols") == nullptr ||
932 171 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
933 171 : "nrows") == nullptr ||
934 171 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
935 171 : "xmin") == nullptr ||
936 171 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
937 171 : "ymin") == nullptr ||
938 171 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
939 171 : "xmax") == nullptr ||
940 171 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
941 171 : "ymax") == nullptr ||
942 171 : strstr(reinterpret_cast<const char *>(poOpenInfo->pabyHeader),
943 : "datatype") == nullptr)
944 : {
945 19 : return FALSE;
946 : }
947 :
948 171 : return TRUE;
949 : }
950 :
951 : /************************************************************************/
952 : /* ComputeSpacing() */
953 : /************************************************************************/
954 :
955 145 : bool RRASTERDataset::ComputeSpacings(const CPLString &osBandOrder, int nCols,
956 : int nRows, int l_nBands, GDALDataType eDT,
957 : int &nPixelOffset, int &nLineOffset,
958 : vsi_l_offset &nBandOffset)
959 : {
960 145 : nPixelOffset = 0;
961 145 : nLineOffset = 0;
962 145 : nBandOffset = 0;
963 145 : const int nPixelSize = GDALGetDataTypeSizeBytes(eDT);
964 145 : if (l_nBands == 1 || EQUAL(osBandOrder, "BIL"))
965 : {
966 139 : nPixelOffset = nPixelSize;
967 139 : if (l_nBands != 0 && nPixelSize != 0 &&
968 139 : nCols > INT_MAX / (l_nBands * nPixelSize))
969 : {
970 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many columns");
971 0 : return false;
972 : }
973 139 : nLineOffset = nPixelSize * nCols * l_nBands;
974 139 : nBandOffset = static_cast<vsi_l_offset>(nPixelSize) * nCols;
975 : }
976 6 : else if (EQUAL(osBandOrder, "BIP"))
977 : {
978 3 : if (l_nBands != 0 && nPixelSize != 0 &&
979 3 : nCols > INT_MAX / (l_nBands * nPixelSize))
980 : {
981 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many columns");
982 0 : return false;
983 : }
984 3 : nPixelOffset = nPixelSize * l_nBands;
985 3 : nLineOffset = nPixelSize * nCols * l_nBands;
986 3 : nBandOffset = nPixelSize;
987 : }
988 3 : else if (EQUAL(osBandOrder, "BSQ"))
989 : {
990 3 : if (nPixelSize != 0 && nCols > INT_MAX / nPixelSize)
991 : {
992 0 : CPLError(CE_Failure, CPLE_AppDefined, "Too many columns");
993 0 : return false;
994 : }
995 3 : nPixelOffset = nPixelSize;
996 3 : nLineOffset = nPixelSize * nCols;
997 3 : nBandOffset = static_cast<vsi_l_offset>(nLineOffset) * nRows;
998 : }
999 0 : else if (l_nBands > 1)
1000 : {
1001 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unknown bandorder");
1002 0 : return false;
1003 : }
1004 145 : return true;
1005 : }
1006 :
1007 : /************************************************************************/
1008 : /* CastToFloat() */
1009 : /************************************************************************/
1010 :
1011 0 : static float CastToFloat(double dfVal)
1012 : {
1013 0 : if (std::isinf(dfVal) || std::isnan(dfVal) ||
1014 0 : (dfVal >= -std::numeric_limits<float>::max() &&
1015 0 : dfVal <= std::numeric_limits<float>::max()))
1016 : {
1017 0 : return static_cast<float>(dfVal);
1018 : }
1019 0 : return std::numeric_limits<float>::quiet_NaN();
1020 : }
1021 :
1022 : /************************************************************************/
1023 : /* Open() */
1024 : /************************************************************************/
1025 :
1026 83 : GDALDataset *RRASTERDataset::Open(GDALOpenInfo *poOpenInfo)
1027 : {
1028 83 : if (!Identify(poOpenInfo))
1029 0 : return nullptr;
1030 :
1031 166 : auto poDS = std::make_unique<RRASTERDataset>();
1032 :
1033 83 : const char *pszLine = nullptr;
1034 83 : int nRows = 0;
1035 83 : int nCols = 0;
1036 83 : double dfXMin = 0.0;
1037 83 : double dfYMin = 0.0;
1038 83 : double dfXMax = 0.0;
1039 83 : double dfYMax = 0.0;
1040 83 : int l_nBands = 1;
1041 166 : CPLString osDataType;
1042 166 : CPLString osProjection;
1043 166 : CPLString osWkt;
1044 166 : CPLString osByteOrder;
1045 166 : CPLString osNoDataValue("NA");
1046 166 : CPLString osMinValue;
1047 166 : CPLString osMaxValue;
1048 166 : CPLString osLayerName;
1049 166 : CPLString osRatNames;
1050 166 : CPLString osRatTypes;
1051 166 : CPLString osRatValues;
1052 83 : bool bInLegend = false;
1053 83 : VSIRewindL(poOpenInfo->fpL);
1054 1574 : while ((pszLine = CPLReadLine2L(poOpenInfo->fpL, 1024 * 1024, nullptr)) !=
1055 : nullptr)
1056 : {
1057 1491 : if (pszLine[0] == '[')
1058 : {
1059 273 : bInLegend = EQUAL(pszLine, "[legend]");
1060 273 : continue;
1061 : }
1062 1218 : if (bInLegend)
1063 : {
1064 6 : poDS->m_osLegend += pszLine;
1065 6 : poDS->m_osLegend += "\n";
1066 : }
1067 1218 : char *pszKey = nullptr;
1068 1218 : const char *pszValue = CPLParseNameValue(pszLine, &pszKey);
1069 1218 : if (pszKey && pszValue)
1070 : {
1071 1218 : if (EQUAL(pszKey, "creator"))
1072 12 : poDS->m_osCreator = pszValue;
1073 1218 : if (EQUAL(pszKey, "created"))
1074 12 : poDS->m_osCreated = pszValue;
1075 1206 : else if (EQUAL(pszKey, "ncols"))
1076 83 : nCols = atoi(pszValue);
1077 1123 : else if (EQUAL(pszKey, "nrows"))
1078 83 : nRows = atoi(pszValue);
1079 1040 : else if (EQUAL(pszKey, "xmin"))
1080 83 : dfXMin = CPLAtof(pszValue);
1081 957 : else if (EQUAL(pszKey, "ymin"))
1082 83 : dfYMin = CPLAtof(pszValue);
1083 874 : else if (EQUAL(pszKey, "xmax"))
1084 83 : dfXMax = CPLAtof(pszValue);
1085 791 : else if (EQUAL(pszKey, "ymax"))
1086 83 : dfYMax = CPLAtof(pszValue);
1087 708 : else if (EQUAL(pszKey, "projection"))
1088 79 : osProjection = pszValue;
1089 629 : else if (EQUAL(pszKey, "wkt"))
1090 66 : osWkt = pszValue;
1091 563 : else if (EQUAL(pszKey, "nbands"))
1092 83 : l_nBands = atoi(pszValue);
1093 480 : else if (EQUAL(pszKey, "bandorder"))
1094 34 : poDS->m_osBandOrder = pszValue;
1095 446 : else if (EQUAL(pszKey, "datatype"))
1096 83 : osDataType = pszValue;
1097 363 : else if (EQUAL(pszKey, "byteorder"))
1098 83 : osByteOrder = pszValue;
1099 280 : else if (EQUAL(pszKey, "nodatavalue"))
1100 12 : osNoDataValue = pszValue;
1101 268 : else if (EQUAL(pszKey, "minvalue"))
1102 49 : osMinValue = pszValue;
1103 219 : else if (EQUAL(pszKey, "maxvalue"))
1104 49 : osMaxValue = pszValue;
1105 170 : else if (EQUAL(pszKey, "ratnames"))
1106 12 : osRatNames = pszValue;
1107 158 : else if (EQUAL(pszKey, "rattypes"))
1108 12 : osRatTypes = pszValue;
1109 146 : else if (EQUAL(pszKey, "ratvalues"))
1110 12 : osRatValues = pszValue;
1111 134 : else if (EQUAL(pszKey, "layername"))
1112 33 : osLayerName = pszValue;
1113 : }
1114 1218 : CPLFree(pszKey);
1115 : }
1116 83 : if (!GDALCheckDatasetDimensions(nCols, nRows))
1117 0 : return nullptr;
1118 83 : if (!GDALCheckBandCount(l_nBands, FALSE))
1119 0 : return nullptr;
1120 :
1121 83 : GDALDataType eDT = GDT_Unknown;
1122 83 : if (EQUAL(osDataType, "LOG1S"))
1123 0 : eDT = GDT_Byte; // mapping TBC
1124 83 : else if (EQUAL(osDataType, "INT1S"))
1125 6 : eDT = GDT_Int8;
1126 77 : else if (EQUAL(osDataType, "INT2S"))
1127 5 : eDT = GDT_Int16;
1128 72 : else if (EQUAL(osDataType, "INT4S"))
1129 5 : eDT = GDT_Int32;
1130 : // else if( EQUAL(osDataType, "INT8S") )
1131 : // eDT = GDT_UInt64; // unhandled
1132 67 : else if (EQUAL(osDataType, "INT1U"))
1133 48 : eDT = GDT_Byte;
1134 19 : else if (EQUAL(osDataType, "INT2U"))
1135 5 : eDT = GDT_UInt16;
1136 14 : else if (EQUAL(osDataType, "INT4U")) // Not documented
1137 5 : eDT = GDT_UInt32;
1138 9 : else if (EQUAL(osDataType, "FLT4S"))
1139 5 : eDT = GDT_Float32;
1140 4 : else if (EQUAL(osDataType, "FLT8S"))
1141 4 : eDT = GDT_Float64;
1142 : else
1143 : {
1144 0 : CPLError(CE_Failure, CPLE_AppDefined, "Unhandled datatype=%s",
1145 : osDataType.c_str());
1146 0 : return nullptr;
1147 : }
1148 83 : if (l_nBands > 1 && poDS->m_osBandOrder.empty())
1149 : {
1150 0 : CPLError(CE_Failure, CPLE_AppDefined, "Missing 'bandorder'");
1151 0 : return nullptr;
1152 : }
1153 :
1154 83 : bool bNativeOrder = true;
1155 83 : if (EQUAL(osByteOrder, "little"))
1156 : {
1157 83 : bNativeOrder = CPL_TO_BOOL(CPL_IS_LSB);
1158 : }
1159 0 : else if (EQUAL(osByteOrder, "big"))
1160 : {
1161 0 : bNativeOrder = CPL_TO_BOOL(!CPL_IS_LSB);
1162 : }
1163 0 : else if (!EQUAL(osByteOrder, ""))
1164 : {
1165 0 : CPLError(CE_Warning, CPLE_AppDefined,
1166 : "Unhandled byteorder=%s. Assuming native order",
1167 : osByteOrder.c_str());
1168 : }
1169 :
1170 83 : int nPixelOffset = 0;
1171 83 : int nLineOffset = 0;
1172 83 : vsi_l_offset nBandOffset = 0;
1173 83 : if (!ComputeSpacings(poDS->m_osBandOrder, nCols, nRows, l_nBands, eDT,
1174 : nPixelOffset, nLineOffset, nBandOffset))
1175 : {
1176 0 : return nullptr;
1177 : }
1178 :
1179 166 : CPLString osDirname(CPLGetDirname(poOpenInfo->pszFilename));
1180 166 : CPLString osBasename(CPLGetBasename(poOpenInfo->pszFilename));
1181 166 : CPLString osGRDExtension(CPLGetExtension(poOpenInfo->pszFilename));
1182 166 : CPLString osGRIExtension((osGRDExtension[0] == 'g') ? "gri" : "GRI");
1183 83 : char **papszSiblings = poOpenInfo->GetSiblingFiles();
1184 83 : if (papszSiblings)
1185 : {
1186 : int iFile =
1187 83 : CSLFindString(papszSiblings,
1188 : CPLFormFilename(nullptr, osBasename, osGRIExtension));
1189 83 : if (iFile < 0)
1190 0 : return nullptr;
1191 83 : poDS->m_osGriFilename =
1192 166 : CPLFormFilename(osDirname, papszSiblings[iFile], nullptr);
1193 : }
1194 : else
1195 : {
1196 0 : poDS->m_osGriFilename =
1197 0 : CPLFormFilename(osDirname, osBasename, osGRIExtension);
1198 : }
1199 :
1200 : VSILFILE *fpImage =
1201 83 : VSIFOpenL(poDS->m_osGriFilename,
1202 83 : (poOpenInfo->eAccess == GA_Update) ? "rb+" : "rb");
1203 83 : if (fpImage == nullptr)
1204 : {
1205 0 : CPLError(CE_Failure, CPLE_FileIO, "Cannot open %s",
1206 0 : poDS->m_osGriFilename.c_str());
1207 0 : return nullptr;
1208 : }
1209 :
1210 83 : if (!RAWDatasetCheckMemoryUsage(nCols, nRows, l_nBands,
1211 : GDALGetDataTypeSizeBytes(eDT), nPixelOffset,
1212 : nLineOffset, 0, nBandOffset, fpImage))
1213 : {
1214 0 : VSIFCloseL(fpImage);
1215 0 : return nullptr;
1216 : }
1217 :
1218 83 : poDS->eAccess = poOpenInfo->eAccess;
1219 83 : poDS->nRasterXSize = nCols;
1220 83 : poDS->nRasterYSize = nRows;
1221 83 : poDS->m_bGeoTransformValid = true;
1222 83 : poDS->m_adfGeoTransform[0] = dfXMin;
1223 83 : poDS->m_adfGeoTransform[1] = (dfXMax - dfXMin) / nCols;
1224 83 : poDS->m_adfGeoTransform[2] = 0.0;
1225 83 : poDS->m_adfGeoTransform[3] = dfYMax;
1226 83 : poDS->m_adfGeoTransform[4] = 0.0;
1227 83 : poDS->m_adfGeoTransform[5] = -(dfYMax - dfYMin) / nRows;
1228 83 : poDS->m_fpImage = fpImage;
1229 83 : poDS->m_bNativeOrder = bNativeOrder;
1230 :
1231 83 : if (osWkt.empty())
1232 : {
1233 17 : if (!osProjection.empty())
1234 : {
1235 13 : poDS->m_oSRS.importFromProj4(osProjection.c_str());
1236 : }
1237 : }
1238 : else
1239 : {
1240 66 : poDS->m_oSRS.importFromWkt(osWkt.c_str());
1241 : }
1242 :
1243 83 : if (!poDS->m_osCreator.empty())
1244 12 : poDS->GDALDataset::SetMetadataItem("CREATOR", poDS->m_osCreator);
1245 :
1246 83 : if (!poDS->m_osCreated.empty())
1247 12 : poDS->GDALDataset::SetMetadataItem("CREATED", poDS->m_osCreated);
1248 :
1249 : // Instantiate RAT
1250 83 : if (!osRatNames.empty() && !osRatTypes.empty() && !osRatValues.empty())
1251 : {
1252 24 : CPLStringList aosRatNames(CSLTokenizeString2(osRatNames, ":", 0));
1253 24 : CPLStringList aosRatTypes(CSLTokenizeString2(osRatTypes, ":", 0));
1254 24 : CPLStringList aosRatValues(CSLTokenizeString2(osRatValues, ":", 0));
1255 12 : if (aosRatNames.size() >= 1 &&
1256 24 : aosRatNames.size() == aosRatTypes.size() &&
1257 12 : (aosRatValues.size() % aosRatNames.size()) == 0)
1258 : {
1259 12 : bool bIsCompatibleOfCT = false;
1260 12 : const int nValues = aosRatValues.size() / aosRatNames.size();
1261 20 : if ((aosRatNames.size() == 4 || aosRatNames.size() == 5) &&
1262 8 : EQUAL(aosRatNames[1], "red") &&
1263 8 : EQUAL(aosRatNames[2], "green") &&
1264 8 : EQUAL(aosRatNames[3], "blue") &&
1265 8 : (aosRatNames.size() == 4 || EQUAL(aosRatNames[4], "alpha")) &&
1266 8 : EQUAL(aosRatTypes[0], "integer") &&
1267 8 : EQUAL(aosRatTypes[1], "integer") &&
1268 8 : EQUAL(aosRatTypes[2], "integer") &&
1269 32 : EQUAL(aosRatTypes[3], "integer") &&
1270 8 : (aosRatTypes.size() == 4 || EQUAL(aosRatTypes[4], "integer")))
1271 : {
1272 8 : bIsCompatibleOfCT = true;
1273 8 : poDS->m_poCT.reset(new GDALColorTable());
1274 24 : for (int i = 0; i < nValues; i++)
1275 : {
1276 16 : const int nIndex = atoi(aosRatValues[i]);
1277 16 : if (nIndex >= 0 && nIndex < 65536)
1278 : {
1279 16 : const int nRed = atoi(aosRatValues[nValues + i]);
1280 16 : const int nGreen = atoi(aosRatValues[2 * nValues + i]);
1281 16 : const int nBlue = atoi(aosRatValues[3 * nValues + i]);
1282 : const int nAlpha =
1283 16 : aosRatTypes.size() == 4
1284 16 : ? 255
1285 8 : : atoi(aosRatValues[4 * nValues + i]);
1286 : const GDALColorEntry oEntry = {
1287 : static_cast<short>(nRed),
1288 : static_cast<short>(nGreen),
1289 : static_cast<short>(nBlue),
1290 16 : static_cast<short>(nAlpha)};
1291 :
1292 16 : poDS->m_poCT->SetColorEntry(nIndex, &oEntry);
1293 : }
1294 : else
1295 : {
1296 0 : bIsCompatibleOfCT = false;
1297 0 : poDS->m_poCT.reset();
1298 0 : break;
1299 : }
1300 : }
1301 : }
1302 :
1303 : // cppcheck-suppress knownConditionTrueFalse
1304 12 : if (!bIsCompatibleOfCT)
1305 : {
1306 4 : poDS->m_poRAT.reset(new GDALDefaultRasterAttributeTable());
1307 44 : for (int i = 0; i < aosRatNames.size(); i++)
1308 : {
1309 120 : poDS->m_poRAT->CreateColumn(
1310 40 : aosRatNames[i],
1311 40 : EQUAL(aosRatTypes[i], "integer") ? GFT_Integer
1312 12 : : EQUAL(aosRatTypes[i], "numeric") ? GFT_Real
1313 : : GFT_String,
1314 40 : EQUAL(aosRatNames[i], "red") ? GFU_Red
1315 68 : : EQUAL(aosRatNames[i], "green") ? GFU_Green
1316 60 : : EQUAL(aosRatNames[i], "blue") ? GFU_Blue
1317 52 : : EQUAL(aosRatNames[i], "alpha") ? GFU_Alpha
1318 44 : : EQUAL(aosRatNames[i], "name") ? GFU_Name
1319 20 : : EQUAL(aosRatNames[i], "pixelcount") ? GFU_PixelCount
1320 40 : : GFU_Generic);
1321 : }
1322 12 : for (int i = 0; i < nValues; i++)
1323 : {
1324 88 : for (int j = 0; j < aosRatTypes.size(); j++)
1325 : {
1326 80 : if (poDS->m_poRAT->GetTypeOfCol(j) == GFT_Integer)
1327 : {
1328 56 : poDS->m_poRAT->SetValue(
1329 56 : i, j, atoi(aosRatValues[j * nValues + i]));
1330 : }
1331 24 : else if (poDS->m_poRAT->GetTypeOfCol(j) == GFT_Real)
1332 : {
1333 16 : poDS->m_poRAT->SetValue(
1334 8 : i, j, CPLAtof(aosRatValues[j * nValues + i]));
1335 : }
1336 : else
1337 : {
1338 16 : poDS->m_poRAT->SetValue(
1339 16 : i, j, aosRatValues[j * nValues + i]);
1340 : }
1341 : }
1342 : }
1343 : }
1344 : }
1345 : }
1346 :
1347 166 : CPLStringList aosMinValues(CSLTokenizeString2(osMinValue, ":", 0));
1348 166 : CPLStringList aosMaxValues(CSLTokenizeString2(osMaxValue, ":", 0));
1349 :
1350 166 : CPLStringList aosLayerNames(CSLTokenizeString2(osLayerName, ":", 0));
1351 235 : for (int i = 1; i <= l_nBands; i++)
1352 : {
1353 : auto poBand = std::make_unique<RRASTERRasterBand>(
1354 152 : poDS.get(), i, fpImage, nBandOffset * (i - 1), nPixelOffset,
1355 152 : nLineOffset, eDT, bNativeOrder);
1356 152 : if (!poBand->IsValid())
1357 : {
1358 0 : return nullptr;
1359 : }
1360 152 : if (!osNoDataValue.empty() && !EQUAL(osNoDataValue, "NA"))
1361 : {
1362 10 : double dfNoDataValue = CPLAtof(osNoDataValue);
1363 10 : if (eDT == GDT_Float32)
1364 0 : dfNoDataValue = CastToFloat(dfNoDataValue);
1365 10 : poBand->m_bHasNoDataValue = true;
1366 10 : poBand->m_dfNoDataValue = dfNoDataValue;
1367 : }
1368 225 : if (i - 1 < static_cast<int>(aosMinValues.size()) &&
1369 73 : i - 1 < static_cast<int>(aosMaxValues.size()))
1370 : {
1371 146 : poBand->SetMinMax(CPLAtof(aosMinValues[i - 1]),
1372 73 : CPLAtof(aosMaxValues[i - 1]));
1373 : }
1374 152 : if (i - 1 < static_cast<int>(aosLayerNames.size()))
1375 : {
1376 156 : const CPLString osName(aosLayerNames[i - 1]);
1377 78 : poBand->GDALRasterBand::SetDescription(osName);
1378 78 : if (EQUAL(osName, "red"))
1379 15 : poBand->SetColorInterpretation(GCI_RedBand);
1380 63 : else if (EQUAL(osName, "green"))
1381 15 : poBand->SetColorInterpretation(GCI_GreenBand);
1382 48 : else if (EQUAL(osName, "blue"))
1383 15 : poBand->SetColorInterpretation(GCI_BlueBand);
1384 33 : else if (EQUAL(osName, "alpha"))
1385 15 : poBand->SetColorInterpretation(GCI_AlphaBand);
1386 : }
1387 152 : poBand->m_poRAT = poDS->m_poRAT;
1388 152 : poBand->m_poCT = poDS->m_poCT;
1389 152 : if (poBand->m_poCT)
1390 8 : poBand->SetColorInterpretation(GCI_PaletteIndex);
1391 152 : poDS->SetBand(i, std::move(poBand));
1392 : }
1393 :
1394 83 : return poDS.release();
1395 : }
1396 :
1397 : /************************************************************************/
1398 : /* Create() */
1399 : /************************************************************************/
1400 :
1401 79 : GDALDataset *RRASTERDataset::Create(const char *pszFilename, int nXSize,
1402 : int nYSize, int nBandsIn,
1403 : GDALDataType eType, char **papszOptions)
1404 :
1405 : {
1406 : // Verify input options.
1407 79 : if (nBandsIn <= 0)
1408 : {
1409 1 : CPLError(CE_Failure, CPLE_NotSupported,
1410 : "RRASTER driver does not support %d bands.", nBandsIn);
1411 1 : return nullptr;
1412 : }
1413 :
1414 78 : if (eType != GDT_Byte && eType != GDT_Int8 && eType != GDT_UInt16 &&
1415 32 : eType != GDT_Int16 && eType != GDT_Int32 && eType != GDT_UInt32 &&
1416 20 : eType != GDT_Float32 && eType != GDT_Float64)
1417 : {
1418 16 : CPLError(CE_Failure, CPLE_NotSupported, "Unsupported data type (%s).",
1419 : GDALGetDataTypeName(eType));
1420 16 : return nullptr;
1421 : }
1422 :
1423 124 : CPLString osGRDExtension(CPLGetExtension(pszFilename));
1424 62 : if (!EQUAL(osGRDExtension, "grd"))
1425 : {
1426 0 : CPLError(CE_Failure, CPLE_NotSupported,
1427 : "RRASTER driver only supports grd extension");
1428 0 : return nullptr;
1429 : }
1430 :
1431 62 : int nPixelOffset = 0;
1432 62 : int nLineOffset = 0;
1433 62 : vsi_l_offset nBandOffset = 0;
1434 : CPLString osBandOrder(
1435 124 : CSLFetchNameValueDef(papszOptions, "INTERLEAVE", "BIL"));
1436 62 : if (!ComputeSpacings(osBandOrder, nXSize, nYSize, nBandsIn, eType,
1437 : nPixelOffset, nLineOffset, nBandOffset))
1438 : {
1439 0 : return nullptr;
1440 : }
1441 :
1442 62 : const std::string osGRIExtension((osGRDExtension[0] == 'g') ? "gri"
1443 186 : : "GRI");
1444 : const std::string osGriFilename(
1445 124 : CPLResetExtension(pszFilename, osGRIExtension.c_str()));
1446 :
1447 : // Try to create the file.
1448 62 : VSILFILE *fpImage = VSIFOpenL(osGriFilename.c_str(), "wb+");
1449 :
1450 62 : if (fpImage == nullptr)
1451 : {
1452 3 : CPLError(CE_Failure, CPLE_OpenFailed,
1453 : "Attempt to create file `%s' failed.", osGriFilename.c_str());
1454 3 : return nullptr;
1455 : }
1456 :
1457 59 : RRASTERDataset *poDS = new RRASTERDataset;
1458 59 : poDS->eAccess = GA_Update;
1459 59 : poDS->m_bHeaderDirty = true;
1460 59 : poDS->m_osGriFilename = osGriFilename;
1461 59 : poDS->nRasterXSize = nXSize;
1462 59 : poDS->nRasterYSize = nYSize;
1463 59 : poDS->m_fpImage = fpImage;
1464 59 : poDS->m_bNativeOrder = true;
1465 59 : poDS->m_osBandOrder = osBandOrder.toupper();
1466 59 : poDS->m_bInitRaster = CPLFetchBool(papszOptions, "@INIT_RASTER", true);
1467 :
1468 59 : const char *pszPixelType = CSLFetchNameValue(papszOptions, "PIXELTYPE");
1469 60 : const bool bByteSigned = (eType == GDT_Byte && pszPixelType &&
1470 1 : EQUAL(pszPixelType, "SIGNEDBYTE"));
1471 59 : if (bByteSigned)
1472 1 : poDS->m_bSignedByte = true;
1473 :
1474 167 : for (int i = 1; i <= nBandsIn; i++)
1475 : {
1476 : RRASTERRasterBand *poBand =
1477 108 : new RRASTERRasterBand(poDS, i, fpImage, nBandOffset * (i - 1),
1478 108 : nPixelOffset, nLineOffset, eType, true);
1479 108 : poDS->SetBand(i, poBand);
1480 : }
1481 :
1482 59 : return poDS;
1483 : }
1484 :
1485 : /************************************************************************/
1486 : /* CreateCopy() */
1487 : /************************************************************************/
1488 :
1489 46 : GDALDataset *RRASTERDataset::CreateCopy(const char *pszFilename,
1490 : GDALDataset *poSrcDS, int bStrict,
1491 : char **papszOptions,
1492 : GDALProgressFunc pfnProgress,
1493 : void *pProgressData)
1494 :
1495 : {
1496 : // Proceed with normal copying using the default createcopy operators.
1497 : GDALDriver *poDriver =
1498 46 : reinterpret_cast<GDALDriver *>(GDALGetDriverByName("RRASTER"));
1499 :
1500 46 : char **papszAdjustedOptions = CSLDuplicate(papszOptions);
1501 : papszAdjustedOptions =
1502 46 : CSLSetNameValue(papszAdjustedOptions, "@INIT_RASTER", "NO");
1503 46 : GDALDataset *poOutDS = poDriver->DefaultCreateCopy(
1504 : pszFilename, poSrcDS, bStrict, papszAdjustedOptions, pfnProgress,
1505 : pProgressData);
1506 46 : CSLDestroy(papszAdjustedOptions);
1507 :
1508 46 : if (poOutDS != nullptr)
1509 36 : poOutDS->FlushCache(false);
1510 :
1511 46 : return poOutDS;
1512 : }
1513 :
1514 : /************************************************************************/
1515 : /* GDALRegister_RRASTER() */
1516 : /************************************************************************/
1517 :
1518 1595 : void GDALRegister_RRASTER()
1519 :
1520 : {
1521 1595 : if (GDALGetDriverByName("RRASTER") != nullptr)
1522 302 : return;
1523 :
1524 1293 : GDALDriver *poDriver = new GDALDriver();
1525 :
1526 1293 : poDriver->SetDescription("RRASTER");
1527 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1528 1293 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "grd");
1529 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "R Raster");
1530 1293 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC,
1531 1293 : "drivers/raster/rraster.html");
1532 1293 : poDriver->SetMetadataItem(
1533 : GDAL_DMD_CREATIONDATATYPES,
1534 1293 : "Byte Int8 Int16 UInt16 Int32 UInt32 Float32 Float64");
1535 :
1536 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1537 :
1538 1293 : poDriver->SetMetadataItem(
1539 : GDAL_DMD_CREATIONOPTIONLIST,
1540 : "<CreationOptionList>"
1541 : " <Option name='PIXELTYPE' type='string' description='(deprecated, "
1542 : "use Int8) "
1543 : "By setting this to SIGNEDBYTE, a new Byte file can be forced to be "
1544 : "written "
1545 : "as signed byte'/>"
1546 : " <Option name='INTERLEAVE' type='string-select' default='BIL'>"
1547 : " <Value>BIP</Value>"
1548 : " <Value>BIL</Value>"
1549 : " <Value>BSQ</Value>"
1550 : " </Option>"
1551 1293 : "</CreationOptionList>");
1552 :
1553 1293 : poDriver->pfnOpen = RRASTERDataset::Open;
1554 1293 : poDriver->pfnIdentify = RRASTERDataset::Identify;
1555 1293 : poDriver->pfnCreate = RRASTERDataset::Create;
1556 1293 : poDriver->pfnCreateCopy = RRASTERDataset::CreateCopy;
1557 :
1558 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
1559 : }
|
