Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GRD Reader
       4             :  * Purpose:  GDAL driver for Northwood Grid Format
       5             :  * Author:   Perry Casson
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2006, Waypoint Information Technology
       9             :  * Copyright (c) 2009-2012, Even Rouault <even dot rouault at spatialys.com>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include <string>
      15             : #include <cstring>
      16             : #include <cstdio>
      17             : #include <climits>
      18             : #include "gdal_frmts.h"
      19             : #include "gdal_pam.h"
      20             : #include "northwood.h"
      21             : #include "ogrmitabspatialref.h"
      22             : 
      23             : #ifndef NO_MITAB_SUPPORT
      24             : #ifdef MSVC
      25             : #include "..\..\ogr\ogrsf_frmts\mitab\mitab.h"
      26             : #else
      27             : #include "../../ogr/ogrsf_frmts/mitab/mitab.h"
      28             : #endif
      29             : #endif
      30             : 
      31             : constexpr float NODATA = -1.e37f;
      32             : constexpr double SCALE16BIT = 65534.0;
      33             : constexpr double SCALE32BIT = 4294967294.0;
      34             : 
      35             : void replaceExt(std::string &s, const std::string &newExt);
      36             : 
      37             : /************************************************************************/
      38             : /* Replace the extension on a filepath with an alternative extension    */
      39             : /************************************************************************/
      40           0 : void replaceExt(std::string &s, const std::string &newExt)
      41             : {
      42             : 
      43           0 :     std::string::size_type i = s.rfind('.', s.length());
      44             : 
      45           0 :     if (i != std::string::npos)
      46             :     {
      47           0 :         s.replace(i + 1, newExt.length(), newExt);
      48             :     }
      49           0 : }
      50             : 
      51             : /************************************************************************/
      52             : /* ==================================================================== */
      53             : /*                      NWT_GRDDataset                                  */
      54             : /* ==================================================================== */
      55             : /************************************************************************/
      56             : class NWT_GRDRasterBand;
      57             : 
      58             : class NWT_GRDDataset final : public GDALPamDataset
      59             : {
      60             :     friend class NWT_GRDRasterBand;
      61             : 
      62             :     VSILFILE *fp;
      63             :     GByte abyHeader[1024];
      64             :     NWT_GRID *pGrd;
      65             :     NWT_RGB ColorMap[4096];
      66             :     bool bUpdateHeader;
      67             :     mutable OGRSpatialReference *m_poSRS = nullptr;
      68             : 
      69             : #ifndef NO_MITAB_SUPPORT
      70             :     // Update the header data with latest changes
      71             :     int UpdateHeader();
      72             :     int WriteTab();
      73             : #endif
      74             : 
      75             :     NWT_GRDDataset(const NWT_GRDDataset &) = delete;
      76             :     NWT_GRDDataset &operator=(const NWT_GRDDataset &) = delete;
      77             : 
      78             :   public:
      79             :     NWT_GRDDataset();
      80             :     ~NWT_GRDDataset();
      81             : 
      82             :     static GDALDataset *Open(GDALOpenInfo *);
      83             :     static int Identify(GDALOpenInfo *);
      84             : 
      85             : #ifndef NO_MITAB_SUPPORT
      86             :     static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
      87             :                                int nBandsIn, GDALDataType eType,
      88             :                                char **papszParamList);
      89             :     static GDALDataset *CreateCopy(const char *pszFilename,
      90             :                                    GDALDataset *poSrcDS, int bStrict,
      91             :                                    char **papszOptions,
      92             :                                    GDALProgressFunc pfnProgress,
      93             :                                    void *pProgressData);
      94             : #endif
      95             : 
      96             :     CPLErr GetGeoTransform(double *padfTransform) override;
      97             :     CPLErr SetGeoTransform(double *padfTransform) override;
      98             :     CPLErr FlushCache(bool bAtClosing) override;
      99             : 
     100             :     const OGRSpatialReference *GetSpatialRef() const override;
     101             : 
     102             : #ifndef NO_MITAB_SUPPORT
     103             :     CPLErr SetSpatialRef(const OGRSpatialReference *poSRS) override;
     104             : #endif
     105             : };
     106             : 
     107             : /************************************************************************/
     108             : /* ==================================================================== */
     109             : /*                            NWT_GRDRasterBand                         */
     110             : /* ==================================================================== */
     111             : /************************************************************************/
     112             : 
     113             : class NWT_GRDRasterBand final : public GDALPamRasterBand
     114             : {
     115             :     friend class NWT_GRDDataset;
     116             : 
     117             :     int bHaveOffsetScale;
     118             :     double dfOffset;
     119             :     double dfScale;
     120             :     double dfNoData;
     121             : 
     122             :   public:
     123             :     NWT_GRDRasterBand(NWT_GRDDataset *, int, int);
     124             : 
     125             :     virtual CPLErr IReadBlock(int, int, void *) override;
     126             :     virtual CPLErr IWriteBlock(int, int, void *) override;
     127             :     virtual double GetNoDataValue(int *pbSuccess) override;
     128             :     virtual CPLErr SetNoDataValue(double dfNoData) override;
     129             : 
     130             :     virtual GDALColorInterp GetColorInterpretation() override;
     131             : };
     132             : 
     133             : /************************************************************************/
     134             : /*                           NWT_GRDRasterBand()                        */
     135             : /************************************************************************/
     136          45 : NWT_GRDRasterBand::NWT_GRDRasterBand(NWT_GRDDataset *poDSIn, int nBandIn,
     137          45 :                                      int nBands)
     138          45 :     : bHaveOffsetScale(FALSE), dfOffset(0.0), dfScale(1.0), dfNoData(0.0)
     139             : {
     140          45 :     poDS = poDSIn;
     141          45 :     nBand = nBandIn;
     142             : 
     143             :     /*
     144             :      * If nBand = 4 we have opened in read mode and have created the 3 'virtual'
     145             :      * RGB bands. so the 4th band is the actual data Otherwise, if we have
     146             :      * opened in update mode, there is only 1 band, which is the actual data
     147             :      */
     148          45 :     if (nBand == 4 || nBands == 1)
     149             :     {
     150          15 :         bHaveOffsetScale = TRUE;
     151          15 :         dfOffset = poDSIn->pGrd->fZMin;
     152             : 
     153          15 :         if (poDSIn->pGrd->cFormat == 0x00)
     154             :         {
     155          15 :             eDataType = GDT_Float32;
     156          15 :             dfScale = (poDSIn->pGrd->fZMax - poDSIn->pGrd->fZMin) / SCALE16BIT;
     157             :         }
     158             :         else
     159             :         {
     160           0 :             eDataType = GDT_Float32;
     161           0 :             dfScale = (poDSIn->pGrd->fZMax - poDSIn->pGrd->fZMin) / SCALE32BIT;
     162             :         }
     163             :     }
     164             :     else
     165             :     {
     166          30 :         bHaveOffsetScale = FALSE;
     167          30 :         dfOffset = 0;
     168          30 :         dfScale = 1.0;
     169          30 :         eDataType = GDT_Byte;
     170             :     }
     171          45 :     nBlockXSize = poDS->GetRasterXSize();
     172          45 :     nBlockYSize = 1;
     173          45 : }
     174             : 
     175         146 : double NWT_GRDRasterBand::GetNoDataValue(int *pbSuccess)
     176             : {
     177         146 :     NWT_GRDDataset *poGDS = cpl::down_cast<NWT_GRDDataset *>(poDS);
     178             :     double dRetval;
     179         146 :     if ((nBand == 4) || (poGDS->nBands == 1))
     180             :     {
     181         146 :         if (pbSuccess != nullptr)
     182         145 :             *pbSuccess = TRUE;
     183         146 :         if (dfNoData != 0.0)
     184             :         {
     185           0 :             dRetval = dfNoData;
     186             :         }
     187             :         else
     188             :         {
     189         146 :             dRetval = NODATA;
     190             :         }
     191             : 
     192         146 :         return dRetval;
     193             :     }
     194             : 
     195           0 :     if (pbSuccess != nullptr)
     196           0 :         *pbSuccess = FALSE;
     197             : 
     198           0 :     return 0;
     199             : }
     200             : 
     201           1 : CPLErr NWT_GRDRasterBand::SetNoDataValue(double dfNoDataIn)
     202             : {
     203             :     // This is essentially a 'virtual' no data value.
     204             :     // Once set, when writing an value == dfNoData will
     205             :     // be converted to the no data value (0) on disk.
     206             :     // If opened again; the no data value will always be the
     207             :     // default (-1.e37f)
     208           1 :     dfNoData = dfNoDataIn;
     209           1 :     return CE_None;
     210             : }
     211             : 
     212           3 : GDALColorInterp NWT_GRDRasterBand::GetColorInterpretation()
     213             : {
     214           3 :     NWT_GRDDataset *poGDS = cpl::down_cast<NWT_GRDDataset *>(poDS);
     215             :     // return GCI_RGB;
     216           3 :     if ((nBand == 4) || (poGDS->nBands == 1))
     217           3 :         return GCI_GrayIndex;
     218           0 :     else if (nBand == 1)
     219           0 :         return GCI_RedBand;
     220           0 :     else if (nBand == 2)
     221           0 :         return GCI_GreenBand;
     222           0 :     else if (nBand == 3)
     223           0 :         return GCI_BlueBand;
     224             : 
     225           0 :     return GCI_Undefined;
     226             : }
     227             : 
     228             : /************************************************************************/
     229             : /*                             IWriteBlock()                            */
     230             : /************************************************************************/
     231          48 : CPLErr NWT_GRDRasterBand::IWriteBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
     232             :                                       void *pImage)
     233             : {
     234             : 
     235             :     // Each block is an entire row of the dataset, so the x offset should always
     236             :     // be 0
     237          48 :     CPLAssert(nBlockXOff == 0);
     238          48 :     NWT_GRDDataset *poGDS = cpl::down_cast<NWT_GRDDataset *>(poDS);
     239             : 
     240          48 :     if (dfScale == 0.0)
     241           1 :         return CE_Failure;
     242             : 
     243             :     // Ensure the blocksize is not beyond the system limits and
     244             :     // initialize the size of the record
     245          47 :     if (nBlockXSize > INT_MAX / 2)
     246             :     {
     247           0 :         return CE_Failure;
     248             :     }
     249          47 :     const int nRecordSize = nBlockXSize * 2;
     250             : 
     251             :     // Seek to the write position in the GRD file
     252          47 :     VSIFSeekL(poGDS->fp,
     253          47 :               1024 + nRecordSize * static_cast<vsi_l_offset>(nBlockYOff),
     254             :               SEEK_SET);
     255             : 
     256             :     // Cast pImage to float
     257          47 :     const float *pfImage = static_cast<const float *>(pImage);
     258             : 
     259             :     // Initialize output array
     260          47 :     GByte *pabyRecord = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nRecordSize));
     261          47 :     if (pabyRecord == nullptr)
     262           0 :         return CE_Failure;
     263             : 
     264             :     // We only ever write to band 4; RGB bands are basically 'virtual'
     265             :     // (i.e. the RGB colour is computed from the raw data).
     266             :     // For all intents and purposes, there is essentially 1 band on disk.
     267          47 :     if (nBand == 1)
     268             :     {
     269        3008 :         for (int i = 0; i < nBlockXSize; i++)
     270             :         {
     271        2961 :             const float fValue = pfImage[i];
     272             :             unsigned short nWrite;  // The stretched value to be written
     273             : 
     274             :             // We allow data to be interpreted by a user-defined null value
     275             :             // (a 'virtual' value, since it is always 0 on disk) or
     276             :             // if not defined we default to the GRD standard of -1E37.
     277             :             // We allow a little bit of flexibility in that if it is below -1E37
     278             :             // it is in all probability still intended as a null value.
     279        2961 :             if ((fValue == dfNoData) || (fValue <= NODATA))
     280             :             {
     281         381 :                 nWrite = 0;
     282             :             }
     283             :             else
     284             :             {
     285        2580 :                 if (fValue < poGDS->pGrd->fZMin)
     286             :                 {
     287           0 :                     poGDS->pGrd->fZMin = fValue;
     288             :                 }
     289        2580 :                 else if (fValue > poGDS->pGrd->fZMax)
     290             :                 {
     291           0 :                     poGDS->pGrd->fZMax = fValue;
     292             :                 }
     293             :                 // Data on disk is stretched within the unsigned short range so
     294             :                 // we must convert (the inverse of what is done in IReadBlock),
     295             :                 // based on the Z value range
     296        2580 :                 nWrite = static_cast<unsigned short>(
     297        2580 :                     ((fValue - dfOffset) / dfScale) + 1);
     298             :             }
     299        2961 :             CPL_LSBPTR16(&nWrite);
     300             :             // Copy the result to the byte array (2 bytes per value)
     301        2961 :             memcpy(pabyRecord + 2 * i, &nWrite, 2);
     302             :         }
     303             : 
     304             :         // Write the buffer to disk
     305          47 :         if (VSIFWriteL(pabyRecord, 1, nRecordSize, poGDS->fp) !=
     306          47 :             static_cast<size_t>(nRecordSize))
     307             :         {
     308           0 :             CPLError(CE_Failure, CPLE_FileIO,
     309             :                      "Failed to write scanline %d to file.\n", nBlockYOff);
     310           0 :             CPLFree(pabyRecord);
     311           0 :             return CE_Failure;
     312             :         }
     313             :     }
     314             :     else
     315             :     {
     316           0 :         CPLError(CE_Failure, CPLE_IllegalArg, "Writing to band %d is not valid",
     317             :                  nBand);
     318           0 :         CPLFree(pabyRecord);
     319           0 :         return CE_Failure;
     320             :     }
     321          47 :     CPLFree(pabyRecord);
     322          47 :     return CE_None;
     323             : }
     324             : 
     325             : /************************************************************************/
     326             : /*                             IReadBlock()                             */
     327             : /************************************************************************/
     328         282 : CPLErr NWT_GRDRasterBand::IReadBlock(CPL_UNUSED int nBlockXOff, int nBlockYOff,
     329             :                                      void *pImage)
     330             : {
     331         282 :     NWT_GRDDataset *poGDS = cpl::down_cast<NWT_GRDDataset *>(poDS);
     332         282 :     if (nBlockXSize > INT_MAX / 2)
     333           0 :         return CE_Failure;
     334         282 :     const int nRecordSize = nBlockXSize * 2;
     335             : 
     336             :     // Seek to the data position
     337         282 :     VSIFSeekL(poGDS->fp,
     338         282 :               1024 + nRecordSize * static_cast<vsi_l_offset>(nBlockYOff),
     339             :               SEEK_SET);
     340             : 
     341         282 :     GByte *pabyRecord = static_cast<GByte *>(VSI_MALLOC_VERBOSE(nRecordSize));
     342         282 :     if (pabyRecord == nullptr)
     343           0 :         return CE_Failure;
     344             : 
     345             :     // Read the data
     346         282 :     if (static_cast<int>(VSIFReadL(pabyRecord, 1, nRecordSize, poGDS->fp)) !=
     347             :         nRecordSize)
     348             :     {
     349           0 :         CPLFree(pabyRecord);
     350           0 :         return CE_Failure;
     351             :     }
     352             : 
     353         282 :     if ((nBand == 4) || (poGDS->nBands == 1))  // Z values
     354             :     {
     355             :         int bSuccess;
     356         141 :         const float fNoData = static_cast<float>(GetNoDataValue(&bSuccess));
     357        9024 :         for (int i = 0; i < nBlockXSize; i++)
     358             :         {
     359             :             unsigned short raw1;
     360        8883 :             memcpy(&raw1, pabyRecord + 2 * i, 2);
     361        8883 :             CPL_LSBPTR16(&raw1);
     362        8883 :             if (raw1 == 0)
     363             :             {
     364        1143 :                 static_cast<float *>(pImage)[i] = fNoData;  // null value
     365             :             }
     366             :             else
     367             :             {
     368        7740 :                 static_cast<float *>(pImage)[i] =
     369        7740 :                     static_cast<float>(dfOffset + ((raw1 - 1) * dfScale));
     370             :             }
     371         141 :         }
     372             :     }
     373         141 :     else if (nBand == 1)  // red values
     374             :     {
     375        3008 :         for (int i = 0; i < nBlockXSize; i++)
     376             :         {
     377             :             unsigned short raw1;
     378        2961 :             memcpy(&raw1, pabyRecord + 2 * i, 2);
     379        2961 :             CPL_LSBPTR16(&raw1);
     380        2961 :             static_cast<unsigned char *>(pImage)[i] =
     381        2961 :                 poGDS->ColorMap[raw1 / 16].r;
     382             :         }
     383             :     }
     384          94 :     else if (nBand == 2)  // green
     385             :     {
     386        3008 :         for (int i = 0; i < nBlockXSize; i++)
     387             :         {
     388             :             unsigned short raw1;
     389        2961 :             memcpy(&raw1, pabyRecord + 2 * i, 2);
     390        2961 :             CPL_LSBPTR16(&raw1);
     391        2961 :             static_cast<unsigned char *>(pImage)[i] =
     392        2961 :                 poGDS->ColorMap[raw1 / 16].g;
     393             :         }
     394             :     }
     395          47 :     else if (nBand == 3)  // blue
     396             :     {
     397        3008 :         for (int i = 0; i < nBlockXSize; i++)
     398             :         {
     399             :             unsigned short raw1;
     400        2961 :             memcpy(&raw1, pabyRecord + 2 * i, 2);
     401        2961 :             CPL_LSBPTR16(&raw1);
     402        2961 :             static_cast<unsigned char *>(pImage)[i] =
     403        2961 :                 poGDS->ColorMap[raw1 / 16].b;
     404             :         }
     405             :     }
     406             :     else
     407             :     {
     408           0 :         CPLError(CE_Failure, CPLE_IllegalArg, "No band number %d", nBand);
     409           0 :         CPLFree(pabyRecord);
     410           0 :         return CE_Failure;
     411             :     }
     412             : 
     413         282 :     CPLFree(pabyRecord);
     414             : 
     415         282 :     return CE_None;
     416             : }
     417             : 
     418             : /************************************************************************/
     419             : /* ==================================================================== */
     420             : /*                             NWT_GRDDataset                           */
     421             : /* ==================================================================== */
     422             : /************************************************************************/
     423             : 
     424          26 : NWT_GRDDataset::NWT_GRDDataset()
     425          26 :     : fp(nullptr), pGrd(nullptr), bUpdateHeader(false)
     426             : {
     427             :     // poCT = NULL;
     428      106522 :     for (size_t i = 0; i < CPL_ARRAYSIZE(ColorMap); ++i)
     429             :     {
     430      106496 :         ColorMap[i].r = 0;
     431      106496 :         ColorMap[i].g = 0;
     432      106496 :         ColorMap[i].b = 0;
     433             :     }
     434          26 : }
     435             : 
     436             : /************************************************************************/
     437             : /*                            ~NWT_GRDDataset()                         */
     438             : /************************************************************************/
     439             : 
     440          52 : NWT_GRDDataset::~NWT_GRDDataset()
     441             : {
     442             : 
     443             :     // Make sure any changes to the header etc are written
     444             :     // if we are in update mode.
     445          26 :     if (eAccess == GA_Update)
     446             :     {
     447          14 :         NWT_GRDDataset::FlushCache(true);
     448             :     }
     449          26 :     pGrd->fp = nullptr;  // this prevents nwtCloseGrid from closing the fp
     450          26 :     nwtCloseGrid(pGrd);
     451          26 :     if (m_poSRS)
     452           1 :         m_poSRS->Release();
     453             : 
     454          26 :     if (fp != nullptr)
     455          25 :         VSIFCloseL(fp);
     456          52 : }
     457             : 
     458             : /************************************************************************/
     459             : /*                 ~FlushCache(bool bAtClosing)                         */
     460             : /************************************************************************/
     461          17 : CPLErr NWT_GRDDataset::FlushCache(bool bAtClosing)
     462             : {
     463             :     // Ensure the header and TAB file are up to date
     464          17 :     if (bUpdateHeader)
     465             :     {
     466             : #ifndef NO_MITAB_SUPPORT
     467           3 :         UpdateHeader();
     468             : #endif
     469             :     }
     470             : 
     471             :     // Call the parent method
     472          17 :     return GDALPamDataset::FlushCache(bAtClosing);
     473             : }
     474             : 
     475             : /************************************************************************/
     476             : /*                          GetGeoTransform()                           */
     477             : /************************************************************************/
     478             : 
     479           3 : CPLErr NWT_GRDDataset::GetGeoTransform(double *padfTransform)
     480             : {
     481           3 :     padfTransform[0] = pGrd->dfMinX - (pGrd->dfStepSize * 0.5);
     482           3 :     padfTransform[3] = pGrd->dfMaxY + (pGrd->dfStepSize * 0.5);
     483           3 :     padfTransform[1] = pGrd->dfStepSize;
     484           3 :     padfTransform[2] = 0.0;
     485             : 
     486           3 :     padfTransform[4] = 0.0;
     487           3 :     padfTransform[5] = -1 * pGrd->dfStepSize;
     488             : 
     489           3 :     return CE_None;
     490             : }
     491             : 
     492             : /************************************************************************/
     493             : /*                          SetGeoTransform()                           */
     494             : /************************************************************************/
     495             : 
     496           3 : CPLErr NWT_GRDDataset::SetGeoTransform(double *padfTransform)
     497             : {
     498           3 :     if (padfTransform[2] != 0.0 || padfTransform[4] != 0.0)
     499             :     {
     500             : 
     501           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     502             :                  "GRD datasets do not support skew/rotation");
     503           0 :         return CE_Failure;
     504             :     }
     505           3 :     pGrd->dfStepSize = padfTransform[1];
     506             : 
     507             :     // GRD format sets the min/max coordinates to the centre of the
     508             :     // cell; We must account for this when copying the GDAL geotransform
     509             :     // which references the top left corner
     510           3 :     pGrd->dfMinX = padfTransform[0] + (pGrd->dfStepSize * 0.5);
     511           3 :     pGrd->dfMaxY = padfTransform[3] - (pGrd->dfStepSize * 0.5);
     512             : 
     513             :     // Now set the miny and maxx
     514           3 :     pGrd->dfMaxX = pGrd->dfMinX + (pGrd->dfStepSize * (nRasterXSize - 1));
     515           3 :     pGrd->dfMinY = pGrd->dfMaxY - (pGrd->dfStepSize * (nRasterYSize - 1));
     516           3 :     bUpdateHeader = true;
     517             : 
     518           3 :     return CE_None;
     519             : }
     520             : 
     521             : /************************************************************************/
     522             : /*                          GetSpatialRef()                             */
     523             : /************************************************************************/
     524           1 : const OGRSpatialReference *NWT_GRDDataset::GetSpatialRef() const
     525             : {
     526             : 
     527             :     // First try getting it from the PAM dataset
     528           1 :     const OGRSpatialReference *poSRS = GDALPamDataset::GetSpatialRef();
     529           1 :     if (poSRS)
     530           0 :         return poSRS;
     531             : 
     532           1 :     if (m_poSRS)
     533           0 :         return m_poSRS;
     534             : 
     535             :     // If that isn't possible, read it from the GRD file. This may be a less
     536             :     //  complete projection string.
     537             :     OGRSpatialReference *poSpatialRef =
     538           1 :         MITABCoordSys2SpatialRef(pGrd->cMICoordSys);
     539           1 :     m_poSRS = poSpatialRef;
     540             : 
     541           1 :     return m_poSRS;
     542             : }
     543             : 
     544             : #ifndef NO_MITAB_SUPPORT
     545             : /************************************************************************/
     546             : /*                            SetSpatialRef()                           */
     547             : /************************************************************************/
     548             : 
     549           3 : CPLErr NWT_GRDDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
     550             : {
     551             : 
     552           3 :     char *psTABProj = MITABSpatialRef2CoordSys(poSRS);
     553           3 :     strncpy(pGrd->cMICoordSys, psTABProj, sizeof(pGrd->cMICoordSys) - 1);
     554           3 :     pGrd->cMICoordSys[255] = '\0';
     555             : 
     556             :     // Free temp projection.
     557           3 :     CPLFree(psTABProj);
     558             :     // Set projection in PAM dataset, so that
     559             :     // GDAL can always retrieve the complete projection.
     560           3 :     GDALPamDataset::SetSpatialRef(poSRS);
     561           3 :     bUpdateHeader = true;
     562             : 
     563           3 :     return CE_None;
     564             : }
     565             : #endif
     566             : 
     567             : /************************************************************************/
     568             : /*                              Identify()                              */
     569             : /************************************************************************/
     570             : 
     571       50898 : int NWT_GRDDataset::Identify(GDALOpenInfo *poOpenInfo)
     572             : {
     573             :     /* -------------------------------------------------------------------- */
     574             :     /*  Look for the header                                                 */
     575             :     /* -------------------------------------------------------------------- */
     576       50898 :     if (poOpenInfo->nHeaderBytes < 1024)
     577       49168 :         return FALSE;
     578             : 
     579        1730 :     if (poOpenInfo->pabyHeader[0] != 'H' || poOpenInfo->pabyHeader[1] != 'G' ||
     580          27 :         poOpenInfo->pabyHeader[2] != 'P' || poOpenInfo->pabyHeader[3] != 'C' ||
     581          27 :         poOpenInfo->pabyHeader[4] != '1')
     582        1706 :         return FALSE;
     583             : 
     584          24 :     return TRUE;
     585             : }
     586             : 
     587             : /************************************************************************/
     588             : /*                                Open()                                */
     589             : /************************************************************************/
     590          12 : GDALDataset *NWT_GRDDataset::Open(GDALOpenInfo *poOpenInfo)
     591             : {
     592          12 :     if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
     593           0 :         return nullptr;
     594             : 
     595             :     /* -------------------------------------------------------------------- */
     596             :     /*      Create a corresponding GDALDataset.                             */
     597             :     /* -------------------------------------------------------------------- */
     598          12 :     int nBandsToCreate = 0;
     599             : 
     600          12 :     NWT_GRDDataset *poDS = new NWT_GRDDataset();
     601             : 
     602          12 :     poDS->fp = poOpenInfo->fpL;
     603          12 :     poOpenInfo->fpL = nullptr;
     604             : 
     605          12 :     if (poOpenInfo->eAccess == GA_Update)
     606             :     {
     607           0 :         nBandsToCreate = 1;
     608             :     }
     609             :     else
     610             :     {
     611          12 :         nBandsToCreate = atoi(CSLFetchNameValueDef(poOpenInfo->papszOpenOptions,
     612             :                                                    "BAND_COUNT", "4"));
     613          12 :         if (nBandsToCreate != 1 && nBandsToCreate != 4)
     614             :         {
     615           0 :             CPLError(CE_Failure, CPLE_AppDefined, "Wrong value for BAND_COUNT");
     616           0 :             delete poDS;
     617           0 :             return nullptr;
     618             :         }
     619             :     }
     620          12 :     poDS->eAccess = poOpenInfo->eAccess;
     621             : 
     622             :     /* -------------------------------------------------------------------- */
     623             :     /*      Read the header.                                                */
     624             :     /* -------------------------------------------------------------------- */
     625          12 :     VSIFSeekL(poDS->fp, 0, SEEK_SET);
     626          12 :     VSIFReadL(poDS->abyHeader, 1, 1024, poDS->fp);
     627          12 :     poDS->pGrd = reinterpret_cast<NWT_GRID *>(calloc(1, sizeof(NWT_GRID)));
     628             : 
     629          12 :     poDS->pGrd->fp = poDS->fp;
     630             : 
     631          24 :     if (!nwt_ParseHeader(poDS->pGrd, poDS->abyHeader) ||
     632          12 :         !GDALCheckDatasetDimensions(poDS->pGrd->nXSide, poDS->pGrd->nYSide))
     633             :     {
     634           0 :         delete poDS;
     635           0 :         return nullptr;
     636             :     }
     637             : 
     638          12 :     poDS->nRasterXSize = poDS->pGrd->nXSide;
     639          12 :     poDS->nRasterYSize = poDS->pGrd->nYSide;
     640             : 
     641             :     // create a colorTable
     642             :     // if( poDS->pGrd->iNumColorInflections > 0 )
     643             :     //   poDS->CreateColorTable();
     644          12 :     nwt_LoadColors(poDS->ColorMap, 4096, poDS->pGrd);
     645             :     /* -------------------------------------------------------------------- */
     646             :     /*      Create band information objects.                                */
     647             :     /* If opening in read-only mode, then we create 4 bands (RGBZ)          */
     648             :     /* with data values being available in band 4. If opening in update mode*/
     649             :     /* we create 1 band (the data values). This is because in reality, there*/
     650             :     /* is only 1 band stored on disk. The RGB bands are 'virtual' - derived */
     651             :     /* from the data values on the fly                                      */
     652             :     /* -------------------------------------------------------------------- */
     653          54 :     for (int i = 0; i < nBandsToCreate; ++i)
     654             :     {
     655          42 :         poDS->SetBand(i + 1,
     656          42 :                       new NWT_GRDRasterBand(poDS, i + 1, nBandsToCreate));
     657             :     }
     658             : 
     659             :     /* -------------------------------------------------------------------- */
     660             :     /*      Initialize any PAM information.                                 */
     661             :     /* -------------------------------------------------------------------- */
     662          12 :     poDS->SetDescription(poOpenInfo->pszFilename);
     663          12 :     poDS->TryLoadXML();
     664             : 
     665             :     /* -------------------------------------------------------------------- */
     666             :     /*      Check for external overviews.                                   */
     667             :     /* -------------------------------------------------------------------- */
     668          24 :     poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename,
     669          12 :                                 poOpenInfo->GetSiblingFiles());
     670             : 
     671          12 :     return poDS;
     672             : }
     673             : 
     674             : #ifndef NO_MITAB_SUPPORT
     675             : /************************************************************************/
     676             : /*                                UpdateHeader()                        */
     677             : /************************************************************************/
     678          16 : int NWT_GRDDataset::UpdateHeader()
     679             : {
     680          16 :     int iStatus = 0;
     681          16 :     TABRawBinBlock *poHeaderBlock = new TABRawBinBlock(TABReadWrite, TRUE);
     682          16 :     poHeaderBlock->InitNewBlock(fp, 1024);
     683             : 
     684             :     // Write the header string
     685          16 :     poHeaderBlock->WriteBytes(5, reinterpret_cast<const GByte *>("HGPC1\0"));
     686             : 
     687             :     // Version number
     688          16 :     poHeaderBlock->WriteFloat(pGrd->fVersion);
     689             : 
     690             :     // Dimensions
     691          16 :     poHeaderBlock->WriteInt16(static_cast<GInt16>(pGrd->nXSide));
     692          16 :     poHeaderBlock->WriteInt16(static_cast<GInt16>(pGrd->nYSide));
     693             : 
     694             :     // Extents
     695          16 :     poHeaderBlock->WriteDouble(pGrd->dfMinX);
     696          16 :     poHeaderBlock->WriteDouble(pGrd->dfMaxX);
     697          16 :     poHeaderBlock->WriteDouble(pGrd->dfMinY);
     698          16 :     poHeaderBlock->WriteDouble(pGrd->dfMaxY);
     699             : 
     700             :     // Z value range
     701          16 :     poHeaderBlock->WriteFloat(pGrd->fZMin);
     702          16 :     poHeaderBlock->WriteFloat(pGrd->fZMax);
     703          16 :     poHeaderBlock->WriteFloat(pGrd->fZMinScale);
     704          16 :     poHeaderBlock->WriteFloat(pGrd->fZMaxScale);
     705             : 
     706             :     // Description String
     707          16 :     int nChar = static_cast<int>(strlen(pGrd->cDescription));
     708          16 :     poHeaderBlock->WriteBytes(
     709          16 :         nChar, reinterpret_cast<const GByte *>(pGrd->cDescription));
     710          16 :     poHeaderBlock->WriteZeros(32 - nChar);
     711             : 
     712             :     // Unit Name String
     713          16 :     nChar = static_cast<int>(strlen(pGrd->cZUnits));
     714          16 :     poHeaderBlock->WriteBytes(nChar,
     715          16 :                               reinterpret_cast<const GByte *>(pGrd->cZUnits));
     716          16 :     poHeaderBlock->WriteZeros(32 - nChar);
     717             : 
     718             :     // Ignore 126 - 141 as unknown usage
     719          16 :     poHeaderBlock->WriteZeros(15);
     720             : 
     721             :     // Hill shading
     722          16 :     poHeaderBlock->WriteInt16(pGrd->bHillShadeExists ? 1 : 0);
     723          16 :     poHeaderBlock->WriteInt16(0);
     724             : 
     725          16 :     poHeaderBlock->WriteByte(pGrd->cHillShadeBrightness);
     726          16 :     poHeaderBlock->WriteByte(pGrd->cHillShadeContrast);
     727             : 
     728             :     // Ignore 147 - 257 as unknown usage
     729          16 :     poHeaderBlock->WriteZeros(110);
     730             : 
     731             :     // Write spatial reference
     732          16 :     poHeaderBlock->WriteBytes(
     733          16 :         static_cast<int>(strlen(pGrd->cMICoordSys)),
     734          16 :         reinterpret_cast<const GByte *>(pGrd->cMICoordSys));
     735          16 :     poHeaderBlock->WriteZeros(256 -
     736          16 :                               static_cast<int>(strlen(pGrd->cMICoordSys)));
     737             : 
     738             :     // Unit code
     739          16 :     poHeaderBlock->WriteByte(static_cast<GByte>(pGrd->iZUnits));
     740             : 
     741             :     // Info on shading
     742          16 :     GByte byDisplayStatus = 0;
     743          16 :     if (pGrd->bShowHillShade)
     744             :     {
     745           0 :         byDisplayStatus |= 1 << 6;
     746             :     }
     747          16 :     if (pGrd->bShowGradient)
     748             :     {
     749           0 :         byDisplayStatus |= 1 << 7;
     750             :     }
     751             : 
     752          16 :     poHeaderBlock->WriteByte(byDisplayStatus);
     753          16 :     poHeaderBlock->WriteInt16(0);  // Data Type?
     754             : 
     755             :     // Colour inflections
     756          16 :     poHeaderBlock->WriteInt16(pGrd->iNumColorInflections);
     757          64 :     for (int i = 0; i < pGrd->iNumColorInflections; i++)
     758             :     {
     759          48 :         poHeaderBlock->WriteFloat(pGrd->stInflection[i].zVal);
     760          48 :         poHeaderBlock->WriteByte(pGrd->stInflection[i].r);
     761          48 :         poHeaderBlock->WriteByte(pGrd->stInflection[i].g);
     762          48 :         poHeaderBlock->WriteByte(pGrd->stInflection[i].b);
     763             :     }
     764             : 
     765             :     // Fill in unused blanks
     766          16 :     poHeaderBlock->WriteZeros((966 - poHeaderBlock->GetCurAddress()));
     767             : 
     768             :     // Azimuth and Inclination
     769          16 :     poHeaderBlock->WriteFloat(pGrd->fHillShadeAzimuth);
     770          16 :     poHeaderBlock->WriteFloat(pGrd->fHillShadeAngle);
     771             : 
     772             :     // Write to disk
     773          16 :     iStatus = poHeaderBlock->CommitToFile();
     774             : 
     775          16 :     delete poHeaderBlock;
     776             : 
     777             :     // Update the TAB file to catch any changes
     778          16 :     if (WriteTab() != 0)
     779          10 :         iStatus = -1;
     780             : 
     781          16 :     return iStatus;
     782             : }
     783             : 
     784          16 : int NWT_GRDDataset::WriteTab()
     785             : {
     786             :     // Create the filename for the .tab file.
     787          32 :     const std::string sTabFile(CPLResetExtension(pGrd->szFileName, "tab"));
     788             : 
     789          16 :     VSILFILE *tabfp = VSIFOpenL(sTabFile.c_str(), "wt");
     790          16 :     if (tabfp == nullptr)
     791             :     {
     792           0 :         CPLError(CE_Failure, CPLE_FileIO, "Failed to create file `%s'",
     793             :                  sTabFile.c_str());
     794           0 :         return -1;
     795             :     }
     796             : 
     797          16 :     bool bOK = true;
     798          16 :     bOK &= VSIFPrintfL(tabfp, "!table\n") > 0;
     799          16 :     bOK &= VSIFPrintfL(tabfp, "!version 500\n") > 0;
     800          16 :     bOK &= VSIFPrintfL(tabfp, "!charset %s\n", "Neutral") > 0;
     801          16 :     bOK &= VSIFPrintfL(tabfp, "\n") > 0;
     802             : 
     803          16 :     bOK &= VSIFPrintfL(tabfp, "Definition Table\n") > 0;
     804          32 :     const std::string path(pGrd->szFileName);
     805          16 :     const std::string basename = path.substr(path.find_last_of("/\\") + 1);
     806          16 :     bOK &= VSIFPrintfL(tabfp, "  File \"%s\"\n", basename.c_str()) > 0;
     807          16 :     bOK &= VSIFPrintfL(tabfp, "  Type \"RASTER\"\n") > 0;
     808             : 
     809          16 :     double dMapUnitsPerPixel =
     810          16 :         (pGrd->dfMaxX - pGrd->dfMinX) / (static_cast<double>(pGrd->nXSide) - 1);
     811          16 :     double dShift = dMapUnitsPerPixel / 2.0;
     812             : 
     813          32 :     bOK &= VSIFPrintfL(tabfp, "  (%f,%f) (%d,%d) Label \"Pt 1\",\n",
     814          16 :                        pGrd->dfMinX - dShift, pGrd->dfMaxY + dShift, 0, 0) > 0;
     815          32 :     bOK &= VSIFPrintfL(tabfp, "  (%f,%f) (%d,%d) Label \"Pt 2\",\n",
     816          16 :                        pGrd->dfMaxX - dShift, pGrd->dfMinY + dShift,
     817          16 :                        pGrd->nXSide - 1, pGrd->nYSide - 1) > 0;
     818          32 :     bOK &= VSIFPrintfL(tabfp, "  (%f,%f) (%d,%d) Label \"Pt 3\"\n",
     819          16 :                        pGrd->dfMinX - dShift, pGrd->dfMinY + dShift, 0,
     820          16 :                        pGrd->nYSide - 1) > 0;
     821             : 
     822          16 :     bOK &= VSIFPrintfL(tabfp, "  CoordSys %s\n", pGrd->cMICoordSys) > 0;
     823          16 :     bOK &= VSIFPrintfL(tabfp, "  Units \"m\"\n") > 0;
     824             : 
     825             :     // Raster Styles.
     826             : 
     827             :     // Raster is a grid, which is style 6.
     828          16 :     bOK &= VSIFPrintfL(tabfp, "  RasterStyle 6 1\n") > 0;
     829             : 
     830             :     // Brightness - style 1
     831          16 :     if (pGrd->style.iBrightness > 0)
     832             :     {
     833          32 :         bOK &= VSIFPrintfL(tabfp, "  RasterStyle 1 %d\n",
     834          16 :                            pGrd->style.iBrightness) > 0;
     835             :     }
     836             : 
     837             :     // Contrast - style 2
     838          16 :     if (pGrd->style.iContrast > 0)
     839             :     {
     840          32 :         bOK &= VSIFPrintfL(tabfp, "  RasterStyle 2 %d\n",
     841          16 :                            pGrd->style.iContrast) > 0;
     842             :     }
     843             : 
     844             :     // Greyscale - style 3; only need to write if TRUE
     845          16 :     if (pGrd->style.bGreyscale == TRUE)
     846             :     {
     847           0 :         bOK &= VSIFPrintfL(tabfp, "  RasterStyle 3 1\n") > 0;
     848             :     }
     849             : 
     850             :     // Flag to render one colour transparent - style 4
     851          16 :     if (pGrd->style.bTransparent == TRUE)
     852             :     {
     853           0 :         bOK &= VSIFPrintfL(tabfp, "  RasterStyle 4 1\n") > 0;
     854           0 :         if (pGrd->style.iTransColour > 0)
     855             :         {
     856           0 :             bOK &= VSIFPrintfL(tabfp, "  RasterStyle 7 %d\n",
     857           0 :                                pGrd->style.iTransColour) > 0;
     858             :         }
     859             :     }
     860             : 
     861             :     // Transparency of immage
     862          16 :     if (pGrd->style.iTranslucency > 0)
     863             :     {
     864           0 :         bOK &= VSIFPrintfL(tabfp, "  RasterStyle 8 %d\n",
     865           0 :                            pGrd->style.iTranslucency) > 0;
     866             :     }
     867             : 
     868          16 :     bOK &= VSIFPrintfL(tabfp, "begin_metadata\n") > 0;
     869          16 :     bOK &= VSIFPrintfL(tabfp, "\"\\MapInfo\" = \"\"\n") > 0;
     870          16 :     bOK &= VSIFPrintfL(tabfp, "\"\\Vm\" = \"\"\n") > 0;
     871          16 :     bOK &= VSIFPrintfL(tabfp, "\"\\Vm\\Grid\" = \"Numeric\"\n") > 0;
     872          16 :     bOK &= VSIFPrintfL(tabfp, "\"\\Vm\\GridName\" = \"%s\"\n",
     873          16 :                        basename.c_str()) > 0;
     874          16 :     bOK &= VSIFPrintfL(tabfp, "\"\\IsReadOnly\" = \"FALSE\"\n") > 0;
     875          16 :     bOK &= VSIFPrintfL(tabfp, "end_metadata\n") > 0;
     876             : 
     877          16 :     if (VSIFCloseL(tabfp) != 0)
     878           0 :         bOK = false;
     879             : 
     880          16 :     return (bOK) ? 0 : -1;
     881             : }
     882             : 
     883             : /************************************************************************/
     884             : /*                                Create()                              */
     885             : /************************************************************************/
     886          57 : GDALDataset *NWT_GRDDataset::Create(const char *pszFilename, int nXSize,
     887             :                                     int nYSize, int nBandsIn,
     888             :                                     GDALDataType eType, char **papszParamList)
     889             : {
     890          57 :     if (nBandsIn != 1)
     891             :     {
     892          18 :         CPLError(CE_Failure, CPLE_FileIO,
     893             :                  "Only single band datasets are supported for writing");
     894          18 :         return nullptr;
     895             :     }
     896          39 :     if (eType != GDT_Float32)
     897             :     {
     898          25 :         CPLError(CE_Failure, CPLE_FileIO,
     899             :                  "Float32 is the only supported data type");
     900          25 :         return nullptr;
     901             :     }
     902          14 :     NWT_GRDDataset *poDS = new NWT_GRDDataset();
     903          14 :     poDS->eAccess = GA_Update;
     904          14 :     poDS->pGrd = static_cast<NWT_GRID *>(calloc(1, sizeof(NWT_GRID)));
     905             : 
     906             :     // We currently only support GRD grid types (could potentially support GRC
     907             :     // in the papszParamList). Also only support GDT_Float32 as the data type.
     908             :     // GRD format allows for data to be stretched to 32bit or 16bit integers on
     909             :     // disk, so it would be feasible to support other data types
     910          14 :     poDS->pGrd->cFormat = 0x00;
     911             : 
     912             :     // File version
     913          14 :     poDS->pGrd->fVersion = 2.0;
     914             : 
     915             :     // Dimensions
     916          14 :     poDS->pGrd->nXSide = nXSize;
     917          14 :     poDS->pGrd->nYSide = nYSize;
     918          14 :     poDS->nRasterXSize = nXSize;
     919          14 :     poDS->nRasterYSize = nYSize;
     920             : 
     921             :     // Some default values to get started with. These will
     922             :     // be altered when SetGeoTransform is called.
     923          14 :     poDS->pGrd->dfMinX = -2E+307;
     924          14 :     poDS->pGrd->dfMinY = -2E+307;
     925          14 :     poDS->pGrd->dfMaxX = 2E+307;
     926          14 :     poDS->pGrd->dfMaxY = 2E+307;
     927             : 
     928             :     float fZMin, fZMax;
     929             :     // See if the user passed the min/max values
     930          14 :     if (CSLFetchNameValue(papszParamList, "ZMIN") == nullptr)
     931             :     {
     932           1 :         fZMin = static_cast<float>(-2E+37);
     933             :     }
     934             :     else
     935             :     {
     936          13 :         fZMin = static_cast<float>(
     937          13 :             CPLAtof(CSLFetchNameValue(papszParamList, "ZMIN")));
     938             :     }
     939             : 
     940          14 :     if (CSLFetchNameValue(papszParamList, "ZMAX") == nullptr)
     941             :     {
     942           1 :         fZMax = static_cast<float>(2E+38);
     943             :     }
     944             :     else
     945             :     {
     946          13 :         fZMax = static_cast<float>(
     947          13 :             CPLAtof(CSLFetchNameValue(papszParamList, "ZMAX")));
     948             :     }
     949             : 
     950          14 :     poDS->pGrd->fZMin = fZMin;
     951          14 :     poDS->pGrd->fZMax = fZMax;
     952             :     // pGrd->dfStepSize = (pGrd->dfMaxX - pGrd->dfMinX) / (pGrd->nXSide - 1);
     953          14 :     poDS->pGrd->fZMinScale = fZMin;
     954          14 :     poDS->pGrd->fZMaxScale = fZMax;
     955             :     // poDS->pGrd->iZUnits
     956          14 :     memset(poDS->pGrd->cZUnits, 0, 32);
     957          14 :     memset(poDS->pGrd->cMICoordSys, 0, 256);
     958             : 
     959             :     // Some default colour inflections; Basic scale from blue to red
     960          14 :     poDS->pGrd->iNumColorInflections = 3;
     961             : 
     962             :     // Lowest inflection
     963          14 :     poDS->pGrd->stInflection[0].zVal = poDS->pGrd->fZMin;
     964          14 :     poDS->pGrd->stInflection[0].r = 0;
     965          14 :     poDS->pGrd->stInflection[0].g = 0;
     966          14 :     poDS->pGrd->stInflection[0].b = 255;
     967             : 
     968             :     // Mean inflection
     969          14 :     poDS->pGrd->stInflection[1].zVal =
     970          14 :         (poDS->pGrd->fZMax - poDS->pGrd->fZMin) / 2;
     971          14 :     poDS->pGrd->stInflection[1].r = 255;
     972          14 :     poDS->pGrd->stInflection[1].g = 255;
     973          14 :     poDS->pGrd->stInflection[1].b = 0;
     974             : 
     975             :     // Highest inflection
     976          14 :     poDS->pGrd->stInflection[2].zVal = poDS->pGrd->fZMax;
     977          14 :     poDS->pGrd->stInflection[2].r = 255;
     978          14 :     poDS->pGrd->stInflection[2].g = 0;
     979          14 :     poDS->pGrd->stInflection[2].b = 0;
     980             : 
     981          14 :     poDS->pGrd->bHillShadeExists = FALSE;
     982          14 :     poDS->pGrd->bShowGradient = FALSE;
     983          14 :     poDS->pGrd->bShowHillShade = FALSE;
     984          14 :     poDS->pGrd->cHillShadeBrightness = 0;
     985          14 :     poDS->pGrd->cHillShadeContrast = 0;
     986          14 :     poDS->pGrd->fHillShadeAzimuth = 0;
     987          14 :     poDS->pGrd->fHillShadeAngle = 0;
     988             : 
     989             :     // Set the raster style settings. These aren't used anywhere other than to
     990             :     // write the TAB file
     991          14 :     if (CSLFetchNameValue(papszParamList, "BRIGHTNESS") == nullptr)
     992             :     {
     993          14 :         poDS->pGrd->style.iBrightness = 50;
     994             :     }
     995             :     else
     996             :     {
     997           0 :         poDS->pGrd->style.iBrightness =
     998           0 :             atoi(CSLFetchNameValue(papszParamList, "BRIGHTNESS"));
     999             :     }
    1000             : 
    1001          14 :     if (CSLFetchNameValue(papszParamList, "CONTRAST") == nullptr)
    1002             :     {
    1003          14 :         poDS->pGrd->style.iContrast = 50;
    1004             :     }
    1005             :     else
    1006             :     {
    1007           0 :         poDS->pGrd->style.iContrast =
    1008           0 :             atoi(CSLFetchNameValue(papszParamList, "CONTRAST"));
    1009             :     }
    1010             : 
    1011          14 :     if (CSLFetchNameValue(papszParamList, "TRANSCOLOR") == nullptr)
    1012             :     {
    1013          14 :         poDS->pGrd->style.iTransColour = 0;
    1014             :     }
    1015             :     else
    1016             :     {
    1017           0 :         poDS->pGrd->style.iTransColour =
    1018           0 :             atoi(CSLFetchNameValue(papszParamList, "TRANSCOLOR"));
    1019             :     }
    1020             : 
    1021          14 :     if (CSLFetchNameValue(papszParamList, "TRANSLUCENCY") == nullptr)
    1022             :     {
    1023          14 :         poDS->pGrd->style.iTranslucency = 0;
    1024             :     }
    1025             :     else
    1026             :     {
    1027           0 :         poDS->pGrd->style.iTranslucency =
    1028           0 :             atoi(CSLFetchNameValue(papszParamList, "TRANSLUCENCY"));
    1029             :     }
    1030             : 
    1031          14 :     poDS->pGrd->style.bGreyscale = FALSE;
    1032          14 :     poDS->pGrd->style.bGrey = FALSE;
    1033          14 :     poDS->pGrd->style.bColour = FALSE;
    1034          14 :     poDS->pGrd->style.bTransparent = FALSE;
    1035             : 
    1036             :     // Open the grid file
    1037          14 :     poDS->fp = VSIFOpenL(pszFilename, "wb");
    1038          14 :     if (poDS->fp == nullptr)
    1039             :     {
    1040           1 :         CPLError(CE_Failure, CPLE_FileIO, "Failed to create GRD file");
    1041           1 :         delete poDS;
    1042           1 :         return nullptr;
    1043             :     }
    1044             : 
    1045          13 :     poDS->pGrd->fp = poDS->fp;
    1046          13 :     strncpy(poDS->pGrd->szFileName, pszFilename,
    1047             :             sizeof(poDS->pGrd->szFileName) - 1);
    1048          13 :     poDS->pGrd->szFileName[sizeof(poDS->pGrd->szFileName) - 1] = '\0';
    1049             : 
    1050             :     // Seek to the start of the file and enter the default header info
    1051          13 :     VSIFSeekL(poDS->fp, 0, SEEK_SET);
    1052          13 :     if (poDS->UpdateHeader() != 0)
    1053             :     {
    1054          10 :         CPLError(CE_Failure, CPLE_FileIO, "Failed to create GRD file");
    1055          10 :         delete poDS;
    1056          10 :         return nullptr;
    1057             :     }
    1058             : 
    1059             :     /* -------------------------------------------------------------------- */
    1060             :     /*      Create band information objects;                                */
    1061             :     /*      Only 1 band is allowed                                          */
    1062             :     /* -------------------------------------------------------------------- */
    1063           3 :     poDS->SetBand(1, new NWT_GRDRasterBand(poDS, 1, 1));  // z
    1064             : 
    1065           3 :     poDS->oOvManager.Initialize(poDS, pszFilename);
    1066           3 :     poDS->FlushCache(false);  // Write the header to disk.
    1067             : 
    1068           3 :     return poDS;
    1069             : }
    1070             : 
    1071             : /************************************************************************/
    1072             : /*                                CreateCopy()                          */
    1073             : /************************************************************************/
    1074          30 : GDALDataset *NWT_GRDDataset::CreateCopy(const char *pszFilename,
    1075             :                                         GDALDataset *poSrcDS, int bStrict,
    1076             :                                         char **papszOptions,
    1077             :                                         GDALProgressFunc pfnProgress,
    1078             :                                         void *pProgressData)
    1079             : {
    1080             : 
    1081          30 :     if (poSrcDS->GetRasterCount() != 1)
    1082             :     {
    1083           5 :         CPLError(CE_Failure, CPLE_FileIO,
    1084             :                  "Only single band datasets are supported for writing");
    1085           5 :         return nullptr;
    1086             :     }
    1087             : 
    1088          25 :     char **tmpOptions = CSLDuplicate(papszOptions);
    1089             : 
    1090             :     /*
    1091             :      * Compute the statistics if ZMAX and ZMIN are not provided
    1092             :      */
    1093          25 :     double dfMin = 0.0;
    1094          25 :     double dfMax = 0.0;
    1095          25 :     double dfMean = 0.0;
    1096          25 :     double dfStdDev = 0.0;
    1097          25 :     GDALRasterBand *pBand = poSrcDS->GetRasterBand(1);
    1098          25 :     char sMax[10] = {};
    1099          25 :     char sMin[10] = {};
    1100             : 
    1101          25 :     if ((CSLFetchNameValue(papszOptions, "ZMAX") == nullptr) ||
    1102           0 :         (CSLFetchNameValue(papszOptions, "ZMIN") == nullptr))
    1103             :     {
    1104          25 :         CPL_IGNORE_RET_VAL(pBand->GetStatistics(FALSE, TRUE, &dfMin, &dfMax,
    1105          25 :                                                 &dfMean, &dfStdDev));
    1106             :     }
    1107             : 
    1108          25 :     if (CSLFetchNameValue(papszOptions, "ZMAX") == nullptr)
    1109             :     {
    1110          25 :         CPLsnprintf(sMax, sizeof(sMax), "%f", dfMax);
    1111          25 :         tmpOptions = CSLSetNameValue(tmpOptions, "ZMAX", sMax);
    1112             :     }
    1113          25 :     if (CSLFetchNameValue(papszOptions, "ZMIN") == nullptr)
    1114             :     {
    1115          25 :         CPLsnprintf(sMin, sizeof(sMin), "%f", dfMin);
    1116          25 :         tmpOptions = CSLSetNameValue(tmpOptions, "ZMIN", sMin);
    1117             :     }
    1118             : 
    1119             :     GDALDriver *poDriver =
    1120          25 :         GDALDriver::FromHandle(GDALGetDriverByName("NWT_GRD"));
    1121          25 :     GDALDataset *poDstDS = poDriver->DefaultCreateCopy(
    1122             :         pszFilename, poSrcDS, bStrict, tmpOptions, pfnProgress, pProgressData);
    1123             : 
    1124          25 :     CSLDestroy(tmpOptions);
    1125             : 
    1126          25 :     return poDstDS;
    1127             : }
    1128             : #endif  // NO_MITAB_SUPPORT
    1129             : 
    1130             : /************************************************************************/
    1131             : /*                          GDALRegister_GRD()                          */
    1132             : /************************************************************************/
    1133        1595 : void GDALRegister_NWT_GRD()
    1134             : {
    1135        1595 :     if (GDALGetDriverByName("NWT_GRD") != nullptr)
    1136         302 :         return;
    1137             : 
    1138        1293 :     GDALDriver *poDriver = new GDALDriver();
    1139             : 
    1140        1293 :     poDriver->SetDescription("NWT_GRD");
    1141        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
    1142        1293 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
    1143        1293 :                               "Northwood Numeric Grid Format .grd/.tab");
    1144        1293 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/nwtgrd.html");
    1145        1293 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "grd");
    1146        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
    1147             : 
    1148        1293 :     poDriver->SetMetadataItem(
    1149             :         GDAL_DMD_OPENOPTIONLIST,
    1150             :         "<OpenOptionList>"
    1151             :         "    <Option name='BAND_COUNT' type='int' description='1 (Z) or 4 "
    1152             :         "(RGBZ). Only used in read-only mode' default='4'/>"
    1153        1293 :         "</OpenOptionList>");
    1154             : 
    1155        1293 :     poDriver->pfnOpen = NWT_GRDDataset::Open;
    1156        1293 :     poDriver->pfnIdentify = NWT_GRDDataset::Identify;
    1157             : 
    1158             : #ifndef NO_MITAB_SUPPORT
    1159             : 
    1160        1293 :     poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Float32");
    1161        1293 :     poDriver->SetMetadataItem(
    1162             :         GDAL_DMD_CREATIONOPTIONLIST,
    1163             :         "<CreationOptionList>"
    1164             :         "    <Option name='ZMIN' type='float' description='Minimum cell value "
    1165             :         "of raster for defining RGB scaling' default='-2E+37'/>"
    1166             :         "    <Option name='ZMAX' type='float' description='Maximum cell value "
    1167             :         "of raster for defining RGB scaling' default='2E+38'/>"
    1168             :         "    <Option name='BRIGHTNESS' type='int' description='Brightness to "
    1169             :         "be recorded in TAB file. Only affects reading with MapInfo' "
    1170             :         "default='50'/>"
    1171             :         "    <Option name='CONTRAST' type='int' description='Contrast to be "
    1172             :         "recorded in TAB file. Only affects reading with MapInfo' "
    1173             :         "default='50'/>"
    1174             :         "    <Option name='TRANSCOLOR' type='int' description='Transparent "
    1175             :         "color to be recorded in TAB file. Only affects reading with MapInfo' "
    1176             :         "default='0'/>"
    1177             :         "    <Option name='TRANSLUCENCY' type='int' description='Level of "
    1178             :         "translucency to be recorded in TAB file. Only affects reading with "
    1179             :         "MapInfo' default='0'/>"
    1180        1293 :         "</CreationOptionList>");
    1181             : 
    1182        1293 :     poDriver->pfnCreate = NWT_GRDDataset::Create;
    1183        1293 :     poDriver->pfnCreateCopy = NWT_GRDDataset::CreateCopy;
    1184             : #endif
    1185             : 
    1186        1293 :     GetGDALDriverManager()->RegisterDriver(poDriver);
    1187             : }