Set CRS in a pipe - idea

[Nowosad]

Below you can find my broad idea how it should work. Basically, I think there should be a function to set a crs in a pipe, such as st_crs(sf_object, crs).

library('sf')
library('dplyr')
filepath <- system.file("shape/nc.shp", package="sf")
nc <- st_read(filepath) %>% st_crs(., "4326")
nc <- st_read(filepath) %>% st_crs(., "+init=epsg:4326")

[edzer]

This works:

nc <- st_read(filepath, quiet = TRUE) %>% `st_crs<-`(4326)

[Nowosad]

Yes, I know - but this way is very unintuitive. I think there will be a large number of people using sf functions in pipes. For example just to read data:

library('sf')
library('dplyr')
filepath <- system.file("shape/nc.shp", package="sf")
nc <- st_read(filepath) %>% 
          st_crs(., "+init=epsg:4326") %>% 
          st_transform(., "+init=epsg:2180")

[edzer]

We also don't want to break type stability as discussed here. I added st_set_crs for use in pipes, e.g. as in

x <- st_sfc(st_point(0:1)) %>% st_set_crs(4326) %>% st_transform(3857)

[Nowosad]

library('sf')
library('dplyr')
filepath <- system.file("shape/nc.shp", package="sf")
x <- st_read(filepath) %>%
        st_set_crs(4326) %>% 
        st_transform(3857)
head(x)

Looks great - thanks!

[rsbivand]

No, this is really not OK:

>filepath <- system.file("shape/nc.shp", package="sf")
> x <- st_read(filepath) %>%
+         st_set_crs(4326) %>% 
+         st_transform(3857)
Reading layer nc from data source /home/rsb/lib/r_libs/sf/shape/nc.shp using driver ESRI Shapefile
features:       100
fields:         14
converted into: MULTIPOLYGON
proj4string:    +proj=longlat +datum=NAD27 +no_defs 
Warning message:
In st_crs<-.sfc(*tmp*, value = list(list(list(c(-81.4727554321289,  :
  st_crs<- : replacing crs does not reproject data; use st_transform for that

Note:

> head(st_read(filepath))
Reading layer nc from data source /home/rsb/lib/r_libs/sf/shape/nc.shp using driver ESRI Shapefile
features:       100
fields:         14
converted into: MULTIPOLYGON
proj4string:    +proj=longlat +datum=NAD27 +no_defs 
Simple feature collection with 6 features and 14 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
epsg (SRID):    NA
proj4string:    +proj=longlat +datum=NAD27 +no_defs

so the pipeline makes it really easy to ignore the fact that the input shapefile has a CRS, which is correct because it is pre-WGS84, and imposing WGS84 is a major blunder. In a pipeline this should hard-fail by default, and only be permitted if an override is set actively. This is hard to do without option setting in a set method which typically only takes the opbject and value.

Where is st_set_crs() defined? Is it autogenerated, it isn't in any file in R/?

By the way, EPSG-guessing is not going to work - the OGR SRS mostly fails miserably.

[edzer]

Hence the warning -- the "correct" example in this case would be

library('sf')
library('dplyr')
filepath <- system.file("shape/nc.shp", package="sf")
x <- st_read(filepath) %>%
    st_set_crs(4267) %>% 
    st_transform(3857)

where, obviously, the st_set_crs() step is obsolete since the shapefile knows its crs, so to speak. It does only add the EPSG code, which st_read cannot guess.

st_set_crs() is in crs.R.

EPSG guessing is only implemented for 4326, or if proj4string contains +init=epsg:, iirc, no OGR code being called here.

[rsbivand]

Just picked up changed crs.R after a struggle. I guess this is non-standard evaluation, so no way for the function to know that it is in a pipeline and prevent the CRS being changed on an object for which it wasn't NA. In sp this is done with local options (warn or no warn, never error unless the warn option is set to throw an error.

OGRSpatialReference.AutoIdentifyEPSG() used in ogrAutoIdentifyEPSG() in src/ogr_proj.cpp, but it very often misses and returns nothing.

[edzer]

Think of the first argument of a function following %>% as stdin - yes, NSE.

ogrAutoIdentifyEPSG() looks useful!