/****************************************************************************** * * Project: GDAL algorithms * Purpose: Apply vertical shift grid * Author: Even Rouault, even.rouault at spatialys.com * ****************************************************************************** * Copyright (c) 2017, Even Rouault * * Permission is hereby granted, free of charge, to any person obtaining a * copy of this software and associated documentation files (the "Software"), * to deal in the Software without restriction, including without limitation * the rights to use, copy, modify, merge, publish, distribute, sublicense, * and/or sell copies of the Software, and to permit persons to whom the * Software is furnished to do so, subject to the following conditions: * * The above copyright notice and this permission notice shall be included * in all copies or substantial portions of the Software. * * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER * DEALINGS IN THE SOFTWARE. ****************************************************************************/ #include "cpl_string.h" #include "gdal.h" #include "gdal_alg.h" #include "gdal_priv.h" #include "gdal_utils.h" #include "gdalwarper.h" #include "vrtdataset.h" #include "ogr_spatialref.h" #ifdef PROJ_STATIC #if defined(PROJ_VERSION) && PROJ_VERSION >= 5 #include "proj.h" #else #include "proj_api.h" #endif #endif #include CPL_CVSID("$Id: gdalapplyverticalshiftgrid.cpp fe2d81c8819bf9794bce0210098e637565728350 2018-05-06 00:49:51 +0200 Even Rouault $") /************************************************************************/ /* GDALApplyVSGDataset */ /************************************************************************/ class GDALApplyVSGDataset final: public GDALDataset { friend class GDALApplyVSGRasterBand; GDALDataset* m_poSrcDataset = nullptr; GDALDataset* m_poReprojectedGrid = nullptr; bool m_bInverse = false; double m_dfSrcUnitToMeter = 0.0; double m_dfDstUnitToMeter = 0.0; CPL_DISALLOW_COPY_ASSIGN(GDALApplyVSGDataset) public: GDALApplyVSGDataset( GDALDataset* poSrcDataset, GDALDataset* poReprojectedGrid, GDALDataType eDT, bool bInverse, double dfSrcUnitToMeter, double dfDstUnitToMeter, int nBlockSize ); virtual ~GDALApplyVSGDataset(); virtual int CloseDependentDatasets() override; virtual CPLErr GetGeoTransform(double* padfGeoTransform) override; virtual const char* GetProjectionRef() 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 ); virtual ~GDALApplyVSGRasterBand(); virtual CPLErr IReadBlock( int nBlockXOff, int nBlockYOff, void * pData ) override; virtual double GetNoDataValue( int* pbSuccess ) override; }; /************************************************************************/ /* GDALApplyVSGDataset() */ /************************************************************************/ GDALApplyVSGDataset::GDALApplyVSGDataset( GDALDataset* poSrcDataset, GDALDataset* poReprojectedGrid, GDALDataType eDT, bool bInverse, double dfSrcUnitToMeter, double dfDstUnitToMeter, int nBlockSize ) : m_poSrcDataset(poSrcDataset), m_poReprojectedGrid(poReprojectedGrid), m_bInverse(bInverse), m_dfSrcUnitToMeter(dfSrcUnitToMeter), m_dfDstUnitToMeter(dfDstUnitToMeter) { 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(double* padfGeoTransform) { return m_poSrcDataset->GetGeoTransform(padfGeoTransform); } /************************************************************************/ /* GetProjectionRef() */ /************************************************************************/ const char* GDALApplyVSGDataset::GetProjectionRef() { return m_poSrcDataset->GetProjectionRef(); } /************************************************************************/ /* IsInitOK() */ /************************************************************************/ bool GDALApplyVSGDataset::IsInitOK() { GDALApplyVSGRasterBand* poBand = reinterpret_cast(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( VSI_MALLOC3_VERBOSE(nBlockXSize, nBlockYSize, sizeof(float))); m_pafGridData = static_cast( VSI_MALLOC3_VERBOSE(nBlockXSize, nBlockYSize, sizeof(float))); } /************************************************************************/ /* ~GDALApplyVSGRasterBand() */ /************************************************************************/ GDALApplyVSGRasterBand::~GDALApplyVSGRasterBand() { VSIFree(m_pafSrcData); VSIFree(m_pafGridData); } /************************************************************************/ /* GetNoDataValue() */ /************************************************************************/ double GDALApplyVSGRasterBand::GetNoDataValue( int* pbSuccess ) { GDALApplyVSGDataset* poGDS = reinterpret_cast(poDS); return poGDS->m_poSrcDataset->GetRasterBand(1)->GetNoDataValue(pbSuccess); } /************************************************************************/ /* IReadBlock() */ /************************************************************************/ CPLErr GDALApplyVSGRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff, void * pData ) { GDALApplyVSGDataset* poGDS = reinterpret_cast(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(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( CPLIsInf(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] = static_cast( (fSrcVal * poGDS->m_dfSrcUnitToMeter - fGridVal) / poGDS->m_dfDstUnitToMeter); } else { m_pafSrcData[iY * nBlockXSize + iX] = static_cast( (fSrcVal * poGDS->m_dfSrcUnitToMeter + fGridVal) / poGDS->m_dfDstUnitToMeter); } } GDALCopyWords( m_pafSrcData + iY * nBlockXSize, GDT_Float32, sizeof(float), static_cast(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: *
    *
  • RESAMPLING=NEAREST/BILINEAR/CUBIC. Defaults to BILINEAR.
    • *
  • 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
    • *
  • DATATYPE=Byte/UInt16/Int16/Float32/Float64. Output data type. If not * specified will be the same as the one of hSrcDataset. *
  • 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) *
  • SRC_SRS=srs_def. Override projection on hSrcDataset; *
* * @return a new dataset corresponding to hSrcDataset adjusted with * hGridDataset, or NULL. If not NULL, it must be closed with GDALClose(). * * @since GDAL 2.2 */ 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 ); double adfSrcGT[6]; if( GDALGetGeoTransform(hSrcDataset, adfSrcGT) != CE_None ) { CPLError(CE_Failure, CPLE_NotSupported, "Source dataset has no geotransform."); return nullptr; } const char* pszSrcProjection = CSLFetchNameValueDef(papszOptions, "SRC_SRS", GDALGetProjectionRef(hSrcDataset)); if( pszSrcProjection == nullptr || pszSrcProjection[0] == '\0' ) { 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; } double adfGridGT[6]; if( GDALGetGeoTransform(hGridDataset, adfGridGT) != CE_None ) { CPLError(CE_Failure, CPLE_NotSupported, "Grid dataset has no geotransform."); return nullptr; } const char* pszGridProjection = GDALGetProjectionRef(hGridDataset); if( pszGridProjection == nullptr || pszGridProjection[0] == '\0' ) { 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); OGRSpatialReference oSRS; CPLString osSrcProjection(pszSrcProjection); oSRS.SetFromUserInput(osSrcProjection); if( oSRS.IsCompound() ) { OGR_SRSNode* poNode = oSRS.GetRoot()->GetChild(1); if( poNode != nullptr ) { char* pszWKT = nullptr; poNode->exportToWkt(&pszWKT); osSrcProjection = pszWKT; CPLFree(pszWKT); } } void* hTransform = GDALCreateGenImgProjTransformer3( pszGridProjection, adfGridGT, osSrcProjection, adfSrcGT ); 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(CPLMalloc(sizeof(double))); psWO->padfSrcNoDataReal[0] = dfSrcNoData; } psWO->padfDstNoDataReal = static_cast(CPLMalloc(sizeof(double))); const bool bErrorOnMissingShift = CPLFetchBool( papszOptions, "ERROR_ON_MISSING_VERT_SHIFT", false ); psWO->padfDstNoDataReal[0] = (bErrorOnMissingShift) ? -std::numeric_limits::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(CPLMalloc(sizeof(int))); psWO->panSrcBands[0] = 1; psWO->panDstBands = static_cast(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(adfSrcGT); poReprojectedGrid->AddBand(GDT_Float32, nullptr); GDALApplyVSGDataset* poOutDS = new GDALApplyVSGDataset( reinterpret_cast(hSrcDataset), poReprojectedGrid, eDT, CPL_TO_BOOL(bInverse), dfSrcUnitToMeter, 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(poOutDS); } /************************************************************************/ /* my_proj4_logger() */ /************************************************************************/ #if defined(PROJ_STATIC) && PJ_VERSION >= 490 && PJ_VERSION <= 493 static void my_proj4_logger(void * user_data, int /*level*/, const char * msg) { CPLString* posMsg = static_cast(user_data); *posMsg += msg; } #endif /************************************************************************/ /* GetProj4Filename() */ /************************************************************************/ static CPLString GetProj4Filename(const char* pszFilename) { CPLString osFilename; /* or fixed path: /name, ./name or ../name */ if ( !CPLIsFilenameRelative(pszFilename) || *pszFilename == '.' ) { return pszFilename; } #if defined(PROJ_STATIC) && PROJ_VERSION >= 5 PJ_GRID_INFO info = proj_grid_info(pszFilename); if( info.filename[0] ) { osFilename = info.filename; } #elif defined(PROJ_STATIC) && PJ_VERSION > 493 osFilename.resize(2048); projCtx ctx = pj_ctx_alloc(); if( pj_find_file(ctx, pszFilename, &osFilename[0], osFilename.size()) ) { osFilename.resize( strlen(osFilename) ); } else { osFilename.clear(); } pj_ctx_free(ctx); #else // Transpose some of the proj.4 pj_open_lib() logic... /* check if ~/name */ char* pszSysname; if (*pszFilename == '~' && (pszFilename[1] == '/' || pszFilename[1] == '\\') ) { if ((pszSysname = getenv("HOME")) != nullptr) { osFilename = CPLFormFilename(pszSysname, pszFilename + 1, nullptr); } return osFilename; } /* or is environment PROJ_LIB defined */ else if ((pszSysname = getenv("PROJ_LIB")) != nullptr) { osFilename = CPLFormFilename(pszSysname, pszFilename, nullptr); VSIStatBufL sStat; if( VSIStatL(osFilename, &sStat) == 0 ) return osFilename; osFilename.clear(); } #if defined(PROJ_STATIC) && PJ_VERSION >= 490 // Super messy. proj.4 up to 4.9.3 had no public API to return the full // path to a resource file, so we rely on the fact that it emits a log // message with it... // Basically this is needed in the case where the file is in the // resource installation directory of proj.4, which we have no way to // know otherwise. CPLString osMsg; projCtx ctx = pj_ctx_alloc(); pj_ctx_set_app_data(ctx, &osMsg); pj_ctx_set_debug(ctx, PJ_LOG_DEBUG_MAJOR); pj_ctx_set_logger(ctx, my_proj4_logger); PAFile f = pj_open_lib(ctx, pszFilename, "rb"); if( f ) { pj_ctx_fclose(ctx, f); size_t nPos = osMsg.find("fopen("); if( nPos != std::string::npos ) { osFilename = osMsg.substr(nPos + strlen("fopen(")); nPos = osFilename.find(")"); if( nPos != std::string::npos ) osFilename = osFilename.substr(0, nPos); } } pj_ctx_free(ctx); #endif #endif 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(). * * @since GDAL 2.2 */ 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; }