Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[Backport release/3.4] gdalwarp: fix issue with vertical shift mode when only one of the CRS is 3D, with PROJ 9.1 #5467

Merged
merged 1 commit into from
Mar 18, 2022
Merged

[Backport release/3.4] gdalwarp: fix issue with vertical shift mode when only one of the CRS is 3D, with PROJ 9.1 #5467

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 4 additions & 2 deletions autotest/alg/applyverticalshiftgrid.py
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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