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]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
