从包含相当大数量(约20000)可能部分重叠的多边形的shapefile开始,我需要提取所有通过交叉它们的不同“边界”而产生的子多边形。
在实践中,从一些模拟数据开始:
library(tibble)
library(dplyr)
library(sf)
ncircles <- 9
rmax <- 120
x_limits <- c(-70,70)
y_limits <- c(-30,30)
set.seed(100)
xy <- data.frame(
id = paste0("id_", 1:ncircles),
x = runif(ncircles, min(x_limits), max(x_limits)),
y = runif(ncircles, min(y_limits), max(y_limits))) %>%
as_tibble()
polys <- st_as_sf(xy, coords = c(2,3)) %>%
st_buffer(runif(ncircles, min = 1, max = 20))
plot(polys[1])
我需要派生一个包含ALL的sf
或sp
多边形,并且只需要交叉点生成的多边形,例如:
(请注意,颜色仅用于举例说明预期结果,其中每个“不同颜色”区域是一个单独的多边形,不覆盖任何其他多边形)
我知道我可以通过一次分析一个多边形,识别并保存所有交叉点然后“擦除”那些区域形成完整的多边形并继续循环来解决问题,但这很慢。
我觉得应该有一个更有效的解决方案,但我无法弄清楚,所以任何帮助将不胜感激! (欢迎基于sf
和sp
的解决方案)
更新:
最后,我发现即使是“一次一个多边形”,任务也远非简单!我真的在努力解决这个显而易见的“简单”问题!任何提示?即使是一个缓慢的解决方案或提示在正确的道路上开始将不胜感激!
更新2:
也许这会澄清一些事情:所需的功能类似于这里描述的功能:
更新3:
我将赏金奖给了@ shuiping-chen(谢谢!),他的回答正确地解决了所提供的示例数据集上的问题。然而,“方法”一般被推广到“四倍”或“n-uple”交叉点可能的情况。如果我管理的话,我会在未来几天尝试解决这个问题并发布一个更通用的解决方案!
现在,在使用单个参数(sf或sfc)调用st_intersection
时,R包sf中已实现此默认结果,请参阅https://r-spatial.github.io/sf/reference/geos_binary_ops.html以获取示例。 (我不确定origins
字段是否包含有用的索引;理想情况下,它们应该只指向x
中的索引,现在它们是自我引用的)。
我稍微修改了模拟数据,以说明处理多个属性的能力。
library(tibble)
library(dplyr)
library(sf)
ncircles <- 9
rmax <- 120
x_limits <- c(-70,70)
y_limits <- c(-30,30)
set.seed(100)
xy <- data.frame(
id = paste0("id_", 1:ncircles),
val = paste0("val_", 1:ncircles),
x = runif(ncircles, min(x_limits), max(x_limits)),
y = runif(ncircles, min(y_limits), max(y_limits)),
stringsAsFactors = FALSE) %>%
as_tibble()
polys <- st_as_sf(xy, coords = c(3,4)) %>%
st_buffer(runif(ncircles, min = 1, max = 20))
plot(polys[1])
然后定义以下两个函数。
cur
:基本多边形的当前索引x
:多边形的索引,与cur
相交input_polys
:多边形的简单特征keep_columns
:几何计算后需要保留的属性名称向量get_difference_region()
得到基本多边形和其他相交多边形之间的差异; get_intersection_region()
得到交叉多边形之间的交叉点。
library(stringr)
get_difference_region <- function(cur, x, input_polys, keep_columns=c("id")){
x <- x[!x==cur] # remove self
len <- length(x)
input_poly_sfc <- st_geometry(input_polys)
input_poly_attr <- as.data.frame(as.data.frame(input_polys)[, keep_columns])
# base poly
res_poly <- input_poly_sfc[[cur]]
res_attr <- input_poly_attr[cur, ]
# substract the intersection parts from base poly
if(len > 0){
for(i in 1:len){
res_poly <- st_difference(res_poly, input_poly_sfc[[x[i]]])
}
}
return(cbind(res_attr, data.frame(geom=st_as_text(res_poly))))
}
get_intersection_region <- function(cur, x, input_polys, keep_columns=c("id"), sep="&"){
x <- x[!x<=cur] # remove self and remove duplicated obj
len <- length(x)
input_poly_sfc <- st_geometry(input_polys)
input_poly_attr <- as.data.frame(as.data.frame(input_polys)[, keep_columns])
res_df <- data.frame()
if(len > 0){
for(i in 1:len){
res_poly <- st_intersection(input_poly_sfc[[cur]], input_poly_sfc[[x[i]]])
res_attr <- list()
for(j in 1:length(keep_columns)){
pred_attr <- str_split(input_poly_attr[cur, j], sep, simplify = TRUE)
next_attr <- str_split(input_poly_attr[x[i], j], sep, simplify = TRUE)
res_attr[[j]] <- paste(sort(unique(c(pred_attr, next_attr))), collapse=sep)
}
res_attr <- as.data.frame(res_attr)
colnames(res_attr) <- keep_columns
res_df <- rbind(res_df, cbind(res_attr, data.frame(geom=st_as_text(res_poly))))
}
}
return(res_df)
}
让我们看看模拟数据的差异函数效应。
flag <- st_intersects(polys, polys)
first_diff <- data.frame()
for(i in 1:length(flag)) {
cur_df <- get_difference_region(i, flag[[i]], polys, keep_column = c("id", "val"))
first_diff <- rbind(first_diff, cur_df)
}
first_diff_sf <- st_as_sf(first_diff, wkt="geom")
first_diff_sf
plot(first_diff_sf[1])
first_inter <- data.frame()
for(i in 1:length(flag)) {
cur_df <- get_intersection_region(i, flag[[i]], polys, keep_column=c("id", "val"))
first_inter <- rbind(first_inter, cur_df)
}
first_inter <- first_inter[row.names(first_inter %>% select(-geom) %>% distinct()),]
first_inter_sf <- st_as_sf(first_inter, wkt="geom")
first_inter_sf
plot(first_inter_sf[1])
使用第一级作为输入,并重复相同的过程。
flag <- st_intersects(first_inter_sf, first_inter_sf)
# Second level difference region
second_diff <- data.frame()
for(i in 1:length(flag)) {
cur_df <- get_difference_region(i, flag[[i]], first_inter_sf, keep_column = c("id", "val"))
second_diff <- rbind(second_diff, cur_df)
}
second_diff_sf <- st_as_sf(second_diff, wkt="geom")
second_diff_sf
plot(second_diff_sf[1])
second_inter <- data.frame()
for(i in 1:length(flag)) {
cur_df <- get_intersection_region(i, flag[[i]], first_inter_sf, keep_column=c("id", "val"))
second_inter <- rbind(second_inter, cur_df)
}
second_inter <- second_inter[row.names(second_inter %>% select(-geom) %>% distinct()),] # remove duplicated shape
second_inter_sf <- st_as_sf(second_inter, wkt="geom")
second_inter_sf
plot(second_inter_sf[1])
获取第二级的不同交叉点,并将它们用作第三级的输入。我们可以得到第三级的交叉结果是NULL
,然后该过程应该结束。
我们将所有差异结果放入关闭列表中,并将所有交集结果放入打开列表中。然后我们有:
因此,我们在这里得到最终的代码(应该声明基本的两个函数):
# init
close_df <- data.frame()
open_sf <- polys
# main loop
while(!is.null(open_sf)) {
flag <- st_intersects(open_sf, open_sf)
for(i in 1:length(flag)) {
cur_df <- get_difference_region(i, flag[[i]], open_sf, keep_column = c("id", "val"))
close_df <- rbind(close_df, cur_df)
}
cur_open <- data.frame()
for(i in 1:length(flag)) {
cur_df <- get_intersection_region(i, flag[[i]], open_sf, keep_column = c("id", "val"))
cur_open <- rbind(cur_open, cur_df)
}
if(nrow(cur_open) != 0) {
cur_open <- cur_open[row.names(cur_open %>% select(-geom) %>% distinct()),]
open_sf <- st_as_sf(cur_open, wkt="geom")
}
else{
open_sf <- NULL
}
}
close_sf <- st_as_sf(close_df, wkt="geom")
close_sf
plot(close_sf[1])
不确定它是否对你有帮助,因为它不在R中,但我认为有一种很好的方法可以使用Python解决这个问题。有一个名为GeoPandas(http://geopandas.org/index.html)的库,它允许您轻松地进行地理操作。您需要做的步骤如下:
完整的示例显示在文档中。
操作前 - 2个多边形
手术后 - 9个多边形
如果有什么不清楚的话随时让我知道!希望能帮助到你!