Coverage Report

Created: 2025-06-09 07:07

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