Search code examples
rrasterr-sfterra

Creating 1m, 5m, 10m, 50m, 100m and 500m grids with OS National Grid codes in R

I previously asked this question on how to add OS National Grid codes to grids 1km and above. However, I now need a method of creating grids below 1km (1m, 5m, 10m, 50m, 100m and 500m) that have OS National Grid codes assigned to them (I will then need to intersect points with these grid layers).

However, I am not sure what is the best way of achieving this in R. I think a raster solution is more likely to work but this is something I am less familiar with. Any ideas on how the solve this are welcome. This is as about as far as I can get:

library(terra)

# Define the extent of the grid
bbox <- c(-100000, 700000, -100000, 1400000)
res <- 50

# Create the grid
grid <- terra::rast(ext=bbox, res=res, crs="+init=epsg:27700")

# Generate unique IDs for each cell
ids <- terra::cellFromRowCol(grid, 1:terra::nrow(grid), 1:terra::ncol(grid))
ids <- sprintf("%03d%03d%s", floor(ids/4000)+1, (ids-1)%%4000+1,
               c("SW", "NW", "NE", "SE")[(floor((ids-1)/4000)+1)%%4+1])

# Set the cell values to the unique IDs
grid[] <- ids

In terms of the logic for the unique cells ID:

• 1 metre grid ID is based on the bottom left coordinates of a cell. The first 2 letters refer to the 100km cell it is in. The next five numbers are the last five numbers of the bottom left x coordinates of that cell. The next five numbers are the last five numbers of the bottom left y coordinates of that cell.

• 5 metre grid ID is based on the bottom left coordinates of a cell. The first 2 letters refer to the 100km cell it is in. The next four numbers are four numbers from the end of the bottom left x coordinates of that cell, except the last number in the coordinate is disregarded. The last four numbers are four numbers from the end of the bottom left y coordinates of that cell, except the last number in the coordinate is disregarded. To distinguish between four cells having the same ID, each cell has either SW, SE, NE or NW added to the end and is assigned in a counterclockwise direction starting in the bottom left.

• 10 metre grid ID is based on the bottom left coordinates of a cell. The first 2 letters refer to the 100km cell it is in. The next four numbers are four numbers from the end of the bottom left x coordinates of that cell, except the last number in the coordinate is disregarded. The last four numbers are four numbers from the end of the bottom left y coordinates of that cell, except the last number in the coordinate is disregarded.

• 50 metre grid ID is based on the bottom left coordinates of a cell. The first 2 letters refer to the 100km cell it is in. The next three numbers are three numbers from the end of the bottom left x coordinates of that cell, except the last two numbers in the coordinate is disregarded. The last three numbers are three numbers from the end of the bottom left y coordinates of that cell, except the last two numbers in the coordinate is disregarded. To distinguish between four cells having the same ID, each cell has either SW, SE, NE or NW added to the end and is assigned in a counterclockwise direction starting in the bottom left.

• 100 metre grid ID is based on the bottom left coordinates of a cell. The first 2 letters refer to the 100km cell it is in. The next three numbers are three numbers from the end of the bottom left x coordinates of that cell, except the last two numbers in the coordinate is disregarded. The last three numbers are three numbers from the end of the bottom left y coordinates of that cell, except the last two numbers in the coordinate is disregarded.

• 500 metre grid ID is based on the bottom left coordinates of a cell. The first 2 letters refer to the 100km cell it is in. The next two numbers are two numbers from the end of the bottom left x coordinates of that cell, except the last three numbers in the coordinate is disregarded. The last two numbers are two numbers from the end of the bottom left y coordinates of that cell, except the last three numbers in the coordinate is disregarded. To distinguish between four cells having the same ID, each cell has either SW, SE, NE or NW added to the end and is assigned in a counterclockwise direction starting in the bottom left.

example of cell ID logic

Solution

  • You example data

    pnts <- rbind(cbind(93555, 256188), 
      cbind(210637, 349798),
      cbind(696457,481704))
colnames(pnts) <- c("x", "y")

    With geosphere v1.5-19 you can use geosphere::OSGB.

    library(geosphere)
OSGB(pnts, "1km")
#[1] "SL9356" "SH1049" "TB9681"
OSGB(pnts, "1m")
#[1] "SL9355556188" "SH1063749798" "TB9645781704"
OSGB(pnts, "5km")
#[1] "SL95NW" "SH14NW" "TB98SE"

    To create your table for the resolutions starting with "1"

    res <- c('1m', '10m', '100m', '1km', '10km', '100km')
data.frame(pnts, sapply(res, \(i) OSGB(pnts, i)), check.names=FALSE)
#       x      y           1m        10m     100m    1km 10km 100km
#1  93555 256188 SL9355556188 SL93555618 SL935561 SL9356 SL95    SL
#2 210637 349798 SH1063749798 SH10634979 SH106497 SH1049 SH14    SH
#3 696457 481704 TB9645781704 TB96458170 TB964817 TB9681 TB98    TB

    And for the resolutions starting with "5"

    res <- c('5m', '50m', '500m', '5km', '50km', '500km')
data.frame(pnts, sapply(res, \(i) OSGB(pnts, i)), check.names=FALSE)
#       x      y           5m        50m     500m    5km 50km 500km
#1  93555 256188 SL93555618NE SL935561NE SL9356SE SL95NW SLNE     S
#2 210637 349798 SH10634979NE SH106497NW SH1049NE SH14NW SHSW     S
#3 696457 481704 TB96458170SE TB964817SE TB9681NW TB98SE TBNE     T

    geosphere v1.5-19 is currently the development version. You can install it with

    install.packages('geosphere', repos='https://rspatial.r-universe.dev')