4

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.

Simon
  • 259
  • 3
  • 11
  • https://www.r-bloggers.com/converting-british-national-grid-references/ I don't know how good it is or if it's really in R, thus a comment. I searched "british national grid" R – mkennedy Jun 26 '19 at 15:46
  • I'll probably have something like that in my https://gitlab.com/b-rowlingson/gridsquares package.... Progress slow at the moment though... – Spacedman Jun 26 '19 at 17:48
  • @mkennedy better to link to the source https://scottishsnow.wordpress.com/2016/10/13/converting-bng/ than an aggregator site – Spacedman Jun 26 '19 at 18:25
  • Hi @mkennedy I did find the link you suggested (and should have included this in my question as a step I had taken myself to resolve the issue) but found that you already need to know the grid letters. It doesn't (as far as I can tell), allow you to work it out going 'backwards' so to speak. Thanks for the suggestions though, appreciated! – Simon Jun 27 '19 at 06:51
  • @Spacedman Nice to hear that you're working on a spatial package! I presume this is focused on the BNG system? I don't have a Gitlab account to follow the project but have bookmarked it and will be sure to check back regularly for progress. Any plans to have this cloned on Github? – Simon Jun 27 '19 at 06:53

3 Answers3

4

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.

Spacedman
  • 63,755
  • 5
  • 81
  • 115
  • Thanks Spacedman. I've now accepted this as the best answer as it does not rely on an external shapefile. I thank @humperderp again though for his contribution to my question. – Simon Jun 28 '19 at 06:53
3

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
humperderp
  • 837
  • 5
  • 14
  • 1
    That's probably overkill. It can be done with a subtract, divide and a lookup! – Spacedman Jun 26 '19 at 17:50
  • 1
    That's certainly a fair point. I'll leave my answer for reference in case anyone else prefer to swat their flies with sledgehammers, or perhaps could benefit from any extended functionality resulting from the approach. – humperderp Jun 26 '19 at 20:51
  • 1
    Thanks humperderp for your answer. Sledgehammer or not, it works as expected so I'm accepting the answer. Thank you! @Spacedman - I'd be keen to learn more about your comment above in relation to this question. Are you able to explain in more detail what you mean? – Simon Jun 27 '19 at 07:03
0

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))
MikeF
  • 1