From 81c4fac41d9aad4a96eacfdf63988360f8366d30 Mon Sep 17 00:00:00 2001 From: Even Rouault Date: Fri, 18 Mar 2022 13:29:59 +0100 Subject: [PATCH] Merge pull request #5466 from rouault/fix_vgridshift gdalwarp: fix issue with vertical shift mode when only one of the CRS is 3D, with PROJ 9.1 --- autotest/alg/applyverticalshiftgrid.py | 6 +- gdal/alg/gdaltransformer.cpp | 6 ++ gdal/apps/gdalwarp_lib.cpp | 79 +++++++++++++++++++------- 3 files changed, 68 insertions(+), 23 deletions(-) diff --git a/autotest/alg/applyverticalshiftgrid.py b/autotest/alg/applyverticalshiftgrid.py index 2c15fb78762e..b95731af19a8 100755 --- a/autotest/alg/applyverticalshiftgrid.py +++ b/autotest/alg/applyverticalshiftgrid.py @@ -295,7 +295,8 @@ def test_applyverticalshiftgrid_6(): grid_ds = None ds = gdal.Warp('', '../gcore/data/byte.tif', format='MEM', - dstSRS='+proj=utm +zone=11 +datum=NAD27 +geoidgrids=./tmp/applyverticalshiftgrid_6.gtx +vunits=m +no_defs') + srcSRS='EPSG:32611', + dstSRS='+proj=utm +zone=11 +datum=WGS84 +geoidgrids=./tmp/applyverticalshiftgrid_6.gtx +vunits=m +no_defs') cs = ds.GetRasterBand(1).Checksum() assert cs == 4783 @@ -314,7 +315,8 @@ def test_applyverticalshiftgrid_7(): grid_ds = None ds = gdal.Warp('', '../gcore/data/byte.tif', format='MEM', - dstSRS='+proj=utm +zone=11 +datum=NAD27 +geoidgrids=./tmp/applyverticalshiftgrid_7.gtx +vunits=m +no_defs') + srcSRS='EPSG:32611', + dstSRS='+proj=utm +zone=11 +datum=WGS84 +geoidgrids=./tmp/applyverticalshiftgrid_7.gtx +vunits=m +no_defs') cs = ds.GetRasterBand(1).Checksum() assert cs == 4783 diff --git a/gdal/alg/gdaltransformer.cpp b/gdal/alg/gdaltransformer.cpp index d57235dad194..33de06df6538 100644 --- a/gdal/alg/gdaltransformer.cpp +++ b/gdal/alg/gdaltransformer.cpp @@ -2151,6 +2151,12 @@ GDALCreateGenImgProjTransformer2( GDALDatasetH hSrcDS, GDALDatasetH hDstDS, InsertCenterLong( hSrcDS, &oSrcSRS, aosOptions ); } + if( CPLFetchBool(papszOptions, "PROMOTE_TO_3D", false) ) + { + oSrcSRS.PromoteTo3D(nullptr); + oDstSRS.PromoteTo3D(nullptr); + } + if( !(dfWestLongitudeDeg == 0.0 && dfSouthLatitudeDeg == 0.0 && dfEastLongitudeDeg == 0.0 && dfNorthLatitudeDeg == 0.0) ) { diff --git a/gdal/apps/gdalwarp_lib.cpp b/gdal/apps/gdalwarp_lib.cpp index 2f8ef9baf788..7728dbf0330e 100644 --- a/gdal/apps/gdalwarp_lib.cpp +++ b/gdal/apps/gdalwarp_lib.cpp @@ -559,27 +559,16 @@ GDALWarpAppOptions* GDALWarpAppOptionsClone(const GDALWarpAppOptions *psOptionsI #ifdef USE_PROJ_BASED_VERTICAL_SHIFT_METHOD -/************************************************************************/ -/* ApplyVerticalShift() */ -/************************************************************************/ - -static bool ApplyVerticalShift( GDALDatasetH hWrkSrcDS, - const GDALWarpAppOptions* psOptions, - GDALWarpOptions *psWO ) +static bool MustApplyVerticalShift( GDALDatasetH hWrkSrcDS, + const GDALWarpAppOptions* psOptions, + OGRSpatialReference& oSRSSrc, + OGRSpatialReference& oSRSDst, + bool& bSrcHasVertAxis, + bool& bDstHasVertAxis ) { - bool bApplyVShift = false; - - if( psOptions->bVShift ) - { - bApplyVShift = true; - psWO->papszWarpOptions = CSLSetNameValue(psWO->papszWarpOptions, - "APPLY_VERTICAL_SHIFT", - "YES"); - } + bool bApplyVShift = psOptions->bVShift; // Check if we must do vertical shift grid transform - OGRSpatialReference oSRSSrc; - OGRSpatialReference oSRSDst; const char* pszSrcWKT = CSLFetchNameValue(psOptions->papszTO, "SRC_SRS"); if( pszSrcWKT ) oSRSSrc.SetFromUserInput( pszSrcWKT ); @@ -594,14 +583,46 @@ static bool ApplyVerticalShift( GDALDatasetH hWrkSrcDS, if( pszDstWKT ) oSRSDst.SetFromUserInput( pszDstWKT ); - const bool bSrcHasVertAxis = + bSrcHasVertAxis = oSRSSrc.IsCompound() || ((oSRSSrc.IsProjected() || oSRSSrc.IsGeographic()) && oSRSSrc.GetAxesCount() == 3); - const bool bDstHasVertAxis = + 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) ) { @@ -2228,6 +2249,22 @@ GDALDatasetH GDALWarpDirect( const char *pszDest, GDALDatasetH hDstDS, psOptions->papszTO = CSLSetNameValue(psOptions->papszTO, "STRIP_VERT_CS", "YES"); } + else if( nSrcCount > 0 ) + { + bool bSrcHasVertAxis = false; + bool bDstHasVertAxis = false; + OGRSpatialReference oSRSSrc; + OGRSpatialReference oSRSDst; + bool bApplyVShift = MustApplyVerticalShift( pahSrcDS[0], psOptions, + oSRSSrc, oSRSDst, + bSrcHasVertAxis, + bDstHasVertAxis ); + if( bApplyVShift ) + { + psOptions->papszTO = CSLSetNameValue(psOptions->papszTO, + "PROMOTE_TO_3D", "YES"); + } + } #else psOptions->papszTO = CSLSetNameValue(psOptions->papszTO, "STRIP_VERT_CS", "YES"); @@ -2715,7 +2752,7 @@ GDALDatasetH GDALWarpDirect( const char *pszDest, GDALDatasetH hDstDS, #ifdef USE_PROJ_BASED_VERTICAL_SHIFT_METHOD if( !psOptions->bNoVShift ) { - // Can modify both psWO->papszWarpOptions + // Can modify psWO->papszWarpOptions if( ApplyVerticalShift( hWrkSrcDS, psOptions, psWO ) ) { bUseApproxTransformer = false;