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

[Roger Bivand]

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