LCOV - code coverage report
Current view: top level - alg - gdalapplyverticalshiftgrid.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 190 254 74.8 %
Date: 2025-07-09 17:50:03 Functions: 13 15 86.7 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL algorithms
       4             :  * Purpose:  Apply vertical shift grid
       5             :  * Author:   Even Rouault, even.rouault at spatialys.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2017, Even Rouault <even.rouault at spatialys.com>
       9             :  *
      10             :  * SPDX-License-Identifier: MIT
      11             :  ****************************************************************************/
      12             : 
      13             : #include "cpl_string.h"
      14             : #include "gdal.h"
      15             : #include "gdal_alg.h"
      16             : #include "gdal_alg_priv.h"
      17             : #include "gdal_priv.h"
      18             : #include "gdal_utils.h"
      19             : #include "gdalwarper.h"
      20             : #include "vrtdataset.h"
      21             : #include "ogr_spatialref.h"
      22             : 
      23             : #include "proj.h"
      24             : 
      25             : #include <cmath>
      26             : #include <limits>
      27             : 
      28             : /************************************************************************/
      29             : /*                        GDALApplyVSGDataset                           */
      30             : /************************************************************************/
      31             : 
      32             : class GDALApplyVSGDataset final : public GDALDataset
      33             : {
      34             :     friend class GDALApplyVSGRasterBand;
      35             : 
      36             :     GDALDataset *m_poSrcDataset = nullptr;
      37             :     GDALDataset *m_poReprojectedGrid = nullptr;
      38             :     bool m_bInverse = false;
      39             :     double m_dfSrcUnitToMeter = 0.0;
      40             :     double m_dfDstUnitToMeter = 0.0;
      41             : 
      42             :     CPL_DISALLOW_COPY_ASSIGN(GDALApplyVSGDataset)
      43             : 
      44             :   public:
      45             :     GDALApplyVSGDataset(GDALDataset *poSrcDataset,
      46             :                         GDALDataset *poReprojectedGrid, GDALDataType eDT,
      47             :                         bool bInverse, double dfSrcUnitToMeter,
      48             :                         double dfDstUnitToMeter, int nBlockSize);
      49             :     virtual ~GDALApplyVSGDataset();
      50             : 
      51             :     virtual int CloseDependentDatasets() override;
      52             : 
      53             :     virtual CPLErr GetGeoTransform(GDALGeoTransform &gt) const override;
      54             :     virtual const OGRSpatialReference *GetSpatialRef() const override;
      55             : 
      56             :     bool IsInitOK();
      57             : };
      58             : 
      59             : /************************************************************************/
      60             : /*                       GDALApplyVSGRasterBand                         */
      61             : /************************************************************************/
      62             : 
      63             : class GDALApplyVSGRasterBand final : public GDALRasterBand
      64             : {
      65             :     friend class GDALApplyVSGDataset;
      66             : 
      67             :     float *m_pafSrcData = nullptr;
      68             :     float *m_pafGridData = nullptr;
      69             : 
      70             :     CPL_DISALLOW_COPY_ASSIGN(GDALApplyVSGRasterBand)
      71             : 
      72             :   public:
      73             :     GDALApplyVSGRasterBand(GDALDataType eDT, int nBlockSize);
      74             :     virtual ~GDALApplyVSGRasterBand();
      75             : 
      76             :     virtual CPLErr IReadBlock(int nBlockXOff, int nBlockYOff,
      77             :                               void *pData) override;
      78             :     virtual double GetNoDataValue(int *pbSuccess) override;
      79             : };
      80             : 
      81             : /************************************************************************/
      82             : /*                        GDALApplyVSGDataset()                         */
      83             : /************************************************************************/
      84             : 
      85          13 : GDALApplyVSGDataset::GDALApplyVSGDataset(GDALDataset *poSrcDataset,
      86             :                                          GDALDataset *poReprojectedGrid,
      87             :                                          GDALDataType eDT, bool bInverse,
      88             :                                          double dfSrcUnitToMeter,
      89             :                                          double dfDstUnitToMeter,
      90          13 :                                          int nBlockSize)
      91             :     : m_poSrcDataset(poSrcDataset), m_poReprojectedGrid(poReprojectedGrid),
      92             :       m_bInverse(bInverse), m_dfSrcUnitToMeter(dfSrcUnitToMeter),
      93          13 :       m_dfDstUnitToMeter(dfDstUnitToMeter)
      94             : {
      95          13 :     m_poSrcDataset->Reference();
      96          13 :     m_poReprojectedGrid->Reference();
      97             : 
      98          13 :     nRasterXSize = poSrcDataset->GetRasterXSize();
      99          13 :     nRasterYSize = poSrcDataset->GetRasterYSize();
     100          13 :     SetBand(1, new GDALApplyVSGRasterBand(eDT, nBlockSize));
     101          13 : }
     102             : 
     103             : /************************************************************************/
     104             : /*                       ~GDALApplyVSGDataset()                         */
     105             : /************************************************************************/
     106             : 
     107          26 : GDALApplyVSGDataset::~GDALApplyVSGDataset()
     108             : {
     109          13 :     GDALApplyVSGDataset::CloseDependentDatasets();
     110          26 : }
     111             : 
     112             : /************************************************************************/
     113             : /*                     CloseDependentDatasets()                         */
     114             : /************************************************************************/
     115             : 
     116          13 : int GDALApplyVSGDataset::CloseDependentDatasets()
     117             : {
     118          13 :     bool bRet = false;
     119          13 :     if (m_poSrcDataset != nullptr)
     120             :     {
     121          13 :         if (m_poSrcDataset->ReleaseRef())
     122             :         {
     123           8 :             bRet = true;
     124             :         }
     125          13 :         m_poSrcDataset = nullptr;
     126             :     }
     127          13 :     if (m_poReprojectedGrid != nullptr)
     128             :     {
     129          13 :         if (m_poReprojectedGrid->ReleaseRef())
     130             :         {
     131          13 :             bRet = true;
     132             :         }
     133          13 :         m_poReprojectedGrid = nullptr;
     134             :     }
     135          13 :     return bRet;
     136             : }
     137             : 
     138             : /************************************************************************/
     139             : /*                          GetGeoTransform()                           */
     140             : /************************************************************************/
     141             : 
     142           2 : CPLErr GDALApplyVSGDataset::GetGeoTransform(GDALGeoTransform &gt) const
     143             : {
     144           2 :     return m_poSrcDataset->GetGeoTransform(gt);
     145             : }
     146             : 
     147             : /************************************************************************/
     148             : /*                          GetSpatialRef()                             */
     149             : /************************************************************************/
     150             : 
     151           2 : const OGRSpatialReference *GDALApplyVSGDataset::GetSpatialRef() const
     152             : {
     153           2 :     return m_poSrcDataset->GetSpatialRef();
     154             : }
     155             : 
     156             : /************************************************************************/
     157             : /*                             IsInitOK()                               */
     158             : /************************************************************************/
     159             : 
     160          13 : bool GDALApplyVSGDataset::IsInitOK()
     161             : {
     162             :     GDALApplyVSGRasterBand *poBand =
     163          13 :         reinterpret_cast<GDALApplyVSGRasterBand *>(GetRasterBand(1));
     164          13 :     return poBand->m_pafSrcData != nullptr && poBand->m_pafGridData != nullptr;
     165             : }
     166             : 
     167             : /************************************************************************/
     168             : /*                       GDALApplyVSGRasterBand()                       */
     169             : /************************************************************************/
     170             : 
     171          13 : GDALApplyVSGRasterBand::GDALApplyVSGRasterBand(GDALDataType eDT, int nBlockSize)
     172             : {
     173          13 :     eDataType = eDT;
     174          13 :     nBlockXSize = nBlockSize;
     175          13 :     nBlockYSize = nBlockSize;
     176          13 :     m_pafSrcData = static_cast<float *>(
     177          13 :         VSI_MALLOC3_VERBOSE(nBlockXSize, nBlockYSize, sizeof(float)));
     178          13 :     m_pafGridData = static_cast<float *>(
     179          13 :         VSI_MALLOC3_VERBOSE(nBlockXSize, nBlockYSize, sizeof(float)));
     180          13 : }
     181             : 
     182             : /************************************************************************/
     183             : /*                      ~GDALApplyVSGRasterBand()                       */
     184             : /************************************************************************/
     185             : 
     186          26 : GDALApplyVSGRasterBand::~GDALApplyVSGRasterBand()
     187             : {
     188          13 :     VSIFree(m_pafSrcData);
     189          13 :     VSIFree(m_pafGridData);
     190          26 : }
     191             : 
     192             : /************************************************************************/
     193             : /*                           GetNoDataValue()                           */
     194             : /************************************************************************/
     195             : 
     196          19 : double GDALApplyVSGRasterBand::GetNoDataValue(int *pbSuccess)
     197             : {
     198          19 :     GDALApplyVSGDataset *poGDS = reinterpret_cast<GDALApplyVSGDataset *>(poDS);
     199          19 :     return poGDS->m_poSrcDataset->GetRasterBand(1)->GetNoDataValue(pbSuccess);
     200             : }
     201             : 
     202             : /************************************************************************/
     203             : /*                              IReadBlock()                            */
     204             : /************************************************************************/
     205             : 
     206          17 : CPLErr GDALApplyVSGRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,
     207             :                                           void *pData)
     208             : {
     209          17 :     GDALApplyVSGDataset *poGDS = reinterpret_cast<GDALApplyVSGDataset *>(poDS);
     210             : 
     211          17 :     const int nXOff = nBlockXOff * nBlockXSize;
     212          34 :     const int nReqXSize = (nXOff > nRasterXSize - nBlockXSize)
     213          17 :                               ? nRasterXSize - nXOff
     214             :                               : nBlockXSize;
     215          17 :     const int nYOff = nBlockYOff * nBlockYSize;
     216          34 :     const int nReqYSize = (nYOff > nRasterYSize - nBlockYSize)
     217          17 :                               ? nRasterYSize - nYOff
     218             :                               : nBlockYSize;
     219             : 
     220          17 :     CPLErr eErr = poGDS->m_poSrcDataset->GetRasterBand(1)->RasterIO(
     221          17 :         GF_Read, nXOff, nYOff, nReqXSize, nReqYSize, m_pafSrcData, nReqXSize,
     222          17 :         nReqYSize, GDT_Float32, sizeof(float), nBlockXSize * sizeof(float),
     223             :         nullptr);
     224          17 :     if (eErr == CE_None)
     225          34 :         eErr = poGDS->m_poReprojectedGrid->GetRasterBand(1)->RasterIO(
     226          17 :             GF_Read, nXOff, nYOff, nReqXSize, nReqYSize, m_pafGridData,
     227             :             nReqXSize, nReqYSize, GDT_Float32, sizeof(float),
     228          17 :             nBlockXSize * sizeof(float), nullptr);
     229          17 :     if (eErr == CE_None)
     230             :     {
     231          17 :         const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);
     232          17 :         int bHasNoData = FALSE;
     233          17 :         float fNoDataValue = static_cast<float>(GetNoDataValue(&bHasNoData));
     234         279 :         for (int iY = 0; iY < nReqYSize; iY++)
     235             :         {
     236        4666 :             for (int iX = 0; iX < nReqXSize; iX++)
     237             :             {
     238        4404 :                 const float fSrcVal = m_pafSrcData[iY * nBlockXSize + iX];
     239        4404 :                 const float fGridVal = m_pafGridData[iY * nBlockXSize + iX];
     240        4404 :                 if (bHasNoData && fSrcVal == fNoDataValue)
     241             :                 {
     242             :                 }
     243        4403 :                 else if (std::isinf(fGridVal))
     244             :                 {
     245           2 :                     CPLError(CE_Failure, CPLE_AppDefined,
     246             :                              "Missing vertical grid value at source (%d,%d)",
     247             :                              nXOff + iX, nYOff + iY);
     248           2 :                     return CE_Failure;
     249             :                 }
     250        4401 :                 else if (poGDS->m_bInverse)
     251             :                 {
     252         800 :                     m_pafSrcData[iY * nBlockXSize + iX] = static_cast<float>(
     253         800 :                         (fSrcVal * poGDS->m_dfSrcUnitToMeter - fGridVal) /
     254         800 :                         poGDS->m_dfDstUnitToMeter);
     255             :                 }
     256             :                 else
     257             :                 {
     258        3601 :                     m_pafSrcData[iY * nBlockXSize + iX] = static_cast<float>(
     259        3601 :                         (fSrcVal * poGDS->m_dfSrcUnitToMeter + fGridVal) /
     260        3601 :                         poGDS->m_dfDstUnitToMeter);
     261             :                 }
     262             :             }
     263         262 :             GDALCopyWords(
     264         262 :                 m_pafSrcData + iY * nBlockXSize, GDT_Float32, sizeof(float),
     265         262 :                 static_cast<GByte *>(pData) + iY * nBlockXSize * nDTSize,
     266             :                 eDataType, nDTSize, nReqXSize);
     267             :         }
     268             :     }
     269             : 
     270          15 :     return eErr;
     271             : }
     272             : 
     273             : /************************************************************************/
     274             : /*                      GDALApplyVerticalShiftGrid()                    */
     275             : /************************************************************************/
     276             : 
     277             : /** Apply a vertical shift grid to a source (DEM typically) dataset.
     278             :  *
     279             :  * hGridDataset will typically use WGS84 as horizontal datum (but this is
     280             :  * not a requirement) and its values are the values to add to go from geoid
     281             :  * elevations to WGS84 ellipsoidal heights.
     282             :  *
     283             :  * hGridDataset will be on-the-fly reprojected and resampled to the projection
     284             :  * and resolution of hSrcDataset, using bilinear resampling by default.
     285             :  *
     286             :  * Both hSrcDataset and hGridDataset must be single band datasets, and have
     287             :  * a valid geotransform and projection.
     288             :  *
     289             :  * On success, a reference will be taken on hSrcDataset and hGridDataset.
     290             :  * Reference counting semantics on the source and grid datasets should be
     291             :  * honoured. That is, don't just GDALClose() it, unless it was opened with
     292             :  * GDALOpenShared(), but rather use GDALReleaseDataset() if wanting to
     293             :  * immediately release the reference(s) and make the returned dataset the
     294             :  * owner of them.
     295             :  *
     296             :  * Valid use cases:
     297             :  *
     298             :  * \code
     299             :  * hSrcDataset = GDALOpen(...)
     300             :  * hGridDataset = GDALOpen(...)
     301             :  * hDstDataset = GDALApplyVerticalShiftGrid(hSrcDataset, hGridDataset, ...)
     302             :  * GDALReleaseDataset(hSrcDataset);
     303             :  * GDALReleaseDataset(hGridDataset);
     304             :  * if( hDstDataset )
     305             :  * {
     306             :  *     // Do things with hDstDataset
     307             :  *     GDALClose(hDstDataset) // will close hSrcDataset and hGridDataset
     308             :  * }
     309             :  * \endcode
     310             : 
     311             :  *
     312             :  * @param hSrcDataset source (DEM) dataset. Must not be NULL.
     313             :  * @param hGridDataset vertical grid shift dataset. Must not be NULL.
     314             :  * @param bInverse if set to FALSE, hGridDataset values will be added to
     315             :  *                 hSrcDataset. If set to TRUE, they will be subtracted.
     316             :  * @param dfSrcUnitToMeter the factor to convert values from hSrcDataset to
     317             :  *                         meters (1.0 if source values are in meter).
     318             :  * @param dfDstUnitToMeter the factor to convert shifted values from meter
     319             :  *                          (1.0 if output values must be in meter).
     320             :  * @param papszOptions list of options, or NULL. Supported options are:
     321             :  * <ul>
     322             :  * <li>RESAMPLING=NEAREST/BILINEAR/CUBIC. Defaults to BILINEAR.</li>
     323             :  * <li>MAX_ERROR=val. Maximum error measured in input pixels that is allowed in
     324             :  * approximating the transformation (0.0 for exact calculations). Defaults
     325             :  * to 0.125</li>
     326             :  * <li>DATATYPE=Byte/UInt16/Int16/Float32/Float64. Output data type. If not
     327             :  * specified will be the same as the one of hSrcDataset.
     328             :  * <li>ERROR_ON_MISSING_VERT_SHIFT=YES/NO. Whether a missing/nodata value in
     329             :  * hGridDataset should cause I/O requests to fail. Default is NO (in which case
     330             :  * 0 will be used)
     331             :  * <li>SRC_SRS=srs_def. Override projection on hSrcDataset;
     332             :  * </ul>
     333             :  *
     334             :  * @return a new dataset corresponding to hSrcDataset adjusted with
     335             :  * hGridDataset, or NULL. If not NULL, it must be closed with GDALClose().
     336             :  *
     337             :  * @since GDAL 2.2
     338             :  * @deprecated GDAL 3.4. Will be removed in GDAL 4.0. This function was used
     339             :  *             by gdalwarp initially, but is no longer needed.
     340             :  */
     341          23 : GDALDatasetH GDALApplyVerticalShiftGrid(GDALDatasetH hSrcDataset,
     342             :                                         GDALDatasetH hGridDataset, int bInverse,
     343             :                                         double dfSrcUnitToMeter,
     344             :                                         double dfDstUnitToMeter,
     345             :                                         const char *const *papszOptions)
     346             : {
     347          23 :     VALIDATE_POINTER1(hSrcDataset, "GDALApplyVerticalShiftGrid", nullptr);
     348          23 :     VALIDATE_POINTER1(hGridDataset, "GDALApplyVerticalShiftGrid", nullptr);
     349             : 
     350          23 :     GDALGeoTransform srcGT;
     351          23 :     if (GDALDataset::FromHandle(hSrcDataset)->GetGeoTransform(srcGT) != CE_None)
     352             :     {
     353           1 :         CPLError(CE_Failure, CPLE_NotSupported,
     354             :                  "Source dataset has no geotransform.");
     355           1 :         return nullptr;
     356             :     }
     357          22 :     const char *pszSrcProjection = CSLFetchNameValue(papszOptions, "SRC_SRS");
     358          44 :     OGRSpatialReference oSrcSRS;
     359          22 :     if (pszSrcProjection != nullptr && pszSrcProjection[0] != '\0')
     360             :     {
     361           0 :         oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
     362           0 :         oSrcSRS.SetFromUserInput(pszSrcProjection);
     363             :     }
     364             :     else
     365             :     {
     366          22 :         auto poSRS = GDALDataset::FromHandle(hSrcDataset)->GetSpatialRef();
     367          22 :         if (poSRS)
     368          21 :             oSrcSRS = *poSRS;
     369             :     }
     370             : 
     371          22 :     if (oSrcSRS.IsCompound())
     372             :     {
     373           0 :         oSrcSRS.StripVertical();
     374             :     }
     375             : 
     376          22 :     if (oSrcSRS.IsEmpty())
     377             :     {
     378           1 :         CPLError(CE_Failure, CPLE_NotSupported,
     379             :                  "Source dataset has no projection.");
     380           1 :         return nullptr;
     381             :     }
     382          21 :     if (GDALGetRasterCount(hSrcDataset) != 1)
     383             :     {
     384           1 :         CPLError(CE_Failure, CPLE_NotSupported,
     385             :                  "Only single band source dataset is supported.");
     386           1 :         return nullptr;
     387             :     }
     388             : 
     389          20 :     GDALGeoTransform gridGT;
     390          20 :     if (GDALDataset::FromHandle(hGridDataset)->GetGeoTransform(gridGT) !=
     391             :         CE_None)
     392             :     {
     393           1 :         CPLError(CE_Failure, CPLE_NotSupported,
     394             :                  "Grid dataset has no geotransform.");
     395           1 :         return nullptr;
     396             :     }
     397             : 
     398          19 :     OGRSpatialReferenceH hGridSRS = GDALGetSpatialRef(hGridDataset);
     399          19 :     if (hGridSRS == nullptr)
     400             :     {
     401           1 :         CPLError(CE_Failure, CPLE_NotSupported,
     402             :                  "Grid dataset has no projection.");
     403           1 :         return nullptr;
     404             :     }
     405          18 :     if (GDALGetRasterCount(hGridDataset) != 1)
     406             :     {
     407           1 :         CPLError(CE_Failure, CPLE_NotSupported,
     408             :                  "Only single band grid dataset is supported.");
     409           1 :         return nullptr;
     410             :     }
     411             : 
     412          17 :     GDALDataType eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDataset, 1));
     413          17 :     const char *pszDataType = CSLFetchNameValue(papszOptions, "DATATYPE");
     414          17 :     if (pszDataType)
     415           2 :         eDT = GDALGetDataTypeByName(pszDataType);
     416          17 :     if (eDT == GDT_Unknown)
     417             :     {
     418           1 :         CPLError(CE_Failure, CPLE_NotSupported, "Invalid DATATYPE=%s",
     419             :                  pszDataType);
     420           1 :         return nullptr;
     421             :     }
     422             : 
     423          16 :     const int nSrcXSize = GDALGetRasterXSize(hSrcDataset);
     424          16 :     const int nSrcYSize = GDALGetRasterYSize(hSrcDataset);
     425             : 
     426          16 :     double dfWestLongitudeDeg = 0.0;
     427          16 :     double dfSouthLatitudeDeg = 0.0;
     428          16 :     double dfEastLongitudeDeg = 0.0;
     429          16 :     double dfNorthLatitudeDeg = 0.0;
     430          16 :     GDALComputeAreaOfInterest(&oSrcSRS, srcGT.data(), nSrcXSize, nSrcYSize,
     431             :                               dfWestLongitudeDeg, dfSouthLatitudeDeg,
     432             :                               dfEastLongitudeDeg, dfNorthLatitudeDeg);
     433             : 
     434          32 :     CPLStringList aosOptions;
     435          16 :     if (!(dfWestLongitudeDeg == 0.0 && dfSouthLatitudeDeg == 0.0 &&
     436           1 :           dfEastLongitudeDeg == 0.0 && dfNorthLatitudeDeg == 0.0))
     437             :     {
     438             :         aosOptions.SetNameValue(
     439             :             "AREA_OF_INTEREST",
     440             :             CPLSPrintf("%.16g,%.16g,%.16g,%.16g", dfWestLongitudeDeg,
     441             :                        dfSouthLatitudeDeg, dfEastLongitudeDeg,
     442          15 :                        dfNorthLatitudeDeg));
     443             :     }
     444          32 :     void *hTransform = GDALCreateGenImgProjTransformer4(
     445          16 :         hGridSRS, gridGT.data(), OGRSpatialReference::ToHandle(&oSrcSRS),
     446          16 :         srcGT.data(), aosOptions.List());
     447          16 :     if (hTransform == nullptr)
     448           3 :         return nullptr;
     449          13 :     GDALWarpOptions *psWO = GDALCreateWarpOptions();
     450          13 :     psWO->hSrcDS = hGridDataset;
     451          13 :     psWO->eResampleAlg = GRA_Bilinear;
     452          13 :     const char *pszResampling = CSLFetchNameValue(papszOptions, "RESAMPLING");
     453          13 :     if (pszResampling)
     454             :     {
     455           3 :         if (EQUAL(pszResampling, "NEAREST"))
     456           1 :             psWO->eResampleAlg = GRA_NearestNeighbour;
     457           2 :         else if (EQUAL(pszResampling, "BILINEAR"))
     458           1 :             psWO->eResampleAlg = GRA_Bilinear;
     459           1 :         else if (EQUAL(pszResampling, "CUBIC"))
     460           1 :             psWO->eResampleAlg = GRA_Cubic;
     461             :     }
     462          13 :     psWO->eWorkingDataType = GDT_Float32;
     463          13 :     int bHasNoData = FALSE;
     464          13 :     const double dfSrcNoData = GDALGetRasterNoDataValue(
     465             :         GDALGetRasterBand(hGridDataset, 1), &bHasNoData);
     466          13 :     if (bHasNoData)
     467             :     {
     468           2 :         psWO->padfSrcNoDataReal =
     469           2 :             static_cast<double *>(CPLMalloc(sizeof(double)));
     470           2 :         psWO->padfSrcNoDataReal[0] = dfSrcNoData;
     471             :     }
     472             : 
     473          13 :     psWO->padfDstNoDataReal = static_cast<double *>(CPLMalloc(sizeof(double)));
     474             :     const bool bErrorOnMissingShift =
     475          13 :         CPLFetchBool(papszOptions, "ERROR_ON_MISSING_VERT_SHIFT", false);
     476          13 :     psWO->padfDstNoDataReal[0] =
     477          13 :         (bErrorOnMissingShift) ? -std::numeric_limits<float>::infinity() : 0.0;
     478          13 :     psWO->papszWarpOptions =
     479          13 :         CSLSetNameValue(psWO->papszWarpOptions, "INIT_DEST", "NO_DATA");
     480             : 
     481          13 :     psWO->pfnTransformer = GDALGenImgProjTransform;
     482          13 :     psWO->pTransformerArg = hTransform;
     483             :     const double dfMaxError =
     484          13 :         CPLAtof(CSLFetchNameValueDef(papszOptions, "MAX_ERROR", "0.125"));
     485          13 :     if (dfMaxError > 0.0)
     486             :     {
     487          13 :         psWO->pTransformerArg = GDALCreateApproxTransformer(
     488             :             psWO->pfnTransformer, psWO->pTransformerArg, dfMaxError);
     489          13 :         psWO->pfnTransformer = GDALApproxTransform;
     490          13 :         GDALApproxTransformerOwnsSubtransformer(psWO->pTransformerArg, TRUE);
     491             :     }
     492          13 :     psWO->nBandCount = 1;
     493          13 :     psWO->panSrcBands = static_cast<int *>(CPLMalloc(sizeof(int)));
     494          13 :     psWO->panSrcBands[0] = 1;
     495          13 :     psWO->panDstBands = static_cast<int *>(CPLMalloc(sizeof(int)));
     496          13 :     psWO->panDstBands[0] = 1;
     497             : 
     498             :     VRTWarpedDataset *poReprojectedGrid =
     499          13 :         new VRTWarpedDataset(nSrcXSize, nSrcYSize);
     500             :     // This takes a reference on hGridDataset
     501          13 :     CPLErr eErr = poReprojectedGrid->Initialize(psWO);
     502          13 :     CPLAssert(eErr == CE_None);
     503          13 :     CPL_IGNORE_RET_VAL(eErr);
     504          13 :     GDALDestroyWarpOptions(psWO);
     505          13 :     poReprojectedGrid->SetGeoTransform(srcGT);
     506          13 :     poReprojectedGrid->AddBand(GDT_Float32, nullptr);
     507             : 
     508             :     GDALApplyVSGDataset *poOutDS = new GDALApplyVSGDataset(
     509          13 :         GDALDataset::FromHandle(hSrcDataset), poReprojectedGrid, eDT,
     510          13 :         CPL_TO_BOOL(bInverse), dfSrcUnitToMeter, dfDstUnitToMeter,
     511             :         // Undocumented option. For testing only
     512          13 :         atoi(CSLFetchNameValueDef(papszOptions, "BLOCKSIZE", "256")));
     513             : 
     514          13 :     poReprojectedGrid->ReleaseRef();
     515             : 
     516          13 :     if (!poOutDS->IsInitOK())
     517             :     {
     518           1 :         delete poOutDS;
     519           1 :         return nullptr;
     520             :     }
     521          12 :     poOutDS->SetDescription(GDALGetDescription(hSrcDataset));
     522          12 :     return reinterpret_cast<GDALDatasetH>(poOutDS);
     523             : }
     524             : 
     525             : /************************************************************************/
     526             : /*                           GetProj4Filename()                         */
     527             : /************************************************************************/
     528             : 
     529           0 : static CPLString GetProj4Filename(const char *pszFilename)
     530             : {
     531           0 :     CPLString osFilename;
     532             : 
     533             :     /* or fixed path: /name, ./name or ../name */
     534           0 :     if (!CPLIsFilenameRelative(pszFilename) || *pszFilename == '.')
     535             :     {
     536           0 :         return pszFilename;
     537             :     }
     538             : 
     539           0 :     PJ_GRID_INFO info = proj_grid_info(pszFilename);
     540           0 :     if (info.filename[0])
     541             :     {
     542           0 :         osFilename = info.filename;
     543             :     }
     544             : 
     545           0 :     return osFilename;
     546             : }
     547             : 
     548             : /************************************************************************/
     549             : /*                       GDALOpenVerticalShiftGrid()                    */
     550             : /************************************************************************/
     551             : 
     552             : /** Load proj.4 geoidgrids as GDAL dataset
     553             :  *
     554             :  * @param pszProj4Geoidgrids Value of proj.4 geoidgrids parameter.
     555             :  * @param pbError If not NULL, the pointed value will be set to TRUE if an
     556             :  *                error occurred.
     557             :  *
     558             :  * @return a dataset. If not NULL, it must be closed with GDALClose().
     559             :  *
     560             :  * @since GDAL 2.2
     561             :  * @deprecated GDAL 3.4. Will be removed in GDAL 4.0. This function was used
     562             :  *             by gdalwarp initially, but is no longer needed.
     563             :  */
     564           0 : GDALDatasetH GDALOpenVerticalShiftGrid(const char *pszProj4Geoidgrids,
     565             :                                        int *pbError)
     566             : {
     567           0 :     char **papszGrids = CSLTokenizeString2(pszProj4Geoidgrids, ",", 0);
     568           0 :     const int nGridCount = CSLCount(papszGrids);
     569           0 :     if (nGridCount == 1)
     570             :     {
     571           0 :         CSLDestroy(papszGrids);
     572             : 
     573           0 :         bool bMissingOk = false;
     574           0 :         if (*pszProj4Geoidgrids == '@')
     575             :         {
     576           0 :             pszProj4Geoidgrids++;
     577           0 :             bMissingOk = true;
     578             :         }
     579           0 :         const CPLString osFilename(GetProj4Filename(pszProj4Geoidgrids));
     580           0 :         const char *const papszOpenOptions[] = {
     581             :             "@SHIFT_ORIGIN_IN_MINUS_180_PLUS_180=YES", nullptr};
     582             :         GDALDatasetH hDS =
     583           0 :             GDALOpenEx(osFilename, 0, nullptr, papszOpenOptions, nullptr);
     584           0 :         if (hDS == nullptr)
     585             :         {
     586           0 :             CPLDebug("GDAL", "Cannot find file corresponding to %s",
     587             :                      pszProj4Geoidgrids);
     588             :         }
     589           0 :         if (pbError)
     590           0 :             *pbError = (!bMissingOk && hDS == nullptr);
     591           0 :         return hDS;
     592             :     }
     593             : 
     594           0 :     CPLStringList aosFilenames;
     595           0 :     for (int i = nGridCount - 1; i >= 0; i--)
     596             :     {
     597           0 :         const char *pszName = papszGrids[i];
     598           0 :         bool bMissingOk = false;
     599           0 :         if (*pszName == '@')
     600             :         {
     601           0 :             pszName++;
     602           0 :             bMissingOk = true;
     603             :         }
     604           0 :         const CPLString osFilename(GetProj4Filename(pszName));
     605             :         VSIStatBufL sStat;
     606           0 :         if (osFilename.empty() || VSIStatL(osFilename, &sStat) != 0)
     607             :         {
     608           0 :             CPLDebug("GDAL", "Cannot find file corresponding to %s", pszName);
     609           0 :             if (!bMissingOk)
     610             :             {
     611           0 :                 if (pbError)
     612           0 :                     *pbError = true;
     613           0 :                 CSLDestroy(papszGrids);
     614           0 :                 return nullptr;
     615             :             }
     616             :         }
     617             :         else
     618             :         {
     619           0 :             aosFilenames.AddString(osFilename);
     620             :         }
     621             :     }
     622             : 
     623           0 :     CSLDestroy(papszGrids);
     624             : 
     625           0 :     if (aosFilenames.empty())
     626             :     {
     627           0 :         if (pbError)
     628           0 :             *pbError = false;
     629           0 :         return nullptr;
     630             :     }
     631             : 
     632           0 :     char **papszArgv = nullptr;
     633           0 :     papszArgv = CSLAddString(papszArgv, "-resolution");
     634           0 :     papszArgv = CSLAddString(papszArgv, "highest");
     635           0 :     papszArgv = CSLAddString(papszArgv, "-vrtnodata");
     636           0 :     papszArgv = CSLAddString(papszArgv, "-inf");
     637           0 :     papszArgv = CSLAddString(papszArgv, "-oo");
     638             :     papszArgv =
     639           0 :         CSLAddString(papszArgv, "@SHIFT_ORIGIN_IN_MINUS_180_PLUS_180=YES");
     640           0 :     GDALBuildVRTOptions *psOptions = GDALBuildVRTOptionsNew(papszArgv, nullptr);
     641           0 :     CSLDestroy(papszArgv);
     642           0 :     GDALDatasetH hDS = GDALBuildVRT("", aosFilenames.size(), nullptr,
     643           0 :                                     aosFilenames.List(), psOptions, nullptr);
     644           0 :     GDALBuildVRTOptionsFree(psOptions);
     645           0 :     if (pbError)
     646           0 :         *pbError = hDS != nullptr;
     647           0 :     return hDS;
     648             : }

Generated by: LCOV version 1.14