Projection differences in R using sf and sp

Question

I have a grid I have converted from GeoTIFFs to a shapefile. I would like to convert and export the shapefile as a GeoPackage and change the projection so it uses the British National Grid as the geographic coordinate system when opened in a GIS. However this only seems to work using sp and not sf (which does not appear to retain aspects like the datum).

This is a problem as I would like to export GeoPackages containing multiple layers which you can only currently do in sf and not sp. Am I doing something wrong?

library(rgdal)
library(sf)

download.file("https://drive.google.com/uc?id=1URbux7Sw25KFTySqRFKXk53DV2UK4lsA&export=download" , destfile="Grid_Shapefile.zip")
unzip("Grid_Shapefile.zip")
Grid_sp <- readOGR(".", "Grid_Shapefile")
Grid_sf <- st_as_sf(Grid_sp)

BNG_Grid_sp <- spTransform(Grid_sp, CRS("+init=epsg:27700"))
BNG_Grid_sf_v1 <- st_transform(Grid_sf, crs=27700)
BNG_Grid_sf_v2 <- st_transform(Grid_sf, crs="+init=epsg:27700 +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs +ellps=airy +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894")

BNG_Grid_sf_v1_geom <- st_geometry(BNG_Grid_sf_v1)
BNG_Grid_sf_v2_geom <- st_geometry(BNG_Grid_sf_v2)

proj4string(BNG_Grid_sp)
attributes(BNG_Grid_sf_v1_geom)
attributes(BNG_Grid_sf_v2_geom)

writeOGR(BNG_Grid_sp, dsn = "BNG_Grid_sp.gpkg", layer = "Grid_sp", driver = "GPKG")
st_write(BNG_Grid_sf_v1, "BNG_Grid_sf_v1.gpkg", "Grid_sf_v1")
st_write(BNG_Grid_sf_v2, "BNG_Grid_sf_v2.gpkg", "Grid_sf_v2")

Answer 1

Please follow up on https://github.com/r-spatial/sf/issues/810, not here.

Answer 2

Using ogrinfo, in particular the command

ogrinfo BNG*v1.gpkg Grid_sf_v1 > info1
ogrinfo BNG*v2.gpkg Grid_sf_v1 > info2

the diff between the two info[1|2] files are, besides the obvious _v1 _v2 naming, this:

13c13
<             TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489]],
---
>             TOWGS84[446.448,-125.157,542.06,0.1502,0.247,0.8421,-20.4894]],

Are the additional digits in _v2 causing your troubles in ArcGIS?

Answer 3

The solution to this (thanks to what Roger has posted here) is using the lwgeom package to do the transformation.The post by Roger on the sf GitHub gives further details.

library(rgdal)
library(sf)

download.file("https://drive.google.com/uc?id=1URbux7Sw25KFTySqRFKXk53DV2UK4lsA&export=download" , destfile="Grid_Shapefile.zip")
unzip("Grid_Shapefile.zip")
Grid_sp <- readOGR(".", "Grid_Shapefile")
Grid_sf <- st_as_sf(Grid_sp)

library(lwgeom)
BNG_Grid_sf_v4 <- st_transform_proj(Grid_sf, crs="+init=epsg:27700 +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs +ellps=airy +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894")
st_crs(BNG_Grid_sf_v4)
st_write(BNG_Grid_sf_v4, "BNG_Grid_sf_v4.gpkg", "Grid_sf_v4")