/****************************************************************************** * * Project: GDAL algorithms * Purpose: Apply vertical shift grid * Author: Even Rouault, even.rouault at spatialys.com * ****************************************************************************** * Copyright (c) 2017, Even Rouault * * 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 #include /************************************************************************/ /* 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 >) 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 >) 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(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 = cpl::down_cast(poDS); return poGDS->m_poSrcDataset->GetRasterBand(1)->GetNoDataValue(pbSuccess); } /************************************************************************/ /* IReadBlock() */ /************************************************************************/ CPLErr GDALApplyVSGRasterBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pData) { GDALApplyVSGDataset *poGDS = cpl::down_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 (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(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(). * * @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(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) ? static_cast(-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(srcGT); poReprojectedGrid->AddBand(GDT_Float32, nullptr); GDALApplyVSGDataset *poOutDS = new GDALApplyVSGDataset( GDALDataset::FromHandle(hSrcDataset), poReprojectedGrid, eDT, CPL_TO_BOOL(bInverse), static_cast(dfSrcUnitToMeter), static_cast(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); } /************************************************************************/ /* 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; }