Supplying proj pipeline to gdalwarp changes output geotransform

[scottyhq]

What is the bug?

I noticed that supplying -ct to gdalwarp can impact the output geotransform. For the case below, there is only one valid pipeline (applying a vertical shift grid per projinfo -s EPSG:9055+3855 -t EPSG:7661 -o proj)

I was expecting outputs to be the same (assuming GDAL behind the scenes must chose the applicable proj pipeline if -ct is not supplied)?

Perhaps the case below is exaggerated or an edge case because of using a VRT with global extent. Despite seeing lots of ERROR 1: PROJ: vgridshift: Invalid latitude in the output, however, both commands succeed and querying the resulting VRTs does give warped (but slightly different values).

Steps to reproduce the issue

%%bash

PROJ_PIPELINE='+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1'
S_SRS='EPSG:9055+3855'
T_SRS='EPSG:7661'

INPUT=vrt:///vsicurl/https://opentopography.s3.sdsc.edu/raster/COP30/COP30_hh.vrt?a_srs=$S_SRS

# NO -ct
CPL_DEBUG=OFF \
 PROJ_NETWORK=ON \
 GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdalwarp -overwrite -wm 500 \
    -s_srs $S_SRS -t_srs $T_SRS \
    ${INPUT} warp_nopipeline.vrt

# With -ct
CPL_DEBUG=OFF \
 PROJ_NETWORK=ON \
 GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR \
 gdalwarp -overwrite -wm 500 \
    -s_srs $S_SRS -t_srs $T_SRS \
    -ct "${PROJ_PIPELINE}" \
    ${INPUT} warp_pipeline.vrt

diff warp_nopipeline.vrt warp_pipeline.vrt

1c1
< <VRTDataset rasterXSize="1296001" rasterYSize="626401" subClass="VRTWarpedDataset">
---
> <VRTDataset rasterXSize="1298450" rasterYSize="621309" subClass="VRTWarpedDataset">
3c3
<   <GeoTransform> -1.8000013890000000e+02,  2.7777777777777778e-04,  0.0000000000000000e+00,  8.4000138899999996e+01,  0.0000000000000000e+00, -2.7777777777777778e-04</GeoTransform>
---
>   <GeoTransform> -1.7999997410495095e+02,  2.7725346045338265e-04,  0.0000000000000000e+00,  8.4000138899999996e+01,  0.0000000000000000e+00, -2.7725367398002905e-04</GeoTransform>
26,27c26,27
<         <DstGeoTransform>-180.0001389,0.00027777777777777778,0,84.000138899999996,0,-0.00027777777777777778</DstGeoTransform>
<         <DstInvGeoTransform>648000.50003999996,3600,0,302400.50004000001,0,-3600</DstInvGeoTransform>
---
>         <DstGeoTransform>-179.99997410495095,0.00027725346045338265,0,84.000138899999996,0,-0.00027725367398002905</DstGeoTransform>
>         <DstInvGeoTransform>649225.34712678951,3606.8080029181092,0,302972.13989687525,0,-3606.8052251384461</DstInvGeoTransform>
33a34
>               <Option key="COORDINATE_OPERATION">+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1 +step +proj=unitconvert +xy_in=rad +xy_out=deg +step +proj=axisswap +order=2,1</Option>

Versions and provenance

gdal 3.10.0 from conda-forge

Additional context

No response

[rouault]

For the case below, there is only one valid pipeline (applying a vertical shift grid per projinfo -s EPSG:9055+3855 -t EPSG:7661 -o proj)

actually, there are 2. The one using the grid which is "only" valid for BBOX[-90,-180,90,180], and the ballpark +proj=noop one which is valid ... "else where".
I haven't investigated but my guess would be that the logic in gdalwarp that guesses the output bounds must overshoot latitudes beyond 90, and when using -ct, then this causes reprojection failures, whereas when not using it, PROJ uses the +proj=noop operation.