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