cannot apply a gridded transformation with sf
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 at nhh.no https://orcid.org/0000-0003-2392-6140 https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en