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

[Jean-Luc Dupouey]

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

-- 
Jean-Luc Dupouey
INRAE
Silva Unit
F-54280 Champenoux
France


	[[alternative HTML version deleted]]