Why does R set CRS inconsistently?

Question

I'm running into a situation (using R on MacOS) where sf() is giving different levels of precision on the CRS. This causes problems for geospatial operations, as indicated below with st_crop().

Can anyone explain what's going on here? I'm admittedly puzzled.

Reproducible example:

# Download data
download.file("http://sedsh127.sedsh.gov.uk/Atom_data/ScotGov/ZippedShapefiles/SG_UrbanRural_2016.zip", 
              destfile = "data/SG_UrbanRural_2016.zip")
unzip("data/SG_UrbanRural_2016.zip", exdir = "data")

urbanrural_sf <- st_read("data/SG_UrbanRural_2016.shp") %>% st_transform(27700)

download.file("https://borders.ukdataservice.ac.uk/ukborders/easy_download/prebuilt/shape/Scotland_ca_2010.zip", 
              destfile = "data/Scotland_ca_2010.zip")
unzip("data/Scotland_ca_2010.zip", exdir = "data")

admin_lev1_sf <- st_read("data/scotland_ca_2010.shp") %>% st_transform(27700)

centralbelt_region <- st_bbox(c(xmin = 224479.2, xmax = 642963.5,
                      ymin = 347475.0, ymax = 711014.5),
                    crs = st_crs("+init=epsg:27700")) %>% 
  st_as_sfc()

Running st_crs shows how an extra decimal of precision is present on the centralbelt_region polygon (scroll to the right to see). Otherwise CRS are exactly the same:

> st_crs(centralbelt_region)
Coordinate Reference System:
  EPSG: 27700 
  proj4string: "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.1502,0.247,0.8421,-20.4894 +units=m +no_defs"
> st_crs(admin_lev1_sf)
Coordinate Reference System:
  EPSG: 27700 
  proj4string: "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs"
> 

This causes problems when I try to run:

> urbanrural_sf_simplified_centralbelt <- st_crop(urbanrural_sf_simplified, centralbelt_region)
Error in geos_op2_geom("intersection", x, y) : 
  st_crs(x) == st_crs(y) is not TRUE

Would be really glad for any insights folk have on this.

Answer 1

I think that sf is ultimately delegating this to two different "authorities"...

The differences boil down to this R code before hitting C code:

> sf:::CPL_crs_from_proj4string("+init=epsg:27700")
Coordinate Reference System:
  EPSG: 27700 
  proj4string: "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy 
   +towgs84=446.448,-125.157,542.06,0.1502,0.247,0.8421,-20.4894 +units=m +no_defs"
> sf:::CPL_crs_from_epsg(as.integer(27700))
Coordinate Reference System:
  EPSG: 27700 
  proj4string: "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy 
  +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs"

The difference is only in the datum shift parameters:

+towgs84=446.448,-125.157,542.06,0.1502,0.247,0.8421,-20.4894
+towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489

These come from proj and gdal. proj can tell you its datum shift parameters with proj -ld:

$ proj -ld 
      [...]
      OSGB36 airy         towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894

and gdal has them in /usr/share/gdal/2.2/datum_shift.csv:

$ grep OSGB  /usr/share/gdal/2.2/datum_shift.csv

530,1314,4277,4326,<etc etc> 1264,49.79,60.94,-8.82,1.92,1,0,9606,
  446.448,-125.157,542.06,0.15,0.247,0.842,-20.489,1

This shows the same difference in the precision, which is why I think lookup by number is doing one thing and lookup by character is doing the other.

Not sure of the solution though - ultimately the two authorities should be made to agree...