cannot apply a gridded transformation with sf

Assuming the Windows CRAN binary. Binary sf (and rgdal) packages bundle PROJ 
and GDAL metadata files taken as fixed for given releases of PROJ and GDAL, 
but from PROJ 7, a content download network (CDN, https://cdn.proj-org) is 
also available. Binary packages set the PROJ_LIB environment variable when 
loaded to the bundled metadata. From PROJ 7, a separate user-writable 
directory is also needed to cache transformation grids downloaded from the 
CDN.

The current difficulty stems both from problems passing through the string 
defining the user-writable directory, and from st_transform() (and 
rgdal::spTransform()) not choosing the most accurate transformation pipeline 
(possibly because the user-writable directory is not properly detected ??). 
The following code works around the problem by coercing the sf object to a 
Spatial object from sp, using spTransform() instead of st_transform, then 
coercing back to sf. spTransform() does not use PROJ through GDAL, but uses 
PROJ directly; it still needs to be told which transformation pipeline to 
use.

#  at the very beginning of the R session, define the user-writable directory
td <- tempfile()
dir.create(td)
Sys.setenv("PROJ_USER_WRITABLE_DIRECTORY"=td)
library(rgdal)
rgdal_extSoftVersion()
#  check the PROJ search paths and check that the user-writable directory is 
#  empty and that the CDN is not enabled
(pths <- get_proj_search_paths())
list.files(pths[1])
is_proj_CDN_enabled()
# define the objects and check that the search paths are still OK
library(sf)
get_proj_search_paths()
pt_27572 <- data.frame(X=55824.4970,Y=2394454.2120)
pt_2154 <- data.frame(X=107242.8310,Y=6832277.1820)
ptSf_27572 <- st_as_sf(pt_27572, coords = c("X", "Y"), crs=27572)
ptSf_2154 <- st_as_sf(pt_2154, coords = c("X", "Y"), crs=2154)
# enable the CDN and that the user-writable directory is still empty
enable_proj_CDN()
is_proj_CDN_enabled()
list.files(pths[1])
# list candidate coordinate operations
(o <- list_coordOps("EPSG:27572", "EPSG:2154"))
# transform using the first returned coordinate operation pipeline
ptSf_2154_grid <- st_as_sf(spTransform(as(ptSf_27572, "Spatial"),
 CRS("EPSG:2154"), coordOp=o[1, "definition"]))
# print the coordinate operation used
get_last_coordOp()
# check that the user-writable directory is populated
list.files(pths[1])
# check distance to defined point in target crs
st_distance(ptSf_2154_grid, ptSf_2154)

With PROJ 7.2.1 this yields a different sub-millimetre distance from PROJ 
8.1.1, because a different grid is preferred, but both do transform your 
input point to your expected output point if half a millimetre is close 
enough.

I've raised https://github.com/r-spatial/sf/issues/1815 with some more 
details. This may take some time to resolve satisfactorily.

Roger

On Mon, 11 Oct 2021, Roger Bivand wrote:

 Dear Jean-Luc,

 Thanks for a very helpful report. I'll answer fully later, but if you
 could provide in addition your platform (I think Windows or macOS, sf
 installed as a CRAN binary), that would help. On these platforms, access
 to the CDN for grids is not what it should be, and your example is helping
 debug the problem.

 Best wishes,

 Roger

 On Sun, 10 Oct 2021, Jean-Luc Dupouey wrote:

  Dear forumites,

  I try to apply a gridded transformation of coordinates in |sf|, but it
  does not work. Is there something else to do apart from calling
|  sf_proj_network(TRUE)|? Here is the 7 lines of code I used:

|  pt_27572 <- data.frame(X=55824.4970,Y=2394454.2120)| # set coordinates
  of one point in EPSG:27572

|  pt_2154 <- data.frame(X=107242.8310,Y=6832277.1820)| # known accurate
  coordinates of the same point in EPSG:2154 (from calculations of the
  French National Geographic Institute)

|  ptSf_27572 <- st_as_sf(pt_27572, coords = c("X", "Y"), crs = 27572)| #
  build sf object

|  ptSf_2154 <- st_as_sf(pt_2154, coords = c("X", "Y"), crs = 2154)| #
  build sf object

|  sf_proj_network(TRUE)| # allow search for online datum grids in the
  PROJ CDN

  [1] "https://cdn.proj.org"

|  sf_proj_pipelines("EPSG:27572", "EPSG:2154")| # grids
  (fr_ign_gr3df97a.tif) seem to be found for accurate transformation from
  EPSG:27572 to EPSG:2154

|  Candidate coordinate operations found: 3 Strict containment: FALSE Axis
  order auth compl: FALSE Source: EPSG:27572 Target: EPSG:2154 Best
  instantiable operation has accuracy: 1 m Description: Inverse of Lambert
  zone II + NTF (Paris) to NTF (1) + NTF to RGF93 (1) + Lambert-93
  Definition: +proj=pipeline +step +inv +proj=lcc +lat_1=46.8 +lat_0=46.8
  +lon_0=0 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +ellps=clrk80ign
  +pm=paris +step +proj=push +v_3 +step +proj=cart +ellps=clrk80ign +step
  +proj=xyzgridshift +grids=fr_ign_gr3df97a.tif +grid_ref=output_crs
  +ellps=GRS80 +step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3
  +step +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000
  +y_0=6600000 +ellps=GRS80 |

|  ptSf_2154_grid <- st_transform(ptSf_27572,crs=2154)| # apply 
|  transformation

|  st_distance(ptSf_2154_grid,ptSf_2154)| # incorrect (ungridded)
  transformation, the distance should be zero. 3.777 m is the known error
  for the ungridded transformation.

|  Units: [m] [,1] [1,] 3.777346 |

  I read in (the so useful) "Spatial Data Science" that by default, the
  most accurate pipeline is chosen by |st_transform|. Why isn't it the
  case here?

  I am running R version 4.1.1, GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1.

  Thanks for your help,

  Jean-Luc