Search code examples
rspatialr-sf

Add OS National Grid names/codes to grid in R

I am looking to recreate the full Ordnance Survey National Grid (as shown here https://upload.wikimedia.org/wikipedia/commons/f/f5/Ordnance_Survey_National_Grid.svg) in R.

I can create the grid for the 4 levels without any issues (the 1km version takes about 20 minutes for me with 36GB of RAM):

library(sf)

OS_National_Grid_BBox <- st_bbox(c(xmin = -1000000, xmax = 1500000, ymax = 2000000, ymin = -500000), crs = st_crs(27700))

OS_National_Grid_500km <- st_make_grid(st_as_sfc(OS_National_Grid_BBox), square = T, cellsize = c(500000, 500000)) %>% st_sf()
OS_National_Grid_100km <- st_make_grid(st_as_sfc(OS_National_Grid_BBox), square = T, cellsize = c(100000, 100000)) %>% st_sf()
OS_National_Grid_10km <- st_make_grid(st_as_sfc(OS_National_Grid_BBox), square = T, cellsize = c(10000, 10000)) %>% st_sf()
OS_National_Grid_1km <- st_make_grid(st_as_sfc(OS_National_Grid_BBox), square = T, cellsize = c(1000, 1000)) %>% st_sf()

However, how can I name all the 500km, 100km, 10km and 1km cells to reflect the naming/coding convention shown in the image?

Solution

  • I guess not too many folks want to tie up their computer for 20 minutes to generate the 6,125,000 polygons en route to testing their answer. As it happens, I had to upgrade my Windows-based cloud server to 32GB just to create the 1km squares...

    Unfortunately when you create your grid squares, they are ordered bottom-to-top, left-to-right, whereas the labels are ordered top-to-bottom, left-to-right. This makes the labeling a little trickier. For the 500km boxes, we would expect to be able to label them with LETTERS[-9], but because of the ordering we need to label them with LETTERS[-9][rep(4:0 * 5, each = 5) + 1:5].

    We can create our named grids by binding together a data frame with the grid names to your grid object, like this:

    gridref500 <- LETTERS[-9][rep(4:0 * 5, each = 5) + 1:5]

OS_National_Grid_500km <- st_as_sfc(OS_National_Grid_BBox)                    %>%
                          st_make_grid(square = TRUE, cellsize = c(5e5, 5e5)) %>% 
                          cbind(data.frame(Grid_Ref = gridref500))            %>%
                          st_sf()

    Now we can plot to ensure we have the correct labels:

    library(ggplot2)

ggplot(OS_National_Grid_500km) + 
  geom_sf(fill = "white")      + 
  geom_sf_text(aes(label = Grid_Ref), size = 5)

    The second level is tougher, since we need to repeat the rows and columns within each of the squares. This requires a bit of modular maths to get the indexing right:

    gridref100 <- rep(gridref500, each = 5) %>%
              split(0:124 %/% 25)       %>%
              lapply(rep, 5)            %>%
              do.call(c, .)             %>%
              paste0(split(gridref500, 0:24 %/% 5) %>%
                     lapply(rep, 5)                %>%
                     do.call(c, .)                 %>%
                     rep(5))

OS_National_Grid_100km <- st_as_sfc(OS_National_Grid_BBox)                    %>%
                          st_make_grid(square = TRUE, cellsize = c(1e5, 1e5)) %>% 
                          cbind(data.frame(Grid_Ref = gridref100))            %>%
                          st_sf()

    But we can see this is also effective:

    ggplot(OS_National_Grid_100km) + 
  geom_sf(fill = "white")      + 
  geom_sf_text(aes(label = Grid_Ref), size = 3)

    enter image description here

    Again, the next layer becomes more complex because of the repetitions, modular maths and subsetting, but can be achieved like this:

    gridref10  <- rep(gridref100, each = 10) %>%
              split(0:6249 %/% 250)      %>%
              lapply(rep, 10)            %>%
              do.call(c, .)              %>%
              paste0(as.character(rep(0:9, 6250))       %>%
                     paste0(rep(rep(0:9, each = 250), 25)))

OS_National_Grid_10km <- st_as_sfc(OS_National_Grid_BBox)                    %>%
                         st_make_grid(square = TRUE, cellsize = c(1e4, 1e4)) %>% 
                         cbind(data.frame(Grid_Ref = gridref10))             %>%
                         st_sf()

    Obviously, I can't plot the entire grid now because it would be far too small to see individual squares (let alone their labels), so I will just pull out TQ to ensure the numbering is correct:

    TQ <- OS_National_Grid_10km[substr(OS_National_Grid_10km$Grid_Ref, 1, 2) == "TQ",]

ggplot(TQ)                + 
  geom_sf(fill = "white") + 
  geom_sf_text(aes(label = Grid_Ref))

    enter image description here

    The finest squares would be possible to label too in a similar way to the 10km boxes, with the added complication that after labeling you need to swap the fourth and fifth digits:

    gridref1  <- rep(gridref10, each = 10)   %>%
  split(0:624999 %/% 2500)    %>%
  lapply(rep, 10)             %>%
  do.call(c, .)               %>%
  paste0(as.character(rep(0:9, 625000))      %>%
           paste0(rep(rep(0:9, each = 2500), 250)))

swapchar               <- substr(gridref1, 4, 4)
substr(gridref1, 4, 4) <- substr(gridref1, 5, 5)
substr(gridref1, 5, 5) <- swapchar

OS_National_Grid_1km  <- st_as_sfc(OS_National_Grid_BBox)                    %>%
  st_make_grid(square = TRUE, cellsize = c(1e3, 1e3)) %>% 
  cbind(data.frame(Grid_Ref = gridref1))              %>%
  st_sf()

    Again, we need to pick out a small subset to show this works:

    ss <- with(OS_National_Grid_1km, 
           which(paste0(substr(Grid_Ref, 1, 3), substr(Grid_Ref, 5, 5)) == "TQ28"))

TQ28 <- OS_National_Grid_1km[ss,]

ggplot(TQ28) + 
  geom_sf(fill = "white") + 
  geom_sf_text(aes(label = Grid_Ref))

    enter image description here