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
              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)

              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")) %>% 

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.

  • Seems to be a difference in the CRS when you specify as numeric or +init - compare st_crs(27700)$proj4string and st_crs("+init=epsg:27700")$proj4string - that shows the problem without any data required. – Spacedman Feb 27 '19 at 13:59
  • At an even lower level its a difference between sf:::CPL_crs_from_proj4string("+init=epsg:27700") and sf:::CPL_crs_from_epsg(as.integer(27700)) - these functions then call C code. – Spacedman Feb 27 '19 at 14:08

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:


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,

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...

  • Thanks for this really detailed breakdown. Makes perfect sense. Though i agree with your conclusion that there shouldn't be inconsistencies in the output. Wondering if I should put up an issue on github for sf() – Jeremy Kidwell Feb 27 '19 at 18:52
  • Its not an sf issue, I suspect Edzer will close it! The problem is only really going to happen if you are inconsistent about assigning CRSs to data. – Spacedman Feb 27 '19 at 19:02
  • So is there an alternate way to build that bbox polygon above which will call sf:::CPL_crs_from_epsg(as.integer()) instead? – Jeremy Kidwell Feb 27 '19 at 19:32
  • In particular, I'm trying to default to using epsg, as this seems to be increasingly preferred to proj4... – Jeremy Kidwell Feb 27 '19 at 19:41
  • Given an integer, you could construct the proj4 string via paste0("+init=epsg:",intepsg) or similar... – Spacedman Feb 27 '19 at 21:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.