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.


# 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]]
  • 1
    It kind of puzzles me why would you want to do that? Package {sf} uses the PROJ backend to reproject coordinates, and the backend is both very efficient in operations and known to be accurate in result – Jindra Lacko Apr 2 at 11:49
  • @JindraLacko it's a fair question. In the same way sqrt((df1_x - df2_x)^2 + (df1_y - df2_y)^2) can be a lot quicker than st_distance() to calculate distances - i'm just seeing if there is a more efficient way of achieving the same result. Efficiency here can be both speed (although less likely with this because as you say {sf} is not exactly sluggish) but also in terms of less resource (not having to create extra objects etc.). If there is a better way of using {sf} than in my example then I would be interested in that as well. – Chris Apr 2 at 15:28

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.

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.