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 GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
58 : int nBands, GDALDataType eType,
59 : char **papszParamList);
60 : static GDALDataset *CreateCopy(const char *pszFilename,
61 : GDALDataset *poSrcDS, int bStrict,
62 : char **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.\n");
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.\n");
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.\n");
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.\n");
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.\n");
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 60208 : int GSBGDataset::Identify(GDALOpenInfo *poOpenInfo)
470 :
471 : {
472 : /* Check for signature */
473 60208 : if (poOpenInfo->nHeaderBytes < 4 ||
474 4666 : !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "DSBB"))
475 : {
476 60118 : return FALSE;
477 : }
478 :
479 90 : return TRUE;
480 : }
481 :
482 : /************************************************************************/
483 : /* Open() */
484 : /************************************************************************/
485 :
486 45 : GDALDataset *GSBGDataset::Open(GDALOpenInfo *poOpenInfo)
487 :
488 : {
489 45 : if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
490 : {
491 0 : return nullptr;
492 : }
493 :
494 : /* -------------------------------------------------------------------- */
495 : /* Create a corresponding GDALDataset. */
496 : /* -------------------------------------------------------------------- */
497 90 : auto poDS = std::make_unique<GSBGDataset>();
498 :
499 45 : poDS->eAccess = poOpenInfo->eAccess;
500 45 : poDS->fp = poOpenInfo->fpL;
501 45 : poOpenInfo->fpL = nullptr;
502 :
503 : /* -------------------------------------------------------------------- */
504 : /* Read the header. */
505 : /* -------------------------------------------------------------------- */
506 45 : if (VSIFSeekL(poDS->fp, 4, SEEK_SET) != 0)
507 : {
508 0 : CPLError(CE_Failure, CPLE_FileIO,
509 : "Unable to seek to start of grid file header.\n");
510 0 : return nullptr;
511 : }
512 :
513 : /* Parse number of X axis grid rows */
514 : GInt16 nTemp;
515 45 : if (VSIFReadL((void *)&nTemp, 2, 1, poDS->fp) != 1)
516 : {
517 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to read raster X size.\n");
518 0 : return nullptr;
519 : }
520 45 : poDS->nRasterXSize = CPL_LSBWORD16(nTemp);
521 :
522 45 : if (VSIFReadL((void *)&nTemp, 2, 1, poDS->fp) != 1)
523 : {
524 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to read raster Y size.\n");
525 0 : return nullptr;
526 : }
527 45 : poDS->nRasterYSize = CPL_LSBWORD16(nTemp);
528 :
529 45 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
530 : {
531 0 : return nullptr;
532 : }
533 :
534 : /* -------------------------------------------------------------------- */
535 : /* Create band information objects. */
536 : /* -------------------------------------------------------------------- */
537 45 : GSBGRasterBand *poBand = new GSBGRasterBand(poDS.get(), 1);
538 45 : poDS->SetBand(1, poBand);
539 :
540 : double dfTemp;
541 45 : if (VSIFReadL((void *)&dfTemp, 8, 1, poDS->fp) != 1)
542 : {
543 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to read minimum X value.\n");
544 0 : return nullptr;
545 : }
546 45 : CPL_LSBPTR64(&dfTemp);
547 45 : poBand->dfMinX = dfTemp;
548 :
549 45 : if (VSIFReadL((void *)&dfTemp, 8, 1, poDS->fp) != 1)
550 : {
551 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to read maximum X value.\n");
552 0 : return nullptr;
553 : }
554 45 : CPL_LSBPTR64(&dfTemp);
555 45 : poBand->dfMaxX = dfTemp;
556 :
557 45 : if (VSIFReadL((void *)&dfTemp, 8, 1, poDS->fp) != 1)
558 : {
559 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to read minimum Y value.\n");
560 0 : return nullptr;
561 : }
562 45 : CPL_LSBPTR64(&dfTemp);
563 45 : poBand->dfMinY = dfTemp;
564 :
565 45 : if (VSIFReadL((void *)&dfTemp, 8, 1, poDS->fp) != 1)
566 : {
567 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to read maximum Y value.\n");
568 0 : return nullptr;
569 : }
570 45 : CPL_LSBPTR64(&dfTemp);
571 45 : poBand->dfMaxY = dfTemp;
572 :
573 45 : if (VSIFReadL((void *)&dfTemp, 8, 1, poDS->fp) != 1)
574 : {
575 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to read minimum Z value.\n");
576 0 : return nullptr;
577 : }
578 45 : CPL_LSBPTR64(&dfTemp);
579 45 : poBand->dfMinZ = dfTemp;
580 :
581 45 : if (VSIFReadL((void *)&dfTemp, 8, 1, poDS->fp) != 1)
582 : {
583 0 : CPLError(CE_Failure, CPLE_FileIO, "Unable to read maximum Z value.\n");
584 0 : return nullptr;
585 : }
586 45 : CPL_LSBPTR64(&dfTemp);
587 45 : poBand->dfMaxZ = dfTemp;
588 :
589 : /* -------------------------------------------------------------------- */
590 : /* Initialize any PAM information. */
591 : /* -------------------------------------------------------------------- */
592 45 : poDS->SetDescription(poOpenInfo->pszFilename);
593 45 : poDS->TryLoadXML();
594 :
595 : /* -------------------------------------------------------------------- */
596 : /* Check for external overviews. */
597 : /* -------------------------------------------------------------------- */
598 90 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename,
599 45 : poOpenInfo->GetSiblingFiles());
600 :
601 45 : return poDS.release();
602 : }
603 :
604 : /************************************************************************/
605 : /* GetGeoTransform() */
606 : /************************************************************************/
607 :
608 29 : CPLErr GSBGDataset::GetGeoTransform(GDALGeoTransform >) const
609 : {
610 : const GSBGRasterBand *poGRB =
611 29 : cpl::down_cast<const GSBGRasterBand *>(GetRasterBand(1));
612 :
613 : /* check if we have a PAM GeoTransform stored */
614 29 : CPLPushErrorHandler(CPLQuietErrorHandler);
615 29 : CPLErr eErr = GDALPamDataset::GetGeoTransform(gt);
616 29 : CPLPopErrorHandler();
617 :
618 29 : if (eErr == CE_None)
619 0 : return CE_None;
620 :
621 29 : if (nRasterXSize == 1 || nRasterYSize == 1)
622 0 : return CE_Failure;
623 :
624 : /* calculate pixel size first */
625 29 : gt[1] = (poGRB->dfMaxX - poGRB->dfMinX) / (nRasterXSize - 1);
626 29 : gt[5] = (poGRB->dfMinY - poGRB->dfMaxY) / (nRasterYSize - 1);
627 :
628 : /* then calculate image origin */
629 29 : gt[0] = poGRB->dfMinX - gt[1] / 2;
630 29 : gt[3] = poGRB->dfMaxY - gt[5] / 2;
631 :
632 : /* tilt/rotation does not supported by the GS grids */
633 29 : gt[4] = 0.0;
634 29 : gt[2] = 0.0;
635 :
636 29 : return CE_None;
637 : }
638 :
639 : /************************************************************************/
640 : /* SetGeoTransform() */
641 : /************************************************************************/
642 :
643 12 : CPLErr GSBGDataset::SetGeoTransform(const GDALGeoTransform >)
644 : {
645 12 : if (eAccess == GA_ReadOnly)
646 : {
647 0 : CPLError(CE_Failure, CPLE_NoWriteAccess,
648 : "Unable to set GeoTransform, dataset opened read only.\n");
649 0 : return CE_Failure;
650 : }
651 :
652 12 : GSBGRasterBand *poGRB = cpl::down_cast<GSBGRasterBand *>(GetRasterBand(1));
653 :
654 : /* non-zero transform 2 or 4 or negative 1 or 5 not supported natively */
655 : // CPLErr eErr = CE_None;
656 : /*if( gt[2] != 0.0 || gt[4] != 0.0
657 : || gt[1] < 0.0 || gt[5] < 0.0 )
658 : eErr = GDALPamDataset::SetGeoTransform( gt );
659 :
660 : if( eErr != CE_None )
661 : return eErr;*/
662 :
663 12 : double dfMinX = gt[0] + gt[1] / 2;
664 12 : double dfMaxX = gt[1] * (nRasterXSize - 0.5) + gt[0];
665 12 : double dfMinY = gt[5] * (nRasterYSize - 0.5) + gt[3];
666 12 : double dfMaxY = gt[3] + gt[5] / 2;
667 :
668 : CPLErr eErr =
669 12 : WriteHeader(fp, poGRB->nRasterXSize, poGRB->nRasterYSize, dfMinX,
670 : dfMaxX, dfMinY, dfMaxY, poGRB->dfMinZ, poGRB->dfMaxZ);
671 :
672 12 : if (eErr == CE_None)
673 : {
674 12 : poGRB->dfMinX = dfMinX;
675 12 : poGRB->dfMaxX = dfMaxX;
676 12 : poGRB->dfMinY = dfMinY;
677 12 : poGRB->dfMaxY = dfMaxY;
678 : }
679 :
680 12 : return eErr;
681 : }
682 :
683 : /************************************************************************/
684 : /* WriteHeader() */
685 : /************************************************************************/
686 :
687 65 : CPLErr GSBGDataset::WriteHeader(VSILFILE *fp, int nXSize, int nYSize,
688 : double dfMinX, double dfMaxX, double dfMinY,
689 : double dfMaxY, double dfMinZ, double dfMaxZ)
690 :
691 : {
692 65 : if (VSIFSeekL(fp, 0, SEEK_SET) != 0)
693 : {
694 0 : CPLError(CE_Failure, CPLE_FileIO,
695 : "Unable to seek to start of grid file.\n");
696 0 : return CE_Failure;
697 : }
698 :
699 65 : if (VSIFWriteL((void *)"DSBB", 1, 4, fp) != 4)
700 : {
701 1 : CPLError(CE_Failure, CPLE_FileIO,
702 : "Unable to write signature to grid file.\n");
703 1 : return CE_Failure;
704 : }
705 :
706 64 : assert(nXSize >= 0 && nXSize <= std::numeric_limits<int16_t>::max());
707 64 : GInt16 nTemp = CPL_LSBWORD16(static_cast<int16_t>(nXSize));
708 64 : if (VSIFWriteL((void *)&nTemp, 2, 1, fp) != 1)
709 : {
710 0 : CPLError(CE_Failure, CPLE_FileIO,
711 : "Unable to write raster X size to grid file.\n");
712 0 : return CE_Failure;
713 : }
714 :
715 64 : assert(nYSize >= 0 && nYSize <= std::numeric_limits<int16_t>::max());
716 64 : nTemp = CPL_LSBWORD16(static_cast<int16_t>(nYSize));
717 64 : if (VSIFWriteL((void *)&nTemp, 2, 1, fp) != 1)
718 : {
719 0 : CPLError(CE_Failure, CPLE_FileIO,
720 : "Unable to write raster Y size to grid file.\n");
721 0 : return CE_Failure;
722 : }
723 :
724 64 : double dfTemp = dfMinX;
725 64 : CPL_LSBPTR64(&dfTemp);
726 64 : if (VSIFWriteL((void *)&dfTemp, 8, 1, fp) != 1)
727 : {
728 0 : CPLError(CE_Failure, CPLE_FileIO,
729 : "Unable to write minimum X value to grid file.\n");
730 0 : return CE_Failure;
731 : }
732 :
733 64 : dfTemp = dfMaxX;
734 64 : CPL_LSBPTR64(&dfTemp);
735 64 : if (VSIFWriteL((void *)&dfTemp, 8, 1, fp) != 1)
736 : {
737 0 : CPLError(CE_Failure, CPLE_FileIO,
738 : "Unable to write maximum X value to grid file.\n");
739 0 : return CE_Failure;
740 : }
741 :
742 64 : dfTemp = dfMinY;
743 64 : CPL_LSBPTR64(&dfTemp);
744 64 : if (VSIFWriteL((void *)&dfTemp, 8, 1, fp) != 1)
745 : {
746 0 : CPLError(CE_Failure, CPLE_FileIO,
747 : "Unable to write minimum Y value to grid file.\n");
748 0 : return CE_Failure;
749 : }
750 :
751 64 : dfTemp = dfMaxY;
752 64 : CPL_LSBPTR64(&dfTemp);
753 64 : if (VSIFWriteL((void *)&dfTemp, 8, 1, fp) != 1)
754 : {
755 0 : CPLError(CE_Failure, CPLE_FileIO,
756 : "Unable to write maximum Y value to grid file.\n");
757 0 : return CE_Failure;
758 : }
759 :
760 64 : dfTemp = dfMinZ;
761 64 : CPL_LSBPTR64(&dfTemp);
762 64 : if (VSIFWriteL((void *)&dfTemp, 8, 1, fp) != 1)
763 : {
764 1 : CPLError(CE_Failure, CPLE_FileIO,
765 : "Unable to write minimum Z value to grid file.\n");
766 1 : return CE_Failure;
767 : }
768 :
769 63 : dfTemp = dfMaxZ;
770 63 : CPL_LSBPTR64(&dfTemp);
771 63 : if (VSIFWriteL((void *)&dfTemp, 8, 1, fp) != 1)
772 : {
773 0 : CPLError(CE_Failure, CPLE_FileIO,
774 : "Unable to write maximum Z value to grid file.\n");
775 0 : return CE_Failure;
776 : }
777 :
778 63 : return CE_None;
779 : }
780 :
781 : /************************************************************************/
782 : /* GSBGCreateCheckDims() */
783 : /************************************************************************/
784 :
785 66 : static bool GSBGCreateCheckDims(int nXSize, int nYSize)
786 : {
787 66 : if (nXSize <= 1 || nYSize <= 1)
788 : {
789 4 : CPLError(CE_Failure, CPLE_IllegalArg,
790 : "Unable to create grid, both X and Y size must be "
791 : "larger or equal to 2.");
792 4 : return false;
793 : }
794 122 : if (nXSize > std::numeric_limits<short>::max() ||
795 60 : nYSize > std::numeric_limits<short>::max())
796 : {
797 4 : CPLError(CE_Failure, CPLE_IllegalArg,
798 : "Unable to create grid, Golden Software Binary Grid format "
799 : "only supports sizes up to %dx%d. %dx%d not supported.",
800 4 : std::numeric_limits<short>::max(),
801 4 : std::numeric_limits<short>::max(), nXSize, nYSize);
802 4 : return false;
803 : }
804 58 : return true;
805 : }
806 :
807 : /************************************************************************/
808 : /* Create() */
809 : /************************************************************************/
810 :
811 37 : GDALDataset *GSBGDataset::Create(const char *pszFilename, int nXSize,
812 : int nYSize, int /* nBands */,
813 : GDALDataType eType,
814 : CPL_UNUSED char **papszParamList)
815 : {
816 37 : if (!GSBGCreateCheckDims(nXSize, nYSize))
817 : {
818 4 : return nullptr;
819 : }
820 :
821 33 : if (eType != GDT_Byte && eType != GDT_Float32 && eType != GDT_UInt16 &&
822 : eType != GDT_Int16)
823 : {
824 20 : CPLError(CE_Failure, CPLE_AppDefined,
825 : "Golden Software Binary Grid only supports Byte, Int16, "
826 : "Uint16, and Float32 datatypes. Unable to create with "
827 : "type %s.\n",
828 : GDALGetDataTypeName(eType));
829 :
830 20 : return nullptr;
831 : }
832 :
833 13 : VSILFILE *fp = VSIFOpenL(pszFilename, "w+b");
834 :
835 13 : if (fp == nullptr)
836 : {
837 0 : CPLError(CE_Failure, CPLE_OpenFailed,
838 : "Attempt to create file '%s' failed.\n", pszFilename);
839 0 : return nullptr;
840 : }
841 :
842 : CPLErr eErr =
843 13 : WriteHeader(fp, nXSize, nYSize, 0.0, nXSize, 0.0, nYSize, 0.0, 0.0);
844 13 : if (eErr != CE_None)
845 : {
846 0 : VSIFCloseL(fp);
847 0 : return nullptr;
848 : }
849 :
850 13 : float fVal = fNODATA_VALUE;
851 13 : CPL_LSBPTR32(&fVal);
852 1233 : for (int iRow = 0; iRow < nYSize; iRow++)
853 : {
854 121620 : for (int iCol = 0; iCol < nXSize; iCol++)
855 : {
856 120400 : if (VSIFWriteL((void *)&fVal, 4, 1, fp) != 1)
857 : {
858 0 : VSIFCloseL(fp);
859 0 : CPLError(CE_Failure, CPLE_FileIO,
860 : "Unable to write grid cell. Disk full?\n");
861 0 : return nullptr;
862 : }
863 : }
864 : }
865 :
866 13 : VSIFCloseL(fp);
867 :
868 13 : return (GDALDataset *)GDALOpen(pszFilename, GA_Update);
869 : }
870 :
871 : /************************************************************************/
872 : /* CreateCopy() */
873 : /************************************************************************/
874 :
875 34 : GDALDataset *GSBGDataset::CreateCopy(const char *pszFilename,
876 : GDALDataset *poSrcDS, int bStrict,
877 : CPL_UNUSED char **papszOptions,
878 : GDALProgressFunc pfnProgress,
879 : void *pProgressData)
880 : {
881 34 : if (pfnProgress == nullptr)
882 0 : pfnProgress = GDALDummyProgress;
883 :
884 34 : int nBands = poSrcDS->GetRasterCount();
885 34 : if (nBands == 0)
886 : {
887 1 : CPLError(
888 : CE_Failure, CPLE_NotSupported,
889 : "GSBG driver does not support source dataset with zero band.\n");
890 1 : return nullptr;
891 : }
892 33 : else if (nBands > 1)
893 : {
894 4 : if (bStrict)
895 : {
896 4 : CPLError(CE_Failure, CPLE_NotSupported,
897 : "Unable to create copy, Golden Software Binary Grid "
898 : "format only supports one raster band.\n");
899 4 : return nullptr;
900 : }
901 : else
902 0 : CPLError(CE_Warning, CPLE_NotSupported,
903 : "Golden Software Binary Grid format only supports one "
904 : "raster band, first band will be copied.\n");
905 : }
906 :
907 29 : GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(1);
908 29 : if (!GSBGCreateCheckDims(poSrcBand->GetXSize(), poSrcBand->GetYSize()))
909 : {
910 4 : return nullptr;
911 : }
912 :
913 25 : if (!pfnProgress(0.0, nullptr, pProgressData))
914 : {
915 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated\n");
916 0 : return nullptr;
917 : }
918 :
919 25 : VSILFILE *fp = VSIFOpenL(pszFilename, "w+b");
920 :
921 25 : if (fp == nullptr)
922 : {
923 3 : CPLError(CE_Failure, CPLE_OpenFailed,
924 : "Attempt to create file '%s' failed.\n", pszFilename);
925 3 : return nullptr;
926 : }
927 :
928 22 : const int nXSize = poSrcBand->GetXSize();
929 22 : const int nYSize = poSrcBand->GetYSize();
930 22 : GDALGeoTransform gt;
931 22 : poSrcDS->GetGeoTransform(gt);
932 :
933 22 : double dfMinX = gt[0] + gt[1] / 2;
934 22 : double dfMaxX = gt[1] * (nXSize - 0.5) + gt[0];
935 22 : double dfMinY = gt[5] * (nYSize - 0.5) + gt[3];
936 22 : double dfMaxY = gt[3] + gt[5] / 2;
937 22 : CPLErr eErr = WriteHeader(fp, nXSize, nYSize, dfMinX, dfMaxX, dfMinY,
938 : dfMaxY, 0.0, 0.0);
939 :
940 22 : if (eErr != CE_None)
941 : {
942 2 : VSIFCloseL(fp);
943 2 : return nullptr;
944 : }
945 :
946 : /* -------------------------------------------------------------------- */
947 : /* Copy band data. */
948 : /* -------------------------------------------------------------------- */
949 20 : float *pfData = (float *)VSI_MALLOC2_VERBOSE(nXSize, sizeof(float));
950 20 : if (pfData == nullptr)
951 : {
952 0 : VSIFCloseL(fp);
953 0 : return nullptr;
954 : }
955 :
956 : int bSrcHasNDValue;
957 20 : float fSrcNoDataValue = (float)poSrcBand->GetNoDataValue(&bSrcHasNDValue);
958 20 : double dfMinZ = std::numeric_limits<double>::max();
959 20 : double dfMaxZ = std::numeric_limits<double>::lowest();
960 185 : for (int iRow = nYSize - 1; iRow >= 0; iRow--)
961 : {
962 173 : eErr = poSrcBand->RasterIO(GF_Read, 0, iRow, nXSize, 1, pfData, nXSize,
963 : 1, GDT_Float32, 0, 0, nullptr);
964 :
965 173 : if (eErr != CE_None)
966 : {
967 0 : VSIFCloseL(fp);
968 0 : VSIFree(pfData);
969 0 : return nullptr;
970 : }
971 :
972 2103 : for (int iCol = 0; iCol < nXSize; iCol++)
973 : {
974 1930 : if (bSrcHasNDValue && pfData[iCol] == fSrcNoDataValue)
975 : {
976 0 : pfData[iCol] = fNODATA_VALUE;
977 : }
978 : else
979 : {
980 1930 : if (pfData[iCol] > dfMaxZ)
981 22 : dfMaxZ = pfData[iCol];
982 :
983 1930 : if (pfData[iCol] < dfMinZ)
984 27 : dfMinZ = pfData[iCol];
985 : }
986 :
987 1930 : CPL_LSBPTR32(pfData + iCol);
988 : }
989 :
990 173 : if (VSIFWriteL((void *)pfData, 4, nXSize, fp) !=
991 173 : static_cast<unsigned>(nXSize))
992 : {
993 8 : VSIFCloseL(fp);
994 8 : VSIFree(pfData);
995 8 : CPLError(CE_Failure, CPLE_FileIO,
996 : "Unable to write grid row. Disk full?\n");
997 8 : return nullptr;
998 : }
999 :
1000 165 : if (!pfnProgress(static_cast<double>(nYSize - iRow) / nYSize, nullptr,
1001 : pProgressData))
1002 : {
1003 0 : VSIFCloseL(fp);
1004 0 : VSIFree(pfData);
1005 0 : CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
1006 0 : return nullptr;
1007 : }
1008 : }
1009 :
1010 12 : VSIFree(pfData);
1011 :
1012 : /* write out the min and max values */
1013 12 : eErr = WriteHeader(fp, nXSize, nYSize, dfMinX, dfMaxX, dfMinY, dfMaxY,
1014 : dfMinZ, dfMaxZ);
1015 :
1016 12 : if (eErr != CE_None)
1017 : {
1018 0 : VSIFCloseL(fp);
1019 0 : return nullptr;
1020 : }
1021 :
1022 12 : VSIFCloseL(fp);
1023 :
1024 12 : GDALPamDataset *poDS = (GDALPamDataset *)GDALOpen(pszFilename, GA_Update);
1025 12 : if (poDS)
1026 : {
1027 12 : poDS->CloneInfo(poSrcDS, GCIF_PAM_DEFAULT);
1028 : }
1029 12 : return poDS;
1030 : }
1031 :
1032 : /************************************************************************/
1033 : /* GDALRegister_GSBG() */
1034 : /************************************************************************/
1035 :
1036 2033 : void GDALRegister_GSBG()
1037 :
1038 : {
1039 2033 : if (GDALGetDriverByName("GSBG") != nullptr)
1040 283 : return;
1041 :
1042 1750 : GDALDriver *poDriver = new GDALDriver();
1043 :
1044 1750 : poDriver->SetDescription("GSBG");
1045 1750 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1046 1750 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
1047 1750 : "Golden Software Binary Grid (.grd)");
1048 1750 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/gsbg.html");
1049 1750 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "grd");
1050 1750 : poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
1051 1750 : "Byte Int16 UInt16 Float32");
1052 1750 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1053 :
1054 1750 : poDriver->pfnIdentify = GSBGDataset::Identify;
1055 1750 : poDriver->pfnOpen = GSBGDataset::Open;
1056 1750 : poDriver->pfnCreate = GSBGDataset::Create;
1057 1750 : poDriver->pfnCreateCopy = GSBGDataset::CreateCopy;
1058 :
1059 1750 : GetGDALDriverManager()->RegisterDriver(poDriver);
1060 : }
|
