Look up OS grid letters from easting/northing in R

Question

Suppose I have an easting/northing point:

data<-c("423851", "300384")

I would like to find out the grid letters for this point using the British National Grid system (Ordnance Survey). If you visit https://gridreferencefinder.com/ and put those coordinates into the 'x (easting)' and 'y (northing)' boxes, it will tell you the two-letter grid - in this case SK.

Is there a package/function within r that will allow me to pass multiple coordinates and return the grid letters? I've searched for some time now and found nothing. I suspect that I'm just using the wrong terms but would like to check for a pre-existing function before considering writing my own.

Answer 1

Here's a quick code I just wrote to compute the grid square given the coordinates:

grid_from_ref <- function(x,y){
### OSGB uses A-Z in a 5x5 grid except "I" (9th LETTER)
### Arrange the matrix correctly...
    m = matrix(LETTERS[-9],ncol=5,byrow=TRUE)[5:1,]
### 500k square letter. There's a false origin so divide by 500k,
### truncate, and add the false origin offset:
    xo = trunc(x/500000) + 3
    yo = trunc(y/500000) + 2
### lookup in the matrix (Y first)
    s1 = m[yo,xo]
### 100k square letter. Take what's left after 500k, divide, trunc, and add 1
    xo = trunc((x %% 500000)/100000) + 1
    yo = trunc((y %% 500000)/100000) + 1
### lookup
    s2 = m[yo,xo]
### stick together
    paste0(s1,s2)
}

Usage:

> grid_from_ref(429509, 479171)
[1] "SE"

It's not sensible vectorised so don't try and do more than one point:

> grid_from_ref(c(429509,45663), c(479171,534231))
[1] "SE" "NZ" "SA" "NV"

and that is mostly nonsense.

A quick test with a CSV from https://gridreferencefinder.com/ seems to work:

8  NY 29127 21621      329127       521621  NY
9  NW 69910 22642      169910       522642  NW
10 TL 78765 60570      578765       260570  TL
11 SJ 48283 23078      348283       323078  SJ
12 SO 11450 52146      311450       252146  SO
13 SP 47205 00337      447205       200337  SP
14 NJ 49721 15950      349721       815950  NJ

The first columns are the original data, the last is the grid square from my code computed from the coordinates, and it matches the grid square in the first column.

Answer 2

One option could be be to download the National Grid as a polygon feature and perform a query to find which tile each point falls on. Following answers to this question I found a link to National Grid shapefiles here: https://github.com/charlesroper/OSGB_Grids

This then allowed me to run this code which returns the correct result for two test points:

library(sp)
library(raster)

# Establishes a point feature of two points
id <- c('a', 'b')
x <- c(423851, 412731)
y <- c(300384, 289999)
pts <- cbind.data.frame(id, x, y)
coordinates(pts) <- ~x+y
proj4string(pts) <- CRS("+init=epsg:27700")

# Loads the downloaded 100 km grid, provided it is located in the working directory
osGrid <- shapefile('OSGB_Grid_100km.shp')
# Makes sure the coordinate systems match
osGrid <- spTransform(osGrid, CRS("+init=epsg:27700"))

# Finds which grid polygons the points fall on
over(pts, osGrid)

This returns the below results:

  TILE_NAME X250K
1        SK     t
2        SP     t

Answer 3

To vectorise the accepted 'grid_from_ref' answer, change the matrix lookups to:

 s1 = m[cbind(yo,xo)]

and

s2 = m[cbind(yo,xo)]

The code should then work with vector lists of eastings and northings.

You can get the whole National Grid Reference by changing the final 'paste0' to:

paste0(s1,s2,substr(x,2,6),substr(y,2,6))