/******************************************************************************

 *

 * Project:  High Performance Image Reprojector

 * Purpose:  Test program for high performance warper API.

 * Author:   Frank Warmerdam <warmerdam@pobox.com>

 *

 ******************************************************************************

 * Copyright (c) 2002, i3 - information integration and imaging

 *                          Fort Collins, CO

 * Copyright (c) 2007-2015, Even Rouault <even dot rouault at spatialys.com>

 * Copyright (c) 2015, Faza Mahamood

 *

 * SPDX-License-Identifier: MIT

 ****************************************************************************/

#include "cpl_port.h"

#include "gdal_utils.h"

#include "gdal_utils_priv.h"

#include "gdalargumentparser.h"

#include <cctype>

#include <cmath>

#include <cstdio>

#include <cstdlib>

#include <cstring>

#include <algorithm>

#include <array>

#include <limits>

#include <set>

#include <utility>

#include <vector>

// Suppress deprecation warning for GDALOpenVerticalShiftGrid and

// GDALApplyVerticalShiftGrid

#ifndef CPL_WARN_DEPRECATED_GDALOpenVerticalShiftGrid

#define CPL_WARN_DEPRECATED_GDALOpenVerticalShiftGrid(x)

#define CPL_WARN_DEPRECATED_GDALApplyVerticalShiftGrid(x)

#endif

#include "commonutils.h"

#include "cpl_conv.h"

#include "cpl_error.h"

#include "cpl_progress.h"

#include "cpl_string.h"

#include "gdal.h"

#include "gdal_alg.h"

#include "gdal_alg_priv.h"

#include "gdal_priv.h"

#include "gdalwarper.h"

#include "ogr_api.h"

#include "ogr_core.h"

#include "ogr_geometry.h"

#include "ogr_spatialref.h"

#include "ogr_srs_api.h"

#include "ogr_proj_p.h"

#include "ogrct_priv.h"

#include "ogrsf_frmts.h"

#include "vrtdataset.h"

#include "../frmts/gtiff/cogdriver.h"

#if PROJ_VERSION_MAJOR > 6 || PROJ_VERSION_MINOR >= 3

#define USE_PROJ_BASED_VERTICAL_SHIFT_METHOD

#endif

/************************************************************************/

/*                          GDALWarpAppOptions                          */

/************************************************************************/

/** Options for use with GDALWarp(). GDALWarpAppOptions* must be allocated and

 * freed with GDALWarpAppOptionsNew() and GDALWarpAppOptionsFree() respectively.

 */

struct GDALWarpAppOptions

{

    /*! Raw program arguments */

    CPLStringList aosArgs{};

    /*! set georeferenced extents of output file to be created (in target SRS by

       default, or in the SRS specified with pszTE_SRS) */

    double dfMinX = 0;

    double dfMinY = 0;

    double dfMaxX = 0;

    double dfMaxY = 0;

    /*! the SRS in which to interpret the coordinates given in

       GDALWarpAppOptions::dfMinX, GDALWarpAppOptions::dfMinY,

       GDALWarpAppOptions::dfMaxX and GDALWarpAppOptions::dfMaxY. The SRS may be

       any of the usual GDAL/OGR forms, complete WKT, PROJ.4, EPSG:n or a file

       containing the WKT. It is a convenience e.g. when knowing the output

       coordinates in a geodetic long/lat SRS, but still wanting a result in a

       projected coordinate system. */

    std::string osTE_SRS{};

    /*! set output file resolution (in target georeferenced units) */

    double dfXRes = 0;

    double dfYRes = 0;

    /*! whether target pixels should have dfXRes == dfYRes */

    bool bSquarePixels = false;

    /*! align the coordinates of the extent of the output file to the values of

       the GDALWarpAppOptions::dfXRes and GDALWarpAppOptions::dfYRes, such that

       the aligned extent includes the minimum extent. */

    bool bTargetAlignedPixels = false;

    /*! set output file size in pixels and lines. If

       GDALWarpAppOptions::nForcePixels or GDALWarpAppOptions::nForceLines is

       set to 0, the other dimension will be guessed from the computed

       resolution. Note that GDALWarpAppOptions::nForcePixels and

        GDALWarpAppOptions::nForceLines cannot be used with

       GDALWarpAppOptions::dfXRes and GDALWarpAppOptions::dfYRes. */

    int nForcePixels = 0;

    int nForceLines = 0;

    /*! allow or suppress progress monitor and other non-error output */

    bool bQuiet = true;

    /*! the progress function to use */

    GDALProgressFunc pfnProgress = GDALDummyProgress;

    /*! pointer to the progress data variable */

    void *pProgressData = nullptr;

    /*! creates an output alpha band to identify nodata (unset/transparent)

       pixels when set to true */

    bool bEnableDstAlpha = false;

    /*! forces the last band of an input file to be considered as alpha band. */

    bool bEnableSrcAlpha = false;

    /*! Prevent a source alpha band from being considered as such */

    bool bDisableSrcAlpha = false;

    /*! output format. Use the short format name. */

    std::string osFormat{};

    bool bCreateOutput = false;

    /*! list of warp options. ("NAME1=VALUE1","NAME2=VALUE2",...). The

        GDALWarpOptions::aosWarpOptions docs show all options. */

    CPLStringList aosWarpOptions{};

    double dfErrorThreshold = -1;

    /*! the amount of memory (in megabytes) that the warp API is allowed

        to use for caching. */

    double dfWarpMemoryLimit = 0;

    /*! list of create options for the output format driver. See format

        specific documentation for legal creation options for each format. */

    CPLStringList aosCreateOptions{};

    /*! the data type of the output bands */

    GDALDataType eOutputType = GDT_Unknown;

    /*! working pixel data type. The data type of pixels in the source

        image and destination image buffers. */

    GDALDataType eWorkingType = GDT_Unknown;

    /*! the resampling method. Available methods are: near, bilinear,

        cubic, cubicspline, lanczos, average, mode, max, min, med,

        q1, q3, sum */

    GDALResampleAlg eResampleAlg = GRA_NearestNeighbour;

    /*! whether -r was specified */

    bool bResampleAlgSpecifiedByUser = false;

    /*! nodata masking values for input bands (different values can be supplied

        for each band). ("value1 value2 ..."). Masked values will not be used

        in interpolation. Use a value of "None" to ignore intrinsic nodata

        settings on the source dataset. */

    std::string osSrcNodata{};

    /*! nodata values for output bands (different values can be supplied for

        each band). ("value1 value2 ..."). New files will be initialized to

        this value and if possible the nodata value will be recorded in the

        output file. Use a value of "None" to ensure that nodata is not defined.

        If this argument is not used then nodata values will be copied from

        the source dataset. */

    std::string osDstNodata{};

    /*! use multithreaded warping implementation. Multiple threads will be used

        to process chunks of image and perform input/output operation

       simultaneously. */

    bool bMulti = false;

    /*! list of transformer options suitable to pass to

       GDALCreateGenImgProjTransformer2().

        ("NAME1=VALUE1","NAME2=VALUE2",...) */

    CPLStringList aosTransformerOptions{};

    /*! enable use of a blend cutline from a vector dataset name or a WKT

     * geometry

     */

    std::string osCutlineDSNameOrWKT{};

    /*! cutline SRS */

    std::string osCutlineSRS{};

    /*! the named layer to be selected from the cutline datasource */

    std::string osCLayer{};

    /*! restrict desired cutline features based on attribute query */

    std::string osCWHERE{};

    /*! SQL query to select the cutline features instead of from a layer

        with osCLayer */

    std::string osCSQL{};

    /*! crop the extent of the target dataset to the extent of the cutline */

    bool bCropToCutline = false;

    /*! copy dataset and band metadata will be copied from the first source

       dataset. Items that differ between source datasets will be set "*" (see

       GDALWarpAppOptions::pszMDConflictValue) */

    bool bCopyMetadata = true;

    /*! copy band information from the first source dataset */

    bool bCopyBandInfo = true;

    /*! value to set metadata items that conflict between source datasets

       (default is "*"). Use "" to remove conflicting items. */

    std::string osMDConflictValue = "*";

    /*! set the color interpretation of the bands of the target dataset from the

     * source dataset */

    bool bSetColorInterpretation = false;

    /*! overview level of source files to be used */

    int nOvLevel = OVR_LEVEL_AUTO;

    /*! Whether to enable vertical shift adjustment */

    bool bVShift = false;

    /*! Whether to disable vertical shift adjustment */

    bool bNoVShift = false;

    /*! Source bands */

    std::vector<int> anSrcBands{};

    /*! Destination bands */

    std::vector<int> anDstBands{};

    /*! Used when using a temporary TIFF file while warping */

    bool bDeleteOutputFileOnceCreated = false;

};

static CPLErr

LoadCutline(const std::string &osCutlineDSNameOrWKT, const std::string &osSRS,

            const std::string &oszCLayer, const std::string &osCWHERE,

            const std::string &osCSQL, OGRGeometryH *phCutlineRet);

static CPLErr TransformCutlineToSource(GDALDataset *poSrcDS,

                                       OGRGeometry *poCutline,

                                       char ***ppapszWarpOptions,

                                       CSLConstList papszTO);

static GDALDatasetH GDALWarpCreateOutput(

    int nSrcCount, GDALDatasetH *pahSrcDS, const char *pszFilename,

    const char *pszFormat, char **papszTO, CSLConstList papszCreateOptions,

    GDALDataType eDT, GDALTransformerArgUniquePtr &hTransformArg,

    bool bSetColorInterpretation, GDALWarpAppOptions *psOptions,

    bool bUpdateTransformerWithDestGT);

static void RemoveConflictingMetadata(GDALMajorObjectH hObj,

                                      CSLConstList papszMetadata,

                                      const char *pszValueConflict);

static double GetAverageSegmentLength(const OGRGeometry *poGeom)

{

    if (!poGeom)

        return 0;

    switch (wkbFlatten(poGeom->getGeometryType()))

    {

        case wkbLineString:

        {

            const auto *poLS = poGeom->toLineString();

            double dfSum = 0;

            const int nPoints = poLS->getNumPoints();

            if (nPoints == 0)

                return 0;

            for (int i = 0; i < nPoints - 1; i++)

            {

                double dfX1 = poLS->getX(i);

                double dfY1 = poLS->getY(i);

                double dfX2 = poLS->getX(i + 1);

                double dfY2 = poLS->getY(i + 1);

                double dfDX = dfX2 - dfX1;

                double dfDY = dfY2 - dfY1;

                dfSum += sqrt(dfDX * dfDX + dfDY * dfDY);

            }

            return dfSum / nPoints;

        }

        case wkbPolygon:

        {

            if (poGeom->IsEmpty())

                return 0;

            double dfSum = 0;

            for (const auto *poLS : poGeom->toPolygon())

            {

                dfSum += GetAverageSegmentLength(poLS);

            }

            return dfSum / (1 + poGeom->toPolygon()->getNumInteriorRings());

        }

        case wkbMultiPolygon:

        case wkbMultiLineString:

        case wkbGeometryCollection:

        {

            if (poGeom->IsEmpty())

                return 0;

            double dfSum = 0;

            for (const auto *poSubGeom : poGeom->toGeometryCollection())

            {

                dfSum += GetAverageSegmentLength(poSubGeom);

            }

            return dfSum / poGeom->toGeometryCollection()->getNumGeometries();

        }

        default:

            return 0;

    }

}

/************************************************************************/

/*                           FetchSrcMethod()                           */

/************************************************************************/

static const char *FetchSrcMethod(CSLConstList papszTO,

                                  const char *pszDefault = nullptr)

{

    const char *pszMethod = CSLFetchNameValue(papszTO, "SRC_METHOD");

    if (!pszMethod)

        pszMethod = CSLFetchNameValueDef(papszTO, "METHOD", pszDefault);

    return pszMethod;

}

static const char *FetchSrcMethod(const CPLStringList &aosTO,

                                  const char *pszDefault = nullptr)

{

    const char *pszMethod = aosTO.FetchNameValue("SRC_METHOD");

    if (!pszMethod)

        pszMethod = aosTO.FetchNameValueDef("METHOD", pszDefault);

    return pszMethod;

}

/************************************************************************/

/*                          GetSrcDSProjection()                        */

/*                                                                      */

/* Takes into account SRC_SRS transformer option in priority, and then  */

/* dataset characteristics as well as the METHOD transformer            */

/* option to determine the source SRS.                                  */

/************************************************************************/

static CPLString GetSrcDSProjection(GDALDatasetH hDS, CSLConstList papszTO)

{

    const char *pszProjection = CSLFetchNameValue(papszTO, "SRC_SRS");

    if (pszProjection != nullptr || hDS == nullptr)

    {

        return pszProjection ? pszProjection : "";

    }

    const char *pszMethod = FetchSrcMethod(papszTO);

    CSLConstList papszMD = nullptr;

    const OGRSpatialReferenceH hSRS = GDALGetSpatialRef(hDS);

    const char *pszGeolocationDataset =

        CSLFetchNameValueDef(papszTO, "SRC_GEOLOC_ARRAY",

                             CSLFetchNameValue(papszTO, "GEOLOC_ARRAY"));

    if (pszGeolocationDataset != nullptr &&

        (pszMethod == nullptr || EQUAL(pszMethod, "GEOLOC_ARRAY")))

    {

        auto aosMD =

            GDALCreateGeolocationMetadata(hDS, pszGeolocationDataset, true);

        pszProjection = aosMD.FetchNameValue("SRS");

        if (pszProjection)

            return pszProjection;  // return in this scope so that aosMD is

                                   // still valid

    }

    else if (hSRS && (pszMethod == nullptr || EQUAL(pszMethod, "GEOTRANSFORM")))

    {

        char *pszWKT = nullptr;

        {

            CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);

            if (OSRExportToWkt(hSRS, &pszWKT) != OGRERR_NONE)

            {

                CPLFree(pszWKT);

                pszWKT = nullptr;

                const char *const apszOptions[] = {"FORMAT=WKT2", nullptr};

                OSRExportToWktEx(hSRS, &pszWKT, apszOptions);

            }

        }

        CPLString osWKT = pszWKT ? pszWKT : "";

        CPLFree(pszWKT);

        return osWKT;

    }

    else if (GDALGetGCPProjection(hDS) != nullptr &&

             strlen(GDALGetGCPProjection(hDS)) > 0 &&

             GDALGetGCPCount(hDS) > 1 &&

             (pszMethod == nullptr || STARTS_WITH_CI(pszMethod, "GCP_")))

    {

        pszProjection = GDALGetGCPProjection(hDS);

    }

    else if (GDALGetMetadata(hDS, "RPC") != nullptr &&

             (pszMethod == nullptr || EQUAL(pszMethod, "RPC")))

    {

        pszProjection = SRS_WKT_WGS84_LAT_LONG;

    }

    else if ((papszMD = GDALGetMetadata(hDS, "GEOLOCATION")) != nullptr &&

             (pszMethod == nullptr || EQUAL(pszMethod, "GEOLOC_ARRAY")))

    {

        pszProjection = CSLFetchNameValue(papszMD, "SRS");

    }

    return pszProjection ? pszProjection : "";

}

/************************************************************************/

/*                        CreateCTCutlineToSrc()                        */

/************************************************************************/

static std::unique_ptr<OGRCoordinateTransformation> CreateCTCutlineToSrc(

    const OGRSpatialReference *poRasterSRS, const OGRSpatialReference *poDstSRS,

    const OGRSpatialReference *poCutlineSRS, CSLConstList papszTO)

{

    const OGRSpatialReference *poCutlineOrTargetSRS =

        poCutlineSRS ? poCutlineSRS : poDstSRS;

    std::unique_ptr<OGRCoordinateTransformation> poCTCutlineToSrc;

    if (poCutlineOrTargetSRS && poRasterSRS &&

        !poCutlineOrTargetSRS->IsSame(poRasterSRS))

    {

        OGRCoordinateTransformationOptions oOptions;

        // If the cutline SRS is the same as the target SRS and there is

        // an explicit -ct between the source SRS and the target SRS, then

        // use it in the reverse way to transform from the cutline SRS to

        // the source SRS.

        if (poDstSRS && poCutlineOrTargetSRS->IsSame(poDstSRS))

        {

            const char *pszCT =

                CSLFetchNameValue(papszTO, "COORDINATE_OPERATION");

            if (pszCT)

            {

                oOptions.SetCoordinateOperation(pszCT, /* bInverse = */ true);

            }

        }

        poCTCutlineToSrc.reset(OGRCreateCoordinateTransformation(

            poCutlineOrTargetSRS, poRasterSRS, oOptions));

    }

    return poCTCutlineToSrc;

}

/************************************************************************/

/*                           CropToCutline()                            */

/************************************************************************/

static CPLErr CropToCutline(const OGRGeometry *poCutline, CSLConstList papszTO,

                            CSLConstList papszWarpOptions, int nSrcCount,

                            GDALDatasetH *pahSrcDS, double &dfMinX,

                            double &dfMinY, double &dfMaxX, double &dfMaxY,

                            const GDALWarpAppOptions *psOptions)

{

    // We could possibly directly reproject from cutline SRS to target SRS,

    // but when applying the cutline, it is reprojected to source raster image

    // space using the source SRS. To be consistent, we reproject

    // the cutline from cutline SRS to source SRS and then from source SRS to

    // target SRS.

    const OGRSpatialReference *poCutlineSRS = poCutline->getSpatialReference();

    const char *pszThisTargetSRS = CSLFetchNameValue(papszTO, "DST_SRS");

    std::unique_ptr<OGRSpatialReference> poSrcSRS;

    std::unique_ptr<OGRSpatialReference> poDstSRS;

    const CPLString osThisSourceSRS =

        GetSrcDSProjection(nSrcCount > 0 ? pahSrcDS[0] : nullptr, papszTO);

    if (!osThisSourceSRS.empty())

    {

        poSrcSRS = std::make_unique<OGRSpatialReference>();

        poSrcSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);

        if (poSrcSRS->SetFromUserInput(osThisSourceSRS) != OGRERR_NONE)

        {

            CPLError(CE_Failure, CPLE_AppDefined,

                     "Cannot compute bounding box of cutline.");

            return CE_Failure;

        }

    }

    else if (!pszThisTargetSRS && !poCutlineSRS)

    {

        OGREnvelope sEnvelope;

        poCutline->getEnvelope(&sEnvelope);

        dfMinX = sEnvelope.MinX;

        dfMinY = sEnvelope.MinY;

        dfMaxX = sEnvelope.MaxX;

        dfMaxY = sEnvelope.MaxY;

        return CE_None;

    }

    else

    {

        CPLError(CE_Failure, CPLE_AppDefined,

                 "Cannot compute bounding box of cutline. Cannot find "

                 "source SRS");

        return CE_Failure;

    }

    if (pszThisTargetSRS)

    {

        poDstSRS = std::make_unique<OGRSpatialReference>();

        poDstSRS->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);

        if (poDstSRS->SetFromUserInput(pszThisTargetSRS) != OGRERR_NONE)

        {

            CPLError(CE_Failure, CPLE_AppDefined,

                     "Cannot compute bounding box of cutline.");

            return CE_Failure;

        }

    }

    else

    {

        poDstSRS.reset(poSrcSRS->Clone());

    }

    auto poCutlineGeom = std::unique_ptr<OGRGeometry>(poCutline->clone());

    auto poCTCutlineToSrc = CreateCTCutlineToSrc(poSrcSRS.get(), poDstSRS.get(),

                                                 poCutlineSRS, papszTO);

    std::unique_ptr<OGRCoordinateTransformation> poCTSrcToDst;

    if (!poSrcSRS->IsSame(poDstSRS.get()))

    {

        poCTSrcToDst.reset(

            OGRCreateCoordinateTransformation(poSrcSRS.get(), poDstSRS.get()));

    }

    // Reproject cutline to target SRS, by doing intermediate vertex

    // densification in source SRS.

    if (poCTSrcToDst || poCTCutlineToSrc)

    {

        OGREnvelope sLastEnvelope, sCurEnvelope;

        std::unique_ptr<OGRGeometry> poTransformedGeom;

        auto poGeomInSrcSRS =

            std::unique_ptr<OGRGeometry>(poCutlineGeom->clone());

        if (poCTCutlineToSrc)

        {

            poGeomInSrcSRS.reset(OGRGeometryFactory::transformWithOptions(

                poGeomInSrcSRS.get(), poCTCutlineToSrc.get(), nullptr));

            if (!poGeomInSrcSRS)

                return CE_Failure;

        }

        // Do not use a smaller epsilon, otherwise it could cause useless

        // segmentization (https://github.com/OSGeo/gdal/issues/4826)

        constexpr double epsilon = 1e-10;

        for (int nIter = 0; nIter < 10; nIter++)

        {

            poTransformedGeom.reset(poGeomInSrcSRS->clone());

            if (poCTSrcToDst)

            {

                poTransformedGeom.reset(

                    OGRGeometryFactory::transformWithOptions(

                        poTransformedGeom.get(), poCTSrcToDst.get(), nullptr));

                if (!poTransformedGeom)

                    return CE_Failure;

            }

            poTransformedGeom->getEnvelope(&sCurEnvelope);

            if (nIter > 0 || !poCTSrcToDst)

            {

                if (std::abs(sCurEnvelope.MinX - sLastEnvelope.MinX) <=

                        epsilon *

                            std::abs(sCurEnvelope.MinX + sLastEnvelope.MinX) &&

                    std::abs(sCurEnvelope.MinY - sLastEnvelope.MinY) <=

                        epsilon *

                            std::abs(sCurEnvelope.MinY + sLastEnvelope.MinY) &&

                    std::abs(sCurEnvelope.MaxX - sLastEnvelope.MaxX) <=

                        epsilon *

                            std::abs(sCurEnvelope.MaxX + sLastEnvelope.MaxX) &&

                    std::abs(sCurEnvelope.MaxY - sLastEnvelope.MaxY) <=

                        epsilon *

                            std::abs(sCurEnvelope.MaxY + sLastEnvelope.MaxY))

                {

                    break;

                }

            }

            double dfAverageSegmentLength =

                GetAverageSegmentLength(poGeomInSrcSRS.get());

            poGeomInSrcSRS->segmentize(dfAverageSegmentLength / 4);

            sLastEnvelope = sCurEnvelope;

        }

        poCutlineGeom = std::move(poTransformedGeom);

    }

    OGREnvelope sEnvelope;

    poCutlineGeom->getEnvelope(&sEnvelope);

    dfMinX = sEnvelope.MinX;

    dfMinY = sEnvelope.MinY;

    dfMaxX = sEnvelope.MaxX;

    dfMaxY = sEnvelope.MaxY;

    if (!poCTSrcToDst && nSrcCount > 0 && psOptions->dfXRes == 0.0 &&

        psOptions->dfYRes == 0.0)

    {

        // No raster reprojection: stick on exact pixel boundaries of the source

        // to preserve resolution and avoid resampling

        double adfGT[6];

        if (GDALGetGeoTransform(pahSrcDS[0], adfGT) == CE_None)

        {

            // We allow for a relative error in coordinates up to 0.1% of the

            // pixel size for rounding purposes.

            constexpr double REL_EPS_PIXEL = 1e-3;

            if (CPLFetchBool(papszWarpOptions, "CUTLINE_ALL_TOUCHED", false))

            {

                // All touched ? Then make the extent a bit larger than the

                // cutline envelope

                dfMinX = adfGT[0] +

                         floor((dfMinX - adfGT[0]) / adfGT[1] + REL_EPS_PIXEL) *

                             adfGT[1];

                dfMinY = adfGT[3] +

                         ceil((dfMinY - adfGT[3]) / adfGT[5] - REL_EPS_PIXEL) *

                             adfGT[5];

                dfMaxX = adfGT[0] +

                         ceil((dfMaxX - adfGT[0]) / adfGT[1] - REL_EPS_PIXEL) *

                             adfGT[1];

                dfMaxY = adfGT[3] +

                         floor((dfMaxY - adfGT[3]) / adfGT[5] + REL_EPS_PIXEL) *

                             adfGT[5];

            }

            else

            {

                // Otherwise, make it a bit smaller

                dfMinX = adfGT[0] +

                         ceil((dfMinX - adfGT[0]) / adfGT[1] - REL_EPS_PIXEL) *

                             adfGT[1];

                dfMinY = adfGT[3] +

                         floor((dfMinY - adfGT[3]) / adfGT[5] + REL_EPS_PIXEL) *

                             adfGT[5];

                dfMaxX = adfGT[0] +

                         floor((dfMaxX - adfGT[0]) / adfGT[1] + REL_EPS_PIXEL) *

                             adfGT[1];

                dfMaxY = adfGT[3] +

                         ceil((dfMaxY - adfGT[3]) / adfGT[5] - REL_EPS_PIXEL) *

                             adfGT[5];

            }

        }

    }

    return CE_None;

}

#ifdef USE_PROJ_BASED_VERTICAL_SHIFT_METHOD

static bool MustApplyVerticalShift(GDALDatasetH hWrkSrcDS,

                                   const GDALWarpAppOptions *psOptions,

                                   OGRSpatialReference &oSRSSrc,

                                   OGRSpatialReference &oSRSDst,

                                   bool &bSrcHasVertAxis, bool &bDstHasVertAxis)

{

    bool bApplyVShift = psOptions->bVShift;

    // Check if we must do vertical shift grid transform

    const char *pszSrcWKT =

        psOptions->aosTransformerOptions.FetchNameValue("SRC_SRS");

    if (pszSrcWKT)

        oSRSSrc.SetFromUserInput(pszSrcWKT);

    else

    {

        auto hSRS = GDALGetSpatialRef(hWrkSrcDS);

        if (hSRS)

            oSRSSrc = *(OGRSpatialReference::FromHandle(hSRS));

        else

            return false;

    }

    const char *pszDstWKT =

        psOptions->aosTransformerOptions.FetchNameValue("DST_SRS");

    if (pszDstWKT)

        oSRSDst.SetFromUserInput(pszDstWKT);

    else

        return false;

    if (oSRSSrc.IsSame(&oSRSDst))

        return false;

    bSrcHasVertAxis = oSRSSrc.IsCompound() ||

                      ((oSRSSrc.IsProjected() || oSRSSrc.IsGeographic()) &&

                       oSRSSrc.GetAxesCount() == 3);

    bDstHasVertAxis = oSRSDst.IsCompound() ||

                      ((oSRSDst.IsProjected() || oSRSDst.IsGeographic()) &&

                       oSRSDst.GetAxesCount() == 3);

    if ((GDALGetRasterCount(hWrkSrcDS) == 1 || psOptions->bVShift) &&

        (bSrcHasVertAxis || bDstHasVertAxis))

    {

        bApplyVShift = true;

    }

    return bApplyVShift;

}

/************************************************************************/

/*                         ApplyVerticalShift()                         */

/************************************************************************/

static bool ApplyVerticalShift(GDALDatasetH hWrkSrcDS,

                               const GDALWarpAppOptions *psOptions,

                               GDALWarpOptions *psWO)

{

    if (psOptions->bVShift)

    {

        psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions,

                                                 "APPLY_VERTICAL_SHIFT", "YES");

    }

    OGRSpatialReference oSRSSrc;

    OGRSpatialReference oSRSDst;

    bool bSrcHasVertAxis = false;

    bool bDstHasVertAxis = false;

    bool bApplyVShift =

        MustApplyVerticalShift(hWrkSrcDS, psOptions, oSRSSrc, oSRSDst,

                               bSrcHasVertAxis, bDstHasVertAxis);

    if ((GDALGetRasterCount(hWrkSrcDS) == 1 || psOptions->bVShift) &&

        (bSrcHasVertAxis || bDstHasVertAxis))

    {

        bApplyVShift = true;

        psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions,

                                                 "APPLY_VERTICAL_SHIFT", "YES");

        if (CSLFetchNameValue(psWO->papszWarpOptions,

                              "MULT_FACTOR_VERTICAL_SHIFT") == nullptr)

        {

            // Select how to go from input dataset units to meters

            double dfToMeterSrc = 1.0;

            const char *pszUnit =

                GDALGetRasterUnitType(GDALGetRasterBand(hWrkSrcDS, 1));

            double dfToMeterSrcAxis = 1.0;

            if (bSrcHasVertAxis)

            {

                oSRSSrc.GetAxis(nullptr, 2, nullptr, &dfToMeterSrcAxis);

            }

            if (pszUnit && (EQUAL(pszUnit, "m") || EQUAL(pszUnit, "meter") ||

                            EQUAL(pszUnit, "metre")))

            {

            }

            else if (pszUnit &&

                     (EQUAL(pszUnit, "ft") || EQUAL(pszUnit, "foot")))

            {

                dfToMeterSrc = CPLAtof(SRS_UL_FOOT_CONV);

            }

            else if (pszUnit && (EQUAL(pszUnit, "US survey foot")))

            {

                dfToMeterSrc = CPLAtof(SRS_UL_US_FOOT_CONV);

            }

            else if (pszUnit && !EQUAL(pszUnit, ""))

            {

                if (bSrcHasVertAxis)

                {

                    dfToMeterSrc = dfToMeterSrcAxis;

                }

                else

                {

                    CPLError(CE_Warning, CPLE_AppDefined,

                             "Unknown units=%s. Assuming metre.", pszUnit);

                }

            }

            else

            {

                if (bSrcHasVertAxis)

                    oSRSSrc.GetAxis(nullptr, 2, nullptr, &dfToMeterSrc);

            }

            double dfToMeterDst = 1.0;

            if (bDstHasVertAxis)

                oSRSDst.GetAxis(nullptr, 2, nullptr, &dfToMeterDst);

            if (dfToMeterSrc > 0 && dfToMeterDst > 0)

            {

                const double dfMultFactorVerticalShift =

                    dfToMeterSrc / dfToMeterDst;

                CPLDebug("WARP", "Applying MULT_FACTOR_VERTICAL_SHIFT=%.18g",

                         dfMultFactorVerticalShift);

                psWO->papszWarpOptions = CSLSetNameValue(

                    psWO->papszWarpOptions, "MULT_FACTOR_VERTICAL_SHIFT",

                    CPLSPrintf("%.18g", dfMultFactorVerticalShift));

                const double dfMultFactorVerticalShiftPipeline =

                    dfToMeterSrcAxis / dfToMeterDst;

                CPLDebug("WARP",

                         "Applying MULT_FACTOR_VERTICAL_SHIFT_PIPELINE=%.18g",

                         dfMultFactorVerticalShiftPipeline);

                psWO->papszWarpOptions = CSLSetNameValue(

                    psWO->papszWarpOptions,

                    "MULT_FACTOR_VERTICAL_SHIFT_PIPELINE",

                    CPLSPrintf("%.18g", dfMultFactorVerticalShiftPipeline));

            }

        }

    }

    return bApplyVShift;

}

#else

/************************************************************************/

/*                       ApplyVerticalShiftGrid()                       */

/************************************************************************/

static GDALDatasetH ApplyVerticalShiftGrid(GDALDatasetH hWrkSrcDS,

                                           const GDALWarpAppOptions *psOptions,

                                           GDALDatasetH hVRTDS,

                                           bool &bErrorOccurredOut)

{

    bErrorOccurredOut = false;

    // Check if we must do vertical shift grid transform

    OGRSpatialReference oSRSSrc;

    OGRSpatialReference oSRSDst;

    const char *pszSrcWKT =

        psOptions->aosTransformerOptions.FetchNameValue("SRC_SRS");

    if (pszSrcWKT)

        oSRSSrc.SetFromUserInput(pszSrcWKT);

    else

    {

        auto hSRS = GDALGetSpatialRef(hWrkSrcDS);

        if (hSRS)

            oSRSSrc = *(OGRSpatialReference::FromHandle(hSRS));

    }

    const char *pszDstWKT =

        psOptions->aosTransformerOptions.FetchNameValue("DST_SRS");

    if (pszDstWKT)

        oSRSDst.SetFromUserInput(pszDstWKT);

    double adfGT[6] = {};

    if (GDALGetRasterCount(hWrkSrcDS) == 1 &&

        GDALGetGeoTransform(hWrkSrcDS, adfGT) == CE_None &&

        !oSRSSrc.IsEmpty() && !oSRSDst.IsEmpty())

    {

        if ((oSRSSrc.IsCompound() ||

             (oSRSSrc.IsGeographic() && oSRSSrc.GetAxesCount() == 3)) ||

            (oSRSDst.IsCompound() ||

             (oSRSDst.IsGeographic() && oSRSDst.GetAxesCount() == 3)))

        {

            const char *pszSrcProj4Geoids =

                oSRSSrc.GetExtension("VERT_DATUM", "PROJ4_GRIDS");

            const char *pszDstProj4Geoids =

                oSRSDst.GetExtension("VERT_DATUM", "PROJ4_GRIDS");

            if (oSRSSrc.IsCompound() && pszSrcProj4Geoids == nullptr)

            {

                CPLDebug("GDALWARP", "Source SRS is a compound CRS but lacks "

                                     "+geoidgrids");

            }

            if (oSRSDst.IsCompound() && pszDstProj4Geoids == nullptr)

            {

                CPLDebug("GDALWARP", "Target SRS is a compound CRS but lacks "

                                     "+geoidgrids");

            }

            if (pszSrcProj4Geoids != nullptr && pszDstProj4Geoids != nullptr &&

                EQUAL(pszSrcProj4Geoids, pszDstProj4Geoids))

            {

                pszSrcProj4Geoids = nullptr;

                pszDstProj4Geoids = nullptr;

            }

            // Select how to go from input dataset units to meters

            const char *pszUnit =

                GDALGetRasterUnitType(GDALGetRasterBand(hWrkSrcDS, 1));

            double dfToMeterSrc = 1.0;

            if (pszUnit && (EQUAL(pszUnit, "m") || EQUAL(pszUnit, "meter") ||

                            EQUAL(pszUnit, "metre")))

            {

            }

            else if (pszUnit &&

                     (EQUAL(pszUnit, "ft") || EQUAL(pszUnit, "foot")))

            {

                dfToMeterSrc = CPLAtof(SRS_UL_FOOT_CONV);

            }

            else if (pszUnit && (EQUAL(pszUnit, "US survey foot")))

            {

                dfToMeterSrc = CPLAtof(SRS_UL_US_FOOT_CONV);

            }

            else

            {

                if (pszUnit && !EQUAL(pszUnit, ""))

                {

                    CPLError(CE_Warning, CPLE_AppDefined, "Unknown units=%s",

                             pszUnit);

                }

                if (oSRSSrc.IsCompound())

                {

                    dfToMeterSrc = oSRSSrc.GetTargetLinearUnits("VERT_CS");

                }

                else if (oSRSSrc.IsProjected())

                {

                    dfToMeterSrc = oSRSSrc.GetLinearUnits();

                }

            }

            double dfToMeterDst = 1.0;

            if (oSRSDst.IsCompound())

            {

                dfToMeterDst = oSRSDst.GetTargetLinearUnits("VERT_CS");

            }

            else if (oSRSDst.IsProjected())

            {

                dfToMeterDst = oSRSDst.GetLinearUnits();

            }

            char **papszOptions = nullptr;

            if (psOptions->eOutputType != GDT_Unknown)

            {

                papszOptions = CSLSetNameValue(

                    papszOptions, "DATATYPE",

                    GDALGetDataTypeName(psOptions->eOutputType));

            }

            papszOptions =

                CSLSetNameValue(papszOptions, "ERROR_ON_MISSING_VERT_SHIFT",

                                psOptions->aosTransformerOptions.FetchNameValue(

                                    "ERROR_ON_MISSING_VERT_SHIFT"));

            papszOptions = CSLSetNameValue(papszOptions, "SRC_SRS", pszSrcWKT);

            if (pszSrcProj4Geoids != nullptr)

            {

                int bError = FALSE;

                GDALDatasetH hGridDataset =

                    GDALOpenVerticalShiftGrid(pszSrcProj4Geoids, &bError);

                if (bError && hGridDataset == nullptr)

                {

                    CPLError(CE_Failure, CPLE_AppDefined, "Cannot open %s.",

                             pszSrcProj4Geoids);

                    bErrorOccurredOut = true;

                    CSLDestroy(papszOptions);

                    return hWrkSrcDS;

                }

                else if (hGridDataset != nullptr)

                {

                    // Transform from source vertical datum to WGS84

                    GDALDatasetH hTmpDS = GDALApplyVerticalShiftGrid(

                        hWrkSrcDS, hGridDataset, FALSE, dfToMeterSrc, 1.0,

                        papszOptions);

                    GDALReleaseDataset(hGridDataset);

                    if (hTmpDS == nullptr)

                    {

                        bErrorOccurredOut = true;

                        CSLDestroy(papszOptions);

                        return hWrkSrcDS;

                    }

                    else

                    {

                        if (hVRTDS)

                        {

                            CPLError(

                                CE_Failure, CPLE_NotSupported,

                                "Warping to VRT with vertical transformation "

                                "not supported with PROJ < 6.3");

                            bErrorOccurredOut = true;

                            CSLDestroy(papszOptions);

                            return hWrkSrcDS;

                        }

                        CPLDebug("GDALWARP",

                                 "Adjusting source dataset "

                                 "with source vertical datum using %s",