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