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