Skip to content

Commit

Permalink
Merge pull request #5467 from OSGeo/backport-5466-to-release/3.4
Browse files Browse the repository at this point in the history 
[Backport release/3.4] gdalwarp: fix issue with vertical shift mode when only one of the CRS is 3D, with PROJ 9.1
  • Loading branch information
@rouault
rouault authored Mar 18, 2022
2 parents a22aa06 + 81c4fac commit 2633fe3
Show file tree
Hide file tree
Showing 3 changed files with 68 additions and 23 deletions.
6 changes: 4 additions & 2 deletions autotest/alg/applyverticalshiftgrid.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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

Expand Down
6 changes: 6 additions & 0 deletions gdal/alg/gdaltransformer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) )
{
Expand Down
79 changes: 58 additions & 21 deletions gdal/apps/gdalwarp_lib.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 );
Expand All @@ -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) )
{
Expand Down Expand Up @@ -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");
Expand Down Expand Up @@ -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;
Expand Down

0 comments on commit 2633fe3

Please sign in to comment.