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

[Roger Bivand]

On Thu, 7 Nov 2019, Even Rouault wrote:

> On jeudi 7 novembre 2019 11:10:26 CET Roger Bivand wrote:
>> On Thu, 7 Nov 2019, Even Rouault wrote:
>>> On jeudi 7 novembre 2019 15:24:44 CET Nyall Dawson wrote:
>>>> On Wed, 30 Oct 2019 at 06:16, Even Rouault <even.rouault at spatialys.com>
>>>
>>> wrote:
>>>>>> 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
>>>>>
>>>>> You get a ~1.5 metre shift to the north east as expected.
>>>>
>>>> This thread has been bouncing around in my mind for the last week.
>>>> Forgive me if I'm misinterpreting the situation, but wouldn't you
>>>> explicitly want legacy results NOT to be reproducible in cases where
>>>> proj 6's more accurate transformation capabilities have affected these
>>>> results?
>>>>
>>>> If a particular study isn't affected by relatively small shifts such
>>>> as the 1.5m shift in the GDA2020 situation, then the results should
>>>> remain reproducible even if this shift is corrected. But if a
>>>> particular study IS affected by distances of this magnitude then I'd
>>>> argue that it would be vital that the older proj behavior can't be
>>>> reproduced exactly, and accordingly the previous study/results are
>>>> invalidated as a direct result..!
>>>
>>> Are you saying that the new behaviour of ExportToProj4() shouldn't cause
>>> issues ? Well, that's more subtle than that. As +towgs84 terms might be
>>> omitted more often now, you can get differences of ~ 100 metres when
>>> dealing with datums dating before the GPS era. So when using GDAL 3
>>> ExportToProj4() you can effectively get "wrong" results in some cases.
>>>
>>> The action I took for now is to document the issue with ExportToProj4():
>>> https://gdal.org/doxygen/
>>> classOGRSpatialReference.html#a271b3de4caf844135b0c61e634860f2b
>>
>> The red-lined Warning is very explicit:
>>
>> "Use of this function is discouraged. Its behaviour in GDAL >= 3 / PROJ >=
>> 6 is significantly different from earlier versions. In particular +datum
>> will only encode WGS84, NAD27 and NAD83, and +towgs84/+nadgrids terms will
>> be missing most of the time. PROJ strings to encode CRS should be
>> considered as a a legacy solution. Using a AUTHORITY:CODE or WKT
>> representation is the recommended way."
>>
>> I think importFromProj4() is more forgiving - is this likely to continue,
>> or should we move away from Proj4 anyway?
>
> importFromProj4() is likely to continue to exist in the foreseable future.
> There are various raster/vector formats that use PROJ.4 strings to encode the
> CRS, so we want to still be able to decode them.
> What is discouraged is the generation of new PROJ.4 CRS strings from now.

Thanks for the clarification. So, as per:

https://proj.org/usage/quickstart.html

and in the development version of R's rgdal package, I can populate a WKT2 
from, for example:

> cat(showSRID("+proj=merc +lat_ts=56.5 +ellps=GRS80", format="WKT2", 
+ multiline="YES"), "\n")
PROJCRS["unknown",
     BASEGEOGCRS["unknown",
         DATUM["Unknown based on GRS80 ellipsoid",
             ELLIPSOID["GRS 1980",6378137,298.257222101,
                 LENGTHUNIT["metre",1],
                 ID["EPSG",7019]]],
         PRIMEM["Greenwich",0,
             ANGLEUNIT["degree",0.0174532925199433],
             ID["EPSG",8901]]],
     CONVERSION["unknown",
         METHOD["Mercator (variant B)",
             ID["EPSG",9805]],
         PARAMETER["Latitude of 1st standard parallel",56.5,
             ANGLEUNIT["degree",0.0174532925199433],
             ID["EPSG",8823]],
         PARAMETER["Longitude of natural origin",0,
             ANGLEUNIT["degree",0.0174532925199433],
             ID["EPSG",8802]],
         PARAMETER["False easting",0,
             LENGTHUNIT["metre",1],
             ID["EPSG",8806]],
         PARAMETER["False northing",0,
             LENGTHUNIT["metre",1],
             ID["EPSG",8807]]],
     CS[Cartesian,2],
         AXIS["(E)",east,
             ORDER[1],
             LENGTHUNIT["metre",1,
                 ID["EPSG",9001]]],
         AXIS["(N)",north,
             ORDER[2],
             LENGTHUNIT["metre",1,
                 ID["EPSG",9001]]]]

So for example, Mike Sumner should be able to continue to use proj4 
strings to instantiate SRS in a shorthand way provided they can be handled 
by importFromProj4(), and we can expect this usage to continue to work 
because, as you say, data sets themselves use Proj4 strings, and this path 
must remain open.

Roger


>
> Even
>
>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: Roger.Bivand at nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en