PROJ 6 CRS handling changes impact workflows using +datum= (and other fun)
rsbivand opened this issue · 58 comments
Ongoing changes in the representation of coordinate reference systems in the PROJ, external software used by sf, sp and other R packages, will impact R packages and workflows using these packages. In particular, the changes affect transformation between coordinate reference systems, as for example shown in this presentation at the OpenGeoHub Summer School 2019.
-
Invite users/developers to contribute cases
-
Add other examples and links, also from other geospatial software.
GDAL RFC 73: https://gdal.org/development/rfc/rfc73_proj6_wkt2_srsbarn.html
Notes from Spatialite
PROJ wiki proj.h adoption status
GRASS-dev thread and PR
cartography CRS issues
geojsonio geojsonio seems vulnerable to PROJ 6
MODIS getTile() sp subset or raster::crop issues with PROJ 6.1 and GDAL 3.0
MODIStsp PROJ 6 failures in MODIStsp
rangeMapper PROJ 6 failures in rangeMapper
rdefra PROJ 6 failures in rdefra
rnrfa PROJ 6 failures in rnfra
rshapemapper PROJ 6 problems?
Incomplete list of packages with CRAN issues after PROJ 6.2.0 installed on some test machines:
https://cran.r-project.org/web/checks/check_results_foieGras.html
https://cran.r-project.org/web/checks/check_results_plotKML.html
https://cran.r-project.org/web/checks/check_results_rangeMapper.html
https://cran.r-project.org/web/checks/check_results_rnrfa.html
https://cran.r-project.org/web/checks/check_results_sf.html
https://cran.r-project.org/web/checks/check_results_vapour.html
- Provide recommendations for checking/testing PROJ versions in use and expected transformation outcomes.
See also r-spatial/sf#1146 for sf and rgdal reprex with PROJ 4.9.3/GDAL 2.4.2 output from the same file contrasted with PROJ 6.2.0/GDAL 3.0.1.
This issue is based on https://github.com/rsbivand/geostat19_talk; https://rsbivand.github.io/geostat19_talk/Discuss_issue_drafts.html
Relevant CDN for serving grid files discussion on https://twitter.com/howardbutler/status/1171778886646022145?s=20. On the client side related to thread: https://lists.osgeo.org/pipermail/proj/2019-September/008825.html.
Adding files (scripts, output) map_package_analyses.zip for packages with reverse dependencies (depends, imports, suggests) on rgdal using pkgapi map_package() to try to see which packages use spTransform(), project(), and raster projectExtent() or projectRaster(). About 80 packages of almost 300 in total seem to use projection mechanisms on the sp side. The scripts may be used to set up similar reports on the sf side. Other columns report CRS() and other vulnerable points.
See also suggestions on https://github.com/pyproj4/pyproj; https://pyproj4.github.io/pyproj/stable/index.html.
Hi Roger, much appreciate your reaching out and pointing to all the resources.
However, I feel a bit overwhelmed by it all and am struggling to turn this into action for how to get my package (geohashTools) working with PROJ>=6.
I see this entry in map_package_analyses/df_all.csv for my package:
raster_projectExtent raster_projectRaster rgdal_checkCRSargs rgdal_CRSargs rgdal_make_EPSG
1: 0 0 0 0 0
rgdal_project rgdal_spTransform sp_CRS sp_proj4string sp_spTransform sum
1: 0 0 1 2 1 4
How can I use this to proceed?
Also, is it safe to use rocker/geospatial for testing (i.e., does that image have the proper system libraries in place already)?
Seems rocker/geospatial is still using PROJ 4.9.3, can you recommend any Docker image to use for testing?
@MichaelChirico for the purpose of testing PROJ6 in the Geocomputation with R book, I've created a new docker image based on rocker:
FROM rocker/rstudio:3.6.1
RUN apt-get update \
&& apt-get install -y --no-install-recommends \
gdb \
git \
libcairo2-dev \
libcurl4-openssl-dev \
libexpat1-dev \
libpq-dev \
libsqlite3-dev \
libudunits2-dev \
make \
pandoc \
qpdf \
r-base-dev \
sqlite3 \
subversion \
valgrind \
vim \
wget
RUN apt-get install -y --no-install-recommends \
libv8-3.14-dev \
libjq-dev \
libprotobuf-dev \
libxml2-dev \
libprotobuf-dev \
protobuf-compiler \
unixodbc-dev \
libssh2-1-dev \
libgit2-dev \
libnetcdf-dev \
locales \
libssl-dev
RUN locale-gen en_US.UTF-8
ENV PROJ_VERSION=6.2.0
# ENV LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
RUN wget http://download.osgeo.org/proj/proj-${PROJ_VERSION}.tar.gz \
&& tar zxf proj-*tar.gz \
&& cd proj* \
&& ./configure \
&& make \
&& make install \
&& cd .. \
&& ldconfig
# install proj-datumgrid:
RUN cd /usr/local/share/proj \
&& wget http://download.osgeo.org/proj/proj-datumgrid-latest.zip \
&& unzip -o proj-datumgrid*zip \
&& rm proj-datumgrid*zip \
&& wget https://download.osgeo.org/proj/proj-datumgrid-europe-latest.zip \
&& unzip -o proj-datumgrid*zip \
&& rm proj-datumgrid*zip \
&& cd -
# GDAL:
ENV GDAL_VERSION=3.0.1
ENV GDAL_VERSION_NAME=3.0.1
RUN wget http://download.osgeo.org/gdal/${GDAL_VERSION}/gdal-${GDAL_VERSION_NAME}.tar.gz \
&& tar -xf gdal-${GDAL_VERSION_NAME}.tar.gz \
&& cd gdal* \
&& ./configure \
&& make \
&& make install \
&& cd .. \
&& ldconfig
#RUN git clone --depth 1 https://github.com/OSGeo/gdal.git
#RUN cd gdal/gdal \
# && ls -l \
# && ./configure \
# && make \
# && make install \
# && cd .. \
# && ldconfig
# GEOS:
ENV GEOS_VERSION=3.7.2
#
RUN wget http://download.osgeo.org/geos/geos-${GEOS_VERSION}.tar.bz2 \
&& bzip2 -d geos-*bz2 \
&& tar xf geos*tar \
&& cd geos* \
&& ./configure \
&& make \
&& make install \
&& cd .. \
&& ldconfig
RUN install2.r --error \
remotes
RUN R -e "remotes::install_github('geocompr/geocompkg')"
@Nowosad great! Is this hosted on dockerhub?
It was not... until today;)
There are two versions:
- https://hub.docker.com/r/jakubnowosad/rspatial_proj6 - with all external libraries (PROJ6, GDAL3, GEOS)
- https://hub.docker.com/r/jakubnowosad/geocompr_proj6 - the same as above + many R packages for geocomputation (https://github.com/geocompr/geocompkg/blob/master/DESCRIPTION)
This is the WIP-vignette from R-Forge/rgdal 1.5-1 rev 888. Comments welcome!
PROJ6_GDAL3.zip
DO make the WKT representations readable using cat and pretty, as in
> library(sf)
Linking to GEOS 3.7.1, GDAL 2.4.2, PROJ 5.2.0
> st_as_text(st_crs(4326))
[1] "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]"
> st_as_text(st_crs(4326), pretty = TRUE)
[1] "GEOGCS[\"WGS 84\",\n DATUM[\"WGS_1984\",\n SPHEROID[\"WGS 84\",6378137,298.257223563,\n AUTHORITY[\"EPSG\",\"7030\"]],\n AUTHORITY[\"EPSG\",\"6326\"]],\n PRIMEM[\"Greenwich\",0,\n AUTHORITY[\"EPSG\",\"8901\"]],\n UNIT[\"degree\",0.0174532925199433,\n AUTHORITY[\"EPSG\",\"9122\"]],\n AUTHORITY[\"EPSG\",\"4326\"]]"
> cat(st_as_text(st_crs(4326), pretty = TRUE), "\n")
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]]
>
(misses the final newline, but nevertheless)
All I need is an indication if someone else is attacking this problem and can see obviously better resolutions, or indications (thanks @Nowosad !) of checks on the current state of play. I do not need comments about aesthetics, they do not matter for the task in hand. I'll wrap the output, but I will not change to multiline out from exportToWkt(), because the idea is to present a string that is machine readable. Unless it can be shown that WKT needs to be human-readable (worse editable), I'll keep then as single line strings.
Just discovered this thread, great to see proactive work to preempt proj6. Just testing my package stplanr with the following command (in case of use/interest to others):
docker run -d -p 8789:8787 -e DISABLE_AUTH=TRUE -v $(pwd):/home/rstudio/ jakubnowosad/geocompr_proj6
what are we to use for testing? github.com/rsbivand/sp@master and rgdal master on r-forge?
Correct?
#remotes::install_github("rsbivand/sp")
#install.packages("rgdal", repos="http://R-Forge.R-project.org")👍 that's what I did, as documented here: ropensci/stplanr#364 (comment)
cool thanks, I'm failing to install rgdal with PROJ 6.1.0 atm but I consider that tangential, got the geocompr_proj6 running now
Thanks, yes, so far I'm running PROJ 6.2.1 and GDAL 3.0.2 for all this,, but will have to see which functions were in 6.1 and 6.0, and modify the code/configure to match.
One thing I notice that confused me, with PROJ 5.2.0 is showSRID() returning WKT2 rather than proj4 or error:
## PROJ 5.2.0
showSRID("+init=epsg:4326", format="PROJ", multiline="NO")
[1]"GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295]]"causing problems downstream with WKT2 in the slot:
Loading required package: sp
rgdal: version: 1.5-2, (SVN revision 897)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 2.4.2, released 2019/06/28
Path to GDAL shared files: /usr/share/gdal
GDAL binary built with GEOS: TRUE
Loaded PROJ.4 runtime: Rel. 5.2.0, September 15th, 2018, [PJ_VERSION: 520]
Path to PROJ.4 shared files: (autodetected)
Linking to sp version: 1.3-3
> checkCRSArgs_ng("+init=epsg:4326")
[[1]]
[1] TRUE
[[2]]
[1] "GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295]]"
[[3]]
[1] "GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.017453292519943295]]"On 6.2.1 all is well (if I understand correctly, I'm still catching up):
showSRID("+init=epsg:4326", format="PROJ", multiline="NO")
[1] "+proj=longlat +datum=WGS84 +no_defs"Is it fair to say that that showSRID() should allow only 1) PROJ and return the input string if in_format == 1L, otherwise error - and 2) WKT1 then the current showWKT() behaviour ?
Well caught! I think that if PROJ < 6 && GDAL < 3, it should revert to showWKT() In that setting, I don't think comment() is added to "CRS" objects. The need for immediate adaptation because CRS degradation cuts in only for PROJ >=6 && GDAL >=3 but there manifests the silent degradation I was so concerned about in Münster at the OpenGeoHub summer school (the video is terribly rambling, start at 37:14), but I could not then see a way to prevent hundreds of packages depending on sp and rgdal suffering major regression when upstream providers moved to GDAL3+PROJ6.
I'm afraid I'm ill at home, so it is harder to use my office desktop to install PROJ 4 + GDAL 2, PROJ 5 + GDAL 2, and the banned PROJ 6 + GDAL 2 (should not work, GDAL 3 requires PROJ 6, but released GDAL 2 couldn't upward-bound PROJ), to test for frailties. I know that some PROJ 6 functionality is only present >= 6.1 or >= 6.2, I was just struggling to clear a path through the undergrowth, so I'm very grateful for reports like this!
Great thanks, I was pretty lost before but this was good to get my head around it.
I have posted on http://spatial.nhh.no/R/rgdal/PROJ6GDAL3/ the *.Rcheck directories for rgdal reverse dependencies failing under rgdal 1.5-1 and PROJ 6.2.1/GDAL 3.0.2. The NOTES file indicates actions taken (and provides the generating script and failure history for earlier releases). I'm posting this just as documentation of the current status, and in case anyone needs more details than I sent or pasted into issues. The whole 212MB bundle is http://spatial.nhh.no/R/rgdal/PROJ6GDAL3/rgdal_1.5-1_revdeps.zip.
Thanks for the docker: really useful. However, to be able to test my package I would need to add HDF4 support to GDAL, which seems not to be available in the GDAL3 provided by the container
To do that, do I only need to modify the dockerfile by adding --with-hdf4 to configure, like in:
RUN wget http://download.osgeo.org/gdal/${GDAL_VERSION}/gdal-${GDAL_VERSION_NAME}.tar.gz \
&& tar -xf gdal-${GDAL_VERSION_NAME}.tar.gz \
&& cd gdal* \
&& ./configure --with-hdf4\
&& make \
&& make install \
&& cd .. \
&& ldconfig
or should I do something else/more? (sorry, if it is a silly question, but I do not have much experience with docker/building GDAL from source)
thanks in advance,
Lorenzo
@lbusett see https://github.com/rocker-org/geospatial/blob/master/Dockerfile for a more complete install, that includes the hdf4 libs.
@lbusett please create a PR at https://github.com/Nowosad/rspatial_proj6
hdf4 is detected automatically by ./configure.
I finally had the time to test the new docker, and the rspatial_proj6 docker now have proper HDF4 support. However, unless I am doing something wrong on my side (I simply did a docker pull on both images), the geocompr_proj6 one seems to be still lacking the support (I guess the image was not "rebuilt" after changes in rspatial_proj6 ? ). Not a big deal, but I think it could be useful to "update" also that. Is that possible ?
The geocompr_proj6 was outdated (I have no idea how to make an automatic connection between rspatial_proj6 and geocompr_proj6). It should be good now - please check.
Pulled it now - seems to be ok, thanks!
I got soooo much value out of this docker help, thanks!!!
This will not affect the CRS problem. Some GDAL drivers behave differently if GDAL is built with GEOS, but this affects the geometries, not the coordinate reference system representation.
Yes, this is "correct" (though not required).
Ok, thanks to both. By the way: while testing we are currently seeing some weird problems in reprojecting from EPSG::4326 - apparently related to axis swapping. We'll do some additional testing to and report back.
Finally, the first axis swapping hell case in the wild - looking forward to a reprex!
The reprex will be what we need as a Christmas present (?)! Keep us busy right into 2020 ... done by St Louis?
This transect defines the perfect example sites cbind(-90:90, -90:90)
(Just off the West Antarctic peninsula is my favourite, always stuff at 61W,61S. )
No, because as far as I've seen so far, it didn't need it on input or output, and this is addressed in list_coordOps(..., visualization_order=TRUE) and SEXP list_coordinate_ops() for listing coordinate operations, calling proj_normalize_for_visualization() from PROJ if TRUE.
Ah, yes, I forgot that rgdal does not use the GDAL interface to reproject.
Dear @rsbivand ,
thank you for managing the transition to PROJ 6 and coordinating issues affecting R packages.
We noticed an issue in spTransform() when the CRS of the input and/or the output vector is associated to a WKT2 in which the axes order is Y – X instead than X – Y (I saw this note in this regard). It is the issue @lbusett anticipated some hours ago.
Since the issue does not occur using standalone GDAL binary (ogr2ogr), I thought to report it here.
You can find a report of this problem here.
Hoping this can help (and apologising if it was already known),
Thanks @ranghetti for the thorough analysis!
So, in the end the problem seems to be appearing only on spTransform on rgdal-dev, while in st_transform on wkt2 everything seems ok, is that right?
Thanks @ranghetti for the thorough analysis!
So, in the end the problem seems to be appearing only onspTransformonrgdal-dev, while inst_transformonwkt2everything seems ok, is that right?
Yes, as well as using the master branch and the CRAN version of sf (but in this last case, if I understood correctly, the reason is that the PROJ.4 is still used).
I’d say what is somehow particularly worrysome / counterintuive is that, starting for example
from a point in a projected system with X-Y axis order:
library(sp)
library(rgdal)
#> rgdal: version: 1.5-2, (SVN revision (unknown))
#> Geospatial Data Abstraction Library extensions to R successfully loaded
#> Loaded GDAL runtime: GDAL 3.0.2, released 2019/10/28
#> Path to GDAL shared files: /usr/local/share/gdal
#> GDAL binary built with GEOS: FALSE
#> Loaded PROJ.4 runtime: Rel. 6.2.1, November 1st, 2019, [PJ_VERSION: 621]
#> Path to PROJ.4 shared files: /usr/local/share/proj
#> Linking to sp version: 1.3-3
library(magrittr)
mypoint <- sp::SpatialPoints(cbind(300000, 400000), CRS(SRS_string = "EPSG:32632"))
coordinates(mypoint)
#> coords.x1 coords.x2
#> [1,] 3e+05 4e+05, the following instructions yield currently different (opposite) results:
coordinates(sp::spTransform(mypoint, CRS("+init=epsg:4326")))
#> coords.x1 coords.x2
#> [1,] 7.199383 3.617083
coordinates(sp::spTransform(mypoint, CRS(SRS_string = "EPSG:4326")))
#> coords.x1 coords.x2
#> [1,] 3.617083 7.199383due probably to the fact that the output WKTS instantiated are different.
This gives X-Y Axis order
cat(showSRID(comment(CRS("+init=epsg:4326")), multiline = "YES"))
#> GEOGCRS["WGS 84",
#> DATUM["World Geodetic System 1984",
#> ELLIPSOID["WGS 84",6378137,298.257223563,
#> LENGTHUNIT["metre",1]],
#> ID["EPSG",6326]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8901]],
#> CS[ellipsoidal,2],
#> AXIS["longitude",east,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]],
#> AXIS["latitude",north,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]],
#> USAGE[
#> SCOPE["unknown"],
#> AREA["World"],
#> BBOX[-90,-180,90,180]]]but this gives Y-X Axis order
cat(showSRID(comment(CRS(SRS_string = "EPSG:4326")), multiline = "YES"))
#> GEOGCRS["WGS 84",
#> DATUM["World Geodetic System 1984",
#> ELLIPSOID["WGS 84",6378137,298.257223563,
#> LENGTHUNIT["metre",1]]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433]],
#> CS[ellipsoidal,2],
#> AXIS["geodetic latitude (Lat)",north,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> AXIS["geodetic longitude (Lon)",east,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> USAGE[
#> SCOPE["unknown"],
#> AREA["World"],
#> BBOX[-90,-180,90,180]],
#> ID["EPSG",4326]]so, it appears to me that it is not so much that the “reprojection” is wrong, but that when the target SRS is Y-X, then the coordinates of the reprojected point are saved/shown/should be interpretes as Lat/Lon order instead than in Lon/Lat order as we are accustomed (so, in the first transformation, coords.x1 is Longitude, in the second it is Latitude) . According to @ranghetti analysis, this appears to happen every time we go from a XY to a YX proj or viceversa.
Indeed, making a round trip works in both cases:
# UTM --> 4325XY --> UTM
coordinates(mypoint)
#> coords.x1 coords.x2
#> [1,] 3e+05 4e+05
sp::spTransform(mypoint, CRS(SRS_string = "+init=epsg:4326")) %>%
sp::spTransform(CRS(SRS_string = "EPSG:32632")) %>%
sp::coordinates()
#> coords.x1 coords.x2
#> [1,] 3e+05 4e+05
# UTM --> 4325YX --> UTM
sp::spTransform(mypoint, CRS(SRS_string = "EPSG:4326")) %>%
sp::spTransform(CRS(SRS_string = "EPSG:32632")) %>%
sp::coordinates()
#> coords.x1 coords.x2
#> [1,] 3e+05 4e+05This “swap” seems however troublesome to me, because I think that then, for example, all plotting routines would need to be made “Axis-order aware”, there (probably) could be troubles for sp retrocompatibility, and also users accustomed to a XY representation could be taken by surprise by the change.
Created on 2019-12-06 by the reprex package (v0.3.0)
It is good to know that this swapping doesn't happen when a CRS is initialised with a PROJ string. As this is up to now the only acceptable way, this problem shouldn't affect any r-spatial work done so far.
Inserting, in rgdal, pkg/src/proj_info6.cpp, line 358, the line
pj_transform = proj_normalize_for_visualization(ctx, pj_transform); // EJP
gives
coordinates(sp::spTransform(mypoint, CRS(SRS_string = "EPSG:4326")))
# coords.x1 coords.x2
# [1,] 7.199383 3.617083but the more important discussion is whether we want this, or want to be OGC compliant. As of: what is wrong with people wanting to do this to themselves. Or default to one and allow for the other?
Yes, that is the problem. How can we adopt a default (visualization axis order) but still permit users who need OGC compliant output to achieve that?
One option would be a global setting. It feels more elegant though to keep data in authority (wkt2/ogc) order, and do an axis swap on the fly when it is needed, e.g.
- before plotting
- before (and after) doing geometry operations
I will run some experiments in sf with setting SetAxisMappingStrategy(OAMS_AUTHORITY_COMPLIANT); by default.
Mmm for that to work, datasets like nc would need to be rewritten to have lat long coordinates first, or undergo an axis swap at read time. But based on what? Then, the bounding box logic breaks (what was xmin?), but after that projecting works fine! Feels like a problem that will not easily go away.
Perhaps restrict sp workflows to be GIS/visualization mode, because that is how the classes were designed? We can't rewrite legacy data sets that themselves may not encode/report axis order because we need to keep things running for both GDAL2/PROJ5 and GDAL3/PROJ6 at the same time (think CentOs).
Sounds fair enough. R CMD check sf_*gz has btw a test where proj_normalize_for_visualization as used above returned NULL, which I didn't understand, but which caused a segfault when left unchecked.
In case others want to experiment with sf and handling objects with authority compliant axes order, see r-spatial/sf#1146 : branch wkt2 in the sf github repo.
Video of talk on this problem on: https://www.youtube.com/playlist?list=PLXUoTpMa_9s10NVk4dBQljNOaOXAOhcE0, third in playlist; presentation at https://rsbivand.github.io/ECS530_h19/ECS530_III.html.
@ranghetti , @lbusett : please see rgdal SVN revision >= 903 (905 current) and scroll down in http://rgdal.r-forge.r-project.org/articles/PROJ6_GDAL3.html to Axis swapping to see a version of the Rpubs round-trip issue (when the served version updates). I've added code to address this, but using list_coordOps() already permitted the creation of a transformation object that would output visualization order objects. In PROJ, it is the transformation object that loses/gains a step to handle coordinate swapping, so setting the GDAL SRS axis policy on a CRS may or may not work (when using PROJ code to transform). So far, I hope, so good ... and thanks for the nice present!
@rsbivand Looks good, thanks! A bit short on time right now, but I'll try to have a better look in the next days.
PS: sorry for the present... ;-(
New fun: see https://lists.osgeo.org/pipermail/proj/2020-February/009325.html. Briefly, when a grid is dowloaded and used in a coordinate operation pipeline, coordinates outside the grid's area of interest may be returned as missing.
Video https://www.youtube.com/watch?v=D4-roPsMz48 and slides https://github.com/rsbivand/celebRation20_files from a talk at celebRation 2020 in Copenhagen published. The PROJ part starts about 37:20; enjoy (if that is the right expression).
Now +towgs84= r-spatial/sf#1328