PROJ 6 drops +datum= specification; transform degraded in workflows

[rsbivand]

This is linked to r-spatial/discuss#28.

This zip contains a simple reprex, the Broad Street pump location in a shapefile used in ASDAR (https://asdar-book.org/), includes a +datum= specification. There are also two R scripts each for sf and rgdal, one pair just using +towgs84=, the other pair using +nadgrids=. The first pair were run with PROJ 4.9.3/GDAL 2.4.2 and PROJ 6.2.0/GDAL 3.0.1 (output Rout files with *4* and *6*), the grid*.R files only under PROJ 6.2.0/GDAL 3.0.1 to use the new PROJ search path extension mechanism tested in test.cpp.

reprex.zip

Basic findings: for PROJ 4.9.3/GDAL 2.4.2, the +datum= spec gives a reasonable transformation of the input to EPSG:4326. For PROJ 6.2.0/GDAL 3.0.1, the +datum= spec is dropped on reading in sf and rgdal, and variants of +towgs84= specs need to be added manually to the CRS to get a reasonable transformation (not off by ~100m).

Side finding: writing the object with a +towgs84= with the ESRI Shapefile driver in both sf and rgdal drops the +towgs84=, but preserves it for the GPKG driver.

Only for PROJ 6: using the NTv2 grid in the Europe datumgrid zip archive indicated by projinfo gives a good transformation with both sf and rgdal using a CRS modified by adding +nadgrids= and extending the search path to include the current directory (where the *.gsb file is located, not included in zipfile, see projinfo output for download link).

Conclusion: Moving to PROJ 6 in software is feasible, but typical workflow input may suffer invisible degradation during processing for rgdal and sf workflows and their reverse dependencies. @tim-salabim (re: mapview) @mtennekes (re: tmap) @Nowosad (re. spData).

[Nowosad]

@rsbivand - thank you for your work in this regard. Do you have any advice on how to proceed with those changes? How can I manually check which dataset in spData could be problematic, and then how to decide on the correct +towgs84 values?

[rsbivand]

I'm thinking that sf and rgdal and any other packages with transform functions may need to include a test on startup to see whether current PROJ diverges from PROJ 4.9.3. So far, CRAN Windows binaries are 4., so the effects are seen only where the packages are installed from source against PROJ 6.. If a package only offers files, like spData, the consequences will only occur once the objects have been read in, so *.rda and native GIS format behaviour will probably differ (the *.rda freeze the CRS representation, but we don't know how a CRS with a +datum= string will be handled - yet). One way is to mapview::mapview or tmap in view mode in the two settings and check visually; another to add the equivalent transformation to EPSG:4326 and check for equivalence of a sample point across the two settings. The QGIS approach is to ask the user to choose the transformation path manually, which is not appropriate in workflows and scripts. For geocompr, I guess adding (hidden) code to trap changes is feasible, say on ch 6. Maybe a blog / addendum to ch 6 might be a goal?

[mtennekes]

Trying to link the concepts ellipsoid/spheroid, datum, grid. As far I as understand it:

• the spheroid (+ellps) defines the global shape of the Earth with the distances from the center to the pole and from the center to the equator.
• the datum (+datum) is then just a modification (i.e. transformation, rotation and scaling) of this spheroid
• in PROJ.4, the +towgs84 tag is always needed, since transformations between CRS's always go via wgs84?
• What is a grid exactly? What I can imagine, is that it is regular grid of lat long coordinates, with as values the difference in elevation between the surface defined by the spheroid. What is the granularity of these grids? I haven't had time to download them myself.
• Last week in Muester, I also thought about combining the two different worlds: one of the projections/datums, and one of the grid systems like Google S2 and Uber H3. Maybe I am too naive, but it sounds like when combining one of these grid systems with the delta-elevation grid I had in mind in my previous point, we don't need to have all these different complex projections. Depending on the level of detail that is needed for the application, different versions of the delta-elevation grid are required.

[rsbivand]

@edzer teaches using Datums and Map Projections by Iliffe and Lott https://www.whittlespublishing.com/Datums_and_Map_Projections. Maybe we need his course online?

[edzer]

+towgs84 gives 3 or 7 parameters for a parametric shift, datum grids give the shift for every location up to some resolution, but can then be interpolated. Reading the book Roger mentioned will make you realize that datum transformation will always be there, and will never be simple.

[rsbivand]

I'm begining to suspect that the GDAL3/PROJ6 combination is unkind to the shapefile driver. I'll try later to add a shapefile without a +datum= but with an equivalent +towgs84= copied in manually from GDAL 2.4.2 data/datum:shift.csv (446.448 | -125.157 | 542.06 | 0.15 | 0.247 | 0.842 | -20.489 I think). We'll need to check a range of drivers to see which handle input CRS differently too.

[edzer]

I just submitted sf 0.8-0 to CRAN addressing the PROJ 6.2.0 issues seen on debian testing, and those caused by the changed tidyr 1.0-0 API.

[Nowosad]

Dear @rsbivand, @edzer and the rest of #rspatial developers,
I still try to wrap my head around the future changes related to the PROJ library (PROJ6 and after). What is the future of the proj4string definitions? Are there plans to (1) keep it, (2) drop it from #rspatial, (3) add some other way of defining projections (e.g., WKT), or (4) do something different?
(https://gdal.org/api/ogr_srs_api.html#_CPPv416OSRExportToProj420OGRSpatialReferenceHPPc).

[edzer]

Good question - what do you suggest?

[rsbivand]

As of now, rgdal on R-Forge (1.5.1, very unstable, rev 862) when built with GDAL >= 3 and PROJ >= 6 will issue warnings when a DATUM is in the input vector SRS but discarded in the PROJ string. Next step to extend this to raster (see http://fwarmerdam.blogspot.com/2010/03/in-last-few-weeks-i-believe-i-have-made.html).

See thread on gdal-dev.

New list_coordOps() takes two CRS strings and returns coordinate operation pipelines. The easiest CRS strings are such as "EPSG:4326". For sf, users are going to have to insert EPSG codes, but there is somewhere to put them in a "crs" object.

In sp and dependencies, "CRS" is S4, so no free slot. I suggest the same kludge as for rings and holes, that is to add a comment() to "CRS" objects to contain the CRS string now needed. This needs to update the "projargs" slot, and needs checking too. Updating it would largely be manual, and users would have to be alerted to the need to do so.

The CRS string could also be added to sf "crs" objects, checking for NULL (as with a NULL comment), to provide somewhere to put WKT, URN, PROJJSON or other string representations. There still isn't clarity about where to put the EPOCH if known.

Feedback on rgdal very welcome. Most ideas should propagate to sf without extra difficulty. Comment now if this matters for you, @mdsumner made useful comments on the gdal-dev thread.

Notes on r-spatial/discuss#28 about packages depending on sf and rgdal for coordinate operations entered last week.

[rsbivand]

Maybe also extend proj string degradation testing to towgs84 and nadgrids; datum was the most obvious problem, but towgs84 is also checked in the current code, and the warning qualified if a towgs84 makes it through OK (helps with -DACCEPT_USE_OF_DEPRECATED_PROJ_API_H pj_transform()), but won't help with PROJ 6 transformation I think (these need +type=crs too anyway).

[rsbivand]

@Nowosad : could you see whether current rgdal tickles any bugs with the spData vector files?

[edzer]

Link to the gdal-dev thread.