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 = cpl::down_cast<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 = cpl::down_cast<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       57969 : 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       57969 :     if (poOpenInfo->nHeaderBytes < 4 ||
     487        4556 :         !STARTS_WITH_CI((const char *)poOpenInfo->pabyHeader, "DSRB"))
     488             :     {
     489       57903 :         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             :     GS7BGRasterBand *poGRB =
     780          23 :         cpl::down_cast<GS7BGRasterBand *>(GetRasterBand(1));
     781             : 
     782          23 :     if (poGRB == nullptr)
     783             :     {
     784           0 :         padfGeoTransform[0] = 0;
     785           0 :         padfGeoTransform[1] = 1;
     786           0 :         padfGeoTransform[2] = 0;
     787           0 :         padfGeoTransform[3] = 0;
     788           0 :         padfGeoTransform[4] = 0;
     789           0 :         padfGeoTransform[5] = 1;
     790           0 :         return CE_Failure;
     791             :     }
     792             : 
     793             :     /* check if we have a PAM GeoTransform stored */
     794          23 :     CPLPushErrorHandler(CPLQuietErrorHandler);
     795          23 :     CPLErr eErr = GDALPamDataset::GetGeoTransform(padfGeoTransform);
     796          23 :     CPLPopErrorHandler();
     797             : 
     798          23 :     if (eErr == CE_None)
     799           0 :         return CE_None;
     800             : 
     801          23 :     if (nRasterXSize == 1 || nRasterYSize == 1)
     802           0 :         return CE_Failure;
     803             : 
     804             :     /* calculate pixel size first */
     805          23 :     padfGeoTransform[1] = (poGRB->dfMaxX - poGRB->dfMinX) / (nRasterXSize - 1);
     806          23 :     padfGeoTransform[5] = (poGRB->dfMinY - poGRB->dfMaxY) / (nRasterYSize - 1);
     807             : 
     808             :     /* then calculate image origin */
     809          23 :     padfGeoTransform[0] = poGRB->dfMinX - padfGeoTransform[1] / 2;
     810          23 :     padfGeoTransform[3] = poGRB->dfMaxY - padfGeoTransform[5] / 2;
     811             : 
     812             :     /* tilt/rotation does not supported by the GS grids */
     813          23 :     padfGeoTransform[4] = 0.0;
     814          23 :     padfGeoTransform[2] = 0.0;
     815             : 
     816          23 :     return CE_None;
     817             : }
     818             : 
     819             : /************************************************************************/
     820             : /*                          SetGeoTransform()                           */
     821             : /************************************************************************/
     822             : 
     823           6 : CPLErr GS7BGDataset::SetGeoTransform(double *padfGeoTransform)
     824             : {
     825           6 :     if (eAccess == GA_ReadOnly)
     826             :     {
     827           0 :         CPLError(CE_Failure, CPLE_NoWriteAccess,
     828             :                  "Unable to set GeoTransform, dataset opened read only.\n");
     829           0 :         return CE_Failure;
     830             :     }
     831             : 
     832             :     GS7BGRasterBand *poGRB =
     833           6 :         cpl::down_cast<GS7BGRasterBand *>(GetRasterBand(1));
     834             : 
     835           6 :     if (padfGeoTransform == nullptr)
     836           0 :         return CE_Failure;
     837             : 
     838             :     /* non-zero transform 2 or 4 or negative 1 or 5 not supported natively */
     839             :     /*if( padfGeoTransform[2] != 0.0 || padfGeoTransform[4] != 0.0
     840             :     || padfGeoTransform[1] < 0.0 || padfGeoTransform[5] < 0.0 )
     841             :     eErr = GDALPamDataset::SetGeoTransform( padfGeoTransform );
     842             : 
     843             :     if( eErr != CE_None )
     844             :     return eErr;*/
     845             : 
     846           6 :     double dfMinX = padfGeoTransform[0] + padfGeoTransform[1] / 2;
     847           6 :     double dfMaxX =
     848           6 :         padfGeoTransform[1] * (nRasterXSize - 0.5) + padfGeoTransform[0];
     849           6 :     double dfMinY =
     850           6 :         padfGeoTransform[5] * (nRasterYSize - 0.5) + padfGeoTransform[3];
     851           6 :     double dfMaxY = padfGeoTransform[3] + padfGeoTransform[5] / 2;
     852             : 
     853             :     CPLErr eErr =
     854           6 :         WriteHeader(fp, poGRB->nRasterXSize, poGRB->nRasterYSize, dfMinX,
     855             :                     dfMaxX, dfMinY, dfMaxY, poGRB->dfMinZ, poGRB->dfMaxZ);
     856             : 
     857           6 :     if (eErr == CE_None)
     858             :     {
     859           6 :         poGRB->dfMinX = dfMinX;
     860           6 :         poGRB->dfMaxX = dfMaxX;
     861           6 :         poGRB->dfMinY = dfMinY;
     862           6 :         poGRB->dfMaxY = dfMaxY;
     863             :     }
     864             : 
     865           6 :     return eErr;
     866             : }
     867             : 
     868             : /************************************************************************/
     869             : /*                             WriteHeader()                            */
     870             : /************************************************************************/
     871             : 
     872          53 : CPLErr GS7BGDataset::WriteHeader(VSILFILE *fp, GInt32 nXSize, GInt32 nYSize,
     873             :                                  double dfMinX, double dfMaxX, double dfMinY,
     874             :                                  double dfMaxY, double dfMinZ, double dfMaxZ)
     875             : 
     876             : {
     877          53 :     if (VSIFSeekL(fp, 0, SEEK_SET) != 0)
     878             :     {
     879           0 :         CPLError(CE_Failure, CPLE_FileIO,
     880             :                  "Unable to seek to start of grid file.\n");
     881           0 :         return CE_Failure;
     882             :     }
     883             : 
     884          53 :     GInt32 nTemp = CPL_LSBWORD32(nHEADER_TAG);
     885          53 :     if (VSIFWriteL((void *)&nTemp, sizeof(GInt32), 1, fp) != 1)
     886             :     {
     887           1 :         CPLError(CE_Failure, CPLE_FileIO,
     888             :                  "Unable to write header tag to grid file.\n");
     889           1 :         return CE_Failure;
     890             :     }
     891             : 
     892          52 :     nTemp = CPL_LSBWORD32(sizeof(GInt32));  // Size of version section.
     893          52 :     if (VSIFWriteL((void *)&nTemp, sizeof(GInt32), 1, fp) != 1)
     894             :     {
     895           0 :         CPLError(CE_Failure, CPLE_FileIO,
     896             :                  "Unable to write size to grid file.\n");
     897           0 :         return CE_Failure;
     898             :     }
     899             : 
     900          52 :     nTemp = CPL_LSBWORD32(1);  // Version
     901          52 :     if (VSIFWriteL((void *)&nTemp, sizeof(GInt32), 1, fp) != 1)
     902             :     {
     903           0 :         CPLError(CE_Failure, CPLE_FileIO,
     904             :                  "Unable to write size to grid file.\n");
     905           0 :         return CE_Failure;
     906             :     }
     907             : 
     908          52 :     nTemp = CPL_LSBWORD32(nGRID_TAG);  // Mark start of grid
     909          52 :     if (VSIFWriteL((void *)&nTemp, sizeof(GInt32), 1, fp) != 1)
     910             :     {
     911           0 :         CPLError(CE_Failure, CPLE_FileIO,
     912             :                  "Unable to write size to grid file.\n");
     913           0 :         return CE_Failure;
     914             :     }
     915             : 
     916          52 :     nTemp = CPL_LSBWORD32(72);  // Grid info size (the remainder of the header)
     917          52 :     if (VSIFWriteL((void *)&nTemp, sizeof(GInt32), 1, fp) != 1)
     918             :     {
     919           0 :         CPLError(CE_Failure, CPLE_FileIO,
     920             :                  "Unable to write size to grid file.\n");
     921           0 :         return CE_Failure;
     922             :     }
     923             : 
     924          52 :     nTemp = CPL_LSBWORD32(nYSize);
     925          52 :     if (VSIFWriteL((void *)&nTemp, sizeof(GInt32), 1, fp) != 1)
     926             :     {
     927           0 :         CPLError(CE_Failure, CPLE_FileIO,
     928             :                  "Unable to write Y size to grid file.\n");
     929           0 :         return CE_Failure;
     930             :     }
     931             : 
     932          52 :     nTemp = CPL_LSBWORD32(nXSize);
     933          52 :     if (VSIFWriteL((void *)&nTemp, sizeof(GInt32), 1, fp) != 1)
     934             :     {
     935           0 :         CPLError(CE_Failure, CPLE_FileIO,
     936             :                  "Unable to write X size to grid file.\n");
     937           0 :         return CE_Failure;
     938             :     }
     939             : 
     940          52 :     double dfTemp = dfMinX;
     941          52 :     CPL_LSBPTR64(&dfTemp);
     942          52 :     if (VSIFWriteL((void *)&dfTemp, sizeof(double), 1, fp) != 1)
     943             :     {
     944           0 :         CPLError(CE_Failure, CPLE_FileIO,
     945             :                  "Unable to write minimum X value to grid file.\n");
     946           0 :         return CE_Failure;
     947             :     }
     948             : 
     949          52 :     dfTemp = dfMinY;
     950          52 :     CPL_LSBPTR64(&dfTemp);
     951          52 :     if (VSIFWriteL((void *)&dfTemp, sizeof(double), 1, fp) != 1)
     952             :     {
     953           0 :         CPLError(CE_Failure, CPLE_FileIO,
     954             :                  "Unable to write minimum Y value to grid file.\n");
     955           0 :         return CE_Failure;
     956             :     }
     957             : 
     958             :     // Write node spacing in x direction
     959          52 :     dfTemp = (dfMaxX - dfMinX) / (nXSize - 1);
     960          52 :     CPL_LSBPTR64(&dfTemp);
     961          52 :     if (VSIFWriteL((void *)&dfTemp, sizeof(double), 1, fp) != 1)
     962             :     {
     963           0 :         CPLError(CE_Failure, CPLE_FileIO,
     964             :                  "Unable to write spacing in X value.\n");
     965           0 :         return CE_Failure;
     966             :     }
     967             : 
     968             :     // Write node spacing in y direction
     969          52 :     dfTemp = (dfMaxY - dfMinY) / (nYSize - 1);
     970          52 :     CPL_LSBPTR64(&dfTemp);
     971          52 :     if (VSIFWriteL((void *)&dfTemp, sizeof(double), 1, fp) != 1)
     972             :     {
     973           0 :         CPLError(CE_Failure, CPLE_FileIO,
     974             :                  "Unable to write spacing in Y value.\n");
     975           0 :         return CE_Failure;
     976             :     }
     977             : 
     978          52 :     dfTemp = dfMinZ;
     979          52 :     CPL_LSBPTR64(&dfTemp);
     980          52 :     if (VSIFWriteL((void *)&dfTemp, sizeof(double), 1, fp) != 1)
     981             :     {
     982           0 :         CPLError(CE_Failure, CPLE_FileIO,
     983             :                  "Unable to write minimum Z value to grid file.\n");
     984           0 :         return CE_Failure;
     985             :     }
     986             : 
     987          52 :     dfTemp = dfMaxZ;
     988          52 :     CPL_LSBPTR64(&dfTemp);
     989          52 :     if (VSIFWriteL((void *)&dfTemp, sizeof(double), 1, fp) != 1)
     990             :     {
     991           0 :         CPLError(CE_Failure, CPLE_FileIO,
     992             :                  "Unable to write maximum Z value to grid file.\n");
     993           0 :         return CE_Failure;
     994             :     }
     995             : 
     996          52 :     dfTemp = 0;  // Rotation value is zero
     997          52 :     CPL_LSBPTR64(&dfTemp);
     998          52 :     if (VSIFWriteL((void *)&dfTemp, sizeof(double), 1, fp) != 1)
     999             :     {
    1000           0 :         CPLError(CE_Failure, CPLE_FileIO,
    1001             :                  "Unable to write rotation value to grid file.\n");
    1002           0 :         return CE_Failure;
    1003             :     }
    1004             : 
    1005          52 :     dfTemp = dfDefaultNoDataValue;
    1006          52 :     CPL_LSBPTR64(&dfTemp);
    1007          52 :     if (VSIFWriteL((void *)&dfTemp, sizeof(double), 1, fp) != 1)
    1008             :     {
    1009           1 :         CPLError(CE_Failure, CPLE_FileIO,
    1010             :                  "Unable to write cell blank value to grid file.\n");
    1011           1 :         return CE_Failure;
    1012             :     }
    1013             : 
    1014             :     // Only supports 1 band so go ahead and write band info here
    1015          51 :     nTemp = CPL_LSBWORD32(nDATA_TAG);  // Mark start of data
    1016          51 :     if (VSIFWriteL((void *)&nTemp, sizeof(GInt32), 1, fp) != 1)
    1017             :     {
    1018           0 :         CPLError(CE_Failure, CPLE_FileIO, "Unable to data tag to grid file.\n");
    1019           0 :         return CE_Failure;
    1020             :     }
    1021             : 
    1022          51 :     int nSize = nXSize * nYSize * (int)sizeof(double);
    1023          51 :     nTemp = CPL_LSBWORD32(nSize);  // Mark size of data
    1024          51 :     if (VSIFWriteL((void *)&nTemp, sizeof(GInt32), 1, fp) != 1)
    1025             :     {
    1026           0 :         CPLError(CE_Failure, CPLE_FileIO,
    1027             :                  "Unable to write data size to grid file.\n");
    1028           0 :         return CE_Failure;
    1029             :     }
    1030             : 
    1031          51 :     return CE_None;
    1032             : }
    1033             : 
    1034             : /************************************************************************/
    1035             : /*                               Create()                               */
    1036             : /************************************************************************/
    1037             : 
    1038          33 : GDALDataset *GS7BGDataset::Create(const char *pszFilename, int nXSize,
    1039             :                                   int nYSize, int nBandsIn, GDALDataType eType,
    1040             :                                   CPL_UNUSED char **papszParamList)
    1041             : 
    1042             : {
    1043          33 :     if (nXSize <= 0 || nYSize <= 0)
    1044             :     {
    1045           0 :         CPLError(CE_Failure, CPLE_IllegalArg,
    1046             :                  "Unable to create grid, both X and Y size must be "
    1047             :                  "non-negative.\n");
    1048             : 
    1049           0 :         return nullptr;
    1050             :     }
    1051             : 
    1052          33 :     if (eType != GDT_Byte && eType != GDT_Float32 && eType != GDT_UInt16 &&
    1053          21 :         eType != GDT_Int16 && eType != GDT_Float64)
    1054             :     {
    1055          18 :         CPLError(
    1056             :             CE_Failure, CPLE_AppDefined,
    1057             :             "GS7BG Grid only supports Byte, Int16, "
    1058             :             "Uint16, Float32, and Float64 datatypes.  Unable to create with "
    1059             :             "type %s.\n",
    1060             :             GDALGetDataTypeName(eType));
    1061             : 
    1062          18 :         return nullptr;
    1063             :     }
    1064             : 
    1065          15 :     if (nBandsIn > 1)
    1066             :     {
    1067           8 :         CPLError(CE_Failure, CPLE_NotSupported,
    1068             :                  "Unable to create copy, "
    1069             :                  "format only supports one raster band.\n");
    1070           8 :         return nullptr;
    1071             :     }
    1072             : 
    1073           7 :     VSILFILE *fp = VSIFOpenL(pszFilename, "w+b");
    1074             : 
    1075           7 :     if (fp == nullptr)
    1076             :     {
    1077           0 :         CPLError(CE_Failure, CPLE_OpenFailed,
    1078             :                  "Attempt to create file '%s' failed.\n", pszFilename);
    1079           0 :         return nullptr;
    1080             :     }
    1081             : 
    1082             :     CPLErr eErr =
    1083           7 :         WriteHeader(fp, nXSize, nYSize, 0.0, nXSize, 0.0, nYSize, 0.0, 0.0);
    1084           7 :     if (eErr != CE_None)
    1085             :     {
    1086           0 :         VSIFCloseL(fp);
    1087           0 :         return nullptr;
    1088             :     }
    1089             : 
    1090           7 :     double dfVal = dfDefaultNoDataValue;
    1091           7 :     CPL_LSBPTR64(&dfVal);
    1092         627 :     for (int iRow = 0; iRow < nYSize; iRow++)
    1093             :     {
    1094       61020 :         for (int iCol = 0; iCol < nXSize; iCol++)
    1095             :         {
    1096       60400 :             if (VSIFWriteL((void *)&dfVal, sizeof(double), 1, fp) != 1)
    1097             :             {
    1098           0 :                 VSIFCloseL(fp);
    1099           0 :                 CPLError(CE_Failure, CPLE_FileIO,
    1100             :                          "Unable to write grid cell.  Disk full?\n");
    1101           0 :                 return nullptr;
    1102             :             }
    1103             :         }
    1104             :     }
    1105             : 
    1106           7 :     VSIFCloseL(fp);
    1107             : 
    1108           7 :     return (GDALDataset *)GDALOpen(pszFilename, GA_Update);
    1109             : }
    1110             : 
    1111             : /************************************************************************/
    1112             : /*                             CreateCopy()                             */
    1113             : /************************************************************************/
    1114             : 
    1115          30 : GDALDataset *GS7BGDataset::CreateCopy(const char *pszFilename,
    1116             :                                       GDALDataset *poSrcDS, int bStrict,
    1117             :                                       CPL_UNUSED char **papszOptions,
    1118             :                                       GDALProgressFunc pfnProgress,
    1119             :                                       void *pProgressData)
    1120             : {
    1121          30 :     if (pfnProgress == nullptr)
    1122           0 :         pfnProgress = GDALDummyProgress;
    1123             : 
    1124          30 :     int nBands = poSrcDS->GetRasterCount();
    1125          30 :     if (nBands == 0)
    1126             :     {
    1127           1 :         CPLError(CE_Failure, CPLE_NotSupported,
    1128             :                  "Driver does not support source dataset with zero band.\n");
    1129           1 :         return nullptr;
    1130             :     }
    1131          29 :     else if (nBands > 1)
    1132             :     {
    1133           4 :         if (bStrict)
    1134             :         {
    1135           4 :             CPLError(CE_Failure, CPLE_NotSupported,
    1136             :                      "Unable to create copy, "
    1137             :                      "format only supports one raster band.\n");
    1138           4 :             return nullptr;
    1139             :         }
    1140             :         else
    1141           0 :             CPLError(CE_Warning, CPLE_NotSupported,
    1142             :                      "Format only supports one "
    1143             :                      "raster band, first band will be copied.\n");
    1144             :     }
    1145             : 
    1146          25 :     GDALRasterBand *poSrcBand = poSrcDS->GetRasterBand(1);
    1147             : 
    1148          25 :     if (!pfnProgress(0.0, nullptr, pProgressData))
    1149             :     {
    1150           0 :         CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated\n");
    1151           0 :         return nullptr;
    1152             :     }
    1153             : 
    1154          25 :     VSILFILE *fp = VSIFOpenL(pszFilename, "w+b");
    1155             : 
    1156          25 :     if (fp == nullptr)
    1157             :     {
    1158           3 :         CPLError(CE_Failure, CPLE_OpenFailed,
    1159             :                  "Attempt to create file '%s' failed.\n", pszFilename);
    1160           3 :         return nullptr;
    1161             :     }
    1162             : 
    1163          22 :     GInt32 nXSize = poSrcBand->GetXSize();
    1164          22 :     GInt32 nYSize = poSrcBand->GetYSize();
    1165             :     double adfGeoTransform[6];
    1166             : 
    1167          22 :     poSrcDS->GetGeoTransform(adfGeoTransform);
    1168             : 
    1169          22 :     double dfMinX = adfGeoTransform[0] + adfGeoTransform[1] / 2;
    1170          22 :     double dfMaxX = adfGeoTransform[1] * (nXSize - 0.5) + adfGeoTransform[0];
    1171          22 :     double dfMinY = adfGeoTransform[5] * (nYSize - 0.5) + adfGeoTransform[3];
    1172          22 :     double dfMaxY = adfGeoTransform[3] + adfGeoTransform[5] / 2;
    1173          22 :     CPLErr eErr = WriteHeader(fp, nXSize, nYSize, dfMinX, dfMaxX, dfMinY,
    1174             :                               dfMaxY, 0.0, 0.0);
    1175             : 
    1176          22 :     if (eErr != CE_None)
    1177             :     {
    1178           2 :         VSIFCloseL(fp);
    1179           2 :         return nullptr;
    1180             :     }
    1181             : 
    1182             :     /* -------------------------------------------------------------------- */
    1183             :     /*      Copy band data.                                                 */
    1184             :     /* -------------------------------------------------------------------- */
    1185          20 :     double *pfData = (double *)VSI_MALLOC2_VERBOSE(nXSize, sizeof(double));
    1186          20 :     if (pfData == nullptr)
    1187             :     {
    1188           0 :         VSIFCloseL(fp);
    1189           0 :         return nullptr;
    1190             :     }
    1191             : 
    1192             :     int bSrcHasNDValue;
    1193          20 :     double dfSrcNoDataValue = poSrcBand->GetNoDataValue(&bSrcHasNDValue);
    1194          20 :     double dfMinZ = std::numeric_limits<double>::max();
    1195          20 :     double dfMaxZ = std::numeric_limits<double>::lowest();
    1196         186 :     for (GInt32 iRow = nYSize - 1; iRow >= 0; iRow--)
    1197             :     {
    1198         174 :         eErr = poSrcBand->RasterIO(GF_Read, 0, iRow, nXSize, 1, pfData, nXSize,
    1199             :                                    1, GDT_Float64, 0, 0, nullptr);
    1200             : 
    1201         174 :         if (eErr != CE_None)
    1202             :         {
    1203           0 :             VSIFCloseL(fp);
    1204           0 :             VSIFree(pfData);
    1205           0 :             return nullptr;
    1206             :         }
    1207             : 
    1208        2114 :         for (int iCol = 0; iCol < nXSize; iCol++)
    1209             :         {
    1210        1940 :             if (bSrcHasNDValue && pfData[iCol] == dfSrcNoDataValue)
    1211             :             {
    1212           0 :                 pfData[iCol] = dfDefaultNoDataValue;
    1213             :             }
    1214             :             else
    1215             :             {
    1216        1940 :                 if (pfData[iCol] > dfMaxZ)
    1217          22 :                     dfMaxZ = pfData[iCol];
    1218             : 
    1219        1940 :                 if (pfData[iCol] < dfMinZ)
    1220          27 :                     dfMinZ = pfData[iCol];
    1221             :             }
    1222             : 
    1223        1940 :             CPL_LSBPTR64(pfData + iCol);
    1224             :         }
    1225             : 
    1226         174 :         if (VSIFWriteL((void *)pfData, sizeof(double), nXSize, fp) !=
    1227         174 :             static_cast<unsigned>(nXSize))
    1228             :         {
    1229           8 :             VSIFCloseL(fp);
    1230           8 :             VSIFree(pfData);
    1231           8 :             CPLError(CE_Failure, CPLE_FileIO,
    1232             :                      "Unable to write grid row. Disk full?\n");
    1233           8 :             return nullptr;
    1234             :         }
    1235             : 
    1236         166 :         if (!pfnProgress(static_cast<double>(nYSize - iRow) / nYSize, nullptr,
    1237             :                          pProgressData))
    1238             :         {
    1239           0 :             VSIFCloseL(fp);
    1240           0 :             VSIFree(pfData);
    1241           0 :             CPLError(CE_Failure, CPLE_UserInterrupt, "User terminated");
    1242           0 :             return nullptr;
    1243             :         }
    1244             :     }
    1245             : 
    1246          12 :     VSIFree(pfData);
    1247             : 
    1248             :     /* write out the min and max values */
    1249          12 :     eErr = WriteHeader(fp, nXSize, nYSize, dfMinX, dfMaxX, dfMinY, dfMaxY,
    1250             :                        dfMinZ, dfMaxZ);
    1251             : 
    1252          12 :     if (eErr != CE_None)
    1253             :     {
    1254           0 :         VSIFCloseL(fp);
    1255           0 :         return nullptr;
    1256             :     }
    1257             : 
    1258          12 :     VSIFCloseL(fp);
    1259             : 
    1260          12 :     GDALPamDataset *poDS = (GDALPamDataset *)GDALOpen(pszFilename, GA_Update);
    1261          12 :     if (poDS)
    1262             :     {
    1263          12 :         poDS->CloneInfo(poSrcDS, GCIF_PAM_DEFAULT);
    1264             :     }
    1265             : 
    1266          12 :     return poDS;
    1267             : }
    1268             : 
    1269             : /************************************************************************/
    1270             : /*                          GDALRegister_GS7BG()                        */
    1271             : /************************************************************************/
    1272        1889 : void GDALRegister_GS7BG()
    1273             : 
    1274             : {
    1275        1889 :     if (GDALGetDriverByName("GS7BG") != nullptr)
    1276         282 :         return;
    1277             : 
    1278        1607 :     GDALDriver *poDriver = new GDALDriver();
    1279             : 
    1280        1607 :     poDriver->SetDescription("GS7BG");
    1281        1607 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    1282        1607 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
    1283        1607 :                               "Golden Software 7 Binary Grid (.grd)");
    1284        1607 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/gs7bg.html");
    1285        1607 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "grd");
    1286        1607 :     poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES,
    1287        1607 :                               "Byte Int16 UInt16 Float32 Float64");
    1288        1607 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
    1289             : 
    1290        1607 :     poDriver->pfnIdentify = GS7BGDataset::Identify;
    1291        1607 :     poDriver->pfnOpen = GS7BGDataset::Open;
    1292        1607 :     poDriver->pfnCreate = GS7BGDataset::Create;
    1293        1607 :     poDriver->pfnCreateCopy = GS7BGDataset::CreateCopy;
    1294             : 
    1295        1607 :     GetGDALDriverManager()->RegisterDriver(poDriver);
    1296             : }