将操作系统国家网格名称/代码添加到R中的网格

Question

将操作系统国家网格名称/代码添加到R中的网格

3

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

我可以轻松创建四个级别的网格，没有任何问题（使用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的单元格，以反映图像中所示的命名/编码约定？

- Chris

2个回答

0

如果您的目标是将点与网格代码相关联（根据您的后续问题），您也可以即时计算代码。

这里是一个执行此操作的函数

getG <- function(pnt) {
    numcodes <- function(x, cell) {
        rc <- rowColFromCell(x, cell)
        c(rc[2]-1, 10-rc[1])
    }
    r <- rast(ext=c(-1000000, 1500000, -500000, 2000000), crs="EPSG:27700", res=500000)
    LTS <- LETTERS[-9]

    # first level
    xy <- rbind(pnt)
    c1 <- cellFromXY(r, xy)

    # second level
    r <- disagg(r[c1, drop=FALSE], 5)
    c2 <- cellFromXY(r, xy)

    # third level
    r <- disagg(r[c2, drop=FALSE], 10)
    cxy <- cellFromXY(r, xy)
    cr3 <- numcodes(r, cxy)

    # fourth level (no need for disagg)
    ext(r) <- ext(r[cxy, drop=FALSE])
    cxy <- cellFromXY(r, xy)
    cr4 <- numcodes(r, cxy)

    paste0(LTS[c1], LTS[c2], cr3[1], cr4[1], cr3[2], cr4[2], collapse="")
}

调用函数

library(terra)
pnts <- rbind(cbind(93555, 256188), 
              cbind(210637, 349798),
              cbind(696457,481704))

getG(pnts[1,])
#[1] "SL9356"

apply(pnts, 1, getG)
#[1] "SL9356" "SH1049" "TB9681"

现在，这个功能也可以在geosphere v1.5-19中使用。

library(geosphere)
pnts <- rbind(cbind(93555, 256188), 
      cbind(210637, 349798),
      cbind(696457,481704))

OSGB(pnts, "1km")
#[1] "SL9356" "SH1049" "TB9681"

geosphere v1.5-19目前是开发版本。你可以使用以下命令进行安装

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

- Robert Hijmans

谢谢。我认为相同的逻辑可以应用于分配1公里以下的单元格？ - Chris

这是我的想法，我稍微整理了一下函数，以便更容易地实现。 - Robert Hijmans

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Allan Cameron · Accepted Answer

我想不会有太多人愿意花费20分钟的时间来生成6125000个多边形以测试他们的答案。恰巧我不得不升级我的基于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))