EPSG info lost on writing to shapefile?

[jreades]

We briefly discussed this on Twitter, so my apologies for not following up right away, but it seems to me like there is an issue with how CRSes/SRIDs are handled by sf on writing to shapefiles. Simple demonstration below. Session info also.

Addendum: had a quick check just before hitting submit and the same issue does not occur with GeoJSON, so perhaps it's specific to the OGR ESRI driver and the parameters that it needs to write .prj files properly? In which case it may be nothing to do with sf, but since shapefiles are a common interchange format (and useful for looking at results) it would be nice if this behaviour were more predictable or prominently warned about -- in QGIS all the files I write from R now come out as 'User Specified CRS' and won't play nicely with data download from the Ordnance Survey since they appear to be in a different projection.

Regardless of minor issues, am really excited about the library overall as is much faster than sp and much more intuitive to boot... Cheers!

library(sf)

# Create a polygon
box <- st_polygon(list(rbind(c(0,0),c(1000,0),c(1000,1000),c(0,1000),c(0,0))))

# Set the CRS to EPSG:27700 (BNG)
box <- st_sfc(box) %>% st_set_crs(NA) %>% st_set_crs(27700)
print(st_crs(box))

# Write it as a shapefile
# Would also be nice if delete_dsn didn't throw an 
# error if test.shp doesn't exist (built-in `if exists()`?)
st_write(box, 'test.shp', delete_dsn=TRUE) 

# Read back what was just written
shp <- st_read('test.shp')

# EPSG info is lost even though sf set it
print(st_crs(shp))

# Test with GeoJSON
st_write(box, 'test.geojson', delete_dsn=TRUE) 

# Read back what was just written
shp2 <- st_read('test.geojson')

# EPSG info is still there
print(st_crs(shp2))

Session info

R version 3.3.3 (2017-03-06)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: macOS Sierra 10.12.5

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] sf_0.5-0

loaded via a namespace (and not attached):
[1] magrittr_1.5  DBI_0.6-1     tools_3.3.3   units_0.4-4   Rcpp_0.12.11  udunits2_0.13 grid_3.3.3  

[rsbivand]

The shapefile format puts the WKT SRS in a *.prj file, but does not provide for storage of the EPSG ID. So you should not expect to read something the format cannot handle. Note that GPKG supports storage and retreival of the EPSG ID.

[jreades]

Aaaah, now I properly understand. For some reason I had managed to convince myself that when I was working in Python with geopandas that I had been able to read/write EPSG info with the shape file format. Having gone back and created a test case I can see that this isn't happening. Not sure where I got that idea from but I'll close this and apologies for not groking the original reply on Twitter.