[gdal-dev] GDAL 3, PROJ 6.2 exportToProj4() drops Datum silently

[Even Rouault]

> The underlying problem is that users may face silent degradation in
> workflows as they upgrade sf and rgdal to PROJ >= 6 and GDAL >= 3 if they
> are not prompted to take remedial steps themselves. I realise that users
> have become too comfortable, but user guidance will be essential to avoid
> too much noise on SO, twitter, etc.

I do understand your issues. If contributors are willing to spend time 
bringing back more backward compatibility into PROJ exportToProjString() (for 
the +datum part I mean. I wouldn't be enthousiastic of attempting to recover 
full backward compatibility on the +towgs84 part), we'd probably accept 
reasonable patches . I consciously did a partial implementation on that part 
of the code, so that people have an incentive to embrace other options, more 
future proof, like WKT or AUTHORITY:CODE, to encode their CRS.

> 
> > If you use directly the WKT representation from this .prj file with PROJ
> > functions that deal with coordinate transformations (or the
> > OGRSpatialReference object with GDAL functions), things should work
> > smoothly.
> 
> WKT is not an option really, though users may wish to look at the WKT2 of
> an EPSG code. My attempts to feed WKT into importFromWkt() have failed so
> far, but that could be string handling (I'm treating the whole WKT block
> as a single text string).

Using WKT should definitely be considered as an option, and investigating why 
that fails is in my opinion a more valuable time investment rather than trying 
to revive a ghost from the past. WKT is a standard, PROJ strings are not.

> I realise this too, but to permit legacy results to be replicated,

Some legacy behaviour was just plain wrong/inaccurate. Trying to replicate it 
in PROJ would be heartbreaking, an extra maintaince burden, and a confusion 
for PROJ users because depending on the path they take, they would get 
different results.

Not-so-random example with GDAL 2.4/PROJ 5 with GDA94 to GDA2020

$ gdalsrsinfo EPSG:4283

PROJ.4 : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 

OGC WKT :
GEOGCS["GDA94",
    DATUM["Geocentric_Datum_of_Australia_1994",
        SPHEROID["GRS 1980",6378137,298.257222101,
            AUTHORITY["EPSG","7019"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6283"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4283"]]

$ gdalsrsinfo EPSG:7844

PROJ.4 : +proj=longlat +ellps=GRS80 +no_defs 

OGC WKT :
GEOGCS["GDA2020",
    DATUM["Geocentric_Datum_of_Australia_2020",
        SPHEROID["GRS 1980",6378137,298.257222101,
            AUTHORITY["EPSG","7019"]],
        AUTHORITY["EPSG","1168"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","7844"]]

Transforming between the above using PROJ string results in a null 
transformation:

$ echo "150 -30" | gdaltransform -s_srs EPSG:4283 -t_srs EPSG:7844 --debug on
OGRCT: Source: +proj=longlat +ellps=GRS80 +no_defs
OGRCT: Target: +proj=longlat +ellps=GRS80 +no_defs
[...]
150 -30 0


Same transformation with GDAL + PROJ master:

$ echo "150 -30" | gdaltransform -s_srs EPSG:4283 -t_srs EPSG:7844 --debug on
OGRCT: proj=pipeline step proj=axisswap order=2,1 step proj=unitconvert 
xy_in=deg xy_out=rad step proj=push v_3 step proj=cart ellps=GRS80 step 
proj=helmert x=0.06155 y=-0.01087 z=-0.04019 rx=-0.0394924 ry=-0.0327221 
rz=-0.0328979 s=-0.009994 convention=coordinate_frame step inv proj=cart 
ellps=GRS80 step proj=pop v_3 step proj=unitconvert xy_in=rad xy_out=deg step 
proj=axisswap order=2,1 (GDA94 to GDA2020 (1))
[...]
150.000006074915 -29.9999871787747 0

You get a ~1.5 metre shift to the north east as expected.

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com