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