Using R to convert British National Grid to/from WGS84 without creating spatial object

Question

I have several million coordinates (which are either British National Grid or Long/Lat) and need to convert them. I know using R you can do this in a spatial operation (such as in my example below using sf). I was however wondering if it could be done without creating a spatial object first (and just be a mathematical function?). Is there a package that can do this, or has someone already written a function that does this that they are willing to share?

Pages 49 to 51 in this PDF detail a fairly long equation for the conversions - but before jumping into coding that I wondered if anyone had any alternative solutions.

library(data.table)
library(sf)

# Sample data
DT1 <- data.table(DT1_id = sprintf("%s%0*d", "UID_", 7, 1:3000000),DT1_x =sample(45000:55000, size = 3000000,replace = TRUE),DT1_y =sample(25000:35000, size = 3000000,replace = TRUE))

# BNG to WGS84
DT1_sf <- st_as_sf(DT1, coords = c("DT1_x","DT1_y"),crs=27700)

DT1_sf_4326 <- DT1_sf %>% st_transform(crs = 4326)  

DT1[, LAT := data.frame(st_coordinates(DT1_sf_4326))[,2]]
DT1[, LONG := data.frame(st_coordinates(DT1_sf_4326))[,1]]

# WGS84 to BNG
DT1_sf <- st_as_sf(DT1, coords = c("LONG","LAT"),crs=4326)

DT1_sf_27700 <- DT1_sf %>% st_transform(crs = 27700)  

DT1[, EASTING := data.frame(st_coordinates(DT1_sf_27700))[,1]]
DT1[, NORTHING := data.frame(st_coordinates(DT1_sf_27700))[,2]]

Answer 1

Use sf (and hence proj with its accurate transformation) to create a grid of test points in lat-long (X,Y) and corresponding accurate OSGB coordinates (A,B).

Now try various functions F optimised to minimise the difference in distance between (A,B) and F(A,B).

For example you could do a scaling of (X,Y):

first = function(A,B,a,b,c,d){
          return(matrix(c(a,b,c,d),2,2) %*% c(A,B))
        }

and then find (a,b,c,d) that minimise the total root mean square distance error. That won't be very good because the OSGB grid isn't oriented like lat-long, and that matrix transformation in the function doesn't do curves.

So you could add more polynomial terms, or different functional forms, and introduce more parameters to optimise over. There are lots of techniques from georeferencing you could use - look at the list of ways the QGIS Georeferencer can warp space, for example.

However at some point you'll either end up with something that is the same algorithm as proj uses, but probably not as efficiently written, or something that takes longer because you have overcomplicated it in order to fit a transform that you don't understand.

But if you are lucky you might hit a transformation function that gives acceptable accuracy for your use case and runs quicker than proj.

I would bet a small amount of magic internet points that nobody has done this yet because proj is just so darn good.

Answer 2

If we consider EPSG:4326 (Geodetic CRS WGS84) and EPSG:4258 (Geodetic CRS ETRS89) equivalent then we could use the package sgo that I released a couple of months ago (https://CRAN.R-project.org/package=sgo) whose main goal is the conversion from/to BNG using the OSTN15 transformation model, which is the definitive OSGB36/ETRS89 transformation within Great Britain.

Using similar data as the OP example we should have an accurate result being also much quicker to compute than using sf for this purpose:

library(sgo)

DT1 <- data.frame(DT1_id = sprintf("%s%0*d", "UID_", 7, 1:3000000),
  DT1_x =sample(45000:55000, size = 3000000, replace = TRUE),
  DT1_y =sample(25000:35000, size = 3000000, replace = TRUE),
  stringsAsFactors = FALSE)

# Create an sgo_points object (which is basically a list)
DT1_sgo <- sgo_points(DT1, coords=c("DT1_x","DT1_y"), epsg=27700)

# BNG to WGS84
DT1_sgo_4326 <- sgo_bng_lonlat(DT1_sgo, to=4326)

# WGS84 to BNG
DT1_sgo_27700 <- sgo_lonlat_bng(DT1_sgo_4326)