将 OS 国家电网 names/codes 添加到 R 中的网格

Question

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

我可以毫无问题地创建 4 个级别的网格（对于 36GB 内存的我，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 的单元格以反映图像中显示的 naming/coding 约定？

Answer 1

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

不幸的是，当您创建网格方块时，它们是从下到上、从左到右排列的，而标签是从上到下、左排列的 - 向右。这使得标签有点棘手。对于 500km 框，我们希望能够用 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))

将 OS 国家电网 names/codes 添加到 R 中的网格

Add OS National Grid names/codes to grid in R

r

spatial

sf