gdalwarp: incorrect shift when reprojecting a raster at the edge of the input SRS validity zone

[olivierdalang]

What is the bug?

Hi !

We spotted some issues when using gdalwarp at the edge of EPSG:2056 (Switzerland) to EPSG:3857.

The zone in red is the area of validity of the SRS (as shown on https://epsg.io/2056), yet you can see a clear shear line inside the validity zone.

The issue seems similar to #10510 (related to auto guessing transformation outside the SRS) but in our case the shift happens inside the SRS.

I'm not very knowledgeable about CRS transformations, naively I would expect gdalwarp to honor the source SRS even outside of the validity domain, but I guess there's some good reason to not do that. What's for sure expected is that it works reliably inside the zone of validity (and have shear lines at the limit or outside the domain of validity).

And in all cases, what would be great is to print out warnings when doing such transforms (either saying that an SRS transformation is used outside of it's nominal domain or that two different transformations where used and the potential consequences).

As @jacmendt (author of the referenced issue), we may be able to fund some bugfixing effort on this if it gets flagged as a bug.

Thanks for all the great work.

Steps to reproduce the issue

• gdalwarp -s_srs epsg:2056 -t_srs EPSG:3857 tile-downsampled.tif output.tif

tile-downsampled.zip (400ko)

Versions and provenance

OSGeo4W
Windows 11
GDAL 3.10.2, released 2025/02/11

Additional context

Note that the issue can be worked around in the same way as #10510

[rouault]

This is a bit tricky. The area of validity of EPSG:2056 and the CH1903+ -> WGS84 transformation is defined by EPSG in a unspecified datum, so depending of if you consider that the Helmert transformation must be applied or not, either the area of use is actually on the very left border of your dataset when you apply the Helmert transformation, or it is rather towards the right side of the image if you don't apply it (what GDAL / PROJ do, which explain the output)
This is a tricky topic, to know when to apply or not a transformation, and is something that is at the intersection of GDAL and PROJ.
In this particular case, as there's a single non-ballpark transformation, disabling the use of ballpark transformations forces PROJ to use the CH1903+ -> WGS84 transformation.

I've tested it successfully with the following patch:

diff --git a/ogr/ogrct.cpp b/ogr/ogrct.cpp
index f04c993677..a8122ed07d 100644
--- a/ogr/ogrct.cpp
+++ b/ogr/ogrct.cpp
@@ -1669,7 +1669,7 @@ int OGRProjCT::Initialize(const OGRSpatialReference *poSourceIn,
             if (options.d->dfAccuracy >= 0)
                 aosOptions.SetNameValue(
                     "ACCURACY", CPLSPrintf("%.17g", options.d->dfAccuracy));
-            if (!options.d->bAllowBallpark)
+            //if (!options.d->bAllowBallpark)
                 aosOptions.SetNameValue("ALLOW_BALLPARK", "NO");
 #if PROJ_VERSION_MAJOR > 9 ||                                                  \
     (PROJ_VERSION_MAJOR == 9 && PROJ_VERSION_MINOR >= 2)

A cleaner version would expose a GDAL warp option like -wo ALLOW_BALLPARK=NO to propagate this wish from the GDAL warp logic to the coordinate transformation (that is to cleanly set this bAllowBallpark flag to false). You can contact me if you wish me to implement that

[olivierdalang]

Oh that's interesting.

For the record here's the boundary as known to PROJ (indeed not the same as what's documented on https://epsg.io/2056), and indeed the seam is where it's expected.

Same case (worse) for our compatriots from Ticino:

If there's a way to interpret the EPSG definitinon that covers the whole of national territory, wouldn't it make sense then to also report/fix this in PROJ ?

Otherwise yes I think exposing the option would be nice, if possible with a warning in this scenario. I'll check with my organisation to see how we can move forward with that.

Thanks for the super quick and detailed feedback !

[jratike80]

Please use epsg.org if you need to know for sure how the data are in the EPSG registry. This area of validity seems to be defined by a polygon, not only by the corner points, but perhaps Proj is made to read only the bbox?

https://epsg.org/extent/gml/id/1286