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