[gdal-dev] Dealing with "default" TOWGS84 values

[Even Rouault]

Hi,

I wanted to get opinions on a topic related to TOWGS84 7-parameter transforms. In the GDAL 
1.x/2.x & PROJ 4.x/5.x era, having TOWGS84 values attached to CRS definitions was critical to 
get correct datum transformation.

Now, in the GDAL 3 + PROJ 6.x/7.x era, there are more an annoyance than anything else, 
since they can prevent more accurate transformations to be used.

Some data formats such as GeoTIFF files written by GDAL 1.x/2.x encode the TOWGS84 value 
of the CRS definition that was built when importing the EPSG dataset. E.g 
TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489] for EPSG:27700 OSGB 1936 / 
British National Grid

In GDAL 3, when such file is read, a BoundCRS object is returned, that is a CRS that has a base 
CRS (potentially with a EPSG code attached. e.g EPSG:27700) + the definition of its 
transformation to WGS84 (which corresponds to the transformation with larger area of use). 
Software like QGIS will then use this BoundCRS (since that's what GDAL returns for the 
dataset CRS) when reprojecting to a target CRS, which will cause the TOWGS84 parameters 
from the file to be used (PROJ >= 6 will honour the TOWGS84 parameters of a BoundCRS 
when transforming from/into it), rather than potentially more accurate transformations 
(typically grid based transformations) from the base CRS to the target CRS.

Whether this behaviour is a bug or a feature, and where this bug belongs to, has been a 
debate in a number of tickets:
https://github.com/OSGeo/gdal/issues/2219
https://github.com/qgis/QGIS/issues/34993

I wanted to have broader opinions regarding if we should do something about that, and 
what.

One possibility would be for GDAL to ignore the TOWGS84 value associated with a base CRS 
if we detect that this TOWGS84 value is the default one that was written by GDAL 1.x/2.x for 
this base CRS, and thus to return just the base CRS, letting PROJ figure out all potential 
transformations when transforming to a target CRS. This could potentially be the default 
behaviour, with a configuration option that could be set to opt for the current behaviour 
(that is return the BoundCRS). This would have to be implemented per driver (logic would be 
similar). I guess GeoTIFF, GPKG, Spatialite could be candidates.
The potential downside of this is if the user relied on those exact default TOWGS84 values to 
be used when reprojecting to WGS84 (or when using WGS84 as the pivot)

I'd note that GDAL 3, in the OGRCoordinateTransformation code, implements partially this 
logic, but only in the transformation code (not impacting the CRS actually returned by the 
dataset), and with limitations. That is when transforming from/into a BoundCRS whose which 
has a unique Helmert transformation to WGS84 that matches the transformation of the 
BoundCRS, it uses by default only the BaseCRS (can be disabled with 
OSR_CT_USE_DEFAULT_EPSG_TOWGS84=YES). But this is for example not sufficient for 
EPSG:27700, because it has several Helmert transformations. So as suggested in https://
github.com/OSGeo/gdal/issues/2219, a broader fix would be for the PROJ database to have a 
GDAL2 authority that would store the TOWGS84 parameters historically used by GDAL 2.

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20200331/ffe64abd/attachment.html>