Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: Implements the Golden Software Binary Grid Format.
5 : * Author: Kevin Locke, kwl7@cornell.edu
6 : * (Based largely on aaigriddataset.cpp by Frank Warmerdam)
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2006, Kevin Locke <kwl7@cornell.edu>
10 : * Copyright (c) 2008-2012, Even Rouault <even dot rouault at spatialys.com>
11 : *
12 : * SPDX-License-Identifier: MIT
13 : ****************************************************************************/
14 :
15 : #include "cpl_conv.h"
16 :
17 : #include <assert.h>
18 : #include <float.h>
19 : #include <limits.h>
20 : #include <limits>
21 :
22 : #include "gdal_frmts.h"
23 : #include "gdal_pam.h"
24 : #include "gdal_driver.h"
25 : #include "gdal_drivermanager.h"
26 : #include "gdal_openinfo.h"
27 : #include "gdal_cpp_functions.h"
28 :
29 : /************************************************************************/
30 : /* ==================================================================== */
31 : /* GSBGDataset */
32 : /* ==================================================================== */
33 : /************************************************************************/
34 :
35 : class GSBGRasterBand;
36 :
37 : class GSBGDataset final : public GDALPamDataset
38 : {
39 : friend class GSBGRasterBand;
40 :
41 : static const float fNODATA_VALUE;
42 : static const size_t nHEADER_SIZE;
43 :
44 : static CPLErr WriteHeader(VSILFILE *fp, int nXSize, int nYSize,
45 : double dfMinX, double dfMaxX, double dfMinY,
46 : double dfMaxY, double dfMinZ, double dfMaxZ);
47 :
48 : VSILFILE *fp = nullptr;
49 :
50 : public:
51 45 : GSBGDataset() = default;
52 :
53 : ~GSBGDataset() override;
54 :
55 : static int Identify(GDALOpenInfo *);
56 : static GDALDataset *Open(GDALOpenInfo *);
57 : static GDALPamDataset *OpenPAM(GDALOpenInfo *);
58 : static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
59 : int nBands, GDALDataType eType, CSLConstList);
60 : static GDALDataset *CreateCopy(const char *pszFilename,
61 : GDALDataset *poSrcDS, int bStrict,
62 : CSLConstList papszOptions,
63 : GDALProgressFunc pfnProgress,
64 : void *pProgressData);
65 :
66 : CPLErr GetGeoTransform(GDALGeoTransform >) const override;
67 : CPLErr SetGeoTransform(const GDALGeoTransform >) override;
68 : };
69 :
70 : /* NOTE: This is not mentioned in the spec, but Surfer 8 uses this value */
71 : /* 0x7effffee (Little Endian: eeffff7e) */
72 : const float GSBGDataset::fNODATA_VALUE = 1.701410009187828e+38f;
73 :
74 : const size_t GSBGDataset::nHEADER_SIZE = 56;
75 :
76 : /************************************************************************/
77 : /* ==================================================================== */
78 : /* GSBGRasterBand */
79 : /* ==================================================================== */
80 : /************************************************************************/
81 :
82 : class GSBGRasterBand final : public GDALPamRasterBand
83 : {
84 : friend class GSBGDataset;
85 :
86 : double dfMinX;
87 : double dfMaxX;
88 : double dfMinY;
89 : double dfMaxY;
90 : double dfMinZ;
91 : double dfMaxZ;
92 :
93 : float *pafRowMinZ;
94 : float *pafRowMaxZ;
95 : int nMinZRow;
96 : int nMaxZRow;
97 :
98 : CPLErr ScanForMinMaxZ();
99 :
100 : public:
101 : GSBGRasterBand(GSBGDataset *, int);
102 : ~GSBGRasterBand() override;
103 :
104 : CPLErr IReadBlock(int, int, void *) override;
105 : CPLErr IWriteBlock(int, int, void *) override;
106 :
107 : double GetNoDataValue(int *pbSuccess = nullptr) override;
108 : double GetMinimum(int *pbSuccess = nullptr) override;
109 : double GetMaximum(int *pbSuccess = nullptr) override;
110 : };
111 :
112 : /************************************************************************/
113 : /* GSBGRasterBand() */
114 : /************************************************************************/
115 :
116 45 : GSBGRasterBand::GSBGRasterBand(GSBGDataset *poDSIn, int nBandIn)
117 : : dfMinX(0.0), dfMaxX(0.0), dfMinY(0.0), dfMaxY(0.0), dfMinZ(0.0),
118 : dfMaxZ(0.0), pafRowMinZ(nullptr), pafRowMaxZ(nullptr), nMinZRow(-1),
119 45 : nMaxZRow(-1)
120 : {
121 45 : this->poDS = poDSIn;
122 45 : this->nBand = nBandIn;
123 :
124 45 : eDataType = GDT_Float32;
125 :
126 45 : nBlockXSize = poDS->GetRasterXSize();
127 45 : nBlockYSize = 1;
128 45 : }
129 :
130 : /************************************************************************/
131 : /* ~GSBGRasterBand() */
132 : /************************************************************************/
133 :
134 90 : GSBGRasterBand::~GSBGRasterBand()
135 :
136 : {
137 45 : if (pafRowMinZ != nullptr)
138 1 : CPLFree(pafRowMinZ);
139 45 : if (pafRowMaxZ != nullptr)
140 1 : CPLFree(pafRowMaxZ);
141 90 : }
142 :
143 : /************************************************************************/
144 : /* ScanForMinMaxZ() */
145 : /************************************************************************/
146 :
147 1 : CPLErr GSBGRasterBand::ScanForMinMaxZ()
148 :
149 : {
150 1 : float *pafRowVals = (float *)VSI_MALLOC2_VERBOSE(nRasterXSize, 4);
151 :
152 1 : if (pafRowVals == nullptr)
153 : {
154 0 : return CE_Failure;
155 : }
156 :
157 1 : double dfNewMinZ = std::numeric_limits<double>::max();
158 1 : double dfNewMaxZ = std::numeric_limits<double>::lowest();
159 1 : int nNewMinZRow = 0;
160 1 : int nNewMaxZRow = 0;
161 :
162 : /* Since we have to scan, lets calc. statistics too */
163 1 : double dfSum = 0.0;
164 1 : double dfSum2 = 0.0;
165 1 : unsigned long nValuesRead = 0;
166 21 : for (int iRow = 0; iRow < nRasterYSize; iRow++)
167 : {
168 20 : CPLErr eErr = IReadBlock(0, iRow, pafRowVals);
169 20 : if (eErr != CE_None)
170 : {
171 0 : VSIFree(pafRowVals);
172 0 : return CE_Failure;
173 : }
174 :
175 20 : pafRowMinZ[iRow] = std::numeric_limits<float>::max();
176 20 : pafRowMaxZ[iRow] = std::numeric_limits<float>::lowest();
177 420 : for (int iCol = 0; iCol < nRasterXSize; iCol++)
178 : {
179 400 : if (pafRowVals[iCol] == GSBGDataset::fNODATA_VALUE)
180 400 : continue;
181 :
182 0 : if (pafRowVals[iCol] < pafRowMinZ[iRow])
183 0 : pafRowMinZ[iRow] = pafRowVals[iCol];
184 :
185 0 : if (pafRowVals[iCol] > pafRowMinZ[iRow])
186 0 : pafRowMaxZ[iRow] = pafRowVals[iCol];
187 :
188 0 : dfSum += pafRowVals[iCol];
189 0 : dfSum2 += static_cast<double>(pafRowVals[iCol]) * pafRowVals[iCol];
190 0 : nValuesRead++;
191 : }
192 :
193 20 : if (pafRowMinZ[iRow] < dfNewMinZ)
194 : {
195 1 : dfNewMinZ = pafRowMinZ[iRow];
196 1 : nNewMinZRow = iRow;
197 : }
198 :
199 20 : if (pafRowMaxZ[iRow] > dfNewMaxZ)
200 : {
201 1 : dfNewMaxZ = pafRowMaxZ[iRow];
202 1 : nNewMaxZRow = iRow;
203 : }
204 : }
205 :
206 1 : VSIFree(pafRowVals);
207 :
208 1 : if (nValuesRead == 0)
209 : {
210 1 : dfMinZ = 0.0;
211 1 : dfMaxZ = 0.0;
212 1 : nMinZRow = 0;
213 1 : nMaxZRow = 0;
214 1 : return CE_None;
215 : }
216 :
217 0 : dfMinZ = dfNewMinZ;
218 0 : dfMaxZ = dfNewMaxZ;
219 0 : nMinZRow = nNewMinZRow;
220 0 : nMaxZRow = nNewMaxZRow;
221 :
222 0 : double dfMean = dfSum / nValuesRead;
223 0 : double dfStdDev = sqrt((dfSum2 / nValuesRead) - (dfMean * dfMean));
224 0 : SetStatistics(dfMinZ, dfMaxZ, dfMean, dfStdDev);
225 :
226 0 : return CE_None;
227 : }
228 :
229 : /************************************************************************/
230 : /* IReadBlock() */
231 : /************************************************************************/
232 :
233 140 : CPLErr GSBGRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage)
234 :
235 : {
236 140 : if (nBlockYOff < 0 || nBlockYOff > nRasterYSize - 1 || nBlockXOff != 0)
237 0 : return CE_Failure;
238 :
239 140 : GSBGDataset *poGDS = cpl::down_cast<GSBGDataset *>(poDS);
240 280 : if (VSIFSeekL(poGDS->fp,
241 140 : GSBGDataset::nHEADER_SIZE +
242 140 : 4 * static_cast<vsi_l_offset>(nRasterXSize) *
243 140 : (nRasterYSize - nBlockYOff - 1),
244 140 : SEEK_SET) != 0)
245 : {
246 0 : CPLError(CE_Failure, CPLE_FileIO,
247 : "Unable to seek to beginning of grid row.");
248 0 : return CE_Failure;
249 : }
250 :
251 140 : if (VSIFReadL(pImage, sizeof(float), nBlockXSize, poGDS->fp) !=
252 140 : static_cast<unsigned>(nBlockXSize))
253 : {
254 0 : CPLError(CE_Failure, CPLE_FileIO,
255 : "Unable to read block from grid file.");
256 0 : return CE_Failure;
257 : }
258 :
259 : #ifdef CPL_MSB
260 : float *pfImage = (float *)pImage;
261 : for (int iPixel = 0; iPixel < nBlockXSize; iPixel++)
262 : {
263 : CPL_LSBPTR32(pfImage + iPixel);
264 : }
265 : #endif
266 :
267 140 : return CE_None;
268 : }
269 :
270 : /************************************************************************/
271 : /* IWriteBlock() */
272 : /************************************************************************/
273 :
274 20 : CPLErr GSBGRasterBand::IWriteBlock(int nBlockXOff, int nBlockYOff, void *pImage)
275 :
276 : {
277 20 : if (eAccess == GA_ReadOnly)
278 : {
279 0 : CPLError(CE_Failure, CPLE_NoWriteAccess,
280 : "Unable to write block, dataset opened read only.");
281 0 : return CE_Failure;
282 : }
283 :
284 20 : if (nBlockYOff < 0 || nBlockYOff > nRasterYSize - 1 || nBlockXOff != 0)
285 0 : return CE_Failure;
286 :
287 20 : GSBGDataset *poGDS = cpl::down_cast<GSBGDataset *>(poDS);
288 :
289 20 : if (pafRowMinZ == nullptr || pafRowMaxZ == nullptr || nMinZRow < 0 ||
290 19 : nMaxZRow < 0)
291 : {
292 1 : pafRowMinZ = (float *)VSI_MALLOC2_VERBOSE(nRasterYSize, sizeof(float));
293 1 : if (pafRowMinZ == nullptr)
294 : {
295 0 : return CE_Failure;
296 : }
297 :
298 1 : pafRowMaxZ = (float *)VSI_MALLOC2_VERBOSE(nRasterYSize, sizeof(float));
299 1 : if (pafRowMaxZ == nullptr)
300 : {
301 0 : VSIFree(pafRowMinZ);
302 0 : pafRowMinZ = nullptr;
303 0 : return CE_Failure;
304 : }
305 :
306 1 : CPLErr eErr = ScanForMinMaxZ();
307 1 : if (eErr != CE_None)
308 0 : return eErr;
309 : }
310 :
311 40 : if (VSIFSeekL(poGDS->fp,
312 20 : GSBGDataset::nHEADER_SIZE +
313 20 : static_cast<vsi_l_offset>(4) * nRasterXSize *
314 20 : (nRasterYSize - nBlockYOff - 1),
315 20 : SEEK_SET) != 0)
316 : {
317 0 : CPLError(CE_Failure, CPLE_FileIO,
318 : "Unable to seek to beginning of grid row.");
319 0 : return CE_Failure;
320 : }
321 :
322 20 : float *pfImage = (float *)pImage;
323 20 : pafRowMinZ[nBlockYOff] = std::numeric_limits<float>::max();
324 20 : pafRowMaxZ[nBlockYOff] = std::numeric_limits<float>::lowest();
325 420 : for (int iPixel = 0; iPixel < nBlockXSize; iPixel++)
326 : {
327 400 : if (pfImage[iPixel] != GSBGDataset::fNODATA_VALUE)
328 : {
329 400 : if (pfImage[iPixel] < pafRowMinZ[nBlockYOff])
330 78 : pafRowMinZ[nBlockYOff] = pfImage[iPixel];
331 :
332 400 : if (pfImage[iPixel] > pafRowMaxZ[nBlockYOff])
333 43 : pafRowMaxZ[nBlockYOff] = pfImage[iPixel];
334 : }
335 :
336 400 : CPL_LSBPTR32(pfImage + iPixel);
337 : }
338 :
339 20 : if (VSIFWriteL(pImage, sizeof(float), nBlockXSize, poGDS->fp) !=
340 20 : static_cast<unsigned>(nBlockXSize))
341 : {
342 0 : CPLError(CE_Failure, CPLE_FileIO,
343 : "Unable to write block to grid file.");
344 0 : return CE_Failure;
345 : }
346 :
347 : /* Update min/max Z values as appropriate */
348 20 : bool bHeaderNeedsUpdate = false;
349 20 : if (nMinZRow == nBlockYOff && pafRowMinZ[nBlockYOff] > dfMinZ)
350 : {
351 1 : double dfNewMinZ = std::numeric_limits<double>::max();
352 21 : for (int iRow = 0; iRow < nRasterYSize; iRow++)
353 : {
354 20 : if (pafRowMinZ[iRow] < dfNewMinZ)
355 : {
356 1 : dfNewMinZ = pafRowMinZ[iRow];
357 1 : nMinZRow = iRow;
358 : }
359 : }
360 :
361 1 : if (dfNewMinZ != dfMinZ)
362 : {
363 1 : dfMinZ = dfNewMinZ;
364 1 : bHeaderNeedsUpdate = true;
365 : }
366 : }
367 :
368 20 : if (nMaxZRow == nBlockYOff && pafRowMaxZ[nBlockYOff] < dfMaxZ)
369 : {
370 0 : double dfNewMaxZ = std::numeric_limits<double>::lowest();
371 0 : for (int iRow = 0; iRow < nRasterYSize; iRow++)
372 : {
373 0 : if (pafRowMaxZ[iRow] > dfNewMaxZ)
374 : {
375 0 : dfNewMaxZ = pafRowMaxZ[iRow];
376 0 : nMaxZRow = iRow;
377 : }
378 : }
379 :
380 0 : if (dfNewMaxZ != dfMaxZ)
381 : {
382 0 : dfMaxZ = dfNewMaxZ;
383 0 : bHeaderNeedsUpdate = true;
384 : }
385 : }
386 :
387 20 : if (pafRowMinZ[nBlockYOff] < dfMinZ || pafRowMaxZ[nBlockYOff] > dfMaxZ)
388 : {
389 6 : if (pafRowMinZ[nBlockYOff] < dfMinZ)
390 : {
391 3 : dfMinZ = pafRowMinZ[nBlockYOff];
392 3 : nMinZRow = nBlockYOff;
393 : }
394 :
395 6 : if (pafRowMaxZ[nBlockYOff] > dfMaxZ)
396 : {
397 4 : dfMaxZ = pafRowMaxZ[nBlockYOff];
398 4 : nMaxZRow = nBlockYOff;
399 : }
400 :
401 6 : bHeaderNeedsUpdate = true;
402 : }
403 :
404 20 : if (bHeaderNeedsUpdate && dfMaxZ > dfMinZ)
405 : {
406 : CPLErr eErr =
407 6 : poGDS->WriteHeader(poGDS->fp, nRasterXSize, nRasterYSize, dfMinX,
408 : dfMaxX, dfMinY, dfMaxY, dfMinZ, dfMaxZ);
409 6 : return eErr;
410 : }
411 :
412 14 : return CE_None;
413 : }
414 :
415 : /************************************************************************/
416 : /* GetNoDataValue() */
417 : /************************************************************************/
418 :
419 27 : double GSBGRasterBand::GetNoDataValue(int *pbSuccess)
420 : {
421 27 : if (pbSuccess)
422 21 : *pbSuccess = TRUE;
423 :
424 27 : return GSBGDataset::fNODATA_VALUE;
425 : }
426 :
427 : /************************************************************************/
428 : /* GetMinimum() */
429 : /************************************************************************/
430 :
431 0 : double GSBGRasterBand::GetMinimum(int *pbSuccess)
432 : {
433 0 : if (pbSuccess)
434 0 : *pbSuccess = TRUE;
435 :
436 0 : return dfMinZ;
437 : }
438 :
439 : /************************************************************************/
440 : /* GetMaximum() */
441 : /************************************************************************/
442 :
443 0 : double GSBGRasterBand::GetMaximum(int *pbSuccess)
444 : {
445 0 : if (pbSuccess)
446 0 : *pbSuccess = TRUE;
447 :
448 0 : return dfMaxZ;
449 : }
450 :
451 : /************************************************************************/
452 : /* ==================================================================== */
453 : /* GSBGDataset */
454 : /* ==================================================================== */
455 : /************************************************************************/
456 :
457 90 : GSBGDataset::~GSBGDataset()
458 :
459 : {
460 45 : FlushCache(true);
461 45 : if (fp != nullptr)
462 45 : VSIFCloseL(fp);
463 90 : }
464 :
465 : /************************************************************************/
466 : /* Identify() */
467 : /************************************************************************/
468 :
469 64201 : int GSBGDataset::Identify(GDALOpenInfo *poOpenInfo)
470 :
471 : {
472 : /* Check for signature */
473 64201 : if (poOpenInfo->nHeaderBytes < 4 ||
474 6250 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "DSBB"))
475 : {
476 64136 : return FALSE;
477 : }
478 :
479 65 : return TRUE;
480 : }
481 :
482 : /************************************************************************/
483 : /* Open() */
484 : /************************************************************************/
485 :
486 33 : GDALDataset *GSBGDataset::Open(GDALOpenInfo *poOpenInfo)
487 : {
488 33 : return OpenPAM(poOpenInfo);
489 : }
490 :
491 45 : GDALPamDataset *GSBGDataset::OpenPAM(GDALOpenInfo *poOpenInfo)
492 :
493 : {
494 45 : if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
495 : {
496 0 : return nullptr;
497 : }
498 :
499 : /* -------------------------------------------------------------------- */
500 : /* Create a corresponding GDALDataset. */
501 : /* -------------------------------------------------------------------- */
502 90 : auto poDS = std::make_unique<GSBGDataset>();
503 :
504 45 : poDS->eAccess = poOpenInfo->eAccess;
505 45 : poDS->fp = poOpenInfo->fpL;
506 45 : poOpenInfo->fpL = nullptr;
507 :
508 : /* -------------------------------------------------------------------- */
509 : /* Read the header. */
510 : /* -------------------------------------------------------------------- */
511 45 : if (VSIFSeekL(poDS->fp, 4, SEEK_SET) != 0)
512 : {
513 0 : CPLError(CE_Failure, CPLE_FileIO,
514 : "Unable to seek to start of grid file header.");
515 0 : return nullptr;
516 : }
517 :
518 : /* Parse number of X axis grid rows */
519 : GInt16 nTemp;
520 45 : if (VSIFReadL((void *)&nTemp, 2, 1, poDS->fp) != 1)
521 : {
522 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to read raster X size.");
523 0 : return nullptr;
524 : }
525 45 : poDS->nRasterXSize = CPL_LSBWORD16(nTemp);
526 :
527 45 : if (VSIFReadL((void *)&nTemp, 2, 1, poDS->fp) != 1)
528 : {
529 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to read raster Y size.");
530 0 : return nullptr;
531 : }
532 45 : poDS->nRasterYSize = CPL_LSBWORD16(nTemp);
533 :
534 45 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
535 : {
536 0 : return nullptr;
537 : }
538 :
539 : /* -------------------------------------------------------------------- */
540 : /* Create band information objects. */
541 : /* -------------------------------------------------------------------- */
542 45 : GSBGRasterBand *poBand = new GSBGRasterBand(poDS.get(), 1);
543 45 : poDS->SetBand(1, poBand);
544 :
545 : double dfTemp;
546 45 : if (VSIFReadL((void *)&dfTemp, 8, 1, poDS->fp) != 1)
547 : {
548 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to read minimum X value.");
549 0 : return nullptr;
550 : }
551 45 : CPL_LSBPTR64(&dfTemp);
552 45 : poBand->dfMinX = dfTemp;
553 :
554 45 : if (VSIFReadL((void *)&dfTemp, 8, 1, poDS->fp) != 1)
555 : {
556 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to read maximum X value.");
557 0 : return nullptr;
558 : }
559 45 : CPL_LSBPTR64(&dfTemp);
560 45 : poBand->dfMaxX = dfTemp;
561 :
562 45 : if (VSIFReadL((void *)&dfTemp, 8, 1, poDS->fp) != 1)
563 : {
564 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to read minimum Y value.");
565 0 : return nullptr;
566 : }
567 45 : CPL_LSBPTR64(&dfTemp);
568 45 : poBand->dfMinY = dfTemp;
569 :
570 45 : if (VSIFReadL((void *)&dfTemp, 8, 1, poDS->fp) != 1)
571 : {
572 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to read maximum Y value.");
573 0 : return nullptr;
574 : }
575 45 : CPL_LSBPTR64(&dfTemp);
576 45 : poBand->dfMaxY = dfTemp;
577 :
578 45 : if (VSIFReadL((void *)&dfTemp, 8, 1, poDS->fp) != 1)
579 : {
580 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to read minimum Z value.");
581 0 : return nullptr;
582 : }
583 45 : CPL_LSBPTR64(&dfTemp);
584 45 : poBand->dfMinZ = dfTemp;
585 :
586 45 : if (VSIFReadL((void *)&dfTemp, 8, 1, poDS->fp) != 1)
587 : {
588 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to read maximum Z value.");
589 0 : return nullptr;
590 : }
591 45 : CPL_LSBPTR64(&dfTemp);
592 45 : poBand->dfMaxZ = dfTemp;
593 :
594 : /* -------------------------------------------------------------------- */
595 : /* Initialize any PAM information. */
596 : /* -------------------------------------------------------------------- */
597 45 : poDS->SetDescription(poOpenInfo->pszFilename);
598 45 : poDS->TryLoadXML();
599 :
600 : /* -------------------------------------------------------------------- */
601 : /* Check for external overviews. */
602 : /* -------------------------------------------------------------------- */
603 90 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename,
604 45 : poOpenInfo->GetSiblingFiles());
605 :
606 45 : return poDS.release();
607 : }
608 :
609 : /************************************************************************/
610 : /* GetGeoTransform() */
611 : /************************************************************************/
612 :
613 29 : CPLErr GSBGDataset::GetGeoTransform(GDALGeoTransform >) const
614 : {
615 : const GSBGRasterBand *poGRB =
616 29 : cpl::down_cast<const GSBGRasterBand *>(GetRasterBand(1));
617 :
618 : /* check if we have a PAM GeoTransform stored */
619 29 : CPLPushErrorHandler(CPLQuietErrorHandler);
620 29 : CPLErr eErr = GDALPamDataset::GetGeoTransform(gt);
621 29 : CPLPopErrorHandler();
622 :
623 29 : if (eErr == CE_None)
624 0 : return CE_None;
625 :
626 29 : if (nRasterXSize == 1 || nRasterYSize == 1)
627 0 : return CE_Failure;
628 :
629 : /* calculate pixel size first */
630 29 : gt.xscale = (poGRB->dfMaxX - poGRB->dfMinX) / (nRasterXSize - 1);
631 29 : gt.yscale = (poGRB->dfMinY - poGRB->dfMaxY) / (nRasterYSize - 1);
632 :
633 : /* then calculate image origin */
634 29 : gt.xorig = poGRB->dfMinX - gt.xscale / 2;
635 29 : gt.yorig = poGRB->dfMaxY - gt.yscale / 2;
636 :
637 : /* tilt/rotation does not supported by the GS grids */
638 29 : gt.yrot = 0.0;
639 29 : gt.xrot = 0.0;
640 :
641 29 : return CE_None;
642 : }
643 :
644 : /************************************************************************/
645 : /* SetGeoTransform() */
646 : /************************************************************************/
647 :
648 12 : CPLErr GSBGDataset::SetGeoTransform(const GDALGeoTransform >)
649 : {
650 12 : if (eAccess == GA_ReadOnly)
651 : {
652 0 : CPLError(CE_Failure, CPLE_NoWriteAccess,
653 : "Unable to set GeoTransform, dataset opened read only.");
654 0 : return CE_Failure;
655 : }
656 :
657 12 : GSBGRasterBand *poGRB = cpl::down_cast<GSBGRasterBand *>(GetRasterBand(1));
658 :
659 : /* non-zero transform 2 or 4 or negative 1 or 5 not supported natively */
660 : // CPLErr eErr = CE_None;
661 : /*if( gt.xrot != 0.0 || gt.yrot != 0.0
662 : || gt.xscale < 0.0 || gt.yscale < 0.0 )
663 : eErr = GDALPamDataset::SetGeoTransform( gt );
664 :
665 : if( eErr != CE_None )
666 : return eErr;*/
667 :
668 12 : double dfMinX = gt.xorig + gt.xscale / 2;
669 12 : double dfMaxX = gt.xscale * (nRasterXSize - 0.5) + gt.xorig;
670 12 : double dfMinY = gt.yscale * (nRasterYSize - 0.5) + gt.yorig;
671 12 : double dfMaxY = gt.yorig + gt.yscale / 2;
672 :
673 : CPLErr eErr =
674 12 : WriteHeader(fp, poGRB->nRasterXSize, poGRB->nRasterYSize, dfMinX,
675 : dfMaxX, dfMinY, dfMaxY, poGRB->dfMinZ, poGRB->dfMaxZ);
676 :
677 12 : if (eErr == CE_None)
678 : {
679 12 : poGRB->dfMinX = dfMinX;
680 12 : poGRB->dfMaxX = dfMaxX;
681 12 : poGRB->dfMinY = dfMinY;
682 12 : poGRB->dfMaxY = dfMaxY;
683 : }
684 :
685 12 : return eErr;
686 : }
687 :
688 : /************************************************************************/
689 : /* WriteHeader() */
690 : /************************************************************************/
691 :
692 65 : CPLErr GSBGDataset::WriteHeader(VSILFILE *fp, int nXSize, int nYSize,
693 : double dfMinX, double dfMaxX, double dfMinY,
694 : double dfMaxY, double dfMinZ, double dfMaxZ)
695 :
696 : {
697 65 : if (VSIFSeekL(fp, 0, SEEK_SET) != 0)
698 : {
699 0 : CPLError(CE_Failure, CPLE_FileIO,
700 : "Unable to seek to start of grid file.");
701 0 : return CE_Failure;
702 : }
703 :
704 65 : if (VSIFWriteL((void *)"DSBB", 1, 4, fp) != 4)
705 : {
706 1 : CPLError(CE_Failure, CPLE_FileIO,
707 : "Unable to write signature to grid file.");
708 1 : return CE_Failure;
709 : }
710 :
711 64 : assert(nXSize >= 0 && nXSize <= std::numeric_limits<int16_t>::max());
712 64 : GInt16 nTemp = CPL_LSBWORD16(static_cast<int16_t>(nXSize));
713 64 : if (VSIFWriteL((void *)&nTemp, 2, 1, fp) != 1)
714 : {
715 0 : CPLError(CE_Failure, CPLE_FileIO,
716 : "Unable to write raster X size to grid file.");
717 0 : return CE_Failure;
718 : }
719 :
720 64 : assert(nYSize >= 0 && nYSize <= std::numeric_limits<int16_t>::max());
721 64 : nTemp = CPL_LSBWORD16(static_cast<int16_t>(nYSize));
722 64 : if (VSIFWriteL((void *)&nTemp, 2, 1, fp) != 1)
723 : {
724 0 : CPLError(CE_Failure, CPLE_FileIO,
725 : "Unable to write raster Y size to grid file.");
726 0 : return CE_Failure;
727 : }
728 :
729 64 : double dfTemp = dfMinX;
730 64 : CPL_LSBPTR64(&dfTemp);
731 64 : if (VSIFWriteL((void *)&dfTemp, 8, 1, fp) != 1)
732 : {
733 0 : CPLError(CE_Failure, CPLE_FileIO,
734 : "Unable to write minimum X value to grid file.");
735 0 : return CE_Failure;
736 : }
737 :
738 64 : dfTemp = dfMaxX;
739 64 : CPL_LSBPTR64(&dfTemp);
740 64 : if (VSIFWriteL((void *)&dfTemp, 8, 1, fp) != 1)
741 : {
742 0 : CPLError(CE_Failure, CPLE_FileIO,
743 : "Unable to write maximum X value to grid file.");
744 0 : return CE_Failure;
745 : }
746 :
747 64 : dfTemp = dfMinY;
748 64 : CPL_LSBPTR64(&dfTemp);
749 64 : if (VSIFWriteL((void *)&dfTemp, 8, 1, fp) != 1)
750 : {
751 0 : CPLError(CE_Failure, CPLE_FileIO,
752 : "Unable to write minimum Y value to grid file.");
753 0 : return CE_Failure;
754 : }
755 :
756 64 : dfTemp = dfMaxY;
757 64 : CPL_LSBPTR64(&dfTemp);
758 64 : if (VSIFWriteL((void *)&dfTemp, 8, 1, fp) != 1)
759 : {
760 0 : CPLError(CE_Failure, CPLE_FileIO,
761 : "Unable to write maximum Y value to grid file.");
762 0 : return CE_Failure;
763 : }
764 :
765 64 : dfTemp = dfMinZ;
766 64 : CPL_LSBPTR64(&dfTemp);
767 64 : if (VSIFWriteL((void *)&dfTemp, 8, 1, fp) != 1)
768 : {
769 1 : CPLError(CE_Failure, CPLE_FileIO,
770 : "Unable to write minimum Z value to grid file.");
771 1 : return CE_Failure;
772 : }
773 :
774 63 : dfTemp = dfMaxZ;
775 63 : CPL_LSBPTR64(&dfTemp);
776 63 : if (VSIFWriteL((void *)&dfTemp, 8, 1, fp) != 1)
777 : {
778 0 : CPLError(CE_Failure, CPLE_FileIO,
779 : "Unable to write maximum Z value to grid file.");
780 0 : return CE_Failure;
781 : }
782 :
783 63 : return CE_None;
784 : }
785 :
786 : /************************************************************************/
787 : /* GSBGCreateCheckDims() */
788 : /************************************************************************/
789 :
790 66 : static bool GSBGCreateCheckDims(int nXSize, int nYSize)
791 : {
792 66 : if (nXSize <= 1 || nYSize <= 1)
793 : {
794 4 : CPLError(CE_Failure, CPLE_IllegalArg,
795 : "Unable to create grid, both X and Y size must be "
796 : "larger or equal to 2.");
797 4 : return false;
798 : }
799 122 : if (nXSize > std::numeric_limits<short>::max() ||
800 60 : nYSize > std::numeric_limits<short>::max())
801 : {
802 4 : CPLError(CE_Failure, CPLE_IllegalArg,
803 : "Unable to create grid, Golden Software Binary Grid format "
804 : "only supports sizes up to %dx%d. %dx%d not supported.",
805 4 : std::numeric_limits<short>::max(),
806 4 : std::numeric_limits<short>::max(), nXSize, nYSize);
807 4 : return false;
808 : }
809 58 : return true;
810 : }
811 :
812 : /************************************************************************/
813 : /* Create() */
814 : /************************************************************************/
815 :
816 37 : GDALDataset *GSBGDataset::Create(const char *pszFilename, int nXSize,
817 : int nYSize, int /* nBands */,
818 : GDALDataType eType, CSLConstList)
819 : {
820 37 : if (!GSBGCreateCheckDims(nXSize, nYSize))
821 : {
822 4 : return nullptr;
823 : }
824 :
825 33 : if (eType != GDT_UInt8 && eType != GDT_Float32 && eType != GDT_UInt16 &&
826 : eType != GDT_Int16)
827 : {
828 20 : CPLError(CE_Failure, CPLE_AppDefined,
829 : "Golden Software Binary Grid only supports Byte, Int16, "
830 : "Uint16, and Float32 datatypes. Unable to create with "
831 : "type %s.\n",
832 : GDALGetDataTypeName(eType));
833 :
834 20 : return nullptr;
835 : }
836 :
837 13 : VSILFILE *fp = VSIFOpenL(pszFilename, "w+b");
838 :
839 13 : if (fp == nullptr)
840 : {
841 0 : CPLError(CE_Failure, CPLE_OpenFailed,
842 : "Attempt to create file '%s' failed.\n", pszFilename);
843 0 : return nullptr;
844 : }
845 :
846 : CPLErr eErr =
847 13 : WriteHeader(fp, nXSize, nYSize, 0.0, nXSize, 0.0, nYSize, 0.0, 0.0);
848 13 : if (eErr != CE_None)
849 : {
850 0 : VSIFCloseL(fp);
851 0 : return nullptr;
852 : }
853 :
854 13 : float fVal = fNODATA_VALUE;
855 13 : CPL_LSBPTR32(&fVal);
856 1233 : for (int iRow = 0; iRow < nYSize; iRow++)
857 : {
858 121620 : for (int iCol = 0; iCol < nXSize; iCol++)
859 : {
860 120400 : if (VSIFWriteL((void *)&fVal, 4, 1, fp) != 1)
861 : {
862 0 : VSIFCloseL(fp);
863 0 : CPLError(CE_Failure, CPLE_FileIO,
864 : "Unable to write grid cell. Disk full?");
865 0 : return nullptr;
866 : }
867 : }
868 : }
869 :
870 13 : VSIFCloseL(fp);
871 :
872 26 : GDALOpenInfo oOpenInfo(pszFilename, GA_Update);
873 13 : return Open(&oOpenInfo);
874 : }
875 :
876 : /************************************************************************/
877 : /* CreateCopy() */
878 : /************************************************************************/
879 :
880 34 : GDALDataset *GSBGDataset::CreateCopy(const char *pszFilename,
881 : GDALDataset *poSrcDS, int bStrict,
882 : CPL_UNUSED CSLConstList papszOptions,
883 : GDALProgressFunc pfnProgress,
884 : void *pProgressData)
885 : {
886 34 : if (pfnProgress == nullptr)
887 0 : pfnProgress = GDALDummyProgress;
888 :
889 34 : int nBands = poSrcDS->GetRasterCount();
890 34 : if (nBands == 0)
891 : {
892 1 : CPLError(
893 : CE_Failure, CPLE_NotSupported,
894 : "GSBG driver does not support source datasets with zero bands.");
895 1 : return nullptr;
896 : }
897 33 : else if (nBands > 1)
898 : {
899 4 : if (bStrict)
900 : {
901 4 : CPLError(CE_Failure, CPLE_NotSupported,
902 : "Unable to create copy, Golden Software Binary Grid "
903 : "format only supports one raster band.");
904 4 : return nullptr;
905 : }
906 : else
907 0 : CPLError(CE_Warning, CPLE_NotSupported,
908 : "Golden Software Binary Grid format only supports one "
909 : "raster band, first band will be copied.");
910 : }
911 :
912 29 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(1);
913 29 : if (!GSBGCreateCheckDims(poSrcBand->GetXSize(), poSrcBand->GetYSize()))
914 : {
915 4 : return nullptr;
916 : }
917 :
918 25 : if (!pfnProgress(0.0, nullptr, pProgressData))
919 : {
920 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
921 0 : return nullptr;
922 : }
923 :
924 25 : VSILFILE *fp = VSIFOpenL(pszFilename, "w+b");
925 :
926 25 : if (fp == nullptr)
927 : {
928 3 : CPLError(CE_Failure, CPLE_OpenFailed,
929 : "Attempt to create file '%s' failed.\n", pszFilename);
930 3 : return nullptr;
931 : }
932 :
933 22 : const int nXSize = poSrcBand->GetXSize();
934 22 : const int nYSize = poSrcBand->GetYSize();
935 22 : GDALGeoTransform gt;
936 22 : poSrcDS->GetGeoTransform(gt);
937 :
938 22 : double dfMinX = gt.xorig + gt.xscale / 2;
939 22 : double dfMaxX = gt.xscale * (nXSize - 0.5) + gt.xorig;
940 22 : double dfMinY = gt.yscale * (nYSize - 0.5) + gt.yorig;
941 22 : double dfMaxY = gt.yorig + gt.yscale / 2;
942 22 : CPLErr eErr = WriteHeader(fp, nXSize, nYSize, dfMinX, dfMaxX, dfMinY,
943 : dfMaxY, 0.0, 0.0);
944 :
945 22 : if (eErr != CE_None)
946 : {
947 2 : VSIFCloseL(fp);
948 2 : return nullptr;
949 : }
950 :
951 : /* -------------------------------------------------------------------- */
952 : /* Copy band data. */
953 : /* -------------------------------------------------------------------- */
954 20 : float *pfData = (float *)VSI_MALLOC2_VERBOSE(nXSize, sizeof(float));
955 20 : if (pfData == nullptr)
956 : {
957 0 : VSIFCloseL(fp);
958 0 : return nullptr;
959 : }
960 :
961 : int bSrcHasNDValue;
962 20 : float fSrcNoDataValue = (float)poSrcBand->GetNoDataValue(&bSrcHasNDValue);
963 20 : double dfMinZ = std::numeric_limits<double>::max();
964 20 : double dfMaxZ = std::numeric_limits<double>::lowest();
965 185 : for (int iRow = nYSize - 1; iRow >= 0; iRow--)
966 : {
967 173 : eErr = poSrcBand->RasterIO(GF_Read, 0, iRow, nXSize, 1, pfData, nXSize,
968 : 1, GDT_Float32, 0, 0, nullptr);
969 :
970 173 : if (eErr != CE_None)
971 : {
972 0 : VSIFCloseL(fp);
973 0 : VSIFree(pfData);
974 0 : return nullptr;
975 : }
976 :
977 2103 : for (int iCol = 0; iCol < nXSize; iCol++)
978 : {
979 1930 : if (bSrcHasNDValue && pfData[iCol] == fSrcNoDataValue)
980 : {
981 0 : pfData[iCol] = fNODATA_VALUE;
982 : }
983 : else
984 : {
985 1930 : if (pfData[iCol] > dfMaxZ)
986 22 : dfMaxZ = pfData[iCol];
987 :
988 1930 : if (pfData[iCol] < dfMinZ)
989 27 : dfMinZ = pfData[iCol];
990 : }
991 :
992 1930 : CPL_LSBPTR32(pfData + iCol);
993 : }
994 :
995 173 : if (VSIFWriteL((void *)pfData, 4, nXSize, fp) !=
996 173 : static_cast<unsigned>(nXSize))
997 : {
998 8 : VSIFCloseL(fp);
999 8 : VSIFree(pfData);
1000 8 : CPLError(CE_Failure, CPLE_FileIO,
1001 : "Unable to write grid row. Disk full?");
1002 8 : return nullptr;
1003 : }
1004 :
1005 165 : if (!pfnProgress(static_cast<double>(nYSize - iRow) / nYSize, nullptr,
1006 : pProgressData))
1007 : {
1008 0 : VSIFCloseL(fp);
1009 0 : VSIFree(pfData);
1010 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
1011 0 : return nullptr;
1012 : }
1013 : }
1014 :
1015 12 : VSIFree(pfData);
1016 :
1017 : /* write out the min and max values */
1018 12 : eErr = WriteHeader(fp, nXSize, nYSize, dfMinX, dfMaxX, dfMinY, dfMaxY,
1019 : dfMinZ, dfMaxZ);
1020 :
1021 12 : if (eErr != CE_None)
1022 : {
1023 0 : VSIFCloseL(fp);
1024 0 : return nullptr;
1025 : }
1026 :
1027 12 : VSIFCloseL(fp);
1028 :
1029 12 : GDALOpenInfo oOpenInfo(pszFilename, GA_Update);
1030 12 : GDALPamDataset *poDS = OpenPAM(&oOpenInfo);
1031 12 : if (poDS)
1032 : {
1033 12 : poDS->CloneInfo(poSrcDS, GCIF_PAM_DEFAULT);
1034 : }
1035 12 : return poDS;
1036 : }
1037 :
1038 : /************************************************************************/
1039 : /* GDALRegister_GSBG() */
1040 : /************************************************************************/
1041 :
1042 2135 : void GDALRegister_GSBG()
1043 :
1044 : {
1045 2135 : if (GDALGetDriverByName("GSBG") != nullptr)
1046 263 : return;
1047 :
1048 1872 : GDALDriver *poDriver = new GDALDriver();
1049 :
1050 1872 : poDriver->SetDescription("GSBG");
1051 1872 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1052 1872 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
1053 1872 : "Golden Software Binary Grid (.grd)");
1054 1872 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/gsbg.html");
1055 1872 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "grd");
1056 1872 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
1057 1872 : "Byte Int16 UInt16 Float32");
1058 1872 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1059 :
1060 1872 : poDriver->pfnIdentify = GSBGDataset::Identify;
1061 1872 : poDriver->pfnOpen = GSBGDataset::Open;
1062 1872 : poDriver->pfnCreate = GSBGDataset::Create;
1063 1872 : poDriver->pfnCreateCopy = GSBGDataset::CreateCopy;
1064 :
1065 1872 : GetGDALDriverManager()->RegisterDriver(poDriver);
1066 : }
|
