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(double *padfGeoTransform) override;
79 : CPLErr SetGeoTransform(double *padfGeoTransform) 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 72 : double SAGARasterBand::GetNoDataValue(int *pbSuccess)
267 : {
268 72 : if (pbSuccess)
269 59 : *pbSuccess = TRUE;
270 :
271 72 : 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 32419 : 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 32419 : const auto &osExtension = poOpenInfo->osExtension;
416 :
417 64576 : if (!EQUAL(osExtension.c_str(), "sdat") &&
418 32161 : !EQUAL(osExtension.c_str(), "sg-grd-z"))
419 : {
420 32154 : return nullptr;
421 : }
422 :
423 521 : 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(double *padfGeoTransform)
698 : {
699 11 : if (padfGeoTransform == nullptr)
700 0 : return CE_Failure;
701 :
702 11 : SAGARasterBand *poGRB = static_cast<SAGARasterBand *>(GetRasterBand(1));
703 :
704 11 : if (poGRB == nullptr)
705 : {
706 0 : padfGeoTransform[0] = 0;
707 0 : padfGeoTransform[1] = 1;
708 0 : padfGeoTransform[2] = 0;
709 0 : padfGeoTransform[3] = 0;
710 0 : padfGeoTransform[4] = 0;
711 0 : padfGeoTransform[5] = 1;
712 0 : return CE_Failure;
713 : }
714 :
715 : /* check if we have a PAM GeoTransform stored */
716 11 : CPLPushErrorHandler(CPLQuietErrorHandler);
717 11 : CPLErr eErr = GDALPamDataset::GetGeoTransform(padfGeoTransform);
718 11 : CPLPopErrorHandler();
719 :
720 11 : if (eErr == CE_None)
721 0 : return CE_None;
722 :
723 11 : padfGeoTransform[1] = poGRB->m_Cellsize;
724 11 : padfGeoTransform[5] = poGRB->m_Cellsize * -1.0;
725 11 : padfGeoTransform[0] = poGRB->m_Xmin - poGRB->m_Cellsize / 2;
726 11 : padfGeoTransform[3] = poGRB->m_Ymin +
727 11 : (nRasterYSize - 1) * poGRB->m_Cellsize +
728 11 : poGRB->m_Cellsize / 2;
729 :
730 : /* tilt/rotation is not supported by SAGA grids */
731 11 : padfGeoTransform[4] = 0.0;
732 11 : padfGeoTransform[2] = 0.0;
733 :
734 11 : return CE_None;
735 : }
736 :
737 : /************************************************************************/
738 : /* SetGeoTransform() */
739 : /************************************************************************/
740 :
741 23 : CPLErr SAGADataset::SetGeoTransform(double *padfGeoTransform)
742 : {
743 :
744 23 : if (eAccess == GA_ReadOnly)
745 : {
746 0 : CPLError(CE_Failure, CPLE_NoWriteAccess,
747 : "Unable to set GeoTransform, dataset opened read only.\n");
748 0 : return CE_Failure;
749 : }
750 :
751 23 : SAGARasterBand *poGRB = static_cast<SAGARasterBand *>(GetRasterBand(1));
752 :
753 23 : if (poGRB == nullptr || padfGeoTransform == nullptr)
754 0 : return CE_Failure;
755 :
756 23 : if (padfGeoTransform[1] != padfGeoTransform[5] * -1.0)
757 : {
758 0 : CPLError(CE_Failure, CPLE_NotSupported,
759 : "Unable to set GeoTransform, SAGA binary grids only support "
760 : "the same cellsize in x-y.\n");
761 0 : return CE_Failure;
762 : }
763 :
764 23 : const double dfMinX = padfGeoTransform[0] + padfGeoTransform[1] / 2;
765 23 : const double dfMinY =
766 23 : padfGeoTransform[5] * (nRasterYSize - 0.5) + padfGeoTransform[3];
767 :
768 23 : poGRB->m_Xmin = dfMinX;
769 23 : poGRB->m_Ymin = dfMinY;
770 23 : poGRB->m_Cellsize = padfGeoTransform[1];
771 23 : headerDirty = true;
772 :
773 23 : return CE_None;
774 : }
775 :
776 : /************************************************************************/
777 : /* WriteHeader() */
778 : /************************************************************************/
779 :
780 73 : CPLErr SAGADataset::WriteHeader(const CPLString &osHDRFilename,
781 : GDALDataType eType, int nXSize, int nYSize,
782 : double dfMinX, double dfMinY, double dfCellsize,
783 : double dfNoData, double dfZFactor,
784 : bool bTopToBottom)
785 :
786 : {
787 73 : VSILFILE *fp = VSIFOpenL(osHDRFilename, "wt");
788 :
789 73 : if (fp == nullptr)
790 : {
791 0 : CPLError(CE_Failure, CPLE_OpenFailed, "Failed to write .sgrd file %s.",
792 : osHDRFilename.c_str());
793 0 : return CE_Failure;
794 : }
795 :
796 73 : VSIFPrintfL(fp, "NAME\t= %s\n", CPLGetBasenameSafe(osHDRFilename).c_str());
797 73 : VSIFPrintfL(fp, "DESCRIPTION\t=\n");
798 73 : VSIFPrintfL(fp, "UNIT\t=\n");
799 73 : VSIFPrintfL(fp, "DATAFILE_OFFSET\t= 0\n");
800 :
801 73 : if (eType == GDT_Int32)
802 8 : VSIFPrintfL(fp, "DATAFORMAT\t= INTEGER\n");
803 65 : else if (eType == GDT_UInt32)
804 8 : VSIFPrintfL(fp, "DATAFORMAT\t= INTEGER_UNSIGNED\n");
805 57 : else if (eType == GDT_Int16)
806 8 : VSIFPrintfL(fp, "DATAFORMAT\t= SHORTINT\n");
807 49 : else if (eType == GDT_UInt16)
808 8 : VSIFPrintfL(fp, "DATAFORMAT\t= SHORTINT_UNSIGNED\n");
809 41 : else if (eType == GDT_Byte)
810 18 : VSIFPrintfL(fp, "DATAFORMAT\t= BYTE_UNSIGNED\n");
811 23 : else if (eType == GDT_Float32)
812 13 : VSIFPrintfL(fp, "DATAFORMAT\t= FLOAT\n");
813 : else // if( eType == GDT_Float64 )
814 10 : VSIFPrintfL(fp, "DATAFORMAT\t= DOUBLE\n");
815 : #ifdef CPL_LSB
816 73 : VSIFPrintfL(fp, "BYTEORDER_BIG\t= FALSE\n");
817 : #else
818 : VSIFPrintfL(fp, "BYTEORDER_BIG\t= TRUE\n");
819 : #endif
820 :
821 73 : VSIFPrintfL(fp, "POSITION_XMIN\t= %.10f\n", dfMinX);
822 73 : VSIFPrintfL(fp, "POSITION_YMIN\t= %.10f\n", dfMinY);
823 73 : VSIFPrintfL(fp, "CELLCOUNT_X\t= %d\n", nXSize);
824 73 : VSIFPrintfL(fp, "CELLCOUNT_Y\t= %d\n", nYSize);
825 73 : VSIFPrintfL(fp, "CELLSIZE\t= %.10f\n", dfCellsize);
826 73 : VSIFPrintfL(fp, "Z_FACTOR\t= %f\n", dfZFactor);
827 73 : VSIFPrintfL(fp, "NODATA_VALUE\t= %f\n", dfNoData);
828 73 : if (bTopToBottom)
829 0 : VSIFPrintfL(fp, "TOPTOBOTTOM\t= TRUE\n");
830 : else
831 73 : VSIFPrintfL(fp, "TOPTOBOTTOM\t= FALSE\n");
832 :
833 73 : VSIFCloseL(fp);
834 :
835 73 : return CE_None;
836 : }
837 :
838 : /************************************************************************/
839 : /* Create() */
840 : /************************************************************************/
841 :
842 81 : GDALDataset *SAGADataset::Create(const char *pszFilename, int nXSize,
843 : int nYSize, int nBandsIn, GDALDataType eType,
844 : char **papszParamList)
845 :
846 : {
847 81 : if (nXSize <= 0 || nYSize <= 0)
848 : {
849 0 : CPLError(CE_Failure, CPLE_IllegalArg,
850 : "Unable to create grid, both X and Y size must be "
851 : "non-negative.\n");
852 :
853 0 : return nullptr;
854 : }
855 :
856 81 : if (nBandsIn != 1)
857 : {
858 18 : CPLError(CE_Failure, CPLE_IllegalArg,
859 : "SAGA Binary Grid only supports 1 band");
860 18 : return nullptr;
861 : }
862 :
863 63 : if (eType != GDT_Byte && eType != GDT_UInt16 && eType != GDT_Int16 &&
864 30 : eType != GDT_UInt32 && eType != GDT_Int32 && eType != GDT_Float32 &&
865 : eType != GDT_Float64)
866 : {
867 11 : CPLError(CE_Failure, CPLE_AppDefined,
868 : "SAGA Binary Grid only supports Byte, UInt16, Int16, "
869 : "UInt32, Int32, Float32 and Float64 datatypes. Unable to "
870 : "create with type %s.\n",
871 : GDALGetDataTypeName(eType));
872 :
873 11 : return nullptr;
874 : }
875 :
876 52 : VSILFILE *fp = VSIFOpenL(pszFilename, "w+b");
877 :
878 52 : if (fp == nullptr)
879 : {
880 3 : CPLError(CE_Failure, CPLE_OpenFailed,
881 : "Attempt to create file '%s' failed.\n", pszFilename);
882 3 : return nullptr;
883 : }
884 :
885 49 : double dfNoDataVal = 0.0;
886 :
887 : const char *pszNoDataValue =
888 49 : CSLFetchNameValue(papszParamList, "NODATA_VALUE");
889 49 : if (pszNoDataValue)
890 : {
891 2 : dfNoDataVal = CPLAtofM(pszNoDataValue);
892 : }
893 : else
894 : {
895 47 : switch (eType) /* GDT_Byte, GDT_UInt16, GDT_Int16, GDT_UInt32 */
896 : { /* GDT_Int32, GDT_Float32, GDT_Float64 */
897 15 : case (GDT_Byte):
898 : {
899 15 : dfNoDataVal = SG_NODATA_GDT_Byte;
900 15 : break;
901 : }
902 5 : case (GDT_UInt16):
903 : {
904 5 : dfNoDataVal = SG_NODATA_GDT_UInt16;
905 5 : break;
906 : }
907 5 : case (GDT_Int16):
908 : {
909 5 : dfNoDataVal = SG_NODATA_GDT_Int16;
910 5 : break;
911 : }
912 5 : case (GDT_UInt32):
913 : {
914 5 : dfNoDataVal = SG_NODATA_GDT_UInt32;
915 5 : break;
916 : }
917 5 : case (GDT_Int32):
918 : {
919 5 : dfNoDataVal = SG_NODATA_GDT_Int32;
920 5 : break;
921 : }
922 6 : default:
923 : case (GDT_Float32):
924 : {
925 6 : dfNoDataVal = SG_NODATA_GDT_Float32;
926 6 : break;
927 : }
928 6 : case (GDT_Float64):
929 : {
930 6 : dfNoDataVal = SG_NODATA_GDT_Float64;
931 6 : break;
932 : }
933 : }
934 : }
935 :
936 : double dfNoDataForAlignment;
937 49 : void *abyNoData = &dfNoDataForAlignment;
938 49 : GDALCopyWords(&dfNoDataVal, GDT_Float64, 0, abyNoData, eType, 0, 1);
939 :
940 98 : const CPLString osHdrFilename = CPLResetExtensionSafe(pszFilename, "sgrd");
941 49 : CPLErr eErr = WriteHeader(osHdrFilename, eType, nXSize, nYSize, 0.0, 0.0,
942 : 1.0, dfNoDataVal, 1.0, false);
943 :
944 49 : if (eErr != CE_None)
945 : {
946 0 : VSIFCloseL(fp);
947 0 : return nullptr;
948 : }
949 :
950 49 : if (CPLFetchBool(papszParamList, "FILL_NODATA", true))
951 : {
952 23 : const int nDataTypeSize = GDALGetDataTypeSize(eType) / 8;
953 : GByte *pabyNoDataBuf =
954 23 : reinterpret_cast<GByte *>(VSIMalloc2(nDataTypeSize, nXSize));
955 23 : if (pabyNoDataBuf == nullptr)
956 : {
957 0 : VSIFCloseL(fp);
958 0 : return nullptr;
959 : }
960 :
961 889 : for (int iCol = 0; iCol < nXSize; iCol++)
962 : {
963 866 : memcpy(pabyNoDataBuf + iCol * nDataTypeSize, abyNoData,
964 : nDataTypeSize);
965 : }
966 :
967 889 : for (int iRow = 0; iRow < nYSize; iRow++)
968 : {
969 866 : if (VSIFWriteL(pabyNoDataBuf, nDataTypeSize, nXSize, fp) !=
970 866 : static_cast<unsigned>(nXSize))
971 : {
972 0 : VSIFCloseL(fp);
973 0 : VSIFree(pabyNoDataBuf);
974 0 : CPLError(CE_Failure, CPLE_FileIO,
975 : "Unable to write grid cell. Disk full?\n");
976 0 : return nullptr;
977 : }
978 : }
979 :
980 23 : VSIFree(pabyNoDataBuf);
981 : }
982 :
983 49 : VSIFCloseL(fp);
984 :
985 49 : return GDALDataset::FromHandle(GDALOpen(pszFilename, GA_Update));
986 : }
987 :
988 : /************************************************************************/
989 : /* CreateCopy() */
990 : /************************************************************************/
991 :
992 38 : GDALDataset *SAGADataset::CreateCopy(const char *pszFilename,
993 : GDALDataset *poSrcDS, int bStrict,
994 : CPL_UNUSED char **papszOptions,
995 : GDALProgressFunc pfnProgress,
996 : void *pProgressData)
997 : {
998 38 : if (pfnProgress == nullptr)
999 0 : pfnProgress = GDALDummyProgress;
1000 :
1001 38 : int nBands = poSrcDS->GetRasterCount();
1002 38 : if (nBands == 0)
1003 : {
1004 1 : CPLError(
1005 : CE_Failure, CPLE_NotSupported,
1006 : "SAGA driver does not support source dataset with zero band.\n");
1007 1 : return nullptr;
1008 : }
1009 37 : else if (nBands > 1)
1010 : {
1011 4 : if (bStrict)
1012 : {
1013 4 : CPLError(CE_Failure, CPLE_NotSupported,
1014 : "Unable to create copy, SAGA Binary Grid "
1015 : "format only supports one raster band.\n");
1016 4 : return nullptr;
1017 : }
1018 : else
1019 0 : CPLError(CE_Warning, CPLE_NotSupported,
1020 : "SAGA Binary Grid format only supports one "
1021 : "raster band, first band will be copied.\n");
1022 : }
1023 :
1024 33 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(1);
1025 :
1026 33 : char **papszCreateOptions = CSLSetNameValue(nullptr, "FILL_NODATA", "NO");
1027 :
1028 33 : int bHasNoDataValue = FALSE;
1029 33 : const double dfNoDataValue = poSrcBand->GetNoDataValue(&bHasNoDataValue);
1030 33 : if (bHasNoDataValue)
1031 : papszCreateOptions =
1032 2 : CSLSetNameValue(papszCreateOptions, "NODATA_VALUE",
1033 : CPLSPrintf("%.16g", dfNoDataValue));
1034 :
1035 : GDALDataset *poDstDS =
1036 33 : Create(pszFilename, poSrcBand->GetXSize(), poSrcBand->GetYSize(), 1,
1037 : poSrcBand->GetRasterDataType(), papszCreateOptions);
1038 33 : CSLDestroy(papszCreateOptions);
1039 :
1040 33 : if (poDstDS == nullptr)
1041 17 : return nullptr;
1042 :
1043 : /* -------------------------------------------------------------------- */
1044 : /* Copy band data. */
1045 : /* -------------------------------------------------------------------- */
1046 :
1047 : CPLErr eErr =
1048 16 : GDALDatasetCopyWholeRaster((GDALDatasetH)poSrcDS, (GDALDatasetH)poDstDS,
1049 : nullptr, pfnProgress, pProgressData);
1050 :
1051 16 : if (eErr == CE_Failure)
1052 : {
1053 0 : delete poDstDS;
1054 0 : return nullptr;
1055 : }
1056 :
1057 : double adfGeoTransform[6];
1058 :
1059 16 : poSrcDS->GetGeoTransform(adfGeoTransform);
1060 16 : poDstDS->SetGeoTransform(adfGeoTransform);
1061 :
1062 16 : poDstDS->SetProjection(poSrcDS->GetProjectionRef());
1063 :
1064 16 : return poDstDS;
1065 : }
1066 :
1067 : /************************************************************************/
1068 : /* GDALRegister_SAGA() */
1069 : /************************************************************************/
1070 :
1071 1889 : void GDALRegister_SAGA()
1072 :
1073 : {
1074 1889 : if (GDALGetDriverByName("SAGA") != nullptr)
1075 282 : return;
1076 :
1077 1607 : GDALDriver *poDriver = new GDALDriver();
1078 :
1079 1607 : poDriver->SetDescription("SAGA");
1080 1607 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1081 1607 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
1082 1607 : "SAGA GIS Binary Grid (.sdat, .sg-grd-z)");
1083 1607 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/sdat.html");
1084 1607 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSIONS, "sdat sg-grd-z");
1085 1607 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
1086 : "Byte Int16 "
1087 1607 : "UInt16 Int32 UInt32 Float32 Float64");
1088 :
1089 1607 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1090 :
1091 1607 : poDriver->pfnOpen = SAGADataset::Open;
1092 1607 : poDriver->pfnCreate = SAGADataset::Create;
1093 1607 : poDriver->pfnCreateCopy = SAGADataset::CreateCopy;
1094 :
1095 1607 : GetGDALDriverManager()->RegisterDriver(poDriver);
1096 : }
|
