将 OS National Grid 名称/代码添加到 R 中的网格

我希望重新创建完整的地形测量国家网格（如此处所示https://upload.wikimedia.org/wikipedia/commons/f/f5/Ordnance_Survey_National_Grid.svg https://upload.wikimedia.org/wikipedia/commons/f/f5/Ordnance_Survey_National_Grid.svg) in R.

我可以毫无问题地创建 4 个级别的网格（使用 36GB RAM 的 1 公里版本大约需要 20 分钟）：

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()

但是，如何命名所有 500 公里、100 公里、10 公里和 1 公里单元格以反映图像中显示的命名/编码约定？

我想没有多少人愿意在测试他们的答案的过程中花费 20 分钟来生成 6,125,000 个多边形。碰巧的是，我不得不将基于 Windows 的云服务器升级到 32GB，只是为了创建 1 公里的方块......

不幸的是，当您创建网格方块时，它们是从下到上、从左到右排序的，而标签是按顺序排序的从上到下， 左到右。这使得标签有点棘手。对于 500 公里的盒子，我们希望能够给它们贴上标签LETTERS[-9]，但由于顺序的原因，我们需要将它们标记为LETTERS[-9][rep(4:0 * 5, each = 5) + 1:5].

我们可以通过将数据框与网格名称绑定到网格对象来创建命名网格，如下所示：

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()

现在我们可以绘图以确保我们拥有正确的标签：

library(ggplot2)

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

第二级更难，因为我们需要重复每个方块内的行和列。这需要一些模块化数学才能获得正确的索引：

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() 

但我们可以看到这也是有效的：

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

同样，由于重复、模块化数学和子集设置，下一层变得更加复杂，但可以这样实现：

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()

显然，我现在无法绘制整个网格，因为它太小而无法看到各个方块（更不用说它们的标签），所以我将直接拉出TQ确保编号正确：

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))

最细的方块也可以用与 10 公里长的盒子类似的方式来标记，但更复杂的是，标记后您需要交换第四和第五位数字：

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()

同样，我们需要挑选一个小子集来展示这个作品：

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))