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

[Roger Bivand]

One basic problem is that a regression was introduced in PROJ from 7.2.0, 
resolved last week by https://github.com/OSGeo/PROJ/pull/2884, where grid 
transformations ceased to work if the prime meridian was not Greenwich. So 
for PROJ 7.2.0-8.1.1, the closest you can get is:

ptSf_2154_via_CRS84 <- st_transform(st_transform(ptSf_27572, 
crs="OGC:CRS84"), crs=2154)
st_distance(ptSf_2154_via_CRS84, ptSf_2154)
# Units: [m]
#          [,1]
# [1,] 3.777444

through a CRS84 hub using a Helmert transformation, unless you can force 
the pipeline as shown in the sf issue: 
https://github.com/r-spatial/sf/issues/1815
This forcing so far only works on Linux, as the user-writable directory is 
not properly accessed in sf for Windows or macOS yet. You may also use the 
rgdal work-around setting the user-writable directory as shown earlier in 
this thread.

Roger

On Mon, 11 Oct 2021, Roger Bivand wrote:

> 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