使用 terra 合并来自栅格的对角接触多边形

Question

考虑以下二进制栅格，其中具有各种形状的面片：

使用

terra

为每个面片获取单独的多边形，我通过

as.polygons()

转换为矢量，然后通过

disagg()

分解。删除值为 0 的多边形后，结果是这样的：

我希望通过单个共享顶点对角连接的补丁具有相同的几何形状。所以最终的结果是：

我的工作解决方案是创建一个数据框来存储多边形 ID 以及一个新字段，该字段是最终合并/聚合的值。然后迭代原始分解的向量数据：

1. 选择一个多边形
2. 获取其范围并以最小程度扩展它
3. 根据这个扩展范围裁剪原始矢量数据
4. 查看哪些几何图形接触原始几何图形
5. 将这些几何图形的合并值与原始值更新为相同

我发现裁剪矢量数据 - 然后仅检查这些“闭合”多边形之间的接触 - 比检查所有多边形是否接触当前孤立的多边形运行得更快。

有没有办法加快这个过程？或者我缺少一个可以完成所有这些但速度更快的功能？我确信有很多解决方案，但最终用途是 3k+ 矢量文件，每个矢量文件都有几千个多边形可以合并。

要重现的代码和我的解决方案，在这里：

### setup    
library(terra)

m <- matrix(
  c(1,1,0,1,1,0,1,
    1,0,0,1,1,0,0,
    0,0,0,0,0,1,0,
    0,1,0,0,0,0,0,
    0,1,1,0,1,1,1,
    0,0,0,1,0,0,0,
    0,0,0,0,0,1,1),
  nrow=7,
  ncol=7,
  byrow=TRUE
)

r <- rast(m)
plot(r)

v <- as.polygons(r)

v <- v[which(v$lyr.1 == 1),]
plot(v, "lyr.1")

v <- disagg(v)
v$id <- 1:nrow(v)
plot(v)
text(v, labels="id")


### solution
nv <- nrow(v)

mv_df <- data.frame(id=1:nv, mval=1:nv)

for (i in 1:nv) {
  i_poly <- v[i,]
  i_ext <- extend(ext(i_poly), 1)
  i_crop <- crop(v, i_ext)
  i_touch <- is.related(i_crop, i_poly, relation="touches")
  i_update <- append(i_crop$id[i_touch], i_poly$id)
  mv_df$mval[i_update] <- min(mv_df$mval[i_update])
}

v <- merge(v, mv_df, by.x="id", by.y="id")

v <- aggregate(v, "mval")
plot(v, "mval")
text(v, labels="mval")

Answer 1

使用基于图形的方法，您可以将其速度提高 7 倍：

library(terra)
library(igraph)
library(microbenchmark)

# Create the raster and polygons
m <- matrix(
  c(1,1,0,1,1,0,1,
    1,0,0,1,1,0,0,
    0,0,0,0,0,1,0,
    0,1,0,0,0,0,0,
    0,1,1,0,1,1,1,
    0,0,0,1,0,0,0,
    0,0,0,0,0,1,1),
  nrow = 7,
  ncol = 7,
  byrow = TRUE
)

r <- rast(m)
v <- as.polygons(r)
v <- v[v$lyr.1 == 1, ]
v <- disagg(v)
v$id <- 1:nrow(v)

# Function 1: Iterative approach
iterative_approach <- function(v) {
  nv <- nrow(v)
  mv_df <- data.frame(id = 1:nv, mval = 1:nv)
  
  for (i in 1:nv) {
    i_poly <- v[i, ]
    i_ext <- extend(ext(i_poly), 1)
    i_crop <- crop(v, i_ext)
    i_touch <- is.related(i_crop, i_poly, relation = "touches")
    i_update <- append(i_crop$id[i_touch], i_poly$id)
    mv_df$mval[i_update] <- min(mv_df$mval[i_update])
  }
  
  v <- merge(v, mv_df, by.x = "id", by.y = "id")
  v <- aggregate(v, "mval")
  return(v)
}
v1 <- iterative_approach(v)
plot(v1, "mval")
text(v1, labels="mval")

# Function 2: Graph-based approach
graph_based_approach <- function(v) {
  touch_matrix <- relate(v, relation = "touches")
  graph <- graph_from_adjacency_matrix(touch_matrix, mode = "undirected")
  components <- components(graph)$membership
  v$group <- components
  v_aggregated <- aggregate(v, by = "group")
  return(v_aggregated)
}

v2 <- graph_based_approach(v)
plot(v2,"group", main = "Aggregated Polygons")



results <- microbenchmark(
  iterative = iterative_approach(v),
  graph_based = graph_based_approach(v),
  times = 10  # Adjust the number of repetitions as needed
)

boxplot(results)