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