将OS国家网格名称/代码添加到R中的网格中

Question

bounty到期在3天内。这个问题的答案有资格获得+ 50声望赏金。 Chris想要引起更多关注对此问题。

我正在尝试在R中重新创建完整的军械测量局国家网格（如此处https://upload.wikimedia.org/wikipedia/commons/f/f5/Ordnance_Survey_National_Grid.svg所示。

我可以毫无问题地为4个级别创建网格（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()

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

Answer 1

[我想没有太多的人想在20分钟内将计算机捆绑起来以生成6,125,000个多边形，以测试他们的答案。碰巧的是，我的基于云的Windows服务器在30分钟后就放弃了，抱怨它的4GB内存已用完，所以我只能向您展示在500km，100km和10km网格中什么有效。希望您也可以将这些方法应用于最小的网格。

不幸的是，当您创建网格正方形时，它们是从下到上，从左到右排序的，而标签是由从上到下，从左到右排序的。这样会使标签更加复杂。对于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公里的盒子，但是如果您被卡住了，您可能需要比我大PC的人来帮助您...