LCOV - code coverage report
Current view: top level - frmts/raw - gtxdataset.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 146 172 84.9 %
Date: 2025-07-09 17:50:03 Functions: 14 15 93.3 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Vertical Datum Transformation
       4             :  * Purpose:  Implementation of NOAA .gtx vertical datum shift file format.
       5             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2010, Frank Warmerdam
       9             :  * Copyright (c) 2010-2013, Even Rouault <even dot rouault at spatialys.com>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "cpl_string.h"
      15             : #include "gdal_frmts.h"
      16             : #include "ogr_srs_api.h"
      17             : #include "rawdataset.h"
      18             : 
      19             : #include <algorithm>
      20             : #include <limits>
      21             : 
      22             : /**
      23             : 
      24             : NOAA .GTX Vertical Datum Grid Shift Format
      25             : 
      26             : All values are bigendian
      27             : 
      28             : Header
      29             : ------
      30             : 
      31             : float64  latitude_of_origin
      32             : float64  longitude_of_origin (0-360)
      33             : float64  cell size (x?y?)
      34             : float64  cell size (x?y?)
      35             : int32    length in pixels
      36             : int32    width in pixels
      37             : 
      38             : Data
      39             : ----
      40             : 
      41             : float32  * width in pixels * length in pixels
      42             : 
      43             : Values are an offset in meters between two vertical datums.
      44             : 
      45             : **/
      46             : 
      47             : /************************************************************************/
      48             : /* ==================================================================== */
      49             : /*                              GTXDataset                              */
      50             : /* ==================================================================== */
      51             : /************************************************************************/
      52             : 
      53             : class GTXDataset final : public RawDataset
      54             : {
      55             :     VSILFILE *fpImage = nullptr;  // image data file.
      56             : 
      57             :     OGRSpatialReference m_oSRS{};
      58             :     GDALGeoTransform m_gt{};
      59             : 
      60             :     CPL_DISALLOW_COPY_ASSIGN(GTXDataset)
      61             : 
      62             :     CPLErr Close() override;
      63             : 
      64             :   public:
      65          47 :     GTXDataset()
      66          47 :     {
      67          47 :         m_oSRS.SetFromUserInput(SRS_WKT_WGS84_LAT_LONG);
      68          47 :         m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
      69          47 :     }
      70             : 
      71             :     ~GTXDataset() override;
      72             : 
      73             :     CPLErr GetGeoTransform(GDALGeoTransform &gt) const override;
      74             :     CPLErr SetGeoTransform(const GDALGeoTransform &gt) override;
      75             : 
      76           1 :     const OGRSpatialReference *GetSpatialRef() const override
      77             :     {
      78           1 :         return &m_oSRS;
      79             :     }
      80             : 
      81             :     static GDALDataset *Open(GDALOpenInfo *);
      82             :     static int Identify(GDALOpenInfo *);
      83             :     static GDALDataset *Create(const char *pszFilename, int nXSize, int nYSize,
      84             :                                int nBands, GDALDataType eType,
      85             :                                char **papszOptions);
      86             : };
      87             : 
      88             : /************************************************************************/
      89             : /* ==================================================================== */
      90             : /*                           GTXRasterBand                              */
      91             : /* ==================================================================== */
      92             : /************************************************************************/
      93             : 
      94             : class GTXRasterBand final : public RawRasterBand
      95             : {
      96             :     CPL_DISALLOW_COPY_ASSIGN(GTXRasterBand)
      97             : 
      98             :   public:
      99             :     GTXRasterBand(GDALDataset *poDS, int nBand, VSILFILE *fpRaw,
     100             :                   vsi_l_offset nImgOffset, int nPixelOffset, int nLineOffset,
     101             :                   GDALDataType eDataType, int bNativeOrder);
     102             : 
     103             :     ~GTXRasterBand() override;
     104             : 
     105             :     double GetNoDataValue(int *pbSuccess = nullptr) override;
     106             : };
     107             : 
     108             : /************************************************************************/
     109             : /*                            GTXRasterBand()                           */
     110             : /************************************************************************/
     111             : 
     112          47 : GTXRasterBand::GTXRasterBand(GDALDataset *poDSIn, int nBandIn,
     113             :                              VSILFILE *fpRawIn, vsi_l_offset nImgOffsetIn,
     114             :                              int nPixelOffsetIn, int nLineOffsetIn,
     115          47 :                              GDALDataType eDataTypeIn, int bNativeOrderIn)
     116             :     : RawRasterBand(poDSIn, nBandIn, fpRawIn, nImgOffsetIn, nPixelOffsetIn,
     117             :                     nLineOffsetIn, eDataTypeIn, bNativeOrderIn,
     118          47 :                     RawRasterBand::OwnFP::NO)
     119             : {
     120          47 : }
     121             : 
     122             : /************************************************************************/
     123             : /*                           ~GTXRasterBand()                           */
     124             : /************************************************************************/
     125             : 
     126          94 : GTXRasterBand::~GTXRasterBand()
     127             : {
     128          94 : }
     129             : 
     130             : /************************************************************************/
     131             : /*                           GetNoDataValue()                           */
     132             : /************************************************************************/
     133             : 
     134           0 : double GTXRasterBand::GetNoDataValue(int *pbSuccess)
     135             : {
     136           0 :     if (pbSuccess)
     137           0 :         *pbSuccess = TRUE;
     138           0 :     int bSuccess = FALSE;
     139           0 :     double dfNoData = GDALPamRasterBand::GetNoDataValue(&bSuccess);
     140           0 :     if (bSuccess)
     141             :     {
     142           0 :         return dfNoData;
     143             :     }
     144           0 :     return -88.8888;
     145             : }
     146             : 
     147             : /************************************************************************/
     148             : /* ==================================================================== */
     149             : /*                              GTXDataset                              */
     150             : /* ==================================================================== */
     151             : /************************************************************************/
     152             : 
     153             : /************************************************************************/
     154             : /*                            ~GTXDataset()                             */
     155             : /************************************************************************/
     156             : 
     157          94 : GTXDataset::~GTXDataset()
     158             : 
     159             : {
     160          47 :     GTXDataset::Close();
     161          94 : }
     162             : 
     163             : /************************************************************************/
     164             : /*                              Close()                                 */
     165             : /************************************************************************/
     166             : 
     167          85 : CPLErr GTXDataset::Close()
     168             : {
     169          85 :     CPLErr eErr = CE_None;
     170          85 :     if (nOpenFlags != OPEN_FLAGS_CLOSED)
     171             :     {
     172          47 :         if (GTXDataset::FlushCache(true) != CE_None)
     173           9 :             eErr = CE_Failure;
     174             : 
     175          47 :         if (fpImage)
     176             :         {
     177          47 :             if (VSIFCloseL(fpImage) != 0)
     178             :             {
     179           0 :                 CPLError(CE_Failure, CPLE_FileIO, "I/O error");
     180           0 :                 eErr = CE_Failure;
     181             :             }
     182             :         }
     183             : 
     184          47 :         if (GDALPamDataset::Close() != CE_None)
     185           0 :             eErr = CE_Failure;
     186             :     }
     187          85 :     return eErr;
     188             : }
     189             : 
     190             : /************************************************************************/
     191             : /*                              Identify()                              */
     192             : /************************************************************************/
     193             : 
     194       57716 : int GTXDataset::Identify(GDALOpenInfo *poOpenInfo)
     195             : 
     196             : {
     197       57716 :     if (poOpenInfo->nHeaderBytes < 40)
     198       54041 :         return FALSE;
     199             : 
     200        3675 :     if (!poOpenInfo->IsExtensionEqualToCI("gtx"))
     201        3582 :         return FALSE;
     202             : 
     203          94 :     return TRUE;
     204             : }
     205             : 
     206             : /************************************************************************/
     207             : /*                                Open()                                */
     208             : /************************************************************************/
     209             : 
     210          47 : GDALDataset *GTXDataset::Open(GDALOpenInfo *poOpenInfo)
     211             : 
     212             : {
     213          47 :     if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
     214           0 :         return nullptr;
     215             : 
     216             :     /* -------------------------------------------------------------------- */
     217             :     /*      Create a corresponding GDALDataset.                             */
     218             :     /* -------------------------------------------------------------------- */
     219          94 :     auto poDS = std::make_unique<GTXDataset>();
     220             : 
     221          47 :     poDS->eAccess = poOpenInfo->eAccess;
     222          47 :     std::swap(poDS->fpImage, poOpenInfo->fpL);
     223             : 
     224             :     /* -------------------------------------------------------------------- */
     225             :     /*      Read the header.                                                */
     226             :     /* -------------------------------------------------------------------- */
     227          47 :     double gt[6] = {0};
     228          47 :     CPL_IGNORE_RET_VAL(VSIFReadL(&gt[3], 8, 1, poDS->fpImage));
     229          47 :     CPL_IGNORE_RET_VAL(VSIFReadL(&gt[0], 8, 1, poDS->fpImage));
     230          47 :     CPL_IGNORE_RET_VAL(VSIFReadL(&gt[5], 8, 1, poDS->fpImage));
     231          47 :     CPL_IGNORE_RET_VAL(VSIFReadL(&gt[1], 8, 1, poDS->fpImage));
     232             : 
     233          47 :     CPL_IGNORE_RET_VAL(VSIFReadL(&(poDS->nRasterYSize), 4, 1, poDS->fpImage));
     234          47 :     CPL_IGNORE_RET_VAL(VSIFReadL(&(poDS->nRasterXSize), 4, 1, poDS->fpImage));
     235             : 
     236          47 :     CPL_MSBPTR32(&(poDS->nRasterYSize));
     237          47 :     CPL_MSBPTR32(&(poDS->nRasterXSize));
     238             : 
     239          47 :     CPL_MSBPTR64(&gt[0]);
     240          47 :     CPL_MSBPTR64(&gt[1]);
     241          47 :     CPL_MSBPTR64(&gt[3]);
     242          47 :     CPL_MSBPTR64(&gt[5]);
     243             : 
     244          47 :     poDS->m_gt = GDALGeoTransform(gt);
     245          47 :     poDS->m_gt[3] +=
     246          47 :         poDS->m_gt[5] * (static_cast<double>(poDS->nRasterYSize) - 1);
     247             : 
     248          47 :     poDS->m_gt[0] -= poDS->m_gt[1] * 0.5;
     249          47 :     poDS->m_gt[3] += poDS->m_gt[5] * 0.5;
     250             : 
     251          47 :     poDS->m_gt[5] *= -1;
     252             : 
     253          47 :     if (CPLFetchBool(poOpenInfo->papszOpenOptions,
     254             :                      "SHIFT_ORIGIN_IN_MINUS_180_PLUS_180", false))
     255             :     {
     256           0 :         if (poDS->m_gt[0] < -180.0 - poDS->m_gt[1])
     257           0 :             poDS->m_gt[0] += 360.0;
     258           0 :         else if (poDS->m_gt[0] > 180.0)
     259           0 :             poDS->m_gt[0] -= 360.0;
     260             :     }
     261             : 
     262          94 :     if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize) ||
     263          47 :         static_cast<vsi_l_offset>(poDS->nRasterXSize) * poDS->nRasterYSize >
     264          47 :             std::numeric_limits<vsi_l_offset>::max() / sizeof(double))
     265             :     {
     266           0 :         return nullptr;
     267             :     }
     268             : 
     269             :     /* -------------------------------------------------------------------- */
     270             :     /*      Guess the data type. Since October 1, 2009, it should be        */
     271             :     /*      Float32. Before it was double.                                  */
     272             :     /* -------------------------------------------------------------------- */
     273          47 :     CPL_IGNORE_RET_VAL(VSIFSeekL(poDS->fpImage, 0, SEEK_END));
     274          47 :     const vsi_l_offset nSize = VSIFTellL(poDS->fpImage);
     275             : 
     276          47 :     GDALDataType eDT = GDT_Float32;
     277          94 :     if (nSize - 40 == sizeof(double) *
     278          47 :                           static_cast<vsi_l_offset>(poDS->nRasterXSize) *
     279          47 :                           poDS->nRasterYSize)
     280           0 :         eDT = GDT_Float64;
     281          47 :     const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
     282          47 :     if (nDTSize <= 0 || poDS->nRasterXSize > INT_MAX / nDTSize)
     283             :     {
     284           0 :         return nullptr;
     285             :     }
     286             : 
     287             :     /* -------------------------------------------------------------------- */
     288             :     /*      Create band information object.                                 */
     289             :     /* -------------------------------------------------------------------- */
     290             :     auto poBand = std::make_unique<GTXRasterBand>(
     291          47 :         poDS.get(), 1, poDS->fpImage,
     292          47 :         static_cast<vsi_l_offset>(poDS->nRasterYSize - 1) * poDS->nRasterXSize *
     293          47 :                 nDTSize +
     294             :             40,
     295         141 :         nDTSize, poDS->nRasterXSize * -nDTSize, eDT, !CPL_IS_LSB);
     296          47 :     if (!poBand->IsValid())
     297           0 :         return nullptr;
     298          47 :     poDS->SetBand(1, std::move(poBand));
     299             : 
     300             :     /* -------------------------------------------------------------------- */
     301             :     /*      Initialize any PAM information.                                 */
     302             :     /* -------------------------------------------------------------------- */
     303          47 :     poDS->SetDescription(poOpenInfo->pszFilename);
     304          47 :     poDS->TryLoadXML();
     305             : 
     306             :     /* -------------------------------------------------------------------- */
     307             :     /*      Check for overviews.                                            */
     308             :     /* -------------------------------------------------------------------- */
     309          47 :     poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
     310             : 
     311          47 :     return poDS.release();
     312             : }
     313             : 
     314             : /************************************************************************/
     315             : /*                          GetGeoTransform()                           */
     316             : /************************************************************************/
     317             : 
     318           3 : CPLErr GTXDataset::GetGeoTransform(GDALGeoTransform &gt) const
     319             : 
     320             : {
     321           3 :     gt = m_gt;
     322           3 :     return CE_None;
     323             : }
     324             : 
     325             : /************************************************************************/
     326             : /*                          SetGeoTransform()                           */
     327             : /************************************************************************/
     328             : 
     329          33 : CPLErr GTXDataset::SetGeoTransform(const GDALGeoTransform &gt)
     330             : {
     331          33 :     if (gt[2] != 0.0 || gt[4] != 0.0)
     332             :     {
     333           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     334             :                  "Attempt to write skewed or rotated geotransform to gtx.");
     335           0 :         return CE_Failure;
     336             :     }
     337             : 
     338          33 :     m_gt = gt;
     339             : 
     340          33 :     const double dfXOrigin = m_gt[0] + 0.5 * m_gt[1];
     341          33 :     const double dfYOrigin = m_gt[3] + (nRasterYSize - 0.5) * m_gt[5];
     342          33 :     const double dfWidth = m_gt[1];
     343          33 :     const double dfHeight = -m_gt[5];
     344             : 
     345          33 :     unsigned char header[32] = {'\0'};
     346          33 :     memcpy(header + 0, &dfYOrigin, 8);
     347          33 :     CPL_MSBPTR64(header + 0);
     348             : 
     349          33 :     memcpy(header + 8, &dfXOrigin, 8);
     350          33 :     CPL_MSBPTR64(header + 8);
     351             : 
     352          33 :     memcpy(header + 16, &dfHeight, 8);
     353          33 :     CPL_MSBPTR64(header + 16);
     354             : 
     355          33 :     memcpy(header + 24, &dfWidth, 8);
     356          33 :     CPL_MSBPTR64(header + 24);
     357             : 
     358          66 :     if (VSIFSeekL(fpImage, SEEK_SET, 0) != 0 ||
     359          33 :         VSIFWriteL(header, 32, 1, fpImage) != 1)
     360             :     {
     361           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     362             :                  "Attempt to write geotransform header to GTX failed.");
     363           0 :         return CE_Failure;
     364             :     }
     365             : 
     366          33 :     return CE_None;
     367             : }
     368             : 
     369             : /************************************************************************/
     370             : /*                               Create()                               */
     371             : /************************************************************************/
     372             : 
     373          83 : GDALDataset *GTXDataset::Create(const char *pszFilename, int nXSize, int nYSize,
     374             :                                 int /* nBands */, GDALDataType eType,
     375             :                                 char ** /* papszOptions */)
     376             : {
     377          83 :     if (eType != GDT_Float32)
     378             :     {
     379          46 :         CPLError(CE_Failure, CPLE_AppDefined,
     380             :                  "Attempt to create gtx file with unsupported data type '%s'.",
     381             :                  GDALGetDataTypeName(eType));
     382          46 :         return nullptr;
     383             :     }
     384             : 
     385          37 :     if (!EQUAL(CPLGetExtensionSafe(pszFilename).c_str(), "gtx"))
     386             :     {
     387           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     388             :                  "Attempt to create gtx file with extension other than gtx.");
     389           0 :         return nullptr;
     390             :     }
     391             : 
     392             :     /* -------------------------------------------------------------------- */
     393             :     /*      Try to create the file.                                         */
     394             :     /* -------------------------------------------------------------------- */
     395          37 :     VSILFILE *fp = VSIFOpenL(pszFilename, "wb");
     396          37 :     if (fp == nullptr)
     397             :     {
     398           3 :         CPLError(CE_Failure, CPLE_OpenFailed,
     399             :                  "Attempt to create file `%s' failed.\n", pszFilename);
     400           3 :         return nullptr;
     401             :     }
     402             : 
     403             :     /* -------------------------------------------------------------------- */
     404             :     /*      Write out the header with stub georeferencing.                  */
     405             :     /* -------------------------------------------------------------------- */
     406             : 
     407          34 :     unsigned char header[40] = {'\0'};
     408          34 :     double dfYOrigin = 0.0;
     409          34 :     memcpy(header + 0, &dfYOrigin, 8);
     410          34 :     CPL_MSBPTR64(header + 0);
     411             : 
     412          34 :     double dfXOrigin = 0.0;
     413          34 :     memcpy(header + 8, &dfXOrigin, 8);
     414          34 :     CPL_MSBPTR64(header + 8);
     415             : 
     416          34 :     double dfYSize = 0.01;
     417          34 :     memcpy(header + 16, &dfYSize, 8);
     418          34 :     CPL_MSBPTR64(header + 16);
     419             : 
     420          34 :     double dfXSize = 0.01;
     421          34 :     memcpy(header + 24, &dfXSize, 8);
     422          34 :     CPL_MSBPTR64(header + 24);
     423             : 
     424          34 :     GInt32 nYSize32 = nYSize;
     425          34 :     memcpy(header + 32, &nYSize32, 4);
     426          34 :     CPL_MSBPTR32(header + 32);
     427             : 
     428          34 :     GInt32 nXSize32 = nXSize;
     429          34 :     memcpy(header + 36, &nXSize32, 4);
     430          34 :     CPL_MSBPTR32(header + 36);
     431             : 
     432          34 :     CPL_IGNORE_RET_VAL(VSIFWriteL(header, 40, 1, fp));
     433          34 :     CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
     434             : 
     435          34 :     return GDALDataset::FromHandle(GDALOpen(pszFilename, GA_Update));
     436             : }
     437             : 
     438             : /************************************************************************/
     439             : /*                          GDALRegister_GTX()                          */
     440             : /************************************************************************/
     441             : 
     442        1935 : void GDALRegister_GTX()
     443             : 
     444             : {
     445        1935 :     if (GDALGetDriverByName("GTX") != nullptr)
     446         282 :         return;
     447             : 
     448        1653 :     GDALDriver *poDriver = new GDALDriver();
     449             : 
     450        1653 :     poDriver->SetDescription("GTX");
     451        1653 :     poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
     452        1653 :     poDriver->SetMetadataItem(GDAL_DMD_LONGNAME, "NOAA Vertical Datum .GTX");
     453        1653 :     poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "gtx");
     454        1653 :     poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
     455             :     // poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
     456             :     //                            "frmt_various.html#GTX" );
     457        1653 :     poDriver->SetMetadataItem(
     458             :         GDAL_DMD_OPENOPTIONLIST,
     459             :         "<OpenOptionList>"
     460             :         "   <Option name='SHIFT_ORIGIN_IN_MINUS_180_PLUS_180' type='boolean' "
     461             :         "description='Whether to apply a +/-360 deg shift to the longitude of "
     462             :         "the top left corner so that it is in the [-180,180] range' "
     463             :         "default='NO'/>"
     464        1653 :         "</OpenOptionList>");
     465             : 
     466        1653 :     poDriver->SetMetadataItem(GDAL_DMD_CREATIONDATATYPES, "Float32");
     467             : 
     468        1653 :     poDriver->pfnOpen = GTXDataset::Open;
     469        1653 :     poDriver->pfnIdentify = GTXDataset::Identify;
     470        1653 :     poDriver->pfnCreate = GTXDataset::Create;
     471             : 
     472        1653 :     GetGDALDriverManager()->RegisterDriver(poDriver);
     473             : }

Generated by: LCOV version 1.14