[R-sig-Geo] cannot apply a gridded transformation with sf

[Roger Bivand]

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
>> 
>> 
>
>

-- 
Roger Bivand
Emeritus Professor
Department of Economics, Norwegian School of Economics,
Postboks 3490 Ytre Sandviken, 5045 Bergen, Norway.
e-mail: Roger.Bivand using nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en