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: [email protected]
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
