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

Question

我希望在 R 中重新创建完整的 Ordnance Survey National Grid(如此处所示 https://upload.wikimedia.org/wikipedia/commons/f/f5/Ordnance_Survey_National_Grid.svg)。

我可以毫无问题地创建 4 个级别的网格(对于 36GB RAM 的我，1km 版本大约需要 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()

但是，我如何命名所有 500km、100km、10km 和 1km 单元以反射(reflect)图像中显示的命名/编码约定？

Answer 1

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

不幸的是，当您创建网格正方形时，它们是从下到上、从左到右排序的，而标签是从从上到下、从左到右排序的。这使得标签有点棘手。对于 500 公里的框，我们希望能够用 LETTERS[-9] 标记它们，但由于顺序我们需要用 LETTERS[-9][rep( 4:0 * 5，每个 = 5) + 1:5]。

我们可以通过将带有网格名称的数据框绑定(bind)到您的网格对象来创建命名网格，如下所示:

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)

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

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

显然，我现在无法绘制整个网格，因为它太小了，看不到单个方 block (更不用说它们的标签了)，所以我将拉出 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))

最好的方 block 也可以用与 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))