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