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

 *

 * Project:  GDAL algorithms

 * Purpose:  Apply vertical shift grid

 * Author:   Even Rouault, even.rouault at spatialys.com

 *

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

 * Copyright (c) 2017, Even Rouault <even.rouault at spatialys.com>

 *

 * SPDX-License-Identifier: MIT

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

#include "cpl_string.h"

#include "gdal.h"

#include "gdal_alg.h"

#include "gdal_alg_priv.h"

#include "gdal_priv.h"

#include "gdal_utils.h"

#include "gdalwarper.h"

#include "vrtdataset.h"

#include "ogr_spatialref.h"

#include "proj.h"

#include <cmath>

#include <limits>

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

/*                        GDALApplyVSGDataset                           */

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

class GDALApplyVSGDataset final : public GDALDataset

{

    friend class GDALApplyVSGRasterBand;

    GDALDataset *m_poSrcDataset = nullptr;

    GDALDataset *m_poReprojectedGrid = nullptr;

    bool m_bInverse = false;

    float m_fSrcUnitToMeter = 0.0f;

    float m_fDstUnitToMeter = 0.0f;

    CPL_DISALLOW_COPY_ASSIGN(GDALApplyVSGDataset)

  public:

    GDALApplyVSGDataset(GDALDataset *poSrcDataset,

                        GDALDataset *poReprojectedGrid, GDALDataType eDT,

                        bool bInverse, float fSrcUnitToMeter,

                        float fDstUnitToMeter, int nBlockSize);

    ~GDALApplyVSGDataset() override;

    int CloseDependentDatasets() override;

    CPLErr GetGeoTransform(GDALGeoTransform &gt) const override;

    const OGRSpatialReference *GetSpatialRef() const override;

    bool IsInitOK();

};

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

/*                       GDALApplyVSGRasterBand                         */

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

class GDALApplyVSGRasterBand final : public GDALRasterBand

{

    friend class GDALApplyVSGDataset;

    float *m_pafSrcData = nullptr;

    float *m_pafGridData = nullptr;

    CPL_DISALLOW_COPY_ASSIGN(GDALApplyVSGRasterBand)

  public:

    GDALApplyVSGRasterBand(GDALDataType eDT, int nBlockSize);

    ~GDALApplyVSGRasterBand() override;

    virtual CPLErr IReadBlock(int nBlockXOff, int nBlockYOff,

                              void *pData) override;

    double GetNoDataValue(int *pbSuccess) override;

};

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

/*                        GDALApplyVSGDataset()                         */

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

GDALApplyVSGDataset::GDALApplyVSGDataset(GDALDataset *poSrcDataset,

                                         GDALDataset *poReprojectedGrid,

                                         GDALDataType eDT, bool bInverse,

                                         float fSrcUnitToMeter,

                                         float fDstUnitToMeter, int nBlockSize)

    : m_poSrcDataset(poSrcDataset), m_poReprojectedGrid(poReprojectedGrid),

      m_bInverse(bInverse), m_fSrcUnitToMeter(fSrcUnitToMeter),

      m_fDstUnitToMeter(fDstUnitToMeter)

{

    m_poSrcDataset->Reference();

    m_poReprojectedGrid->Reference();

    nRasterXSize = poSrcDataset->GetRasterXSize();

    nRasterYSize = poSrcDataset->GetRasterYSize();

    SetBand(1, new GDALApplyVSGRasterBand(eDT, nBlockSize));

}

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

/*                       ~GDALApplyVSGDataset()                         */

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

GDALApplyVSGDataset::~GDALApplyVSGDataset()

{

    GDALApplyVSGDataset::CloseDependentDatasets();

}

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

/*                     CloseDependentDatasets()                         */

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

int GDALApplyVSGDataset::CloseDependentDatasets()

{

    bool bRet = false;

    if (m_poSrcDataset != nullptr)

    {

        if (m_poSrcDataset->ReleaseRef())

        {

            bRet = true;

        }

        m_poSrcDataset = nullptr;

    }

    if (m_poReprojectedGrid != nullptr)

    {

        if (m_poReprojectedGrid->ReleaseRef())

        {

            bRet = true;

        }

        m_poReprojectedGrid = nullptr;

    }

    return bRet;

}

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

/*                          GetGeoTransform()                           */

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

CPLErr GDALApplyVSGDataset::GetGeoTransform(GDALGeoTransform &gt) const

{

    return m_poSrcDataset->GetGeoTransform(gt);

}

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

/*                          GetSpatialRef()                             */

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

const OGRSpatialReference *GDALApplyVSGDataset::GetSpatialRef() const

{

    return m_poSrcDataset->GetSpatialRef();

}

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

/*                             IsInitOK()                               */

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

bool GDALApplyVSGDataset::IsInitOK()

{

    GDALApplyVSGRasterBand *poBand =

        cpl::down_cast<GDALApplyVSGRasterBand *>(GetRasterBand(1));

    return poBand->m_pafSrcData != nullptr && poBand->m_pafGridData != nullptr;

}

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

/*                       GDALApplyVSGRasterBand()                       */

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

GDALApplyVSGRasterBand::GDALApplyVSGRasterBand(GDALDataType eDT, int nBlockSize)

{

    eDataType = eDT;

    nBlockXSize = nBlockSize;

    nBlockYSize = nBlockSize;

    m_pafSrcData = static_cast<float *>(

        VSI_MALLOC3_VERBOSE(nBlockXSize, nBlockYSize, sizeof(float)));

    m_pafGridData = static_cast<float *>(

        VSI_MALLOC3_VERBOSE(nBlockXSize, nBlockYSize, sizeof(float)));

}

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

/*                      ~GDALApplyVSGRasterBand()                       */

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

GDALApplyVSGRasterBand::~GDALApplyVSGRasterBand()

{

    VSIFree(m_pafSrcData);

    VSIFree(m_pafGridData);

}

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

/*                           GetNoDataValue()                           */

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

double GDALApplyVSGRasterBand::GetNoDataValue(int *pbSuccess)

{

    GDALApplyVSGDataset *poGDS = cpl::down_cast<GDALApplyVSGDataset *>(poDS);

    return poGDS->m_poSrcDataset->GetRasterBand(1)->GetNoDataValue(pbSuccess);

}

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

/*                              IReadBlock()                            */

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

CPLErr GDALApplyVSGRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff,

                                          void *pData)

{

    GDALApplyVSGDataset *poGDS = cpl::down_cast<GDALApplyVSGDataset *>(poDS);

    const int nXOff = nBlockXOff * nBlockXSize;

    const int nReqXSize = (nXOff > nRasterXSize - nBlockXSize)

                              ? nRasterXSize - nXOff

                              : nBlockXSize;

    const int nYOff = nBlockYOff * nBlockYSize;

    const int nReqYSize = (nYOff > nRasterYSize - nBlockYSize)

                              ? nRasterYSize - nYOff

                              : nBlockYSize;

    CPLErr eErr = poGDS->m_poSrcDataset->GetRasterBand(1)->RasterIO(

        GF_Read, nXOff, nYOff, nReqXSize, nReqYSize, m_pafSrcData, nReqXSize,

        nReqYSize, GDT_Float32, sizeof(float), nBlockXSize * sizeof(float),

        nullptr);

    if (eErr == CE_None)

        eErr = poGDS->m_poReprojectedGrid->GetRasterBand(1)->RasterIO(

            GF_Read, nXOff, nYOff, nReqXSize, nReqYSize, m_pafGridData,

            nReqXSize, nReqYSize, GDT_Float32, sizeof(float),

            nBlockXSize * sizeof(float), nullptr);

    if (eErr == CE_None)

    {

        const int nDTSize = GDALGetDataTypeSizeBytes(eDataType);

        int bHasNoData = FALSE;

        float fNoDataValue = static_cast<float>(GetNoDataValue(&bHasNoData));

        for (int iY = 0; iY < nReqYSize; iY++)

        {

            for (int iX = 0; iX < nReqXSize; iX++)

            {

                const float fSrcVal = m_pafSrcData[iY * nBlockXSize + iX];

                const float fGridVal = m_pafGridData[iY * nBlockXSize + iX];

                if (bHasNoData && fSrcVal == fNoDataValue)

                {

                }

                else if (std::isinf(fGridVal))

                {

                    CPLError(CE_Failure, CPLE_AppDefined,

                             "Missing vertical grid value at source (%d,%d)",

                             nXOff + iX, nYOff + iY);

                    return CE_Failure;

                }

                else if (poGDS->m_bInverse)

                {

                    m_pafSrcData[iY * nBlockXSize + iX] =

                        (fSrcVal * poGDS->m_fSrcUnitToMeter - fGridVal) /

                        poGDS->m_fDstUnitToMeter;

                }

                else

                {

                    m_pafSrcData[iY * nBlockXSize + iX] =

                        (fSrcVal * poGDS->m_fSrcUnitToMeter + fGridVal) /

                        poGDS->m_fDstUnitToMeter;

                }

            }

            GDALCopyWords(

                m_pafSrcData + iY * nBlockXSize, GDT_Float32, sizeof(float),

                static_cast<GByte *>(pData) + iY * nBlockXSize * nDTSize,

                eDataType, nDTSize, nReqXSize);

        }

    }

    return eErr;

}

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

/*                      GDALApplyVerticalShiftGrid()                    */

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

/** Apply a vertical shift grid to a source (DEM typically) dataset.

 *

 * hGridDataset will typically use WGS84 as horizontal datum (but this is

 * not a requirement) and its values are the values to add to go from geoid

 * elevations to WGS84 ellipsoidal heights.

 *

 * hGridDataset will be on-the-fly reprojected and resampled to the projection

 * and resolution of hSrcDataset, using bilinear resampling by default.

 *

 * Both hSrcDataset and hGridDataset must be single band datasets, and have

 * a valid geotransform and projection.

 *

 * On success, a reference will be taken on hSrcDataset and hGridDataset.

 * Reference counting semantics on the source and grid datasets should be

 * honoured. That is, don't just GDALClose() it, unless it was opened with

 * GDALOpenShared(), but rather use GDALReleaseDataset() if wanting to

 * immediately release the reference(s) and make the returned dataset the

 * owner of them.

 *

 * Valid use cases:

 *

 * \code

 * hSrcDataset = GDALOpen(...)

 * hGridDataset = GDALOpen(...)

 * hDstDataset = GDALApplyVerticalShiftGrid(hSrcDataset, hGridDataset, ...)

 * GDALReleaseDataset(hSrcDataset);

 * GDALReleaseDataset(hGridDataset);

 * if( hDstDataset )

 * {

 *     // Do things with hDstDataset

 *     GDALClose(hDstDataset) // will close hSrcDataset and hGridDataset

 * }

 * \endcode

 *

 * @param hSrcDataset source (DEM) dataset. Must not be NULL.

 * @param hGridDataset vertical grid shift dataset. Must not be NULL.

 * @param bInverse if set to FALSE, hGridDataset values will be added to

 *                 hSrcDataset. If set to TRUE, they will be subtracted.

 * @param dfSrcUnitToMeter the factor to convert values from hSrcDataset to

 *                         meters (1.0 if source values are in meter).

 * @param dfDstUnitToMeter the factor to convert shifted values from meter

 *                          (1.0 if output values must be in meter).

 * @param papszOptions list of options, or NULL. Supported options are:

 * <ul>

 * <li>RESAMPLING=NEAREST/BILINEAR/CUBIC. Defaults to BILINEAR.</li>

 * <li>MAX_ERROR=val. Maximum error measured in input pixels that is allowed in

 * approximating the transformation (0.0 for exact calculations). Defaults

 * to 0.125</li>

 * <li>DATATYPE=Byte/UInt16/Int16/Float32/Float64. Output data type. If not

 * specified will be the same as the one of hSrcDataset.

 * <li>ERROR_ON_MISSING_VERT_SHIFT=YES/NO. Whether a missing/nodata value in

 * hGridDataset should cause I/O requests to fail. Default is NO (in which case

 * 0 will be used)

 * <li>SRC_SRS=srs_def. Override projection on hSrcDataset;

 * </ul>

 *

 * @return a new dataset corresponding to hSrcDataset adjusted with

 * hGridDataset, or NULL. If not NULL, it must be closed with GDALClose().

 *

 * @deprecated GDAL 3.4. Will be removed in GDAL 4.0. This function was used

 *             by gdalwarp initially, but is no longer needed.

 */

GDALDatasetH GDALApplyVerticalShiftGrid(GDALDatasetH hSrcDataset,

                                        GDALDatasetH hGridDataset, int bInverse,

                                        double dfSrcUnitToMeter,

                                        double dfDstUnitToMeter,

                                        const char *const *papszOptions)

{

    VALIDATE_POINTER1(hSrcDataset, "GDALApplyVerticalShiftGrid", nullptr);

    VALIDATE_POINTER1(hGridDataset, "GDALApplyVerticalShiftGrid", nullptr);

    GDALGeoTransform srcGT;

    if (GDALDataset::FromHandle(hSrcDataset)->GetGeoTransform(srcGT) != CE_None)

    {

        CPLError(CE_Failure, CPLE_NotSupported,

                 "Source dataset has no geotransform.");

        return nullptr;

    }

    const char *pszSrcProjection = CSLFetchNameValue(papszOptions, "SRC_SRS");

    OGRSpatialReference oSrcSRS;

    if (pszSrcProjection != nullptr && pszSrcProjection[0] != '\0')

    {

        oSrcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);

        oSrcSRS.SetFromUserInput(pszSrcProjection);

    }

    else

    {

        auto poSRS = GDALDataset::FromHandle(hSrcDataset)->GetSpatialRef();

        if (poSRS)

            oSrcSRS = *poSRS;

    }

    if (oSrcSRS.IsCompound())

    {

        oSrcSRS.StripVertical();

    }

    if (oSrcSRS.IsEmpty())

    {

        CPLError(CE_Failure, CPLE_NotSupported,

                 "Source dataset has no projection.");

        return nullptr;

    }

    if (GDALGetRasterCount(hSrcDataset) != 1)

    {

        CPLError(CE_Failure, CPLE_NotSupported,

                 "Only single band source dataset is supported.");

        return nullptr;

    }

    GDALGeoTransform gridGT;

    if (GDALDataset::FromHandle(hGridDataset)->GetGeoTransform(gridGT) !=

        CE_None)

    {

        CPLError(CE_Failure, CPLE_NotSupported,

                 "Grid dataset has no geotransform.");

        return nullptr;

    }

    OGRSpatialReferenceH hGridSRS = GDALGetSpatialRef(hGridDataset);

    if (hGridSRS == nullptr)

    {

        CPLError(CE_Failure, CPLE_NotSupported,

                 "Grid dataset has no projection.");

        return nullptr;

    }

    if (GDALGetRasterCount(hGridDataset) != 1)

    {

        CPLError(CE_Failure, CPLE_NotSupported,

                 "Only single band grid dataset is supported.");

        return nullptr;

    }

    GDALDataType eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDataset, 1));

    const char *pszDataType = CSLFetchNameValue(papszOptions, "DATATYPE");

    if (pszDataType)

        eDT = GDALGetDataTypeByName(pszDataType);

    if (eDT == GDT_Unknown)

    {

        CPLError(CE_Failure, CPLE_NotSupported, "Invalid DATATYPE=%s",

                 pszDataType);

        return nullptr;

    }

    const int nSrcXSize = GDALGetRasterXSize(hSrcDataset);

    const int nSrcYSize = GDALGetRasterYSize(hSrcDataset);

    double dfWestLongitudeDeg = 0.0;

    double dfSouthLatitudeDeg = 0.0;

    double dfEastLongitudeDeg = 0.0;

    double dfNorthLatitudeDeg = 0.0;

    GDALComputeAreaOfInterest(&oSrcSRS, srcGT.data(), nSrcXSize, nSrcYSize,

                              dfWestLongitudeDeg, dfSouthLatitudeDeg,

                              dfEastLongitudeDeg, dfNorthLatitudeDeg);

    CPLStringList aosOptions;

    if (!(dfWestLongitudeDeg == 0.0 && dfSouthLatitudeDeg == 0.0 &&

          dfEastLongitudeDeg == 0.0 && dfNorthLatitudeDeg == 0.0))

    {

        aosOptions.SetNameValue(

            "AREA_OF_INTEREST",

            CPLSPrintf("%.16g,%.16g,%.16g,%.16g", dfWestLongitudeDeg,

                       dfSouthLatitudeDeg, dfEastLongitudeDeg,

                       dfNorthLatitudeDeg));

    }

    void *hTransform = GDALCreateGenImgProjTransformer4(

        hGridSRS, gridGT.data(), OGRSpatialReference::ToHandle(&oSrcSRS),

        srcGT.data(), aosOptions.List());

    if (hTransform == nullptr)

        return nullptr;

    GDALWarpOptions *psWO = GDALCreateWarpOptions();

    psWO->hSrcDS = hGridDataset;

    psWO->eResampleAlg = GRA_Bilinear;

    const char *pszResampling = CSLFetchNameValue(papszOptions, "RESAMPLING");

    if (pszResampling)

    {

        if (EQUAL(pszResampling, "NEAREST"))

            psWO->eResampleAlg = GRA_NearestNeighbour;

        else if (EQUAL(pszResampling, "BILINEAR"))

            psWO->eResampleAlg = GRA_Bilinear;

        else if (EQUAL(pszResampling, "CUBIC"))

            psWO->eResampleAlg = GRA_Cubic;

    }

    psWO->eWorkingDataType = GDT_Float32;

    int bHasNoData = FALSE;

    const double dfSrcNoData = GDALGetRasterNoDataValue(

        GDALGetRasterBand(hGridDataset, 1), &bHasNoData);

    if (bHasNoData)

    {

        psWO->padfSrcNoDataReal =

            static_cast<double *>(CPLMalloc(sizeof(double)));

        psWO->padfSrcNoDataReal[0] = dfSrcNoData;

    }

    psWO->padfDstNoDataReal = static_cast<double *>(CPLMalloc(sizeof(double)));

    const bool bErrorOnMissingShift =

        CPLFetchBool(papszOptions, "ERROR_ON_MISSING_VERT_SHIFT", false);

    psWO->padfDstNoDataReal[0] =

        (bErrorOnMissingShift)

            ? static_cast<double>(-std::numeric_limits<float>::infinity())

            : 0.0;

    psWO->papszWarpOptions =

        CSLSetNameValue(psWO->papszWarpOptions, "INIT_DEST", "NO_DATA");

    psWO->pfnTransformer = GDALGenImgProjTransform;

    psWO->pTransformerArg = hTransform;

    const double dfMaxError =

        CPLAtof(CSLFetchNameValueDef(papszOptions, "MAX_ERROR", "0.125"));

    if (dfMaxError > 0.0)

    {

        psWO->pTransformerArg = GDALCreateApproxTransformer(

            psWO->pfnTransformer, psWO->pTransformerArg, dfMaxError);

        psWO->pfnTransformer = GDALApproxTransform;

        GDALApproxTransformerOwnsSubtransformer(psWO->pTransformerArg, TRUE);

    }

    psWO->nBandCount = 1;

    psWO->panSrcBands = static_cast<int *>(CPLMalloc(sizeof(int)));

    psWO->panSrcBands[0] = 1;

    psWO->panDstBands = static_cast<int *>(CPLMalloc(sizeof(int)));

    psWO->panDstBands[0] = 1;

    VRTWarpedDataset *poReprojectedGrid =

        new VRTWarpedDataset(nSrcXSize, nSrcYSize);

    // This takes a reference on hGridDataset

    CPLErr eErr = poReprojectedGrid->Initialize(psWO);

    CPLAssert(eErr == CE_None);

    CPL_IGNORE_RET_VAL(eErr);

    GDALDestroyWarpOptions(psWO);

    poReprojectedGrid->SetGeoTransform(srcGT);

    poReprojectedGrid->AddBand(GDT_Float32, nullptr);

    GDALApplyVSGDataset *poOutDS = new GDALApplyVSGDataset(

        GDALDataset::FromHandle(hSrcDataset), poReprojectedGrid, eDT,

        CPL_TO_BOOL(bInverse), static_cast<float>(dfSrcUnitToMeter),

        static_cast<float>(dfDstUnitToMeter),

        // Undocumented option. For testing only

        atoi(CSLFetchNameValueDef(papszOptions, "BLOCKSIZE", "256")));

    poReprojectedGrid->ReleaseRef();

    if (!poOutDS->IsInitOK())

    {

        delete poOutDS;

        return nullptr;

    }

    poOutDS->SetDescription(GDALGetDescription(hSrcDataset));

    return reinterpret_cast<GDALDatasetH>(poOutDS);

}

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

/*                           GetProj4Filename()                         */

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

static CPLString GetProj4Filename(const char *pszFilename)

{

    CPLString osFilename;

    /* or fixed path: /name, ./name or ../name */

    if (!CPLIsFilenameRelative(pszFilename) || *pszFilename == '.')

    {

        return pszFilename;

    }

    PJ_GRID_INFO info = proj_grid_info(pszFilename);

    if (info.filename[0])

    {

        osFilename = info.filename;

    }

    return osFilename;

}

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

/*                       GDALOpenVerticalShiftGrid()                    */

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

/** Load proj.4 geoidgrids as GDAL dataset

 *

 * @param pszProj4Geoidgrids Value of proj.4 geoidgrids parameter.

 * @param pbError If not NULL, the pointed value will be set to TRUE if an

 *                error occurred.

 *

 * @return a dataset. If not NULL, it must be closed with GDALClose().

 *

 * @deprecated GDAL 3.4. Will be removed in GDAL 4.0. This function was used

 *             by gdalwarp initially, but is no longer needed.

 */

GDALDatasetH GDALOpenVerticalShiftGrid(const char *pszProj4Geoidgrids,

                                       int *pbError)

{

    char **papszGrids = CSLTokenizeString2(pszProj4Geoidgrids, ",", 0);

    const int nGridCount = CSLCount(papszGrids);

    if (nGridCount == 1)

    {

        CSLDestroy(papszGrids);

        bool bMissingOk = false;

        if (*pszProj4Geoidgrids == '@')

        {

            pszProj4Geoidgrids++;

            bMissingOk = true;

        }

        const CPLString osFilename(GetProj4Filename(pszProj4Geoidgrids));

        const char *const papszOpenOptions[] = {

            "@SHIFT_ORIGIN_IN_MINUS_180_PLUS_180=YES", nullptr};

        GDALDatasetH hDS =

            GDALOpenEx(osFilename, 0, nullptr, papszOpenOptions, nullptr);

        if (hDS == nullptr)

        {

            CPLDebug("GDAL", "Cannot find file corresponding to %s",

                     pszProj4Geoidgrids);

        }

        if (pbError)

            *pbError = (!bMissingOk && hDS == nullptr);

        return hDS;

    }

    CPLStringList aosFilenames;

    for (int i = nGridCount - 1; i >= 0; i--)

    {

        const char *pszName = papszGrids[i];

        bool bMissingOk = false;

        if (*pszName == '@')

        {

            pszName++;

            bMissingOk = true;

        }

        const CPLString osFilename(GetProj4Filename(pszName));

        VSIStatBufL sStat;

        if (osFilename.empty() || VSIStatL(osFilename, &sStat) != 0)

        {

            CPLDebug("GDAL", "Cannot find file corresponding to %s", pszName);

            if (!bMissingOk)

            {

                if (pbError)

                    *pbError = true;

                CSLDestroy(papszGrids);

                return nullptr;

            }

        }

        else

        {

            aosFilenames.AddString(osFilename);

        }

    }

    CSLDestroy(papszGrids);

    if (aosFilenames.empty())

    {

        if (pbError)

            *pbError = false;

        return nullptr;

    }

    char **papszArgv = nullptr;

    papszArgv = CSLAddString(papszArgv, "-resolution");

    papszArgv = CSLAddString(papszArgv, "highest");

    papszArgv = CSLAddString(papszArgv, "-vrtnodata");

    papszArgv = CSLAddString(papszArgv, "-inf");

    papszArgv = CSLAddString(papszArgv, "-oo");

    papszArgv =

        CSLAddString(papszArgv, "@SHIFT_ORIGIN_IN_MINUS_180_PLUS_180=YES");

    GDALBuildVRTOptions *psOptions = GDALBuildVRTOptionsNew(papszArgv, nullptr);

    CSLDestroy(papszArgv);

    GDALDatasetH hDS = GDALBuildVRT("", aosFilenames.size(), nullptr,

                                    aosFilenames.List(), psOptions, nullptr);

    GDALBuildVRTOptionsFree(psOptions);

    if (pbError)

        *pbError = hDS != nullptr;

    return hDS;

}