Line data Source code
1 : /******************************************************************************
2 : * Project: SAGA GIS Binary Driver
3 : * Purpose: Implements the SAGA GIS Binary Grid Format.
4 : * Author: Volker Wichmann, wichmann@laserdata.at
5 : * (Based on gsbgdataset.cpp by Kevin Locke and Frank Warmerdam)
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2009, Volker Wichmann <wichmann@laserdata.at>
9 : * Copyright (c) 2009-2011, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_conv.h"
15 :
16 : #include <cassert>
17 : #include <cfloat>
18 : #include <climits>
19 :
20 : #include "gdal_frmts.h"
21 : #include "gdal_pam.h"
22 : #include "ogr_spatialref.h"
23 :
24 : #ifndef INT_MAX
25 : #define INT_MAX 2147483647
26 : #endif /* INT_MAX */
27 :
28 : /* NODATA Values */
29 : // #define SG_NODATA_GDT_Bit 0.0
30 : constexpr GByte SG_NODATA_GDT_Byte = 255;
31 : #define SG_NODATA_GDT_UInt16 65535
32 : #define SG_NODATA_GDT_Int16 -32767
33 : #define SG_NODATA_GDT_UInt32 4294967295U
34 : #define SG_NODATA_GDT_Int32 -2147483647
35 : #define SG_NODATA_GDT_Float32 -99999.0
36 : #define SG_NODATA_GDT_Float64 -99999.0
37 :
38 : /************************************************************************/
39 : /* ==================================================================== */
40 : /* SAGADataset */
41 : /* ==================================================================== */
42 : /************************************************************************/
43 :
44 : class SAGARasterBand;
45 :
46 : class SAGADataset final : public GDALPamDataset
47 : {
48 : friend class SAGARasterBand;
49 :
50 : static CPLErr WriteHeader(CPLString osHDRFilename, GDALDataType eType,
51 : int nXSize, int nYSize, double dfMinX,
52 : double dfMinY, double dfCellsize, double dfNoData,
53 : double dfZFactor, bool bTopToBottom);
54 : VSILFILE *fp;
55 : OGRSpatialReference m_oSRS{};
56 : bool headerDirty = false;
57 :
58 : public:
59 : SAGADataset();
60 : virtual ~SAGADataset();
61 :
62 : static GDALDataset *Open(GDALOpenInfo *);
63 : static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
64 : int nBandsIn, GDALDataType eType,
65 : char **papszParamList);
66 : static GDALDataset *CreateCopy(const char *pszFilename,
67 : GDALDataset *poSrcDS, int bStrict,
68 : char **papszOptions,
69 : GDALProgressFunc pfnProgress,
70 : void *pProgressData);
71 :
72 : const OGRSpatialReference *GetSpatialRef() const override;
73 : CPLErr SetSpatialRef(const OGRSpatialReference *poSRS) override;
74 :
75 : virtual char **GetFileList() override;
76 :
77 : CPLErr GetGeoTransform(double *padfGeoTransform) override;
78 : CPLErr SetGeoTransform(double *padfGeoTransform) override;
79 : };
80 :
81 : /************************************************************************/
82 : /* ==================================================================== */
83 : /* SAGARasterBand */
84 : /* ==================================================================== */
85 : /************************************************************************/
86 :
87 : class SAGARasterBand final : public GDALPamRasterBand
88 : {
89 : friend class SAGADataset;
90 :
91 : double m_Xmin;
92 : double m_Ymin;
93 : double m_Cellsize;
94 : double m_NoData;
95 : int m_ByteOrder;
96 : int m_nBits;
97 :
98 : void SetDataType(GDALDataType eType);
99 : void SwapBuffer(void *pImage) const;
100 :
101 : public:
102 : SAGARasterBand(SAGADataset *, int);
103 :
104 : CPLErr IReadBlock(int, int, void *) override;
105 : CPLErr IWriteBlock(int, int, void *) override;
106 :
107 : double GetNoDataValue(int *pbSuccess = nullptr) override;
108 : CPLErr SetNoDataValue(double dfNoData) override;
109 : };
110 :
111 : /************************************************************************/
112 : /* SAGARasterBand() */
113 : /************************************************************************/
114 :
115 108 : SAGARasterBand::SAGARasterBand(SAGADataset *poDS_, int nBand_)
116 : : m_Xmin(0.0), m_Ymin(0.0), m_Cellsize(0.0), m_NoData(0.0), m_ByteOrder(0),
117 108 : m_nBits(0)
118 : {
119 108 : poDS = poDS_;
120 108 : nBand = nBand_;
121 :
122 108 : eDataType = GDT_Float32;
123 :
124 108 : nBlockXSize = poDS->GetRasterXSize();
125 108 : nBlockYSize = 1;
126 108 : }
127 :
128 : /************************************************************************/
129 : /* SetDataType() */
130 : /************************************************************************/
131 :
132 108 : void SAGARasterBand::SetDataType(GDALDataType eType)
133 :
134 : {
135 108 : eDataType = eType;
136 108 : return;
137 : }
138 :
139 : /************************************************************************/
140 : /* SwapBuffer() */
141 : /************************************************************************/
142 :
143 1127 : void SAGARasterBand::SwapBuffer(void *pImage) const
144 : {
145 :
146 : #ifdef CPL_LSB
147 1127 : const bool bSwap = (m_ByteOrder == 1);
148 : #else
149 : const bool bSwap = (m_ByteOrder == 0);
150 : #endif
151 :
152 1127 : if (bSwap)
153 : {
154 0 : if (m_nBits == 16)
155 : {
156 0 : short *pImage16 = reinterpret_cast<short *>(pImage);
157 0 : for (int iPixel = 0; iPixel < nBlockXSize; iPixel++)
158 : {
159 0 : CPL_SWAP16PTR(pImage16 + iPixel);
160 : }
161 : }
162 0 : else if (m_nBits == 32)
163 : {
164 0 : int *pImage32 = reinterpret_cast<int *>(pImage);
165 0 : for (int iPixel = 0; iPixel < nBlockXSize; iPixel++)
166 : {
167 0 : CPL_SWAP32PTR(pImage32 + iPixel);
168 : }
169 : }
170 0 : else if (m_nBits == 64)
171 : {
172 0 : double *pImage64 = reinterpret_cast<double *>(pImage);
173 0 : for (int iPixel = 0; iPixel < nBlockXSize; iPixel++)
174 : {
175 0 : CPL_SWAP64PTR(pImage64 + iPixel);
176 : }
177 : }
178 : }
179 1127 : }
180 :
181 : /************************************************************************/
182 : /* IReadBlock() */
183 : /************************************************************************/
184 :
185 367 : CPLErr SAGARasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage)
186 :
187 : {
188 367 : if (nBlockYOff < 0 || nBlockYOff > nRasterYSize - 1 || nBlockXOff != 0)
189 0 : return CE_Failure;
190 :
191 367 : SAGADataset *poGDS = static_cast<SAGADataset *>(poDS);
192 367 : vsi_l_offset offset = static_cast<vsi_l_offset>(m_nBits / 8) *
193 367 : nRasterXSize * (nRasterYSize - nBlockYOff - 1);
194 :
195 367 : if (VSIFSeekL(poGDS->fp, offset, SEEK_SET) != 0)
196 : {
197 0 : CPLError(CE_Failure, CPLE_FileIO,
198 : "Unable to seek to beginning of grid row.\n");
199 0 : return CE_Failure;
200 : }
201 367 : if (VSIFReadL(pImage, m_nBits / 8, nBlockXSize, poGDS->fp) !=
202 367 : static_cast<unsigned>(nBlockXSize))
203 : {
204 0 : CPLError(CE_Failure, CPLE_FileIO,
205 : "Unable to read block from grid file.\n");
206 0 : return CE_Failure;
207 : }
208 :
209 367 : SwapBuffer(pImage);
210 :
211 367 : return CE_None;
212 : }
213 :
214 : /************************************************************************/
215 : /* IWriteBlock() */
216 : /************************************************************************/
217 :
218 380 : CPLErr SAGARasterBand::IWriteBlock(int nBlockXOff, int nBlockYOff, void *pImage)
219 :
220 : {
221 380 : if (eAccess == GA_ReadOnly)
222 : {
223 0 : CPLError(CE_Failure, CPLE_NoWriteAccess,
224 : "Unable to write block, dataset opened read only.\n");
225 0 : return CE_Failure;
226 : }
227 :
228 380 : if (nBlockYOff < 0 || nBlockYOff > nRasterYSize - 1 || nBlockXOff != 0)
229 0 : return CE_Failure;
230 :
231 380 : const vsi_l_offset offset = static_cast<vsi_l_offset>(m_nBits / 8) *
232 380 : nRasterXSize * (nRasterYSize - nBlockYOff - 1);
233 380 : SAGADataset *poGDS = static_cast<SAGADataset *>(poDS);
234 380 : assert(poGDS != nullptr);
235 :
236 380 : if (VSIFSeekL(poGDS->fp, offset, SEEK_SET) != 0)
237 : {
238 0 : CPLError(CE_Failure, CPLE_FileIO,
239 : "Unable to seek to beginning of grid row.\n");
240 0 : return CE_Failure;
241 : }
242 :
243 380 : SwapBuffer(pImage);
244 :
245 : const bool bSuccess =
246 380 : (VSIFWriteL(pImage, m_nBits / 8, nBlockXSize, poGDS->fp) ==
247 380 : static_cast<unsigned>(nBlockXSize));
248 :
249 380 : SwapBuffer(pImage);
250 :
251 380 : if (!bSuccess)
252 : {
253 0 : CPLError(CE_Failure, CPLE_FileIO,
254 : "Unable to write block to grid file.\n");
255 0 : return CE_Failure;
256 : }
257 :
258 380 : return CE_None;
259 : }
260 :
261 : /************************************************************************/
262 : /* GetNoDataValue() */
263 : /************************************************************************/
264 :
265 72 : double SAGARasterBand::GetNoDataValue(int *pbSuccess)
266 : {
267 72 : if (pbSuccess)
268 59 : *pbSuccess = TRUE;
269 :
270 72 : return m_NoData;
271 : }
272 :
273 : /************************************************************************/
274 : /* SetNoDataValue() */
275 : /************************************************************************/
276 :
277 2 : CPLErr SAGARasterBand::SetNoDataValue(double dfNoData)
278 : {
279 2 : if (eAccess == GA_ReadOnly)
280 : {
281 1 : CPLError(CE_Failure, CPLE_NoWriteAccess,
282 : "Unable to set no data value, dataset opened read only.\n");
283 1 : return CE_Failure;
284 : }
285 :
286 1 : m_NoData = dfNoData;
287 1 : SAGADataset *poSAGADS = static_cast<SAGADataset *>(poDS);
288 1 : poSAGADS->headerDirty = true;
289 1 : return CE_None;
290 : }
291 :
292 : /************************************************************************/
293 : /* ==================================================================== */
294 : /* SAGADataset */
295 : /* ==================================================================== */
296 : /************************************************************************/
297 :
298 108 : SAGADataset::SAGADataset() : fp(nullptr)
299 : {
300 108 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
301 108 : }
302 :
303 216 : SAGADataset::~SAGADataset()
304 : {
305 108 : if (headerDirty)
306 : {
307 24 : SAGARasterBand *poGRB = static_cast<SAGARasterBand *>(GetRasterBand(1));
308 48 : const CPLString osPath = CPLGetPathSafe(GetDescription());
309 48 : const CPLString osName = CPLGetBasenameSafe(GetDescription());
310 : const CPLString osFilename =
311 24 : CPLFormCIFilenameSafe(osPath, osName, ".sgrd");
312 24 : WriteHeader(osFilename, poGRB->GetRasterDataType(), poGRB->nRasterXSize,
313 : poGRB->nRasterYSize, poGRB->m_Xmin, poGRB->m_Ymin,
314 : poGRB->m_Cellsize, poGRB->m_NoData, 1.0, false);
315 : }
316 108 : FlushCache(true);
317 108 : if (fp != nullptr)
318 108 : VSIFCloseL(fp);
319 216 : }
320 :
321 : /************************************************************************/
322 : /* GetFileList() */
323 : /************************************************************************/
324 :
325 25 : char **SAGADataset::GetFileList()
326 : {
327 50 : const CPLString osPath = CPLGetPathSafe(GetDescription());
328 25 : const CPLString osName = CPLGetBasenameSafe(GetDescription());
329 :
330 : // Main data file, etc.
331 25 : char **papszFileList = GDALPamDataset::GetFileList();
332 :
333 25 : if (!EQUAL(CPLGetExtensionSafe(GetDescription()).c_str(), "sg-grd-z"))
334 : {
335 : // Header file.
336 48 : CPLString osFilename = CPLFormCIFilenameSafe(osPath, osName, ".sgrd");
337 24 : papszFileList = CSLAddString(papszFileList, osFilename);
338 :
339 : // projections file.
340 24 : osFilename = CPLFormCIFilenameSafe(osPath, osName, "prj");
341 : VSIStatBufL sStatBuf;
342 24 : if (VSIStatExL(osFilename, &sStatBuf, VSI_STAT_EXISTS_FLAG) == 0)
343 10 : papszFileList = CSLAddString(papszFileList, osFilename);
344 : }
345 :
346 50 : return papszFileList;
347 : }
348 :
349 : /************************************************************************/
350 : /* GetSpatialRef() */
351 : /************************************************************************/
352 :
353 6 : const OGRSpatialReference *SAGADataset::GetSpatialRef() const
354 :
355 : {
356 6 : if (!m_oSRS.IsEmpty())
357 6 : return &m_oSRS;
358 :
359 0 : return GDALPamDataset::GetSpatialRef();
360 : }
361 :
362 : /************************************************************************/
363 : /* SetSpatialRef() */
364 : /************************************************************************/
365 :
366 23 : CPLErr SAGADataset::SetSpatialRef(const OGRSpatialReference *poSRS)
367 :
368 : {
369 : /* -------------------------------------------------------------------- */
370 : /* Reset coordinate system on the dataset. */
371 : /* -------------------------------------------------------------------- */
372 23 : m_oSRS.Clear();
373 23 : if (poSRS == nullptr)
374 0 : return CE_None;
375 23 : m_oSRS = *poSRS;
376 :
377 : /* -------------------------------------------------------------------- */
378 : /* Convert to ESRI WKT. */
379 : /* -------------------------------------------------------------------- */
380 23 : char *pszESRI_SRS = nullptr;
381 23 : const char *const apszOptions[] = {"FORMAT=WKT1_ESRI", nullptr};
382 23 : m_oSRS.exportToWkt(&pszESRI_SRS, apszOptions);
383 :
384 : /* -------------------------------------------------------------------- */
385 : /* Write to .prj file. */
386 : /* -------------------------------------------------------------------- */
387 : const CPLString osPrjFilename =
388 23 : CPLResetExtensionSafe(GetDescription(), "prj");
389 23 : VSILFILE *l_fp = VSIFOpenL(osPrjFilename.c_str(), "wt");
390 23 : if (l_fp != nullptr)
391 : {
392 23 : VSIFWriteL(pszESRI_SRS, 1, strlen(pszESRI_SRS), l_fp);
393 23 : VSIFWriteL(reinterpret_cast<void *>(const_cast<char *>("\n")), 1, 1,
394 : l_fp);
395 23 : VSIFCloseL(l_fp);
396 : }
397 :
398 23 : CPLFree(pszESRI_SRS);
399 :
400 23 : return CE_None;
401 : }
402 :
403 : /************************************************************************/
404 : /* Open() */
405 : /************************************************************************/
406 :
407 30197 : GDALDataset *SAGADataset::Open(GDALOpenInfo *poOpenInfo)
408 :
409 : {
410 : /* -------------------------------------------------------------------- */
411 : /* We assume the user is pointing to the binary (i.e. .sdat) file or a */
412 : /* compressed raster (.sg-grd-z) file. */
413 : /* -------------------------------------------------------------------- */
414 30197 : const auto &osExtension = poOpenInfo->osExtension;
415 :
416 60128 : if (!EQUAL(osExtension.c_str(), "sdat") &&
417 29937 : !EQUAL(osExtension.c_str(), "sg-grd-z"))
418 : {
419 29922 : return nullptr;
420 : }
421 :
422 527 : CPLString osPath, osFullname, osName, osHDRFilename;
423 :
424 263 : if (EQUAL(osExtension.c_str(), "sg-grd-z") &&
425 2 : !STARTS_WITH(poOpenInfo->pszFilename, "/vsizip"))
426 : {
427 2 : osPath = "/vsizip/{";
428 2 : osPath += poOpenInfo->pszFilename;
429 2 : osPath += "}/";
430 :
431 2 : char **filesinzip = VSIReadDir(osPath);
432 2 : if (filesinzip == nullptr)
433 0 : return nullptr; // empty zip file
434 :
435 2 : CPLString file;
436 4 : for (int iFile = 0; filesinzip[iFile] != nullptr; iFile++)
437 : {
438 4 : if (EQUAL(CPLGetExtensionSafe(filesinzip[iFile]).c_str(), "sdat"))
439 : {
440 2 : file = filesinzip[iFile];
441 2 : break;
442 : }
443 : }
444 :
445 2 : CSLDestroy(filesinzip);
446 :
447 2 : osFullname = CPLFormFilenameSafe(osPath, file, nullptr);
448 2 : osName = CPLGetBasenameSafe(file);
449 2 : osHDRFilename = CPLFormFilenameSafe(
450 4 : osPath, CPLGetBasenameSafe(file).c_str(), "sgrd");
451 : }
452 : else
453 : {
454 259 : osFullname = poOpenInfo->pszFilename;
455 259 : osPath = CPLGetPathSafe(poOpenInfo->pszFilename);
456 259 : osName = CPLGetBasenameSafe(poOpenInfo->pszFilename);
457 259 : osHDRFilename = CPLFormCIFilenameSafe(
458 518 : osPath, CPLGetBasenameSafe(poOpenInfo->pszFilename).c_str(),
459 259 : "sgrd");
460 : }
461 :
462 261 : VSILFILE *fp = VSIFOpenL(osHDRFilename, "r");
463 261 : if (fp == nullptr)
464 : {
465 143 : return nullptr;
466 : }
467 :
468 : /* -------------------------------------------------------------------- */
469 : /* Is this file a SAGA header file? Read a few lines of text */
470 : /* searching for something starting with nrows or ncols. */
471 : /* -------------------------------------------------------------------- */
472 118 : int nRows = -1;
473 118 : int nCols = -1;
474 118 : double dXmin = 0.0;
475 118 : double dYmin = 0.0;
476 118 : double dCellsize = 0.0;
477 118 : double dNoData = 0.0;
478 118 : double dZFactor = 0.0;
479 118 : int nLineCount = 0;
480 118 : char szDataFormat[20] = "DOUBLE";
481 118 : char szByteOrderBig[10] = "FALSE";
482 118 : char szTopToBottom[10] = "FALSE";
483 :
484 118 : const char *pszLine = nullptr;
485 1656 : while ((pszLine = CPLReadLineL(fp)) != nullptr)
486 : {
487 1538 : nLineCount++;
488 :
489 1538 : if (nLineCount > 50 || strlen(pszLine) > 1000)
490 : break;
491 :
492 : char **papszTokens =
493 1538 : CSLTokenizeStringComplex(pszLine, " =", TRUE, FALSE);
494 1538 : if (CSLCount(papszTokens) < 2)
495 : {
496 229 : CSLDestroy(papszTokens);
497 229 : continue;
498 : }
499 :
500 1309 : char **papszHDR = CSLAddString(nullptr, pszLine);
501 :
502 1309 : if (STARTS_WITH_CI(papszTokens[0], "CELLCOUNT_X"))
503 108 : nCols = atoi(papszTokens[1]);
504 1201 : else if (STARTS_WITH_CI(papszTokens[0], "CELLCOUNT_Y"))
505 108 : nRows = atoi(papszTokens[1]);
506 1093 : else if (STARTS_WITH_CI(papszTokens[0], "POSITION_XMIN"))
507 108 : dXmin = CPLAtofM(papszTokens[1]);
508 985 : else if (STARTS_WITH_CI(papszTokens[0], "POSITION_YMIN"))
509 108 : dYmin = CPLAtofM(papszTokens[1]);
510 877 : else if (STARTS_WITH_CI(papszTokens[0], "CELLSIZE"))
511 108 : dCellsize = CPLAtofM(papszTokens[1]);
512 769 : else if (STARTS_WITH_CI(papszTokens[0], "NODATA_VALUE"))
513 108 : dNoData = CPLAtofM(papszTokens[1]);
514 661 : else if (STARTS_WITH_CI(papszTokens[0], "DATAFORMAT"))
515 108 : strncpy(szDataFormat, papszTokens[1], sizeof(szDataFormat) - 1);
516 553 : else if (STARTS_WITH_CI(papszTokens[0], "BYTEORDER_BIG"))
517 109 : strncpy(szByteOrderBig, papszTokens[1], sizeof(szByteOrderBig) - 1);
518 444 : else if (STARTS_WITH_CI(papszTokens[0], "TOPTOBOTTOM"))
519 108 : strncpy(szTopToBottom, papszTokens[1], sizeof(szTopToBottom) - 1);
520 336 : else if (STARTS_WITH_CI(papszTokens[0], "Z_FACTOR"))
521 108 : dZFactor = CPLAtofM(papszTokens[1]);
522 :
523 1309 : CSLDestroy(papszTokens);
524 1309 : CSLDestroy(papszHDR);
525 : }
526 :
527 118 : VSIFCloseL(fp);
528 :
529 : /* -------------------------------------------------------------------- */
530 : /* Did we get the required keywords? If not we return with */
531 : /* this never having been considered to be a match. This isn't */
532 : /* an error! */
533 : /* -------------------------------------------------------------------- */
534 118 : if (nRows == -1 || nCols == -1)
535 : {
536 10 : return nullptr;
537 : }
538 :
539 108 : if (!GDALCheckDatasetDimensions(nCols, nRows))
540 : {
541 0 : return nullptr;
542 : }
543 :
544 108 : if (STARTS_WITH_CI(szTopToBottom, "TRUE"))
545 : {
546 0 : CPLError(CE_Failure, CPLE_AppDefined,
547 : "Currently the SAGA Binary Grid driver does not support\n"
548 : "SAGA grids written TOPTOBOTTOM.\n");
549 0 : return nullptr;
550 : }
551 108 : if (dZFactor != 1.0)
552 : {
553 0 : CPLError(CE_Warning, CPLE_AppDefined,
554 : "Currently the SAGA Binary Grid driver does not support\n"
555 : "ZFACTORs other than 1.\n");
556 : }
557 :
558 : /* -------------------------------------------------------------------- */
559 : /* Create a corresponding GDALDataset. */
560 : /* -------------------------------------------------------------------- */
561 108 : SAGADataset *poDS = new SAGADataset();
562 :
563 108 : poDS->eAccess = poOpenInfo->eAccess;
564 108 : if (poOpenInfo->eAccess == GA_ReadOnly)
565 68 : poDS->fp = VSIFOpenL(osFullname.c_str(), "rb");
566 : else
567 40 : poDS->fp = VSIFOpenL(osFullname.c_str(), "r+b");
568 :
569 108 : if (poDS->fp == nullptr)
570 : {
571 0 : delete poDS;
572 0 : CPLError(CE_Failure, CPLE_OpenFailed,
573 : "VSIFOpenL(%s) failed unexpectedly.", osFullname.c_str());
574 0 : return nullptr;
575 : }
576 :
577 108 : poDS->nRasterXSize = nCols;
578 108 : poDS->nRasterYSize = nRows;
579 :
580 108 : SAGARasterBand *poBand = new SAGARasterBand(poDS, 1);
581 :
582 : /* -------------------------------------------------------------------- */
583 : /* Figure out the byte order. */
584 : /* -------------------------------------------------------------------- */
585 108 : if (STARTS_WITH_CI(szByteOrderBig, "TRUE"))
586 0 : poBand->m_ByteOrder = 1;
587 108 : else if (STARTS_WITH_CI(szByteOrderBig, "FALSE"))
588 108 : poBand->m_ByteOrder = 0;
589 :
590 : /* -------------------------------------------------------------------- */
591 : /* Figure out the data type. */
592 : /* -------------------------------------------------------------------- */
593 108 : if (EQUAL(szDataFormat, "BIT"))
594 : {
595 0 : poBand->SetDataType(GDT_Byte);
596 0 : poBand->m_nBits = 8;
597 : }
598 108 : else if (EQUAL(szDataFormat, "BYTE_UNSIGNED"))
599 : {
600 13 : poBand->SetDataType(GDT_Byte);
601 13 : poBand->m_nBits = 8;
602 : }
603 95 : else if (EQUAL(szDataFormat, "BYTE"))
604 : {
605 0 : poBand->SetDataType(GDT_Byte);
606 0 : poBand->m_nBits = 8;
607 : }
608 95 : else if (EQUAL(szDataFormat, "SHORTINT_UNSIGNED"))
609 : {
610 13 : poBand->SetDataType(GDT_UInt16);
611 13 : poBand->m_nBits = 16;
612 : }
613 82 : else if (EQUAL(szDataFormat, "SHORTINT"))
614 : {
615 13 : poBand->SetDataType(GDT_Int16);
616 13 : poBand->m_nBits = 16;
617 : }
618 69 : else if (EQUAL(szDataFormat, "INTEGER_UNSIGNED"))
619 : {
620 13 : poBand->SetDataType(GDT_UInt32);
621 13 : poBand->m_nBits = 32;
622 : }
623 56 : else if (EQUAL(szDataFormat, "INTEGER"))
624 : {
625 13 : poBand->SetDataType(GDT_Int32);
626 13 : poBand->m_nBits = 32;
627 : }
628 43 : else if (EQUAL(szDataFormat, "FLOAT"))
629 : {
630 29 : poBand->SetDataType(GDT_Float32);
631 29 : poBand->m_nBits = 32;
632 : }
633 14 : else if (EQUAL(szDataFormat, "DOUBLE"))
634 : {
635 14 : poBand->SetDataType(GDT_Float64);
636 14 : poBand->m_nBits = 64;
637 : }
638 : else
639 : {
640 0 : CPLError(CE_Failure, CPLE_NotSupported,
641 : "SAGA driver does not support the dataformat %s.",
642 : szDataFormat);
643 0 : delete poBand;
644 0 : delete poDS;
645 0 : return nullptr;
646 : }
647 :
648 : /* -------------------------------------------------------------------- */
649 : /* Save band information */
650 : /* -------------------------------------------------------------------- */
651 108 : poBand->m_Xmin = dXmin;
652 108 : poBand->m_Ymin = dYmin;
653 108 : poBand->m_NoData = dNoData;
654 108 : poBand->m_Cellsize = dCellsize;
655 :
656 108 : poDS->SetBand(1, poBand);
657 :
658 : /* -------------------------------------------------------------------- */
659 : /* Initialize any PAM information. */
660 : /* -------------------------------------------------------------------- */
661 108 : poDS->SetDescription(poOpenInfo->pszFilename);
662 108 : poDS->TryLoadXML();
663 :
664 : /* -------------------------------------------------------------------- */
665 : /* Check for a .prj file. */
666 : /* -------------------------------------------------------------------- */
667 : const std::string osPrjFilename =
668 108 : CPLFormCIFilenameSafe(osPath, osName, "prj");
669 :
670 108 : fp = VSIFOpenL(osPrjFilename.c_str(), "r");
671 :
672 108 : if (fp != nullptr)
673 : {
674 32 : VSIFCloseL(fp);
675 :
676 32 : char **papszLines = CSLLoad(osPrjFilename.c_str());
677 :
678 32 : poDS->m_oSRS.importFromESRI(papszLines);
679 :
680 32 : CSLDestroy(papszLines);
681 : }
682 :
683 : /* -------------------------------------------------------------------- */
684 : /* Check for external overviews. */
685 : /* -------------------------------------------------------------------- */
686 216 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename,
687 108 : poOpenInfo->GetSiblingFiles());
688 :
689 108 : return poDS;
690 : }
691 :
692 : /************************************************************************/
693 : /* GetGeoTransform() */
694 : /************************************************************************/
695 :
696 11 : CPLErr SAGADataset::GetGeoTransform(double *padfGeoTransform)
697 : {
698 11 : if (padfGeoTransform == nullptr)
699 0 : return CE_Failure;
700 :
701 11 : SAGARasterBand *poGRB = static_cast<SAGARasterBand *>(GetRasterBand(1));
702 :
703 11 : if (poGRB == nullptr)
704 : {
705 0 : padfGeoTransform[0] = 0;
706 0 : padfGeoTransform[1] = 1;
707 0 : padfGeoTransform[2] = 0;
708 0 : padfGeoTransform[3] = 0;
709 0 : padfGeoTransform[4] = 0;
710 0 : padfGeoTransform[5] = 1;
711 0 : return CE_Failure;
712 : }
713 :
714 : /* check if we have a PAM GeoTransform stored */
715 11 : CPLPushErrorHandler(CPLQuietErrorHandler);
716 11 : CPLErr eErr = GDALPamDataset::GetGeoTransform(padfGeoTransform);
717 11 : CPLPopErrorHandler();
718 :
719 11 : if (eErr == CE_None)
720 0 : return CE_None;
721 :
722 11 : padfGeoTransform[1] = poGRB->m_Cellsize;
723 11 : padfGeoTransform[5] = poGRB->m_Cellsize * -1.0;
724 11 : padfGeoTransform[0] = poGRB->m_Xmin - poGRB->m_Cellsize / 2;
725 11 : padfGeoTransform[3] = poGRB->m_Ymin +
726 11 : (nRasterYSize - 1) * poGRB->m_Cellsize +
727 11 : poGRB->m_Cellsize / 2;
728 :
729 : /* tilt/rotation is not supported by SAGA grids */
730 11 : padfGeoTransform[4] = 0.0;
731 11 : padfGeoTransform[2] = 0.0;
732 :
733 11 : return CE_None;
734 : }
735 :
736 : /************************************************************************/
737 : /* SetGeoTransform() */
738 : /************************************************************************/
739 :
740 23 : CPLErr SAGADataset::SetGeoTransform(double *padfGeoTransform)
741 : {
742 :
743 23 : if (eAccess == GA_ReadOnly)
744 : {
745 0 : CPLError(CE_Failure, CPLE_NoWriteAccess,
746 : "Unable to set GeoTransform, dataset opened read only.\n");
747 0 : return CE_Failure;
748 : }
749 :
750 23 : SAGARasterBand *poGRB = static_cast<SAGARasterBand *>(GetRasterBand(1));
751 :
752 23 : if (poGRB == nullptr || padfGeoTransform == nullptr)
753 0 : return CE_Failure;
754 :
755 23 : if (padfGeoTransform[1] != padfGeoTransform[5] * -1.0)
756 : {
757 0 : CPLError(CE_Failure, CPLE_NotSupported,
758 : "Unable to set GeoTransform, SAGA binary grids only support "
759 : "the same cellsize in x-y.\n");
760 0 : return CE_Failure;
761 : }
762 :
763 23 : const double dfMinX = padfGeoTransform[0] + padfGeoTransform[1] / 2;
764 23 : const double dfMinY =
765 23 : padfGeoTransform[5] * (nRasterYSize - 0.5) + padfGeoTransform[3];
766 :
767 23 : poGRB->m_Xmin = dfMinX;
768 23 : poGRB->m_Ymin = dfMinY;
769 23 : poGRB->m_Cellsize = padfGeoTransform[1];
770 23 : headerDirty = true;
771 :
772 23 : return CE_None;
773 : }
774 :
775 : /************************************************************************/
776 : /* WriteHeader() */
777 : /************************************************************************/
778 :
779 73 : CPLErr SAGADataset::WriteHeader(CPLString osHDRFilename, GDALDataType eType,
780 : int nXSize, int nYSize, double dfMinX,
781 : double dfMinY, double dfCellsize,
782 : double dfNoData, double dfZFactor,
783 : bool bTopToBottom)
784 :
785 : {
786 73 : VSILFILE *fp = VSIFOpenL(osHDRFilename, "wt");
787 :
788 73 : if (fp == nullptr)
789 : {
790 0 : CPLError(CE_Failure, CPLE_OpenFailed, "Failed to write .sgrd file %s.",
791 : osHDRFilename.c_str());
792 0 : return CE_Failure;
793 : }
794 :
795 73 : VSIFPrintfL(fp, "NAME\t= %s\n", CPLGetBasenameSafe(osHDRFilename).c_str());
796 73 : VSIFPrintfL(fp, "DESCRIPTION\t=\n");
797 73 : VSIFPrintfL(fp, "UNIT\t=\n");
798 73 : VSIFPrintfL(fp, "DATAFILE_OFFSET\t= 0\n");
799 :
800 73 : if (eType == GDT_Int32)
801 8 : VSIFPrintfL(fp, "DATAFORMAT\t= INTEGER\n");
802 65 : else if (eType == GDT_UInt32)
803 8 : VSIFPrintfL(fp, "DATAFORMAT\t= INTEGER_UNSIGNED\n");
804 57 : else if (eType == GDT_Int16)
805 8 : VSIFPrintfL(fp, "DATAFORMAT\t= SHORTINT\n");
806 49 : else if (eType == GDT_UInt16)
807 8 : VSIFPrintfL(fp, "DATAFORMAT\t= SHORTINT_UNSIGNED\n");
808 41 : else if (eType == GDT_Byte)
809 18 : VSIFPrintfL(fp, "DATAFORMAT\t= BYTE_UNSIGNED\n");
810 23 : else if (eType == GDT_Float32)
811 13 : VSIFPrintfL(fp, "DATAFORMAT\t= FLOAT\n");
812 : else // if( eType == GDT_Float64 )
813 10 : VSIFPrintfL(fp, "DATAFORMAT\t= DOUBLE\n");
814 : #ifdef CPL_LSB
815 73 : VSIFPrintfL(fp, "BYTEORDER_BIG\t= FALSE\n");
816 : #else
817 : VSIFPrintfL(fp, "BYTEORDER_BIG\t= TRUE\n");
818 : #endif
819 :
820 73 : VSIFPrintfL(fp, "POSITION_XMIN\t= %.10f\n", dfMinX);
821 73 : VSIFPrintfL(fp, "POSITION_YMIN\t= %.10f\n", dfMinY);
822 73 : VSIFPrintfL(fp, "CELLCOUNT_X\t= %d\n", nXSize);
823 73 : VSIFPrintfL(fp, "CELLCOUNT_Y\t= %d\n", nYSize);
824 73 : VSIFPrintfL(fp, "CELLSIZE\t= %.10f\n", dfCellsize);
825 73 : VSIFPrintfL(fp, "Z_FACTOR\t= %f\n", dfZFactor);
826 73 : VSIFPrintfL(fp, "NODATA_VALUE\t= %f\n", dfNoData);
827 73 : if (bTopToBottom)
828 0 : VSIFPrintfL(fp, "TOPTOBOTTOM\t= TRUE\n");
829 : else
830 73 : VSIFPrintfL(fp, "TOPTOBOTTOM\t= FALSE\n");
831 :
832 73 : VSIFCloseL(fp);
833 :
834 73 : return CE_None;
835 : }
836 :
837 : /************************************************************************/
838 : /* Create() */
839 : /************************************************************************/
840 :
841 81 : GDALDataset *SAGADataset::Create(const char *pszFilename, int nXSize,
842 : int nYSize, int nBandsIn, GDALDataType eType,
843 : char **papszParamList)
844 :
845 : {
846 81 : if (nXSize <= 0 || nYSize <= 0)
847 : {
848 0 : CPLError(CE_Failure, CPLE_IllegalArg,
849 : "Unable to create grid, both X and Y size must be "
850 : "non-negative.\n");
851 :
852 0 : return nullptr;
853 : }
854 :
855 81 : if (nBandsIn != 1)
856 : {
857 18 : CPLError(CE_Failure, CPLE_IllegalArg,
858 : "SAGA Binary Grid only supports 1 band");
859 18 : return nullptr;
860 : }
861 :
862 63 : if (eType != GDT_Byte && eType != GDT_UInt16 && eType != GDT_Int16 &&
863 30 : eType != GDT_UInt32 && eType != GDT_Int32 && eType != GDT_Float32 &&
864 : eType != GDT_Float64)
865 : {
866 11 : CPLError(CE_Failure, CPLE_AppDefined,
867 : "SAGA Binary Grid only supports Byte, UInt16, Int16, "
868 : "UInt32, Int32, Float32 and Float64 datatypes. Unable to "
869 : "create with type %s.\n",
870 : GDALGetDataTypeName(eType));
871 :
872 11 : return nullptr;
873 : }
874 :
875 52 : VSILFILE *fp = VSIFOpenL(pszFilename, "w+b");
876 :
877 52 : if (fp == nullptr)
878 : {
879 3 : CPLError(CE_Failure, CPLE_OpenFailed,
880 : "Attempt to create file '%s' failed.\n", pszFilename);
881 3 : return nullptr;
882 : }
883 :
884 49 : double dfNoDataVal = 0.0;
885 :
886 : const char *pszNoDataValue =
887 49 : CSLFetchNameValue(papszParamList, "NODATA_VALUE");
888 49 : if (pszNoDataValue)
889 : {
890 2 : dfNoDataVal = CPLAtofM(pszNoDataValue);
891 : }
892 : else
893 : {
894 47 : switch (eType) /* GDT_Byte, GDT_UInt16, GDT_Int16, GDT_UInt32 */
895 : { /* GDT_Int32, GDT_Float32, GDT_Float64 */
896 15 : case (GDT_Byte):
897 : {
898 15 : dfNoDataVal = SG_NODATA_GDT_Byte;
899 15 : break;
900 : }
901 5 : case (GDT_UInt16):
902 : {
903 5 : dfNoDataVal = SG_NODATA_GDT_UInt16;
904 5 : break;
905 : }
906 5 : case (GDT_Int16):
907 : {
908 5 : dfNoDataVal = SG_NODATA_GDT_Int16;
909 5 : break;
910 : }
911 5 : case (GDT_UInt32):
912 : {
913 5 : dfNoDataVal = SG_NODATA_GDT_UInt32;
914 5 : break;
915 : }
916 5 : case (GDT_Int32):
917 : {
918 5 : dfNoDataVal = SG_NODATA_GDT_Int32;
919 5 : break;
920 : }
921 6 : default:
922 : case (GDT_Float32):
923 : {
924 6 : dfNoDataVal = SG_NODATA_GDT_Float32;
925 6 : break;
926 : }
927 6 : case (GDT_Float64):
928 : {
929 6 : dfNoDataVal = SG_NODATA_GDT_Float64;
930 6 : break;
931 : }
932 : }
933 : }
934 :
935 : double dfNoDataForAlignment;
936 49 : void *abyNoData = &dfNoDataForAlignment;
937 49 : GDALCopyWords(&dfNoDataVal, GDT_Float64, 0, abyNoData, eType, 0, 1);
938 :
939 98 : const CPLString osHdrFilename = CPLResetExtensionSafe(pszFilename, "sgrd");
940 49 : CPLErr eErr = WriteHeader(osHdrFilename, eType, nXSize, nYSize, 0.0, 0.0,
941 : 1.0, dfNoDataVal, 1.0, false);
942 :
943 49 : if (eErr != CE_None)
944 : {
945 0 : VSIFCloseL(fp);
946 0 : return nullptr;
947 : }
948 :
949 49 : if (CPLFetchBool(papszParamList, "FILL_NODATA", true))
950 : {
951 23 : const int nDataTypeSize = GDALGetDataTypeSize(eType) / 8;
952 : GByte *pabyNoDataBuf =
953 23 : reinterpret_cast<GByte *>(VSIMalloc2(nDataTypeSize, nXSize));
954 23 : if (pabyNoDataBuf == nullptr)
955 : {
956 0 : VSIFCloseL(fp);
957 0 : return nullptr;
958 : }
959 :
960 889 : for (int iCol = 0; iCol < nXSize; iCol++)
961 : {
962 866 : memcpy(pabyNoDataBuf + iCol * nDataTypeSize, abyNoData,
963 : nDataTypeSize);
964 : }
965 :
966 889 : for (int iRow = 0; iRow < nYSize; iRow++)
967 : {
968 866 : if (VSIFWriteL(pabyNoDataBuf, nDataTypeSize, nXSize, fp) !=
969 866 : static_cast<unsigned>(nXSize))
970 : {
971 0 : VSIFCloseL(fp);
972 0 : VSIFree(pabyNoDataBuf);
973 0 : CPLError(CE_Failure, CPLE_FileIO,
974 : "Unable to write grid cell. Disk full?\n");
975 0 : return nullptr;
976 : }
977 : }
978 :
979 23 : VSIFree(pabyNoDataBuf);
980 : }
981 :
982 49 : VSIFCloseL(fp);
983 :
984 49 : return GDALDataset::FromHandle(GDALOpen(pszFilename, GA_Update));
985 : }
986 :
987 : /************************************************************************/
988 : /* CreateCopy() */
989 : /************************************************************************/
990 :
991 38 : GDALDataset *SAGADataset::CreateCopy(const char *pszFilename,
992 : GDALDataset *poSrcDS, int bStrict,
993 : CPL_UNUSED char **papszOptions,
994 : GDALProgressFunc pfnProgress,
995 : void *pProgressData)
996 : {
997 38 : if (pfnProgress == nullptr)
998 0 : pfnProgress = GDALDummyProgress;
999 :
1000 38 : int nBands = poSrcDS->GetRasterCount();
1001 38 : if (nBands == 0)
1002 : {
1003 1 : CPLError(
1004 : CE_Failure, CPLE_NotSupported,
1005 : "SAGA driver does not support source dataset with zero band.\n");
1006 1 : return nullptr;
1007 : }
1008 37 : else if (nBands > 1)
1009 : {
1010 4 : if (bStrict)
1011 : {
1012 4 : CPLError(CE_Failure, CPLE_NotSupported,
1013 : "Unable to create copy, SAGA Binary Grid "
1014 : "format only supports one raster band.\n");
1015 4 : return nullptr;
1016 : }
1017 : else
1018 0 : CPLError(CE_Warning, CPLE_NotSupported,
1019 : "SAGA Binary Grid format only supports one "
1020 : "raster band, first band will be copied.\n");
1021 : }
1022 :
1023 33 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(1);
1024 :
1025 33 : char **papszCreateOptions = CSLSetNameValue(nullptr, "FILL_NODATA", "NO");
1026 :
1027 33 : int bHasNoDataValue = FALSE;
1028 33 : const double dfNoDataValue = poSrcBand->GetNoDataValue(&bHasNoDataValue);
1029 33 : if (bHasNoDataValue)
1030 : papszCreateOptions =
1031 2 : CSLSetNameValue(papszCreateOptions, "NODATA_VALUE",
1032 : CPLSPrintf("%.16g", dfNoDataValue));
1033 :
1034 : GDALDataset *poDstDS =
1035 33 : Create(pszFilename, poSrcBand->GetXSize(), poSrcBand->GetYSize(), 1,
1036 : poSrcBand->GetRasterDataType(), papszCreateOptions);
1037 33 : CSLDestroy(papszCreateOptions);
1038 :
1039 33 : if (poDstDS == nullptr)
1040 17 : return nullptr;
1041 :
1042 : /* -------------------------------------------------------------------- */
1043 : /* Copy band data. */
1044 : /* -------------------------------------------------------------------- */
1045 :
1046 : CPLErr eErr =
1047 16 : GDALDatasetCopyWholeRaster((GDALDatasetH)poSrcDS, (GDALDatasetH)poDstDS,
1048 : nullptr, pfnProgress, pProgressData);
1049 :
1050 16 : if (eErr == CE_Failure)
1051 : {
1052 0 : delete poDstDS;
1053 0 : return nullptr;
1054 : }
1055 :
1056 : double adfGeoTransform[6];
1057 :
1058 16 : poSrcDS->GetGeoTransform(adfGeoTransform);
1059 16 : poDstDS->SetGeoTransform(adfGeoTransform);
1060 :
1061 16 : poDstDS->SetProjection(poSrcDS->GetProjectionRef());
1062 :
1063 16 : return poDstDS;
1064 : }
1065 :
1066 : /************************************************************************/
1067 : /* GDALRegister_SAGA() */
1068 : /************************************************************************/
1069 :
1070 1682 : void GDALRegister_SAGA()
1071 :
1072 : {
1073 1682 : if (GDALGetDriverByName("SAGA") != nullptr)
1074 301 : return;
1075 :
1076 1381 : GDALDriver *poDriver = new GDALDriver();
1077 :
1078 1381 : poDriver->SetDescription("SAGA");
1079 1381 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1080 1381 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
1081 1381 : "SAGA GIS Binary Grid (.sdat, .sg-grd-z)");
1082 1381 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/sdat.html");
1083 1381 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSIONS, "sdat sg-grd-z");
1084 1381 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
1085 : "Byte Int16 "
1086 1381 : "UInt16 Int32 UInt32 Float32 Float64");
1087 :
1088 1381 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1089 :
1090 1381 : poDriver->pfnOpen = SAGADataset::Open;
1091 1381 : poDriver->pfnCreate = SAGADataset::Create;
1092 1381 : poDriver->pfnCreateCopy = SAGADataset::CreateCopy;
1093 :
1094 1381 : GetGDALDriverManager()->RegisterDriver(poDriver);
1095 : }
|
