Line data    Source code

       1             : /****************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  Implements the Golden Software Surfer 7 Binary Grid Format.
       5             :  * Author:   Adam Guernsey, adam@ctech.com
       6             :  *           (Based almost entirely on gsbgdataset.cpp by Kevin Locke)
       7             :  *           Create functions added by Russell Jurgensen.
       8             :  *
       9             :  ****************************************************************************
      10             :  * Copyright (c) 2007, Adam Guernsey <adam@ctech.com>
      11             :  * Copyright (c) 2009-2011, Even Rouault <even dot rouault at spatialys.com>
      12             :  *
      13             :  * SPDX-License-Identifier: MIT
      14             :  ****************************************************************************/
      15             : 
      16             : #include <cassert>
      17             : #include <cfloat>
      18             : #include <climits>
      19             : #include <cmath>
      20             : #include <limits>
      21             : 
      22             : #include "gdal_frmts.h"
      23             : #include "gdal_pam.h"
      24             : 
      25             : /************************************************************************/
      26             : /* ==================================================================== */
      27             : /*                GS7BGDataset                */
      28             : /* ==================================================================== */
      29             : /************************************************************************/
      30             : 
      31             : class GS7BGRasterBand;
      32             : 
      33             : constexpr double dfDefaultNoDataValue = 1.701410009187828e+38f;
      34             : 
      35             : class GS7BGDataset final : public GDALPamDataset
      36             : {
      37             :     friend class GS7BGRasterBand;
      38             : 
      39             :     double dfNoData_Value;
      40             :     static const size_t nHEADER_SIZE;
      41             :     size_t nData_Position;
      42             : 
      43             :     static CPLErr WriteHeader(VSILFILE *fp, GInt32 nXSize, GInt32 nYSize,
      44             :                               double dfMinX, double dfMaxX, double dfMinY,
      45             :                               double dfMaxY, double dfMinZ, double dfMaxZ);
      46             : 
      47             :     VSILFILE *fp;
      48             : 
      49             :   public:
      50          33 :     GS7BGDataset()
      51          33 :         : /* NOTE:  This is not mentioned in the spec, but Surfer 8 uses this
      52             :              value */
      53             :           /* 0x7effffee (Little Endian: eeffff7e) */
      54          33 :           dfNoData_Value(dfDefaultNoDataValue), nData_Position(0), fp(nullptr)
      55             :     {
      56          33 :     }
      57             : 
      58             :     ~GS7BGDataset();
      59             : 
      60             :     static int Identify(GDALOpenInfo *);
      61             :     static GDALDataset *Open(GDALOpenInfo *);
      62             :     static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
      63             :                                int nBandsIn, GDALDataType eType,
      64             :                                char **papszParamList);
      65             :     static GDALDataset *CreateCopy(const char *pszFilename,
      66             :                                    GDALDataset *poSrcDS, int bStrict,
      67             :                                    char **papszOptions,
      68             :                                    GDALProgressFunc pfnProgress,
      69             :                                    void *pProgressData);
      70             : 
      71             :     CPLErr GetGeoTransform(double *padfGeoTransform) override;
      72             :     CPLErr SetGeoTransform(double *padfGeoTransform) override;
      73             : };
      74             : 
      75             : const size_t GS7BGDataset::nHEADER_SIZE = 100;
      76             : 
      77             : constexpr long nHEADER_TAG = 0x42525344;
      78             : constexpr long nGRID_TAG = 0x44495247;
      79             : constexpr long nDATA_TAG = 0x41544144;
      80             : #if 0 /* Unused */
      81             : const long  nFAULT_TAG = 0x49544c46;
      82             : #endif
      83             : 
      84             : /************************************************************************/
      85             : /* ==================================================================== */
      86             : /*                            GS7BGRasterBand                           */
      87             : /* ==================================================================== */
      88             : /************************************************************************/
      89             : 
      90             : class GS7BGRasterBand final : public GDALPamRasterBand
      91             : {
      92             :     friend class GS7BGDataset;
      93             : 
      94             :     double dfMinX;
      95             :     double dfMaxX;
      96             :     double dfMinY;
      97             :     double dfMaxY;
      98             :     double dfMinZ;
      99             :     double dfMaxZ;
     100             : 
     101             :     double *pafRowMinZ;
     102             :     double *pafRowMaxZ;
     103             :     int nMinZRow;
     104             :     int nMaxZRow;
     105             : 
     106             :     CPLErr ScanForMinMaxZ();
     107             : 
     108             :   public:
     109             :     GS7BGRasterBand(GS7BGDataset *, int);
     110             :     ~GS7BGRasterBand();
     111             : 
     112             :     CPLErr IReadBlock(int, int, void *) override;
     113             :     CPLErr IWriteBlock(int, int, void *) override;
     114             :     double GetMinimum(int *pbSuccess = nullptr) override;
     115             :     double GetMaximum(int *pbSuccess = nullptr) override;
     116             : 
     117             :     double GetNoDataValue(int *pbSuccess = nullptr) override;
     118             : };
     119             : 
     120             : /************************************************************************/
     121             : /*                           GS7BGRasterBand()                          */
     122             : /************************************************************************/
     123             : 
     124          33 : GS7BGRasterBand::GS7BGRasterBand(GS7BGDataset *poDSIn, int nBandIn)
     125             :     : dfMinX(0.0), dfMaxX(0.0), dfMinY(0.0), dfMaxY(0.0), dfMinZ(0.0),
     126             :       dfMaxZ(0.0), pafRowMinZ(nullptr), pafRowMaxZ(nullptr), nMinZRow(-1),
     127          33 :       nMaxZRow(-1)
     128             : 
     129             : {
     130          33 :     poDS = poDSIn;
     131          33 :     nBand = nBandIn;
     132             : 
     133          33 :     eDataType = GDT_Float64;
     134             : 
     135          33 :     nBlockXSize = poDS->GetRasterXSize();
     136          33 :     nBlockYSize = 1;
     137          33 : }
     138             : 
     139             : /************************************************************************/
     140             : /*                           ~GSBGRasterBand()                          */
     141             : /************************************************************************/
     142             : 
     143          66 : GS7BGRasterBand::~GS7BGRasterBand()
     144             : 
     145             : {
     146          33 :     CPLFree(pafRowMinZ);
     147          33 :     CPLFree(pafRowMaxZ);
     148          66 : }
     149             : 
     150             : /************************************************************************/
     151             : /*                          ScanForMinMaxZ()                            */
     152             : /************************************************************************/
     153             : 
     154           1 : CPLErr GS7BGRasterBand::ScanForMinMaxZ()
     155             : 
     156             : {
     157           1 :     GS7BGDataset *poGDS = reinterpret_cast<GS7BGDataset *>(poDS);
     158             :     double *pafRowVals =
     159           1 :         (double *)VSI_MALLOC2_VERBOSE(nRasterXSize, sizeof(double));
     160             : 
     161           1 :     if (pafRowVals == nullptr)
     162             :     {
     163           0 :         return CE_Failure;
     164             :     }
     165             : 
     166           1 :     double dfNewMinZ = std::numeric_limits<double>::max();
     167           1 :     double dfNewMaxZ = std::numeric_limits<double>::lowest();
     168           1 :     int nNewMinZRow = 0;
     169           1 :     int nNewMaxZRow = 0;
     170             : 
     171             :     /* Since we have to scan, lets calc. statistics too */
     172           1 :     double dfSum = 0.0;
     173           1 :     double dfSum2 = 0.0;
     174           1 :     unsigned long nValuesRead = 0;
     175          21 :     for (int iRow = 0; iRow < nRasterYSize; iRow++)
     176             :     {
     177          20 :         CPLErr eErr = IReadBlock(0, iRow, pafRowVals);
     178          20 :         if (eErr != CE_None)
     179             :         {
     180           0 :             VSIFree(pafRowVals);
     181           0 :             return CE_Failure;
     182             :         }
     183             : 
     184          20 :         pafRowMinZ[iRow] = std::numeric_limits<float>::max();
     185          20 :         pafRowMaxZ[iRow] = std::numeric_limits<float>::lowest();
     186         420 :         for (int iCol = 0; iCol < nRasterXSize; iCol++)
     187             :         {
     188         400 :             if (pafRowVals[iCol] == poGDS->dfNoData_Value)
     189         400 :                 continue;
     190             : 
     191           0 :             if (pafRowVals[iCol] < pafRowMinZ[iRow])
     192           0 :                 pafRowMinZ[iRow] = pafRowVals[iCol];
     193             : 
     194           0 :             if (pafRowVals[iCol] > pafRowMinZ[iRow])
     195           0 :                 pafRowMaxZ[iRow] = pafRowVals[iCol];
     196             : 
     197           0 :             dfSum += pafRowVals[iCol];
     198           0 :             dfSum2 += pafRowVals[iCol] * pafRowVals[iCol];
     199           0 :             nValuesRead++;
     200             :         }
     201             : 
     202          20 :         if (pafRowMinZ[iRow] < dfNewMinZ)
     203             :         {
     204           1 :             dfNewMinZ = pafRowMinZ[iRow];
     205           1 :             nNewMinZRow = iRow;
     206             :         }
     207             : 
     208          20 :         if (pafRowMaxZ[iRow] > dfNewMaxZ)
     209             :         {
     210           1 :             dfNewMaxZ = pafRowMaxZ[iRow];
     211           1 :             nNewMaxZRow = iRow;
     212             :         }
     213             :     }
     214             : 
     215           1 :     VSIFree(pafRowVals);
     216             : 
     217           1 :     if (nValuesRead == 0)
     218             :     {
     219           1 :         dfMinZ = 0.0;
     220           1 :         dfMaxZ = 0.0;
     221           1 :         nMinZRow = 0;
     222           1 :         nMaxZRow = 0;
     223           1 :         return CE_None;
     224             :     }
     225             : 
     226           0 :     dfMinZ = dfNewMinZ;
     227           0 :     dfMaxZ = dfNewMaxZ;
     228           0 :     nMinZRow = nNewMinZRow;
     229           0 :     nMaxZRow = nNewMaxZRow;
     230             : 
     231           0 :     double dfMean = dfSum / nValuesRead;
     232           0 :     double dfStdDev = sqrt((dfSum2 / nValuesRead) - (dfMean * dfMean));
     233           0 :     SetStatistics(dfMinZ, dfMaxZ, dfMean, dfStdDev);
     234             : 
     235           0 :     return CE_None;
     236             : }
     237             : 
     238             : /************************************************************************/
     239             : /*                             IReadBlock()                             */
     240             : /************************************************************************/
     241             : 
     242         140 : CPLErr GS7BGRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage)
     243             : 
     244             : {
     245         140 :     if (nBlockYOff < 0 || nBlockYOff > nRasterYSize - 1 || nBlockXOff != 0)
     246           0 :         return CE_Failure;
     247             : 
     248         140 :     GS7BGDataset *poGDS = (GS7BGDataset *)(poDS);
     249             : 
     250         280 :     if (VSIFSeekL(poGDS->fp,
     251         140 :                   (poGDS->nData_Position +
     252         140 :                    sizeof(double) * static_cast<vsi_l_offset>(nRasterXSize) *
     253         140 :                        (nRasterYSize - nBlockYOff - 1)),
     254         140 :                   SEEK_SET) != 0)
     255             :     {
     256           0 :         CPLError(CE_Failure, CPLE_FileIO,
     257             :                  "Unable to seek to beginning of grid row.\n");
     258           0 :         return CE_Failure;
     259             :     }
     260             : 
     261         140 :     if (VSIFReadL(pImage, sizeof(double), nBlockXSize, poGDS->fp) !=
     262         140 :         static_cast<unsigned>(nBlockXSize))
     263             :     {
     264           0 :         CPLError(CE_Failure, CPLE_FileIO,
     265             :                  "Unable to read block from grid file.\n");
     266           0 :         return CE_Failure;
     267             :     }
     268             : 
     269             : #ifdef CPL_MSB
     270             :     double *pfImage = (double *)pImage;
     271             :     for (int iPixel = 0; iPixel < nBlockXSize; iPixel++)
     272             :         CPL_LSBPTR64(pfImage + iPixel);
     273             : #endif
     274             : 
     275         140 :     return CE_None;
     276             : }
     277             : 
     278             : /************************************************************************/
     279             : /*                            IWriteBlock()                             */
     280             : /************************************************************************/
     281             : 
     282          20 : CPLErr GS7BGRasterBand::IWriteBlock(int nBlockXOff, int nBlockYOff,
     283             :                                     void *pImage)
     284             : 
     285             : {
     286          20 :     if (eAccess == GA_ReadOnly)
     287             :     {
     288           0 :         CPLError(CE_Failure, CPLE_NoWriteAccess,
     289             :                  "Unable to write block, dataset opened read only.\n");
     290           0 :         return CE_Failure;
     291             :     }
     292             : 
     293          20 :     if (nBlockYOff < 0 || nBlockYOff > nRasterYSize - 1 || nBlockXOff != 0)
     294           0 :         return CE_Failure;
     295             : 
     296          20 :     GS7BGDataset *poGDS = (GS7BGDataset *)(poDS);
     297             : 
     298          20 :     if (pafRowMinZ == nullptr || pafRowMaxZ == nullptr || nMinZRow < 0 ||
     299          19 :         nMaxZRow < 0)
     300             :     {
     301           1 :         pafRowMinZ =
     302           1 :             (double *)VSI_MALLOC2_VERBOSE(nRasterYSize, sizeof(double));
     303           1 :         if (pafRowMinZ == nullptr)
     304             :         {
     305           0 :             return CE_Failure;
     306             :         }
     307             : 
     308           1 :         pafRowMaxZ =
     309           1 :             (double *)VSI_MALLOC2_VERBOSE(nRasterYSize, sizeof(double));
     310           1 :         if (pafRowMaxZ == nullptr)
     311             :         {
     312           0 :             VSIFree(pafRowMinZ);
     313           0 :             pafRowMinZ = nullptr;
     314           0 :             return CE_Failure;
     315             :         }
     316             : 
     317           1 :         CPLErr eErr = ScanForMinMaxZ();
     318           1 :         if (eErr != CE_None)
     319           0 :             return eErr;
     320             :     }
     321             : 
     322          40 :     if (VSIFSeekL(poGDS->fp,
     323          20 :                   GS7BGDataset::nHEADER_SIZE +
     324          20 :                       sizeof(double) * nRasterXSize *
     325          20 :                           (nRasterYSize - nBlockYOff - 1),
     326          20 :                   SEEK_SET) != 0)
     327             :     {
     328           0 :         CPLError(CE_Failure, CPLE_FileIO,
     329             :                  "Unable to seek to beginning of grid row.\n");
     330           0 :         return CE_Failure;
     331             :     }
     332             : 
     333          20 :     double *pdfImage = (double *)pImage;
     334          20 :     pafRowMinZ[nBlockYOff] = std::numeric_limits<double>::max();
     335          20 :     pafRowMaxZ[nBlockYOff] = std::numeric_limits<double>::lowest();
     336         420 :     for (int iPixel = 0; iPixel < nBlockXSize; iPixel++)
     337             :     {
     338         400 :         if (pdfImage[iPixel] != poGDS->dfNoData_Value)
     339             :         {
     340         400 :             if (pdfImage[iPixel] < pafRowMinZ[nBlockYOff])
     341          78 :                 pafRowMinZ[nBlockYOff] = pdfImage[iPixel];
     342             : 
     343         400 :             if (pdfImage[iPixel] > pafRowMaxZ[nBlockYOff])
     344          43 :                 pafRowMaxZ[nBlockYOff] = pdfImage[iPixel];
     345             :         }
     346             : 
     347         400 :         CPL_LSBPTR64(pdfImage + iPixel);
     348             :     }
     349             : 
     350          20 :     if (VSIFWriteL(pImage, sizeof(double), nBlockXSize, poGDS->fp) !=
     351          20 :         static_cast<unsigned>(nBlockXSize))
     352             :     {
     353           0 :         CPLError(CE_Failure, CPLE_FileIO,
     354             :                  "Unable to write block to grid file.\n");
     355           0 :         return CE_Failure;
     356             :     }
     357             : 
     358             :     /* Update min/max Z values as appropriate */
     359          20 :     bool bHeaderNeedsUpdate = false;
     360          20 :     if (nMinZRow == nBlockYOff && pafRowMinZ[nBlockYOff] > dfMinZ)
     361             :     {
     362           1 :         double dfNewMinZ = std::numeric_limits<double>::max();
     363          21 :         for (int iRow = 0; iRow < nRasterYSize; iRow++)
     364             :         {
     365          20 :             if (pafRowMinZ[iRow] < dfNewMinZ)
     366             :             {
     367           1 :                 dfNewMinZ = pafRowMinZ[iRow];
     368           1 :                 nMinZRow = iRow;
     369             :             }
     370             :         }
     371             : 
     372           1 :         if (dfNewMinZ != dfMinZ)
     373             :         {
     374           1 :             dfMinZ = dfNewMinZ;
     375           1 :             bHeaderNeedsUpdate = true;
     376             :         }
     377             :     }
     378             : 
     379          20 :     if (nMaxZRow == nBlockYOff && pafRowMaxZ[nBlockYOff] < dfMaxZ)
     380             :     {
     381           0 :         double dfNewMaxZ = std::numeric_limits<double>::lowest();
     382           0 :         for (int iRow = 0; iRow < nRasterYSize; iRow++)
     383             :         {
     384           0 :             if (pafRowMaxZ[iRow] > dfNewMaxZ)
     385             :             {
     386           0 :                 dfNewMaxZ = pafRowMaxZ[iRow];
     387           0 :                 nMaxZRow = iRow;
     388             :             }
     389             :         }
     390             : 
     391           0 :         if (dfNewMaxZ != dfMaxZ)
     392             :         {
     393           0 :             dfMaxZ = dfNewMaxZ;
     394           0 :             bHeaderNeedsUpdate = true;
     395             :         }
     396             :     }
     397             : 
     398          20 :     if (pafRowMinZ[nBlockYOff] < dfMinZ || pafRowMaxZ[nBlockYOff] > dfMaxZ)
     399             :     {
     400           6 :         if (pafRowMinZ[nBlockYOff] < dfMinZ)
     401             :         {
     402           3 :             dfMinZ = pafRowMinZ[nBlockYOff];
     403           3 :             nMinZRow = nBlockYOff;
     404             :         }
     405             : 
     406           6 :         if (pafRowMaxZ[nBlockYOff] > dfMaxZ)
     407             :         {
     408           4 :             dfMaxZ = pafRowMaxZ[nBlockYOff];
     409           4 :             nMaxZRow = nBlockYOff;
     410             :         }
     411             : 
     412           6 :         bHeaderNeedsUpdate = true;
     413             :     }
     414             : 
     415          20 :     if (bHeaderNeedsUpdate && dfMaxZ > dfMinZ)
     416             :     {
     417             :         CPLErr eErr =
     418           6 :             poGDS->WriteHeader(poGDS->fp, nRasterXSize, nRasterYSize, dfMinX,
     419             :                                dfMaxX, dfMinY, dfMaxY, dfMinZ, dfMaxZ);
     420           6 :         return eErr;
     421             :     }
     422             : 
     423          14 :     return CE_None;
     424             : }
     425             : 
     426             : /************************************************************************/
     427             : /*                           GetNoDataValue()                           */
     428             : /************************************************************************/
     429             : 
     430          11 : double GS7BGRasterBand::GetNoDataValue(int *pbSuccess)
     431             : {
     432          11 :     GS7BGDataset *poGDS = reinterpret_cast<GS7BGDataset *>(poDS);
     433          11 :     if (pbSuccess)
     434          10 :         *pbSuccess = TRUE;
     435             : 
     436          11 :     return poGDS->dfNoData_Value;
     437             : }
     438             : 
     439             : /************************************************************************/
     440             : /*                             GetMinimum()                             */
     441             : /************************************************************************/
     442             : 
     443           0 : double GS7BGRasterBand::GetMinimum(int *pbSuccess)
     444             : {
     445           0 :     if (pbSuccess)
     446           0 :         *pbSuccess = TRUE;
     447             : 
     448           0 :     return dfMinZ;
     449             : }
     450             : 
     451             : /************************************************************************/
     452             : /*                             GetMaximum()                             */
     453             : /************************************************************************/
     454             : 
     455           0 : double GS7BGRasterBand::GetMaximum(int *pbSuccess)
     456             : {
     457           0 :     if (pbSuccess)
     458           0 :         *pbSuccess = TRUE;
     459             : 
     460           0 :     return dfMaxZ;
     461             : }
     462             : 
     463             : /************************************************************************/
     464             : /* ==================================================================== */
     465             : /*                            GS7BGDataset                              */
     466             : /* ==================================================================== */
     467             : /************************************************************************/
     468             : 
     469          66 : GS7BGDataset::~GS7BGDataset()
     470             : 
     471             : {
     472          33 :     FlushCache(true);
     473          33 :     if (fp != nullptr)
     474          33 :         VSIFCloseL(fp);
     475          66 : }
     476             : 
     477             : /************************************************************************/
     478             : /*                              Identify()                              */
     479             : /************************************************************************/
     480             : 
     481       52914 : int GS7BGDataset::Identify(GDALOpenInfo *poOpenInfo)
     482             : 
     483             : {
     484             :     /* Check for signature - for GS7BG the signature is the */
     485             :     /* nHEADER_TAG with reverse byte order.                 */
     486       52914 :     if (poOpenInfo->nHeaderBytes < 4 ||
     487        4516 :         !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "DSRB"))
     488             :     {
     489       52848 :         return FALSE;
     490             :     }
     491             : 
     492          66 :     return TRUE;
     493             : }
     494             : 
     495             : /************************************************************************/
     496             : /*                                Open()                                */
     497             : /************************************************************************/
     498             : 
     499          33 : GDALDataset *GS7BGDataset::Open(GDALOpenInfo *poOpenInfo)
     500             : 
     501             : {
     502          33 :     if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
     503             :     {
     504           0 :         return nullptr;
     505             :     }
     506             : 
     507             :     /* ------------------------------------------------------------------- */
     508             :     /*      Create a corresponding GDALDataset.                            */
     509             :     /* ------------------------------------------------------------------- */
     510          33 :     GS7BGDataset *poDS = new GS7BGDataset();
     511          33 :     poDS->eAccess = poOpenInfo->eAccess;
     512          33 :     poDS->fp = poOpenInfo->fpL;
     513          33 :     poOpenInfo->fpL = nullptr;
     514             : 
     515             :     /* ------------------------------------------------------------------- */
     516             :     /*      Read the header. The Header section must be the first section  */
     517             :     /*      in the file.                                                   */
     518             :     /* ------------------------------------------------------------------- */
     519          33 :     if (VSIFSeekL(poDS->fp, 0, SEEK_SET) != 0)
     520             :     {
     521           0 :         delete poDS;
     522           0 :         CPLError(CE_Failure, CPLE_FileIO,
     523             :                  "Unable to seek to start of grid file header.\n");
     524           0 :         return nullptr;
     525             :     }
     526             : 
     527             :     GInt32 nTag;
     528          33 :     if (VSIFReadL((void *)&nTag, sizeof(GInt32), 1, poDS->fp) != 1)
     529             :     {
     530           0 :         delete poDS;
     531           0 :         CPLError(CE_Failure, CPLE_FileIO, "Unable to read Tag.\n");
     532           0 :         return nullptr;
     533             :     }
     534             : 
     535          33 :     CPL_LSBPTR32(&nTag);
     536             : 
     537          33 :     if (nTag != nHEADER_TAG)
     538             :     {
     539           0 :         delete poDS;
     540           0 :         CPLError(CE_Failure, CPLE_FileIO, "Header tag not found.\n");
     541           0 :         return nullptr;
     542             :     }
     543             : 
     544             :     GUInt32 nSize;
     545          33 :     if (VSIFReadL((void *)&nSize, sizeof(GUInt32), 1, poDS->fp) != 1)
     546             :     {
     547           0 :         delete poDS;
     548           0 :         CPLError(CE_Failure, CPLE_FileIO,
     549             :                  "Unable to read file section size.\n");
     550           0 :         return nullptr;
     551             :     }
     552             : 
     553          33 :     CPL_LSBPTR32(&nSize);
     554             : 
     555             :     GInt32 nVersion;
     556          33 :     if (VSIFReadL((void *)&nVersion, sizeof(GInt32), 1, poDS->fp) != 1)
     557             :     {
     558           0 :         delete poDS;
     559           0 :         CPLError(CE_Failure, CPLE_FileIO, "Unable to read file version.\n");
     560           0 :         return nullptr;
     561             :     }
     562             : 
     563          33 :     CPL_LSBPTR32(&nVersion);
     564             : 
     565          33 :     if (nVersion != 1 && nVersion != 2)
     566             :     {
     567           0 :         delete poDS;
     568           0 :         CPLError(CE_Failure, CPLE_FileIO, "Incorrect file version (%d).",
     569             :                  nVersion);
     570           0 :         return nullptr;
     571             :     }
     572             : 
     573             :     // advance until the grid tag is found
     574          66 :     while (nTag != nGRID_TAG)
     575             :     {
     576          33 :         if (VSIFReadL((void *)&nTag, sizeof(GInt32), 1, poDS->fp) != 1)
     577             :         {
     578           0 :             delete poDS;
     579           0 :             CPLError(CE_Failure, CPLE_FileIO, "Unable to read Tag.\n");
     580           0 :             return nullptr;
     581             :         }
     582             : 
     583          33 :         CPL_LSBPTR32(&nTag);
     584             : 
     585          33 :         if (VSIFReadL((void *)&nSize, sizeof(GUInt32), 1, poDS->fp) != 1)
     586             :         {
     587           0 :             delete poDS;
     588           0 :             CPLError(CE_Failure, CPLE_FileIO,
     589             :                      "Unable to read file section size.\n");
     590           0 :             return nullptr;
     591             :         }
     592             : 
     593          33 :         CPL_LSBPTR32(&nSize);
     594             : 
     595          33 :         if (nTag != nGRID_TAG)
     596             :         {
     597           0 :             if (VSIFSeekL(poDS->fp, nSize, SEEK_CUR) != 0)
     598             :             {
     599           0 :                 delete poDS;
     600           0 :                 CPLError(CE_Failure, CPLE_FileIO,
     601             :                          "Unable to seek to end of file section.\n");
     602           0 :                 return nullptr;
     603             :             }
     604             :         }
     605             :     }
     606             : 
     607             :     /* --------------------------------------------------------------------*/
     608             :     /*      Read the grid.                                                 */
     609             :     /* --------------------------------------------------------------------*/
     610             :     /* Parse number of Y axis grid rows */
     611             :     GInt32 nRows;
     612          33 :     if (VSIFReadL((void *)&nRows, sizeof(GInt32), 1, poDS->fp) != 1)
     613             :     {
     614           0 :         delete poDS;
     615           0 :         CPLError(CE_Failure, CPLE_FileIO, "Unable to read raster Y size.\n");
     616           0 :         return nullptr;
     617             :     }
     618          33 :     CPL_LSBPTR32(&nRows);
     619          33 :     poDS->nRasterYSize = nRows;
     620             : 
     621             :     /* Parse number of X axis grid columns */
     622             :     GInt32 nCols;
     623          33 :     if (VSIFReadL((void *)&nCols, sizeof(GInt32), 1, poDS->fp) != 1)
     624             :     {
     625           0 :         delete poDS;
     626           0 :         CPLError(CE_Failure, CPLE_FileIO, "Unable to read raster X size.\n");
     627           0 :         return nullptr;
     628             :     }
     629          33 :     CPL_LSBPTR32(&nCols);
     630          33 :     poDS->nRasterXSize = nCols;
     631             : 
     632          33 :     if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
     633             :     {
     634           0 :         delete poDS;
     635           0 :         return nullptr;
     636             :     }
     637             : 
     638             :     /* --------------------------------------------------------------------*/
     639             :     /*      Create band information objects.                               */
     640             :     /* --------------------------------------------------------------------*/
     641          33 :     GS7BGRasterBand *poBand = new GS7BGRasterBand(poDS, 1);
     642          33 :     poDS->SetBand(1, poBand);
     643             : 
     644             :     // find the min X Value of the grid
     645             :     double dfTemp;
     646          33 :     if (VSIFReadL((void *)&dfTemp, sizeof(double), 1, poDS->fp) != 1)
     647             :     {
     648           0 :         delete poDS;
     649           0 :         CPLError(CE_Failure, CPLE_FileIO, "Unable to read minimum X value.\n");
     650           0 :         return nullptr;
     651             :     }
     652          33 :     CPL_LSBPTR64(&dfTemp);
     653          33 :     poBand->dfMinX = dfTemp;
     654             : 
     655             :     // find the min Y value of the grid
     656          33 :     if (VSIFReadL((void *)&dfTemp, sizeof(double), 1, poDS->fp) != 1)
     657             :     {
     658           0 :         delete poDS;
     659           0 :         CPLError(CE_Failure, CPLE_FileIO, "Unable to read minimum X value.\n");
     660           0 :         return nullptr;
     661             :     }
     662          33 :     CPL_LSBPTR64(&dfTemp);
     663          33 :     poBand->dfMinY = dfTemp;
     664             : 
     665             :     // find the spacing between adjacent nodes in the X direction
     666             :     // (between columns)
     667          33 :     if (VSIFReadL((void *)&dfTemp, sizeof(double), 1, poDS->fp) != 1)
     668             :     {
     669           0 :         delete poDS;
     670           0 :         CPLError(CE_Failure, CPLE_FileIO,
     671             :                  "Unable to read spacing in X value.\n");
     672           0 :         return nullptr;
     673             :     }
     674          33 :     CPL_LSBPTR64(&dfTemp);
     675          33 :     poBand->dfMaxX = poBand->dfMinX + (dfTemp * (nCols - 1));
     676             : 
     677             :     // find the spacing between adjacent nodes in the Y direction
     678             :     // (between rows)
     679          33 :     if (VSIFReadL((void *)&dfTemp, sizeof(double), 1, poDS->fp) != 1)
     680             :     {
     681           0 :         delete poDS;
     682           0 :         CPLError(CE_Failure, CPLE_FileIO,
     683             :                  "Unable to read spacing in Y value.\n");
     684           0 :         return nullptr;
     685             :     }
     686          33 :     CPL_LSBPTR64(&dfTemp);
     687          33 :     poBand->dfMaxY = poBand->dfMinY + (dfTemp * (nRows - 1));
     688             : 
     689             :     // set the z min
     690          33 :     if (VSIFReadL((void *)&dfTemp, sizeof(double), 1, poDS->fp) != 1)
     691             :     {
     692           0 :         delete poDS;
     693           0 :         CPLError(CE_Failure, CPLE_FileIO, "Unable to read Z min value.\n");
     694           0 :         return nullptr;
     695             :     }
     696          33 :     CPL_LSBPTR64(&dfTemp);
     697          33 :     poBand->dfMinZ = dfTemp;
     698             : 
     699             :     // set the z max
     700          33 :     if (VSIFReadL((void *)&dfTemp, sizeof(double), 1, poDS->fp) != 1)
     701             :     {
     702           0 :         delete poDS;
     703           0 :         CPLError(CE_Failure, CPLE_FileIO, "Unable to read Z max value.\n");
     704           0 :         return nullptr;
     705             :     }
     706          33 :     CPL_LSBPTR64(&dfTemp);
     707          33 :     poBand->dfMaxZ = dfTemp;
     708             : 
     709             :     // read and ignore the rotation value
     710             :     //(This is not used in the current version).
     711          33 :     if (VSIFReadL((void *)&dfTemp, sizeof(double), 1, poDS->fp) != 1)
     712             :     {
     713           0 :         delete poDS;
     714           0 :         CPLError(CE_Failure, CPLE_FileIO, "Unable to read rotation value.\n");
     715           0 :         return nullptr;
     716             :     }
     717             : 
     718             :     // read and set the cell blank value
     719          33 :     if (VSIFReadL((void *)&dfTemp, sizeof(double), 1, poDS->fp) != 1)
     720             :     {
     721           0 :         delete poDS;
     722           0 :         CPLError(CE_Failure, CPLE_FileIO, "Unable to Blank value.\n");
     723           0 :         return nullptr;
     724             :     }
     725          33 :     CPL_LSBPTR64(&dfTemp);
     726          33 :     poDS->dfNoData_Value = dfTemp;
     727             : 
     728             :     /* --------------------------------------------------------------------*/
     729             :     /*      Set the current offset of the grid data.                       */
     730             :     /* --------------------------------------------------------------------*/
     731          33 :     if (VSIFReadL((void *)&nTag, sizeof(GInt32), 1, poDS->fp) != 1)
     732             :     {
     733           0 :         delete poDS;
     734           0 :         CPLError(CE_Failure, CPLE_FileIO, "Unable to read Tag.\n");
     735           0 :         return nullptr;
     736             :     }
     737             : 
     738          33 :     CPL_LSBPTR32(&nTag);
     739          33 :     if (nTag != nDATA_TAG)
     740             :     {
     741           0 :         delete poDS;
     742           0 :         CPLError(CE_Failure, CPLE_FileIO, "Data tag not found.\n");
     743           0 :         return nullptr;
     744             :     }
     745             : 
     746          33 :     if (VSIFReadL((void *)&nSize, sizeof(GInt32), 1, poDS->fp) != 1)
     747             :     {
     748           0 :         delete poDS;
     749           0 :         CPLError(CE_Failure, CPLE_FileIO, "Unable to data section size.\n");
     750           0 :         return nullptr;
     751             :     }
     752             : 
     753          33 :     poDS->nData_Position = (size_t)VSIFTellL(poDS->fp);
     754             : 
     755             :     /* --------------------------------------------------------------------*/
     756             :     /*      Initialize any PAM information.                                */
     757             :     /* --------------------------------------------------------------------*/
     758          33 :     poDS->SetDescription(poOpenInfo->pszFilename);
     759          33 :     poDS->TryLoadXML();
     760             : 
     761             :     /* -------------------------------------------------------------------- */
     762             :     /*      Check for external overviews.                                   */
     763             :     /* -------------------------------------------------------------------- */
     764          66 :     poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename,
     765          33 :                                 poOpenInfo->GetSiblingFiles());
     766             : 
     767          33 :     return poDS;
     768             : }
     769             : 
     770             : /************************************************************************/
     771             : /*                          GetGeoTransform()                           */
     772             : /************************************************************************/
     773             : 
     774          23 : CPLErr GS7BGDataset::GetGeoTransform(double *padfGeoTransform)
     775             : {
     776          23 :     if (padfGeoTransform == nullptr)
     777           0 :         return CE_Failure;
     778             : 
     779          23 :     GS7BGRasterBand *poGRB = (GS7BGRasterBand *)GetRasterBand(1);
     780             : 
     781          23 :     if (poGRB == nullptr)
     782             :     {
     783           0 :         padfGeoTransform[0] = 0;
     784           0 :         padfGeoTransform[1] = 1;
     785           0 :         padfGeoTransform[2] = 0;
     786           0 :         padfGeoTransform[3] = 0;
     787           0 :         padfGeoTransform[4] = 0;
     788           0 :         padfGeoTransform[5] = 1;
     789           0 :         return CE_Failure;
     790             :     }
     791             : 
     792             :     /* check if we have a PAM GeoTransform stored */
     793          23 :     CPLPushErrorHandler(CPLQuietErrorHandler);
     794          23 :     CPLErr eErr = GDALPamDataset::GetGeoTransform(padfGeoTransform);
     795          23 :     CPLPopErrorHandler();
     796             : 
     797          23 :     if (eErr == CE_None)
     798           0 :         return CE_None;
     799             : 
     800          23 :     if (nRasterXSize == 1 || nRasterYSize == 1)
     801           0 :         return CE_Failure;
     802             : 
     803             :     /* calculate pixel size first */
     804          23 :     padfGeoTransform[1] = (poGRB->dfMaxX - poGRB->dfMinX) / (nRasterXSize - 1);
     805          23 :     padfGeoTransform[5] = (poGRB->dfMinY - poGRB->dfMaxY) / (nRasterYSize - 1);
     806             : 
     807             :     /* then calculate image origin */
     808          23 :     padfGeoTransform[0] = poGRB->dfMinX - padfGeoTransform[1] / 2;
     809          23 :     padfGeoTransform[3] = poGRB->dfMaxY - padfGeoTransform[5] / 2;
     810             : 
     811             :     /* tilt/rotation does not supported by the GS grids */
     812          23 :     padfGeoTransform[4] = 0.0;
     813          23 :     padfGeoTransform[2] = 0.0;
     814             : 
     815          23 :     return CE_None;
     816             : }
     817             : 
     818             : /************************************************************************/
     819             : /*                          SetGeoTransform()                           */
     820             : /************************************************************************/
     821             : 
     822           6 : CPLErr GS7BGDataset::SetGeoTransform(double *padfGeoTransform)
     823             : {
     824           6 :     if (eAccess == GA_ReadOnly)
     825             :     {
     826           0 :         CPLError(CE_Failure, CPLE_NoWriteAccess,
     827             :                  "Unable to set GeoTransform, dataset opened read only.\n");
     828           0 :         return CE_Failure;
     829             :     }
     830             : 
     831             :     GS7BGRasterBand *poGRB =
     832           6 :         cpl::down_cast<GS7BGRasterBand *>(GetRasterBand(1));
     833             : 
     834           6 :     if (padfGeoTransform == nullptr)
     835           0 :         return CE_Failure;
     836             : 
     837             :     /* non-zero transform 2 or 4 or negative 1 or 5 not supported natively */
     838             :     /*if( padfGeoTransform[2] != 0.0 || padfGeoTransform[4] != 0.0
     839             :     || padfGeoTransform[1] < 0.0 || padfGeoTransform[5] < 0.0 )
     840             :     eErr = GDALPamDataset::SetGeoTransform( padfGeoTransform );
     841             : 
     842             :     if( eErr != CE_None )
     843             :     return eErr;*/
     844             : 
     845           6 :     double dfMinX = padfGeoTransform[0] + padfGeoTransform[1] / 2;
     846           6 :     double dfMaxX =
     847           6 :         padfGeoTransform[1] * (nRasterXSize - 0.5) + padfGeoTransform[0];
     848           6 :     double dfMinY =
     849           6 :         padfGeoTransform[5] * (nRasterYSize - 0.5) + padfGeoTransform[3];
     850           6 :     double dfMaxY = padfGeoTransform[3] + padfGeoTransform[5] / 2;
     851             : 
     852             :     CPLErr eErr =
     853           6 :         WriteHeader(fp, poGRB->nRasterXSize, poGRB->nRasterYSize, dfMinX,
     854             :                     dfMaxX, dfMinY, dfMaxY, poGRB->dfMinZ, poGRB->dfMaxZ);
     855             : 
     856           6 :     if (eErr == CE_None)
     857             :     {
     858           6 :         poGRB->dfMinX = dfMinX;
     859           6 :         poGRB->dfMaxX = dfMaxX;
     860           6 :         poGRB->dfMinY = dfMinY;
     861           6 :         poGRB->dfMaxY = dfMaxY;
     862             :     }
     863             : 
     864           6 :     return eErr;
     865             : }
     866             : 
     867             : /************************************************************************/
     868             : /*                             WriteHeader()                            */
     869             : /************************************************************************/
     870             : 
     871          53 : CPLErr GS7BGDataset::WriteHeader(VSILFILE *fp, GInt32 nXSize, GInt32 nYSize,
     872             :                                  double dfMinX, double dfMaxX, double dfMinY,
     873             :                                  double dfMaxY, double dfMinZ, double dfMaxZ)
     874             : 
     875             : {
     876          53 :     if (VSIFSeekL(fp, 0, SEEK_SET) != 0)
     877             :     {
     878           0 :         CPLError(CE_Failure, CPLE_FileIO,
     879             :                  "Unable to seek to start of grid file.\n");
     880           0 :         return CE_Failure;
     881             :     }
     882             : 
     883          53 :     GInt32 nTemp = CPL_LSBWORD32(nHEADER_TAG);
     884          53 :     if (VSIFWriteL((void *)&nTemp, sizeof(GInt32), 1, fp) != 1)
     885             :     {
     886           1 :         CPLError(CE_Failure, CPLE_FileIO,
     887             :                  "Unable to write header tag to grid file.\n");
     888           1 :         return CE_Failure;
     889             :     }
     890             : 
     891          52 :     nTemp = CPL_LSBWORD32(sizeof(GInt32));  // Size of version section.
     892          52 :     if (VSIFWriteL((void *)&nTemp, sizeof(GInt32), 1, fp) != 1)
     893             :     {
     894           0 :         CPLError(CE_Failure, CPLE_FileIO,
     895             :                  "Unable to write size to grid file.\n");
     896           0 :         return CE_Failure;
     897             :     }
     898             : 
     899          52 :     nTemp = CPL_LSBWORD32(1);  // Version
     900          52 :     if (VSIFWriteL((void *)&nTemp, sizeof(GInt32), 1, fp) != 1)
     901             :     {
     902           0 :         CPLError(CE_Failure, CPLE_FileIO,
     903             :                  "Unable to write size to grid file.\n");
     904           0 :         return CE_Failure;
     905             :     }
     906             : 
     907          52 :     nTemp = CPL_LSBWORD32(nGRID_TAG);  // Mark start of grid
     908          52 :     if (VSIFWriteL((void *)&nTemp, sizeof(GInt32), 1, fp) != 1)
     909             :     {
     910           0 :         CPLError(CE_Failure, CPLE_FileIO,
     911             :                  "Unable to write size to grid file.\n");
     912           0 :         return CE_Failure;
     913             :     }
     914             : 
     915          52 :     nTemp = CPL_LSBWORD32(72);  // Grid info size (the remainder of the header)
     916          52 :     if (VSIFWriteL((void *)&nTemp, sizeof(GInt32), 1, fp) != 1)
     917             :     {
     918           0 :         CPLError(CE_Failure, CPLE_FileIO,
     919             :                  "Unable to write size to grid file.\n");
     920           0 :         return CE_Failure;
     921             :     }
     922             : 
     923          52 :     nTemp = CPL_LSBWORD32(nYSize);
     924          52 :     if (VSIFWriteL((void *)&nTemp, sizeof(GInt32), 1, fp) != 1)
     925             :     {
     926           0 :         CPLError(CE_Failure, CPLE_FileIO,
     927             :                  "Unable to write Y size to grid file.\n");
     928           0 :         return CE_Failure;
     929             :     }
     930             : 
     931          52 :     nTemp = CPL_LSBWORD32(nXSize);
     932          52 :     if (VSIFWriteL((void *)&nTemp, sizeof(GInt32), 1, fp) != 1)
     933             :     {
     934           0 :         CPLError(CE_Failure, CPLE_FileIO,
     935             :                  "Unable to write X size to grid file.\n");
     936           0 :         return CE_Failure;
     937             :     }
     938             : 
     939          52 :     double dfTemp = dfMinX;
     940          52 :     CPL_LSBPTR64(&dfTemp);
     941          52 :     if (VSIFWriteL((void *)&dfTemp, sizeof(double), 1, fp) != 1)
     942             :     {
     943           0 :         CPLError(CE_Failure, CPLE_FileIO,
     944             :                  "Unable to write minimum X value to grid file.\n");
     945           0 :         return CE_Failure;
     946             :     }
     947             : 
     948          52 :     dfTemp = dfMinY;
     949          52 :     CPL_LSBPTR64(&dfTemp);
     950          52 :     if (VSIFWriteL((void *)&dfTemp, sizeof(double), 1, fp) != 1)
     951             :     {
     952           0 :         CPLError(CE_Failure, CPLE_FileIO,
     953             :                  "Unable to write minimum Y value to grid file.\n");
     954           0 :         return CE_Failure;
     955             :     }
     956             : 
     957             :     // Write node spacing in x direction
     958          52 :     dfTemp = (dfMaxX - dfMinX) / (nXSize - 1);
     959          52 :     CPL_LSBPTR64(&dfTemp);
     960          52 :     if (VSIFWriteL((void *)&dfTemp, sizeof(double), 1, fp) != 1)
     961             :     {
     962           0 :         CPLError(CE_Failure, CPLE_FileIO,
     963             :                  "Unable to write spacing in X value.\n");
     964           0 :         return CE_Failure;
     965             :     }
     966             : 
     967             :     // Write node spacing in y direction
     968          52 :     dfTemp = (dfMaxY - dfMinY) / (nYSize - 1);
     969          52 :     CPL_LSBPTR64(&dfTemp);
     970          52 :     if (VSIFWriteL((void *)&dfTemp, sizeof(double), 1, fp) != 1)
     971             :     {
     972           0 :         CPLError(CE_Failure, CPLE_FileIO,
     973             :                  "Unable to write spacing in Y value.\n");
     974           0 :         return CE_Failure;
     975             :     }
     976             : 
     977          52 :     dfTemp = dfMinZ;
     978          52 :     CPL_LSBPTR64(&dfTemp);
     979          52 :     if (VSIFWriteL((void *)&dfTemp, sizeof(double), 1, fp) != 1)
     980             :     {
     981           0 :         CPLError(CE_Failure, CPLE_FileIO,
     982             :                  "Unable to write minimum Z value to grid file.\n");
     983           0 :         return CE_Failure;
     984             :     }
     985             : 
     986          52 :     dfTemp = dfMaxZ;
     987          52 :     CPL_LSBPTR64(&dfTemp);
     988          52 :     if (VSIFWriteL((void *)&dfTemp, sizeof(double), 1, fp) != 1)
     989             :     {
     990           0 :         CPLError(CE_Failure, CPLE_FileIO,
     991             :                  "Unable to write maximum Z value to grid file.\n");
     992           0 :         return CE_Failure;
     993             :     }
     994             : 
     995          52 :     dfTemp = 0;  // Rotation value is zero
     996          52 :     CPL_LSBPTR64(&dfTemp);
     997          52 :     if (VSIFWriteL((void *)&dfTemp, sizeof(double), 1, fp) != 1)
     998             :     {
     999           0 :         CPLError(CE_Failure, CPLE_FileIO,
    1000             :                  "Unable to write rotation value to grid file.\n");
    1001           0 :         return CE_Failure;
    1002             :     }
    1003             : 
    1004          52 :     dfTemp = dfDefaultNoDataValue;
    1005          52 :     CPL_LSBPTR64(&dfTemp);
    1006          52 :     if (VSIFWriteL((void *)&dfTemp, sizeof(double), 1, fp) != 1)
    1007             :     {
    1008           1 :         CPLError(CE_Failure, CPLE_FileIO,
    1009             :                  "Unable to write cell blank value to grid file.\n");
    1010           1 :         return CE_Failure;
    1011             :     }
    1012             : 
    1013             :     // Only supports 1 band so go ahead and write band info here
    1014          51 :     nTemp = CPL_LSBWORD32(nDATA_TAG);  // Mark start of data
    1015          51 :     if (VSIFWriteL((void *)&nTemp, sizeof(GInt32), 1, fp) != 1)
    1016             :     {
    1017           0 :         CPLError(CE_Failure, CPLE_FileIO, "Unable to data tag to grid file.\n");
    1018           0 :         return CE_Failure;
    1019             :     }
    1020             : 
    1021          51 :     int nSize = nXSize * nYSize * (int)sizeof(double);
    1022          51 :     nTemp = CPL_LSBWORD32(nSize);  // Mark size of data
    1023          51 :     if (VSIFWriteL((void *)&nTemp, sizeof(GInt32), 1, fp) != 1)
    1024             :     {
    1025           0 :         CPLError(CE_Failure, CPLE_FileIO,
    1026             :                  "Unable to write data size to grid file.\n");
    1027           0 :         return CE_Failure;
    1028             :     }
    1029             : 
    1030          51 :     return CE_None;
    1031             : }
    1032             : 
    1033             : /************************************************************************/
    1034             : /*                               Create()                               */
    1035             : /************************************************************************/
    1036             : 
    1037          33 : GDALDataset *GS7BGDataset::Create(const char *pszFilename, int nXSize,
    1038             :                                   int nYSize, int nBandsIn, GDALDataType eType,
    1039             :                                   CPL_UNUSED char **papszParamList)
    1040             : 
    1041             : {
    1042          33 :     if (nXSize <= 0 || nYSize <= 0)
    1043             :     {
    1044           0 :         CPLError(CE_Failure, CPLE_IllegalArg,
    1045             :                  "Unable to create grid, both X and Y size must be "
    1046             :                  "non-negative.\n");
    1047             : 
    1048           0 :         return nullptr;
    1049             :     }
    1050             : 
    1051          33 :     if (eType != GDT_Byte && eType != GDT_Float32 && eType != GDT_UInt16 &&
    1052          21 :         eType != GDT_Int16 && eType != GDT_Float64)
    1053             :     {
    1054          18 :         CPLError(
    1055             :             CE_Failure, CPLE_AppDefined,
    1056             :             "GS7BG Grid only supports Byte, Int16, "
    1057             :             "Uint16, Float32, and Float64 datatypes.  Unable to create with "
    1058             :             "type %s.\n",
    1059             :             GDALGetDataTypeName(eType));
    1060             : 
    1061          18 :         return nullptr;
    1062             :     }
    1063             : 
    1064          15 :     if (nBandsIn > 1)
    1065             :     {
    1066           8 :         CPLError(CE_Failure, CPLE_NotSupported,
    1067             :                  "Unable to create copy, "
    1068             :                  "format only supports one raster band.\n");
    1069           8 :         return nullptr;
    1070             :     }
    1071             : 
    1072           7 :     VSILFILE *fp = VSIFOpenL(pszFilename, "w+b");
    1073             : 
    1074           7 :     if (fp == nullptr)
    1075             :     {
    1076           0 :         CPLError(CE_Failure, CPLE_OpenFailed,
    1077             :                  "Attempt to create file '%s' failed.\n", pszFilename);
    1078           0 :         return nullptr;
    1079             :     }
    1080             : 
    1081             :     CPLErr eErr =
    1082           7 :         WriteHeader(fp, nXSize, nYSize, 0.0, nXSize, 0.0, nYSize, 0.0, 0.0);
    1083           7 :     if (eErr != CE_None)
    1084             :     {
    1085           0 :         VSIFCloseL(fp);
    1086           0 :         return nullptr;
    1087             :     }
    1088             : 
    1089           7 :     double dfVal = dfDefaultNoDataValue;
    1090           7 :     CPL_LSBPTR64(&dfVal);
    1091         627 :     for (int iRow = 0; iRow < nYSize; iRow++)
    1092             :     {
    1093       61020 :         for (int iCol = 0; iCol < nXSize; iCol++)
    1094             :         {
    1095       60400 :             if (VSIFWriteL((void *)&dfVal, sizeof(double), 1, fp) != 1)
    1096             :             {
    1097           0 :                 VSIFCloseL(fp);
    1098           0 :                 CPLError(CE_Failure, CPLE_FileIO,
    1099             :                          "Unable to write grid cell.  Disk full?\n");
    1100           0 :                 return nullptr;
    1101             :             }
    1102             :         }
    1103             :     }
    1104             : 
    1105           7 :     VSIFCloseL(fp);
    1106             : 
    1107           7 :     return (GDALDataset *)GDALOpen(pszFilename, GA_Update);
    1108             : }
    1109             : 
    1110             : /************************************************************************/
    1111             : /*                             CreateCopy()                             */
    1112             : /************************************************************************/
    1113             : 
    1114          30 : GDALDataset *GS7BGDataset::CreateCopy(const char *pszFilename,
    1115             :                                       GDALDataset *poSrcDS, int bStrict,
    1116             :                                       CPL_UNUSED char **papszOptions,
    1117             :                                       GDALProgressFunc pfnProgress,
    1118             :                                       void *pProgressData)
    1119             : {
    1120          30 :     if (pfnProgress == nullptr)
    1121           0 :         pfnProgress = GDALDummyProgress;
    1122             : 
    1123          30 :     int nBands = poSrcDS->GetRasterCount();
    1124          30 :     if (nBands == 0)
    1125             :     {
    1126           1 :         CPLError(CE_Failure, CPLE_NotSupported,
    1127             :                  "Driver does not support source dataset with zero band.\n");
    1128           1 :         return nullptr;
    1129             :     }
    1130          29 :     else if (nBands > 1)
    1131             :     {
    1132           4 :         if (bStrict)
    1133             :         {
    1134           4 :             CPLError(CE_Failure, CPLE_NotSupported,
    1135             :                      "Unable to create copy, "
    1136             :                      "format only supports one raster band.\n");
    1137           4 :             return nullptr;
    1138             :         }
    1139             :         else
    1140           0 :             CPLError(CE_Warning, CPLE_NotSupported,
    1141             :                      "Format only supports one "
    1142             :                      "raster band, first band will be copied.\n");
    1143             :     }
    1144             : 
    1145          25 :     GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(1);
    1146             : 
    1147          25 :     if (!pfnProgress(0.0, nullptr, pProgressData))
    1148             :     {
    1149           0 :         CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated\n");
    1150           0 :         return nullptr;
    1151             :     }
    1152             : 
    1153          25 :     VSILFILE *fp = VSIFOpenL(pszFilename, "w+b");
    1154             : 
    1155          25 :     if (fp == nullptr)
    1156             :     {
    1157           3 :         CPLError(CE_Failure, CPLE_OpenFailed,
    1158             :                  "Attempt to create file '%s' failed.\n", pszFilename);
    1159           3 :         return nullptr;
    1160             :     }
    1161             : 
    1162          22 :     GInt32 nXSize = poSrcBand->GetXSize();
    1163          22 :     GInt32 nYSize = poSrcBand->GetYSize();
    1164             :     double adfGeoTransform[6];
    1165             : 
    1166          22 :     poSrcDS->GetGeoTransform(adfGeoTransform);
    1167             : 
    1168          22 :     double dfMinX = adfGeoTransform[0] + adfGeoTransform[1] / 2;
    1169          22 :     double dfMaxX = adfGeoTransform[1] * (nXSize - 0.5) + adfGeoTransform[0];
    1170          22 :     double dfMinY = adfGeoTransform[5] * (nYSize - 0.5) + adfGeoTransform[3];
    1171          22 :     double dfMaxY = adfGeoTransform[3] + adfGeoTransform[5] / 2;
    1172          22 :     CPLErr eErr = WriteHeader(fp, nXSize, nYSize, dfMinX, dfMaxX, dfMinY,
    1173             :                               dfMaxY, 0.0, 0.0);
    1174             : 
    1175          22 :     if (eErr != CE_None)
    1176             :     {
    1177           2 :         VSIFCloseL(fp);
    1178           2 :         return nullptr;
    1179             :     }
    1180             : 
    1181             :     /* -------------------------------------------------------------------- */
    1182             :     /*      Copy band data.                                                 */
    1183             :     /* -------------------------------------------------------------------- */
    1184          20 :     double *pfData = (double *)VSI_MALLOC2_VERBOSE(nXSize, sizeof(double));
    1185          20 :     if (pfData == nullptr)
    1186             :     {
    1187           0 :         VSIFCloseL(fp);
    1188           0 :         return nullptr;
    1189             :     }
    1190             : 
    1191             :     int bSrcHasNDValue;
    1192          20 :     double dfSrcNoDataValue = poSrcBand->GetNoDataValue(&bSrcHasNDValue);
    1193          20 :     double dfMinZ = std::numeric_limits<double>::max();
    1194          20 :     double dfMaxZ = std::numeric_limits<double>::lowest();
    1195         186 :     for (GInt32 iRow = nYSize - 1; iRow >= 0; iRow--)
    1196             :     {
    1197         174 :         eErr = poSrcBand->RasterIO(GF_Read, 0, iRow, nXSize, 1, pfData, nXSize,
    1198             :                                    1, GDT_Float64, 0, 0, nullptr);
    1199             : 
    1200         174 :         if (eErr != CE_None)
    1201             :         {
    1202           0 :             VSIFCloseL(fp);
    1203           0 :             VSIFree(pfData);
    1204           0 :             return nullptr;
    1205             :         }
    1206             : 
    1207        2114 :         for (int iCol = 0; iCol < nXSize; iCol++)
    1208             :         {
    1209        1940 :             if (bSrcHasNDValue && pfData[iCol] == dfSrcNoDataValue)
    1210             :             {
    1211           0 :                 pfData[iCol] = dfDefaultNoDataValue;
    1212             :             }
    1213             :             else
    1214             :             {
    1215        1940 :                 if (pfData[iCol] > dfMaxZ)
    1216          22 :                     dfMaxZ = pfData[iCol];
    1217             : 
    1218        1940 :                 if (pfData[iCol] < dfMinZ)
    1219          27 :                     dfMinZ = pfData[iCol];
    1220             :             }
    1221             : 
    1222        1940 :             CPL_LSBPTR64(pfData + iCol);
    1223             :         }
    1224             : 
    1225         174 :         if (VSIFWriteL((void *)pfData, sizeof(double), nXSize, fp) !=
    1226         174 :             static_cast<unsigned>(nXSize))
    1227             :         {
    1228           8 :             VSIFCloseL(fp);
    1229           8 :             VSIFree(pfData);
    1230           8 :             CPLError(CE_Failure, CPLE_FileIO,
    1231             :                      "Unable to write grid row. Disk full?\n");
    1232           8 :             return nullptr;
    1233             :         }
    1234             : 
    1235         166 :         if (!pfnProgress(static_cast<double>(nYSize - iRow) / nYSize, nullptr,
    1236             :                          pProgressData))
    1237             :         {
    1238           0 :             VSIFCloseL(fp);
    1239           0 :             VSIFree(pfData);
    1240           0 :             CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
    1241           0 :             return nullptr;
    1242             :         }
    1243             :     }
    1244             : 
    1245          12 :     VSIFree(pfData);
    1246             : 
    1247             :     /* write out the min and max values */
    1248          12 :     eErr = WriteHeader(fp, nXSize, nYSize, dfMinX, dfMaxX, dfMinY, dfMaxY,
    1249             :                        dfMinZ, dfMaxZ);
    1250             : 
    1251          12 :     if (eErr != CE_None)
    1252             :     {
    1253           0 :         VSIFCloseL(fp);
    1254           0 :         return nullptr;
    1255             :     }
    1256             : 
    1257          12 :     VSIFCloseL(fp);
    1258             : 
    1259          12 :     GDALPamDataset *poDS = (GDALPamDataset *)GDALOpen(pszFilename, GA_Update);
    1260          12 :     if (poDS)
    1261             :     {
    1262          12 :         poDS->CloneInfo(poSrcDS, GCIF_PAM_DEFAULT);
    1263             :     }
    1264             : 
    1265          12 :     return poDS;
    1266             : }
    1267             : 
    1268             : /************************************************************************/
    1269             : /*                          GDALRegister_GS7BG()                        */
    1270             : /************************************************************************/
    1271        1595 : void GDALRegister_GS7BG()
    1272             : 
    1273             : {
    1274        1595 :     if (GDALGetDriverByName("GS7BG") != nullptr)
    1275         302 :         return;
    1276             : 
    1277        1293 :     GDALDriver *poDriver = new GDALDriver();
    1278             : 
    1279        1293 :     poDriver->SetDescription("GS7BG");
    1280        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    1281        1293 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
    1282        1293 :                               "Golden Software 7 Binary Grid (.grd)");
    1283        1293 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/gs7bg.html");
    1284        1293 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "grd");
    1285        1293 :     poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
    1286        1293 :                               "Byte Int16 UInt16 Float32 Float64");
    1287        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
    1288             : 
    1289        1293 :     poDriver->pfnIdentify = GS7BGDataset::Identify;
    1290        1293 :     poDriver->pfnOpen = GS7BGDataset::Open;
    1291        1293 :     poDriver->pfnCreate = GS7BGDataset::Create;
    1292        1293 :     poDriver->pfnCreateCopy = GS7BGDataset::CreateCopy;
    1293             : 
    1294        1293 :     GetGDALDriverManager()->RegisterDriver(poDriver);
    1295             : }