Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Natural Resources Canada's Geoid BYN file format
       4             :  * Purpose:  Implementation of BYN format
       5             :  * Author:   Ivan Lucena, ivan.lucena@outlook.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2018, Ivan Lucena
       9             :  * Copyright (c) 2018, Even Rouault
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "cpl_port.h"
      15             : #include "byndataset.h"
      16             : #include "rawdataset.h"
      17             : 
      18             : #include "cpl_string.h"
      19             : #include "gdal_frmts.h"
      20             : #include "ogr_spatialref.h"
      21             : #include "ogr_srs_api.h"
      22             : 
      23             : #include <algorithm>
      24             : #include <cstdlib>
      25             : #include <limits>
      26             : 
      27             : // Specification at
      28             : // https://www.nrcan.gc.ca/sites/www.nrcan.gc.ca/files/earthsciences/pdf/gpshgrid_e.pdf
      29             : 
      30             : const static BYNEllipsoids EllipsoidTable[] = {
      31             :     {"GRS80", 6378137.0, 298.257222101},
      32             :     {"WGS84", 6378137.0, 298.257223564},
      33             :     {"ALT1", 6378136.3, 298.256415099},
      34             :     {"GRS67", 6378160.0, 298.247167427},
      35             :     {"ELLIP1", 6378136.46, 298.256415099},
      36             :     {"ALT2", 6378136.3, 298.257},
      37             :     {"ELLIP2", 6378136.0, 298.257},
      38             :     {"CLARKE 1866", 6378206.4, 294.9786982}};
      39             : 
      40             : /************************************************************************/
      41             : /*                            BYNRasterBand()                           */
      42             : /************************************************************************/
      43             : 
      44           8 : BYNRasterBand::BYNRasterBand(GDALDataset *poDSIn, int nBandIn,
      45             :                              VSILFILE *fpRawIn, vsi_l_offset nImgOffsetIn,
      46             :                              int nPixelOffsetIn, int nLineOffsetIn,
      47           8 :                              GDALDataType eDataTypeIn, int bNativeOrderIn)
      48             :     : RawRasterBand(poDSIn, nBandIn, fpRawIn, nImgOffsetIn, nPixelOffsetIn,
      49             :                     nLineOffsetIn, eDataTypeIn, bNativeOrderIn,
      50           8 :                     RawRasterBand::OwnFP::NO)
      51             : {
      52           8 : }
      53             : 
      54             : /************************************************************************/
      55             : /*                           ~BYNRasterBand()                           */
      56             : /************************************************************************/
      57             : 
      58          16 : BYNRasterBand::~BYNRasterBand()
      59             : {
      60          16 : }
      61             : 
      62             : /************************************************************************/
      63             : /*                           GetNoDataValue()                           */
      64             : /************************************************************************/
      65             : 
      66           6 : double BYNRasterBand::GetNoDataValue(int *pbSuccess)
      67             : {
      68           6 :     if (pbSuccess)
      69           5 :         *pbSuccess = TRUE;
      70           6 :     int bSuccess = FALSE;
      71           6 :     double dfNoData = GDALPamRasterBand::GetNoDataValue(&bSuccess);
      72           6 :     if (bSuccess)
      73             :     {
      74           2 :         return dfNoData;
      75             :     }
      76           4 :     const double dfFactor =
      77           4 :         reinterpret_cast<BYNDataset *>(poDS)->hHeader.dfFactor;
      78           4 :     return eDataType == GDT_Int16 ? 32767.0 : 9999.0 * dfFactor;
      79             : }
      80             : 
      81             : /************************************************************************/
      82             : /*                              GetScale()                              */
      83             : /************************************************************************/
      84             : 
      85           1 : double BYNRasterBand::GetScale(int *pbSuccess)
      86             : {
      87           1 :     if (pbSuccess != nullptr)
      88           1 :         *pbSuccess = TRUE;
      89           1 :     const double dfFactor =
      90           1 :         reinterpret_cast<BYNDataset *>(poDS)->hHeader.dfFactor;
      91           1 :     return (dfFactor != 0.0) ? 1.0 / dfFactor : 0.0;
      92             : }
      93             : 
      94             : /************************************************************************/
      95             : /*                              SetScale()                              */
      96             : /************************************************************************/
      97             : 
      98           1 : CPLErr BYNRasterBand::SetScale(double dfNewValue)
      99             : {
     100           1 :     BYNDataset *poIDS = reinterpret_cast<BYNDataset *>(poDS);
     101           1 :     poIDS->hHeader.dfFactor = 1.0 / dfNewValue;
     102           1 :     return CE_None;
     103             : }
     104             : 
     105             : /************************************************************************/
     106             : /* ==================================================================== */
     107             : /*                              BYNDataset                              */
     108             : /* ==================================================================== */
     109             : /************************************************************************/
     110             : 
     111           8 : BYNDataset::BYNDataset()
     112             :     : fpImage(nullptr), hHeader{0, 0, 0, 0, 0, 0,   0,   0, 0.0, 0,   0, 0,
     113           8 :                                 0, 0, 0, 0, 0, 0.0, 0.0, 0, 0,   0.0, 0}
     114             : {
     115           8 :     m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     116           8 :     adfGeoTransform[0] = 0.0;
     117           8 :     adfGeoTransform[1] = 1.0;
     118           8 :     adfGeoTransform[2] = 0.0;
     119           8 :     adfGeoTransform[3] = 0.0;
     120           8 :     adfGeoTransform[4] = 0.0;
     121           8 :     adfGeoTransform[5] = 1.0;
     122           8 : }
     123             : 
     124             : /************************************************************************/
     125             : /*                            ~BYNDataset()                             */
     126             : /************************************************************************/
     127             : 
     128          16 : BYNDataset::~BYNDataset()
     129             : 
     130             : {
     131           8 :     BYNDataset::Close();
     132          16 : }
     133             : 
     134             : /************************************************************************/
     135             : /*                              Close()                                 */
     136             : /************************************************************************/
     137             : 
     138          16 : CPLErr BYNDataset::Close()
     139             : {
     140          16 :     CPLErr eErr = CE_None;
     141          16 :     if (nOpenFlags != OPEN_FLAGS_CLOSED)
     142             :     {
     143           8 :         if (BYNDataset::FlushCache(true) != CE_None)
     144           0 :             eErr = CE_Failure;
     145             : 
     146           8 :         if (GetAccess() == GA_Update)
     147           1 :             UpdateHeader();
     148             : 
     149           8 :         if (fpImage != nullptr)
     150             :         {
     151           8 :             if (VSIFCloseL(fpImage) != 0)
     152             :             {
     153           0 :                 eErr = CE_Failure;
     154           0 :                 CPLError(CE_Failure, CPLE_FileIO, "I/O error");
     155             :             }
     156             :         }
     157             : 
     158           8 :         if (GDALPamDataset::Close() != CE_None)
     159           0 :             eErr = CE_Failure;
     160             :     }
     161          16 :     return eErr;
     162             : }
     163             : 
     164             : /************************************************************************/
     165             : /*                              Identify()                              */
     166             : /************************************************************************/
     167             : 
     168       51299 : int BYNDataset::Identify(GDALOpenInfo *poOpenInfo)
     169             : 
     170             : {
     171       51299 :     if (poOpenInfo->nHeaderBytes < BYN_HDR_SZ)
     172       48089 :         return FALSE;
     173             : 
     174             : /* -------------------------------------------------------------------- */
     175             : /*      Check file extension (.byn/.err)                                */
     176             : /* -------------------------------------------------------------------- */
     177             : #ifndef FUZZING_BUILD_MODE_UNSAFE_FOR_PRODUCTION
     178        3210 :     const char *pszFileExtension = CPLGetExtension(poOpenInfo->pszFilename);
     179             : 
     180        3210 :     if (!EQUAL(pszFileExtension, "byn") && !EQUAL(pszFileExtension, "err"))
     181             :     {
     182        3194 :         return FALSE;
     183             :     }
     184             : #endif
     185             : 
     186             :     /* -------------------------------------------------------------------- */
     187             :     /*      Check some value's ranges on header                             */
     188             :     /* -------------------------------------------------------------------- */
     189             : 
     190          16 :     BYNHeader hHeader = {0, 0, 0, 0, 0, 0,   0,   0, 0.0, 0,   0, 0,
     191             :                          0, 0, 0, 0, 0, 0.0, 0.0, 0, 0,   0.0, 0};
     192             : 
     193          16 :     buffer2header(poOpenInfo->pabyHeader, &hHeader);
     194             : 
     195          16 :     if (hHeader.nGlobal < 0 || hHeader.nGlobal > 1 || hHeader.nType < 0 ||
     196          16 :         hHeader.nType > 9 || (hHeader.nSizeOf != 2 && hHeader.nSizeOf != 4) ||
     197          16 :         hHeader.nVDatum < 0 || hHeader.nVDatum > 3 || hHeader.nDescrip < 0 ||
     198          16 :         hHeader.nDescrip > 3 || hHeader.nSubType < 0 || hHeader.nSubType > 9 ||
     199          16 :         hHeader.nDatum < 0 || hHeader.nDatum > 1 || hHeader.nEllipsoid < 0 ||
     200          16 :         hHeader.nEllipsoid > 7 || hHeader.nByteOrder < 0 ||
     201          16 :         hHeader.nByteOrder > 1 || hHeader.nScale < 0 || hHeader.nScale > 1)
     202           0 :         return FALSE;
     203             : 
     204             : #if 0
     205             :     // We have disabled those checks as invalid values are often found in some
     206             :     // datasets, such as http://s3.microsurvey.com/os/fieldgenius/geoids/Lithuania.zip
     207             :     // We don't use those fields, so we may just ignore them.
     208             :     if((hHeader.nTideSys   < 0 || hHeader.nTideSys   > 2 ||
     209             :         hHeader.nPtType    < 0 || hHeader.nPtType    > 1 ))
     210             :     {
     211             :         // Some datasets use 0xCC as a marker for invalidity for
     212             :         // records starting from Geopotential Wo
     213             :         for( int i = 52; i < 78; i++ )
     214             :         {
     215             :             if( poOpenInfo->pabyHeader[i] != 0xCC )
     216             :                 return FALSE;
     217             :         }
     218             :     }
     219             : #endif
     220             : 
     221          16 :     if (hHeader.nScale == 0)
     222             :     {
     223          32 :         if ((std::abs(static_cast<GIntBig>(hHeader.nSouth) -
     224          32 :                       (hHeader.nDLat / 2)) > BYN_MAX_LAT) ||
     225          16 :             (std::abs(static_cast<GIntBig>(hHeader.nNorth) +
     226          32 :                       (hHeader.nDLat / 2)) > BYN_MAX_LAT) ||
     227          16 :             (std::abs(static_cast<GIntBig>(hHeader.nWest) -
     228          48 :                       (hHeader.nDLon / 2)) > BYN_MAX_LON) ||
     229          16 :             (std::abs(static_cast<GIntBig>(hHeader.nEast) +
     230          16 :                       (hHeader.nDLon / 2)) > BYN_MAX_LON))
     231           0 :             return FALSE;
     232             :     }
     233             :     else
     234             :     {
     235           0 :         if ((std::abs(static_cast<GIntBig>(hHeader.nSouth) -
     236           0 :                       (hHeader.nDLat / 2)) > BYN_MAX_LAT_SCL) ||
     237           0 :             (std::abs(static_cast<GIntBig>(hHeader.nNorth) +
     238           0 :                       (hHeader.nDLat / 2)) > BYN_MAX_LAT_SCL) ||
     239           0 :             (std::abs(static_cast<GIntBig>(hHeader.nWest) -
     240           0 :                       (hHeader.nDLon / 2)) > BYN_MAX_LON_SCL) ||
     241           0 :             (std::abs(static_cast<GIntBig>(hHeader.nEast) +
     242           0 :                       (hHeader.nDLon / 2)) > BYN_MAX_LON_SCL))
     243           0 :             return FALSE;
     244             :     }
     245             : 
     246          16 :     return TRUE;
     247             : }
     248             : 
     249             : /************************************************************************/
     250             : /*                                Open()                                */
     251             : /************************************************************************/
     252             : 
     253           8 : GDALDataset *BYNDataset::Open(GDALOpenInfo *poOpenInfo)
     254             : 
     255             : {
     256           8 :     if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
     257           0 :         return nullptr;
     258             : 
     259             :     /* -------------------------------------------------------------------- */
     260             :     /*      Create a corresponding GDALDataset.                             */
     261             :     /* -------------------------------------------------------------------- */
     262             : 
     263          16 :     auto poDS = std::make_unique<BYNDataset>();
     264             : 
     265           8 :     poDS->eAccess = poOpenInfo->eAccess;
     266           8 :     std::swap(poDS->fpImage, poOpenInfo->fpL);
     267             : 
     268             :     /* -------------------------------------------------------------------- */
     269             :     /*      Read the header.                                                */
     270             :     /* -------------------------------------------------------------------- */
     271             : 
     272           8 :     buffer2header(poOpenInfo->pabyHeader, &poDS->hHeader);
     273             : 
     274             :     /********************************/
     275             :     /* Scale boundaries and spacing */
     276             :     /********************************/
     277             : 
     278           8 :     double dfSouth = poDS->hHeader.nSouth;
     279           8 :     double dfNorth = poDS->hHeader.nNorth;
     280           8 :     double dfWest = poDS->hHeader.nWest;
     281           8 :     double dfEast = poDS->hHeader.nEast;
     282           8 :     double dfDLat = poDS->hHeader.nDLat;
     283           8 :     double dfDLon = poDS->hHeader.nDLon;
     284             : 
     285           8 :     if (poDS->hHeader.nScale == 1)
     286             :     {
     287           0 :         dfSouth *= BYN_SCALE;
     288           0 :         dfNorth *= BYN_SCALE;
     289           0 :         dfWest *= BYN_SCALE;
     290           0 :         dfEast *= BYN_SCALE;
     291           0 :         dfDLat *= BYN_SCALE;
     292           0 :         dfDLon *= BYN_SCALE;
     293             :     }
     294             : 
     295             :     /******************************/
     296             :     /* Calculate rows and columns */
     297             :     /******************************/
     298             : 
     299           8 :     double dfXSize = -1;
     300           8 :     double dfYSize = -1;
     301             : 
     302           8 :     poDS->nRasterXSize = -1;
     303           8 :     poDS->nRasterYSize = -1;
     304             : 
     305           8 :     if (dfDLat != 0.0 && dfDLon != 0.0)
     306             :     {
     307           8 :         dfXSize = ((dfEast - dfWest + 1.0) / dfDLon) + 1.0;
     308           8 :         dfYSize = ((dfNorth - dfSouth + 1.0) / dfDLat) + 1.0;
     309             :     }
     310             : 
     311           8 :     if (dfXSize > 0.0 && dfXSize < std::numeric_limits<double>::max() &&
     312          16 :         dfYSize > 0.0 && dfYSize < std::numeric_limits<double>::max())
     313             :     {
     314           8 :         poDS->nRasterXSize = static_cast<GInt32>(dfXSize);
     315           8 :         poDS->nRasterYSize = static_cast<GInt32>(dfYSize);
     316             :     }
     317             : 
     318           8 :     if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize))
     319             :     {
     320           0 :         return nullptr;
     321             :     }
     322             : 
     323             :     /*****************************/
     324             :     /* Build GeoTransform matrix */
     325             :     /*****************************/
     326             : 
     327           8 :     poDS->adfGeoTransform[0] = (dfWest - (dfDLon / 2.0)) / 3600.0;
     328           8 :     poDS->adfGeoTransform[1] = dfDLon / 3600.0;
     329           8 :     poDS->adfGeoTransform[2] = 0.0;
     330           8 :     poDS->adfGeoTransform[3] = (dfNorth + (dfDLat / 2.0)) / 3600.0;
     331           8 :     poDS->adfGeoTransform[4] = 0.0;
     332           8 :     poDS->adfGeoTransform[5] = -1 * dfDLat / 3600.0;
     333             : 
     334             :     /*********************/
     335             :     /* Set data type     */
     336             :     /*********************/
     337             : 
     338           8 :     GDALDataType eDT = GDT_Unknown;
     339             : 
     340           8 :     if (poDS->hHeader.nSizeOf == 2)
     341           0 :         eDT = GDT_Int16;
     342           8 :     else if (poDS->hHeader.nSizeOf == 4)
     343           8 :         eDT = GDT_Int32;
     344             :     else
     345             :     {
     346           0 :         return nullptr;
     347             :     }
     348             : 
     349             :     /* -------------------------------------------------------------------- */
     350             :     /*      Create band information object.                                 */
     351             :     /* -------------------------------------------------------------------- */
     352             : 
     353           8 :     const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
     354             : 
     355           8 :     int bIsLSB = poDS->hHeader.nByteOrder == 1 ? 1 : 0;
     356             : 
     357             :     auto poBand = std::make_unique<BYNRasterBand>(
     358           8 :         poDS.get(), 1, poDS->fpImage, BYN_HDR_SZ, nDTSize,
     359          24 :         poDS->nRasterXSize * nDTSize, eDT, CPL_IS_LSB == bIsLSB);
     360           8 :     if (!poBand->IsValid())
     361           0 :         return nullptr;
     362           8 :     poDS->SetBand(1, std::move(poBand));
     363             : 
     364             :     /* -------------------------------------------------------------------- */
     365             :     /*      Initialize any PAM information.                                 */
     366             :     /* -------------------------------------------------------------------- */
     367             : 
     368           8 :     poDS->SetDescription(poOpenInfo->pszFilename);
     369           8 :     poDS->TryLoadXML();
     370             : 
     371             :     /* -------------------------------------------------------------------- */
     372             :     /*      Check for overviews.                                            */
     373             :     /* -------------------------------------------------------------------- */
     374             : 
     375           8 :     poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
     376             : 
     377           8 :     return poDS.release();
     378             : }
     379             : 
     380             : /************************************************************************/
     381             : /*                          GetGeoTransform()                           */
     382             : /************************************************************************/
     383             : 
     384           2 : CPLErr BYNDataset::GetGeoTransform(double *padfTransform)
     385             : 
     386             : {
     387           2 :     memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
     388           2 :     return CE_None;
     389             : }
     390             : 
     391             : /************************************************************************/
     392             : /*                          SetGeoTransform()                           */
     393             : /************************************************************************/
     394             : 
     395           1 : CPLErr BYNDataset::SetGeoTransform(double *padfTransform)
     396             : 
     397             : {
     398           1 :     if (padfTransform[2] != 0.0 || padfTransform[4] != 0.0)
     399             :     {
     400           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     401             :                  "Attempt to write skewed or rotated geotransform to byn.");
     402           0 :         return CE_Failure;
     403             :     }
     404           1 :     memcpy(adfGeoTransform, padfTransform, sizeof(double) * 6);
     405           1 :     return CE_None;
     406             : }
     407             : 
     408             : /************************************************************************/
     409             : /*                          GetSpatialRef()                             */
     410             : /************************************************************************/
     411             : 
     412           1 : const OGRSpatialReference *BYNDataset::GetSpatialRef() const
     413             : 
     414             : {
     415           1 :     if (!m_oSRS.IsEmpty())
     416           0 :         return &m_oSRS;
     417             : 
     418             :     /* Try to use a prefefined EPSG compound CS */
     419             : 
     420           1 :     if (hHeader.nDatum == 1 && hHeader.nVDatum == 2)
     421             :     {
     422           0 :         m_oSRS.importFromEPSG(BYN_DATUM_1_VDATUM_2);
     423           0 :         return &m_oSRS;
     424             :     }
     425             : 
     426             :     /* Build the GEOGCS based on Datum ( or Ellipsoid )*/
     427             : 
     428           1 :     bool bNoGeogCS = false;
     429             : 
     430           1 :     if (hHeader.nDatum == 0)
     431           1 :         m_oSRS.importFromEPSG(BYN_DATUM_0);
     432           0 :     else if (hHeader.nDatum == 1)
     433           0 :         m_oSRS.importFromEPSG(BYN_DATUM_1);
     434             :     else
     435             :     {
     436             :         /* Build GEOGCS based on Ellipsoid (Table 3) */
     437             : 
     438           0 :         if (hHeader.nEllipsoid > -1 &&
     439           0 :             hHeader.nEllipsoid <
     440             :                 static_cast<GInt16>(CPL_ARRAYSIZE(EllipsoidTable)))
     441           0 :             m_oSRS.SetGeogCS(
     442           0 :                 CPLSPrintf("BYN Ellipsoid(%d)", hHeader.nEllipsoid),
     443           0 :                 "Unspecified", EllipsoidTable[hHeader.nEllipsoid].pszName,
     444           0 :                 EllipsoidTable[hHeader.nEllipsoid].dfSemiMajor,
     445           0 :                 EllipsoidTable[hHeader.nEllipsoid].dfInvFlattening);
     446             :         else
     447           0 :             bNoGeogCS = true;
     448             :     }
     449             : 
     450             :     /* Build the VERT_CS based on VDatum */
     451             : 
     452           2 :     OGRSpatialReference oSRSComp;
     453           2 :     OGRSpatialReference oSRSVert;
     454             : 
     455           1 :     int nVertCS = 0;
     456             : 
     457           1 :     if (hHeader.nVDatum == 1)
     458           0 :         nVertCS = BYN_VDATUM_1;
     459           1 :     else if (hHeader.nVDatum == 2)
     460           1 :         nVertCS = BYN_VDATUM_2;
     461           0 :     else if (hHeader.nVDatum == 3)
     462           0 :         nVertCS = BYN_VDATUM_3;
     463             :     else
     464             :     {
     465             :         /* Return GEOGCS ( .err files ) */
     466             : 
     467           0 :         if (bNoGeogCS)
     468           0 :             return nullptr;
     469             : 
     470           0 :         return &m_oSRS;
     471             :     }
     472             : 
     473           1 :     oSRSVert.importFromEPSG(nVertCS);
     474             : 
     475             :     /* Create CPMPD_CS with GEOGCS and VERT_CS */
     476             : 
     477           1 :     if (oSRSComp.SetCompoundCS(CPLSPrintf("BYN Datum(%d) & VDatum(%d)",
     478           1 :                                           hHeader.nDatum, hHeader.nDatum),
     479           2 :                                &m_oSRS, &oSRSVert) == CE_None)
     480             :     {
     481             :         /* Return COMPD_CS with GEOGCS and VERT_CS */
     482             : 
     483           1 :         m_oSRS = std::move(oSRSComp);
     484           1 :         m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     485           1 :         return &m_oSRS;
     486             :     }
     487             : 
     488           0 :     return nullptr;
     489             : }
     490             : 
     491             : /************************************************************************/
     492             : /*                          SetSpatialRef()                             */
     493             : /************************************************************************/
     494             : 
     495           1 : CPLErr BYNDataset::SetSpatialRef(const OGRSpatialReference *poSRS)
     496             : 
     497             : {
     498           1 :     if (poSRS == nullptr)
     499           0 :         return CE_Failure;
     500             : 
     501           1 :     m_oSRS = *poSRS;
     502             : 
     503             :     /* Try to recognize prefefined EPSG compound CS */
     504             : 
     505           1 :     if (poSRS->IsCompound())
     506             :     {
     507           1 :         const char *pszAuthName = poSRS->GetAuthorityName("COMPD_CS");
     508           1 :         const char *pszAuthCode = poSRS->GetAuthorityCode("COMPD_CS");
     509             : 
     510           1 :         if (pszAuthName != nullptr && pszAuthCode != nullptr &&
     511           0 :             EQUAL(pszAuthName, "EPSG") &&
     512           0 :             atoi(pszAuthCode) == BYN_DATUM_1_VDATUM_2)
     513             :         {
     514           0 :             hHeader.nVDatum = 2;
     515           0 :             hHeader.nDatum = 1;
     516           0 :             return CE_None;
     517             :         }
     518             :     }
     519             : 
     520           1 :     OGRSpatialReference oSRSTemp;
     521             : 
     522             :     /* Try to match GEOGCS */
     523             : 
     524           1 :     if (poSRS->IsGeographic())
     525             :     {
     526           1 :         oSRSTemp.importFromEPSG(BYN_DATUM_0);
     527           1 :         if (poSRS->IsSameGeogCS(&oSRSTemp))
     528           1 :             hHeader.nDatum = 0;
     529             :         else
     530             :         {
     531           0 :             oSRSTemp.importFromEPSG(BYN_DATUM_1);
     532           0 :             if (poSRS->IsSameGeogCS(&oSRSTemp))
     533           0 :                 hHeader.nDatum = 1;
     534             :         }
     535             :     }
     536             : 
     537             :     /* Try to match VERT_CS */
     538             : 
     539           1 :     if (poSRS->IsVertical())
     540             :     {
     541           1 :         oSRSTemp.importFromEPSG(BYN_VDATUM_1);
     542           1 :         if (poSRS->IsSameVertCS(&oSRSTemp))
     543           0 :             hHeader.nVDatum = 1;
     544             :         else
     545             :         {
     546           1 :             oSRSTemp.importFromEPSG(BYN_VDATUM_2);
     547           1 :             if (poSRS->IsSameVertCS(&oSRSTemp))
     548           1 :                 hHeader.nVDatum = 2;
     549             :             else
     550             :             {
     551           0 :                 oSRSTemp.importFromEPSG(BYN_VDATUM_3);
     552           0 :                 if (poSRS->IsSameVertCS(&oSRSTemp))
     553           0 :                     hHeader.nVDatum = 3;
     554             :             }
     555             :         }
     556             :     }
     557             : 
     558           1 :     return CE_None;
     559             : }
     560             : 
     561             : /*----------------------------------------------------------------------*/
     562             : /*                           buffer2header()                            */
     563             : /*----------------------------------------------------------------------*/
     564             : 
     565          24 : void BYNDataset::buffer2header(const GByte *pabyBuf, BYNHeader *pohHeader)
     566             : 
     567             : {
     568          24 :     memcpy(&pohHeader->nSouth, pabyBuf, 4);
     569          24 :     memcpy(&pohHeader->nNorth, pabyBuf + 4, 4);
     570          24 :     memcpy(&pohHeader->nWest, pabyBuf + 8, 4);
     571          24 :     memcpy(&pohHeader->nEast, pabyBuf + 12, 4);
     572          24 :     memcpy(&pohHeader->nDLat, pabyBuf + 16, 2);
     573          24 :     memcpy(&pohHeader->nDLon, pabyBuf + 18, 2);
     574          24 :     memcpy(&pohHeader->nGlobal, pabyBuf + 20, 2);
     575          24 :     memcpy(&pohHeader->nType, pabyBuf + 22, 2);
     576          24 :     memcpy(&pohHeader->dfFactor, pabyBuf + 24, 8);
     577          24 :     memcpy(&pohHeader->nSizeOf, pabyBuf + 32, 2);
     578          24 :     memcpy(&pohHeader->nVDatum, pabyBuf + 34, 2);
     579          24 :     memcpy(&pohHeader->nDescrip, pabyBuf + 40, 2);
     580          24 :     memcpy(&pohHeader->nSubType, pabyBuf + 42, 2);
     581          24 :     memcpy(&pohHeader->nDatum, pabyBuf + 44, 2);
     582          24 :     memcpy(&pohHeader->nEllipsoid, pabyBuf + 46, 2);
     583          24 :     memcpy(&pohHeader->nByteOrder, pabyBuf + 48, 2);
     584          24 :     memcpy(&pohHeader->nScale, pabyBuf + 50, 2);
     585          24 :     memcpy(&pohHeader->dfWo, pabyBuf + 52, 8);
     586          24 :     memcpy(&pohHeader->dfGM, pabyBuf + 60, 8);
     587          24 :     memcpy(&pohHeader->nTideSys, pabyBuf + 68, 2);
     588          24 :     memcpy(&pohHeader->nRealiz, pabyBuf + 70, 2);
     589          24 :     memcpy(&pohHeader->dEpoch, pabyBuf + 72, 4);
     590          24 :     memcpy(&pohHeader->nPtType, pabyBuf + 76, 2);
     591             : 
     592             : #if defined(CPL_MSB)
     593             :     CPL_LSBPTR32(&pohHeader->nSouth);
     594             :     CPL_LSBPTR32(&pohHeader->nNorth);
     595             :     CPL_LSBPTR32(&pohHeader->nWest);
     596             :     CPL_LSBPTR32(&pohHeader->nEast);
     597             :     CPL_LSBPTR16(&pohHeader->nDLat);
     598             :     CPL_LSBPTR16(&pohHeader->nDLon);
     599             :     CPL_LSBPTR16(&pohHeader->nGlobal);
     600             :     CPL_LSBPTR16(&pohHeader->nType);
     601             :     CPL_LSBPTR64(&pohHeader->dfFactor);
     602             :     CPL_LSBPTR16(&pohHeader->nSizeOf);
     603             :     CPL_LSBPTR16(&pohHeader->nVDatum);
     604             :     CPL_LSBPTR16(&pohHeader->nDescrip);
     605             :     CPL_LSBPTR16(&pohHeader->nSubType);
     606             :     CPL_LSBPTR16(&pohHeader->nDatum);
     607             :     CPL_LSBPTR16(&pohHeader->nEllipsoid);
     608             :     CPL_LSBPTR16(&pohHeader->nByteOrder);
     609             :     CPL_LSBPTR16(&pohHeader->nScale);
     610             :     CPL_LSBPTR64(&pohHeader->dfWo);
     611             :     CPL_LSBPTR64(&pohHeader->dfGM);
     612             :     CPL_LSBPTR16(&pohHeader->nTideSys);
     613             :     CPL_LSBPTR16(&pohHeader->nRealiz);
     614             :     CPL_LSBPTR32(&pohHeader->dEpoch);
     615             :     CPL_LSBPTR16(&pohHeader->nPtType);
     616             : #endif
     617             : 
     618             : #if DEBUG
     619          24 :     CPLDebug("BYN", "South         = %d", pohHeader->nSouth);
     620          24 :     CPLDebug("BYN", "North         = %d", pohHeader->nNorth);
     621          24 :     CPLDebug("BYN", "West          = %d", pohHeader->nWest);
     622          24 :     CPLDebug("BYN", "East          = %d", pohHeader->nEast);
     623          24 :     CPLDebug("BYN", "DLat          = %d", pohHeader->nDLat);
     624          24 :     CPLDebug("BYN", "DLon          = %d", pohHeader->nDLon);
     625          24 :     CPLDebug("BYN", "DGlobal       = %d", pohHeader->nGlobal);
     626          24 :     CPLDebug("BYN", "DType         = %d", pohHeader->nType);
     627          24 :     CPLDebug("BYN", "Factor        = %f", pohHeader->dfFactor);
     628          24 :     CPLDebug("BYN", "SizeOf        = %d", pohHeader->nSizeOf);
     629          24 :     CPLDebug("BYN", "VDatum        = %d", pohHeader->nVDatum);
     630          24 :     CPLDebug("BYN", "Data          = %d", pohHeader->nDescrip);
     631          24 :     CPLDebug("BYN", "SubType       = %d", pohHeader->nSubType);
     632          24 :     CPLDebug("BYN", "Datum         = %d", pohHeader->nDatum);
     633          24 :     CPLDebug("BYN", "Ellipsoid     = %d", pohHeader->nEllipsoid);
     634          24 :     CPLDebug("BYN", "ByteOrder     = %d", pohHeader->nByteOrder);
     635          24 :     CPLDebug("BYN", "Scale         = %d", pohHeader->nScale);
     636          24 :     CPLDebug("BYN", "Wo            = %f", pohHeader->dfWo);
     637          24 :     CPLDebug("BYN", "GM            = %f", pohHeader->dfGM);
     638          24 :     CPLDebug("BYN", "TideSystem    = %d", pohHeader->nTideSys);
     639          24 :     CPLDebug("BYN", "RefRealzation = %d", pohHeader->nRealiz);
     640          24 :     CPLDebug("BYN", "Epoch         = %f", pohHeader->dEpoch);
     641          24 :     CPLDebug("BYN", "PtType        = %d", pohHeader->nPtType);
     642             : #endif
     643          24 : }
     644             : 
     645             : /*----------------------------------------------------------------------*/
     646             : /*                           header2buffer()                            */
     647             : /*----------------------------------------------------------------------*/
     648             : 
     649           2 : void BYNDataset::header2buffer(const BYNHeader *pohHeader, GByte *pabyBuf)
     650             : 
     651             : {
     652           2 :     memcpy(pabyBuf, &pohHeader->nSouth, 4);
     653           2 :     memcpy(pabyBuf + 4, &pohHeader->nNorth, 4);
     654           2 :     memcpy(pabyBuf + 8, &pohHeader->nWest, 4);
     655           2 :     memcpy(pabyBuf + 12, &pohHeader->nEast, 4);
     656           2 :     memcpy(pabyBuf + 16, &pohHeader->nDLat, 2);
     657           2 :     memcpy(pabyBuf + 18, &pohHeader->nDLon, 2);
     658           2 :     memcpy(pabyBuf + 20, &pohHeader->nGlobal, 2);
     659           2 :     memcpy(pabyBuf + 22, &pohHeader->nType, 2);
     660           2 :     memcpy(pabyBuf + 24, &pohHeader->dfFactor, 8);
     661           2 :     memcpy(pabyBuf + 32, &pohHeader->nSizeOf, 2);
     662           2 :     memcpy(pabyBuf + 34, &pohHeader->nVDatum, 2);
     663           2 :     memcpy(pabyBuf + 40, &pohHeader->nDescrip, 2);
     664           2 :     memcpy(pabyBuf + 42, &pohHeader->nSubType, 2);
     665           2 :     memcpy(pabyBuf + 44, &pohHeader->nDatum, 2);
     666           2 :     memcpy(pabyBuf + 46, &pohHeader->nEllipsoid, 2);
     667           2 :     memcpy(pabyBuf + 48, &pohHeader->nByteOrder, 2);
     668           2 :     memcpy(pabyBuf + 50, &pohHeader->nScale, 2);
     669           2 :     memcpy(pabyBuf + 52, &pohHeader->dfWo, 8);
     670           2 :     memcpy(pabyBuf + 60, &pohHeader->dfGM, 8);
     671           2 :     memcpy(pabyBuf + 68, &pohHeader->nTideSys, 2);
     672           2 :     memcpy(pabyBuf + 70, &pohHeader->nRealiz, 2);
     673           2 :     memcpy(pabyBuf + 72, &pohHeader->dEpoch, 4);
     674           2 :     memcpy(pabyBuf + 76, &pohHeader->nPtType, 2);
     675             : 
     676             : #if defined(CPL_MSB)
     677             :     CPL_LSBPTR32(pabyBuf);
     678             :     CPL_LSBPTR32(pabyBuf + 4);
     679             :     CPL_LSBPTR32(pabyBuf + 8);
     680             :     CPL_LSBPTR32(pabyBuf + 12);
     681             :     CPL_LSBPTR16(pabyBuf + 16);
     682             :     CPL_LSBPTR16(pabyBuf + 18);
     683             :     CPL_LSBPTR16(pabyBuf + 20);
     684             :     CPL_LSBPTR16(pabyBuf + 22);
     685             :     CPL_LSBPTR64(pabyBuf + 24);
     686             :     CPL_LSBPTR16(pabyBuf + 32);
     687             :     CPL_LSBPTR16(pabyBuf + 34);
     688             :     CPL_LSBPTR16(pabyBuf + 40);
     689             :     CPL_LSBPTR16(pabyBuf + 42);
     690             :     CPL_LSBPTR16(pabyBuf + 44);
     691             :     CPL_LSBPTR16(pabyBuf + 46);
     692             :     CPL_LSBPTR16(pabyBuf + 48);
     693             :     CPL_LSBPTR16(pabyBuf + 50);
     694             :     CPL_LSBPTR64(pabyBuf + 52);
     695             :     CPL_LSBPTR64(pabyBuf + 60);
     696             :     CPL_LSBPTR16(pabyBuf + 68);
     697             :     CPL_LSBPTR16(pabyBuf + 70);
     698             :     CPL_LSBPTR32(pabyBuf + 72);
     699             :     CPL_LSBPTR16(pabyBuf + 76);
     700             : #endif
     701           2 : }
     702             : 
     703             : /************************************************************************/
     704             : /*                               Create()                               */
     705             : /************************************************************************/
     706             : 
     707          50 : GDALDataset *BYNDataset::Create(const char *pszFilename, int nXSize, int nYSize,
     708             :                                 int /* nBands */, GDALDataType eType,
     709             :                                 char ** /* papszOptions */)
     710             : {
     711          50 :     if (eType != GDT_Int16 && eType != GDT_Int32)
     712             :     {
     713          43 :         CPLError(CE_Failure, CPLE_AppDefined,
     714             :                  "Attempt to create byn file with unsupported data type '%s'.",
     715             :                  GDALGetDataTypeName(eType));
     716          43 :         return nullptr;
     717             :     }
     718             : 
     719             :     /* -------------------------------------------------------------------- */
     720             :     /*      Check file extension (.byn/.err)                                */
     721             :     /* -------------------------------------------------------------------- */
     722             : 
     723           7 :     char *pszFileExtension = CPLStrdup(CPLGetExtension(pszFilename));
     724             : 
     725           7 :     if (!EQUAL(pszFileExtension, "byn") && !EQUAL(pszFileExtension, "err"))
     726             :     {
     727           6 :         CPLError(
     728             :             CE_Failure, CPLE_AppDefined,
     729             :             "Attempt to create byn file with extension other than byn/err.");
     730           6 :         CPLFree(pszFileExtension);
     731           6 :         return nullptr;
     732             :     }
     733             : 
     734           1 :     CPLFree(pszFileExtension);
     735             : 
     736             :     /* -------------------------------------------------------------------- */
     737             :     /*      Try to create the file.                                         */
     738             :     /* -------------------------------------------------------------------- */
     739             : 
     740           1 :     VSILFILE *fp = VSIFOpenL(pszFilename, "wb+");
     741           1 :     if (fp == nullptr)
     742             :     {
     743           0 :         CPLError(CE_Failure, CPLE_OpenFailed,
     744             :                  "Attempt to create file `%s' failed.\n", pszFilename);
     745           0 :         return nullptr;
     746             :     }
     747             : 
     748             :     /* -------------------------------------------------------------------- */
     749             :     /*      Write an empty header.                                          */
     750             :     /* -------------------------------------------------------------------- */
     751             : 
     752           1 :     GByte abyBuf[BYN_HDR_SZ] = {'\0'};
     753             : 
     754             :     /* Load header with some commum values */
     755             : 
     756           1 :     BYNHeader hHeader = {0, 0, 0, 0, 0, 0,   0,   0, 0.0, 0,   0, 0,
     757             :                          0, 0, 0, 0, 0, 0.0, 0.0, 0, 0,   0.0, 0};
     758             : 
     759             :     /* Set temporary values on header */
     760             : 
     761           1 :     hHeader.nSouth = 0;
     762           1 :     hHeader.nNorth = nYSize - 2;
     763           1 :     hHeader.nWest = 0;
     764           1 :     hHeader.nEast = nXSize - 2;
     765           1 :     hHeader.nDLat = 1;
     766           1 :     hHeader.nDLon = 1;
     767           1 :     hHeader.nSizeOf = static_cast<GInt16>(GDALGetDataTypeSizeBytes(eType));
     768             : 
     769             :     /* Prepare buffer for writing */
     770             : 
     771           1 :     header2buffer(&hHeader, abyBuf);
     772             : 
     773             :     /* Write initial header */
     774             : 
     775           1 :     VSIFWriteL(abyBuf, BYN_HDR_SZ, 1, fp);
     776           1 :     VSIFCloseL(fp);
     777             : 
     778           1 :     return GDALDataset::FromHandle(GDALOpen(pszFilename, GA_Update));
     779             : }
     780             : 
     781             : /************************************************************************/
     782             : /*                          UpdateHeader()                              */
     783             : /************************************************************************/
     784             : 
     785           1 : void BYNDataset::UpdateHeader()
     786             : {
     787           1 :     double dfDLon = (adfGeoTransform[1] * 3600.0);
     788           1 :     double dfDLat = (adfGeoTransform[5] * 3600.0 * -1);
     789           1 :     double dfWest = (adfGeoTransform[0] * 3600.0) + (dfDLon / 2);
     790           1 :     double dfNorth = (adfGeoTransform[3] * 3600.0) - (dfDLat / 2);
     791           1 :     double dfSouth = dfNorth - ((nRasterYSize - 1) * dfDLat);
     792           1 :     double dfEast = dfWest + ((nRasterXSize - 1) * dfDLon);
     793             : 
     794           1 :     if (hHeader.nScale == 1)
     795             :     {
     796           0 :         dfSouth /= BYN_SCALE;
     797           0 :         dfNorth /= BYN_SCALE;
     798           0 :         dfWest /= BYN_SCALE;
     799           0 :         dfEast /= BYN_SCALE;
     800           0 :         dfDLat /= BYN_SCALE;
     801           0 :         dfDLon /= BYN_SCALE;
     802             :     }
     803             : 
     804           1 :     hHeader.nSouth = static_cast<GInt32>(dfSouth);
     805           1 :     hHeader.nNorth = static_cast<GInt32>(dfNorth);
     806           1 :     hHeader.nWest = static_cast<GInt32>(dfWest);
     807           1 :     hHeader.nEast = static_cast<GInt32>(dfEast);
     808           1 :     hHeader.nDLat = static_cast<GInt16>(dfDLat);
     809           1 :     hHeader.nDLon = static_cast<GInt16>(dfDLon);
     810             : 
     811             :     GByte abyBuf[BYN_HDR_SZ];
     812             : 
     813           1 :     header2buffer(&hHeader, abyBuf);
     814             : 
     815           1 :     const char *pszValue = GetMetadataItem("GLOBAL");
     816           1 :     if (pszValue != nullptr)
     817           0 :         hHeader.nGlobal = static_cast<GInt16>(atoi(pszValue));
     818             : 
     819           1 :     pszValue = GetMetadataItem("TYPE");
     820           1 :     if (pszValue != nullptr)
     821           0 :         hHeader.nType = static_cast<GInt16>(atoi(pszValue));
     822             : 
     823           1 :     pszValue = GetMetadataItem("DESCRIPTION");
     824           1 :     if (pszValue != nullptr)
     825           0 :         hHeader.nDescrip = static_cast<GInt16>(atoi(pszValue));
     826             : 
     827           1 :     pszValue = GetMetadataItem("SUBTYPE");
     828           1 :     if (pszValue != nullptr)
     829           0 :         hHeader.nSubType = static_cast<GInt16>(atoi(pszValue));
     830             : 
     831           1 :     pszValue = GetMetadataItem("WO");
     832           1 :     if (pszValue != nullptr)
     833           0 :         hHeader.dfWo = CPLAtof(pszValue);
     834             : 
     835           1 :     pszValue = GetMetadataItem("GM");
     836           1 :     if (pszValue != nullptr)
     837           0 :         hHeader.dfGM = CPLAtof(pszValue);
     838             : 
     839           1 :     pszValue = GetMetadataItem("TIDESYSTEM");
     840           1 :     if (pszValue != nullptr)
     841           0 :         hHeader.nTideSys = static_cast<GInt16>(atoi(pszValue));
     842             : 
     843           1 :     pszValue = GetMetadataItem("REALIZATION");
     844           1 :     if (pszValue != nullptr)
     845           0 :         hHeader.nRealiz = static_cast<GInt16>(atoi(pszValue));
     846             : 
     847           1 :     pszValue = GetMetadataItem("EPOCH");
     848           1 :     if (pszValue != nullptr)
     849           0 :         hHeader.dEpoch = static_cast<float>(CPLAtof(pszValue));
     850             : 
     851           1 :     pszValue = GetMetadataItem("PTTYPE");
     852           1 :     if (pszValue != nullptr)
     853           0 :         hHeader.nPtType = static_cast<GInt16>(atoi(pszValue));
     854             : 
     855           1 :     CPL_IGNORE_RET_VAL(VSIFSeekL(fpImage, 0, SEEK_SET));
     856           1 :     CPL_IGNORE_RET_VAL(VSIFWriteL(abyBuf, BYN_HDR_SZ, 1, fpImage));
     857             : 
     858             :     /* GDALPam metadata update */
     859             : 
     860           1 :     SetMetadataItem("GLOBAL", CPLSPrintf("%d", hHeader.nGlobal), "BYN");
     861           1 :     SetMetadataItem("TYPE", CPLSPrintf("%d", hHeader.nType), "BYN");
     862           1 :     SetMetadataItem("DESCRIPTION", CPLSPrintf("%d", hHeader.nDescrip), "BYN");
     863           1 :     SetMetadataItem("SUBTYPE", CPLSPrintf("%d", hHeader.nSubType), "BYN");
     864           1 :     SetMetadataItem("WO", CPLSPrintf("%g", hHeader.dfWo), "BYN");
     865           1 :     SetMetadataItem("GM", CPLSPrintf("%g", hHeader.dfGM), "BYN");
     866           1 :     SetMetadataItem("TIDESYSTEM", CPLSPrintf("%d", hHeader.nTideSys), "BYN");
     867           1 :     SetMetadataItem("REALIZATION", CPLSPrintf("%d", hHeader.nRealiz), "BYN");
     868           1 :     SetMetadataItem("EPOCH", CPLSPrintf("%g", hHeader.dEpoch), "BYN");
     869           1 :     SetMetadataItem("PTTYPE", CPLSPrintf("%d", hHeader.nPtType), "BYN");
     870           1 : }
     871             : 
     872             : /************************************************************************/
     873             : /*                          GDALRegister_BYN()                          */
     874             : /************************************************************************/
     875             : 
     876        1595 : void GDALRegister_BYN()
     877             : 
     878             : {
     879        1595 :     if (GDALGetDriverByName("BYN") != nullptr)
     880         302 :         return;
     881             : 
     882        1293 :     GDALDriver *poDriver = new GDALDriver();
     883             : 
     884        1293 :     poDriver->SetDescription("BYN");
     885        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     886        1293 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
     887        1293 :                               "Natural Resources Canada's Geoid");
     888        1293 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSIONS, "byn err");
     889        1293 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     890        1293 :     poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/byn.html");
     891        1293 :     poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Int16 Int32");
     892             : 
     893        1293 :     poDriver->pfnOpen = BYNDataset::Open;
     894        1293 :     poDriver->pfnIdentify = BYNDataset::Identify;
     895        1293 :     poDriver->pfnCreate = BYNDataset::Create;
     896             : 
     897        1293 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     898             : }