如何计算多边形之间的所有成对交互以及 R 中 sf 的百分比覆盖率？
Posted
技术标签:
【中文标题】如何计算多边形之间的所有成对交互以及 R 中 sf 的百分比覆盖率？【英文标题】：How to compute all pairwise interaction between polygons and the the percentage coverage in R with sf? 【发布时间】：2021-12-28 17:26:49 【问题描述】：我有多边形，我想计算它们之间的重叠百分比。这个想法是，当一个多边形与另一个多边形相交时，可以从一个多边形或另一个多边形的角度计算百分比（即最大值）。因此，我想制作一个脚本，生成多边形之间的百分比覆盖率，获取一个多边形和另一个多边形的重叠百分比，然后将所有结果放在一个数据框中。
这是我目前拥有的代码：
set.seed(131)
library(sf)
library(mapview)
m = rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))
p = st_polygon(list(m))
n = 5
l = vector("list", n)
for (i in 1:n)
l[[i]] = p + 2 * runif(2)
s = st_sfc(l)
s.f = st_sf(s)
s.f$id = c(1,1,2,2,3)
s.f.2 = s.f %>% group_by(id) %>% summarise(geometry = sf::st_union(s)) # %>% summarise(area = st_area(s))
s.f.2$area = st_area(s.f.2)
i = s.f.2 %>%
st_intersection(.) %>%
mutate(intersect_area = st_area(.)) #%>%
st_intersection(s.f.2) %>%
mutate(intersect_area = st_area(.),
# id.int = sapply(i$origins, function(x) paste0(as.character(hr.pol$id)[x], collapse = ", ")),
id1 = sapply(i$origins, function(x) paste0(as.character(s.f.2$id)[x][1])),
id2 = sapply(i$origins, function(x) paste0(as.character(s.f.2$id)[x][2])),
area.id1 = sapply(i$origins, function(x) s.f.2$area[x][1]),
area.id2 = sapply(i$origins, function(x) s.f.2$area[x][2]),
perc1 = as.vector(intersect_area/area.id1),
perc2 = as.vector(intersect_area/area.id2)) %>% # create new column with shape area
filter(n.overlaps ==2) %>%
dplyr::select(id, intersect_area, id1, id2,
# id.int,
perc1,perc2) %>% # only select columns needed to merge
st_drop_geometry() %>% # drop geometry as we don't need it
select(-id) %>%
pivot_longer(#names_prefix = "id",
names_to = "perc",
cols = starts_with("perc"))
这段代码给出了多边形之间重叠的百分比（我只为 2 个重叠做，但如果这可以推广到多个重叠，那就太好了！）
mapview(s.f.2,zcol = "id")
最后，我正在寻找的是这样的：
id `1` `2` `3`
1 100 31.6 0
2 27.0 100 0
3 0 0 100
所以多边形“1”占多边形“2”面积的 31.6%，多边形“2”占多边形“1”面积的 27.0%。
我目前拥有的是（但非常很慢）：
data.sp = s.f.2 %>%
st_as_sf(.) %>%
mutate(area.m = st_area(geometry),
area.ha = units::set_units(area.m, ha)) %>%
select(-c(area,area.m))
id.sort = sort(unique(data.sp$id)) # used to reorder columns based on ID
df.fill =data.frame(id1 = NULL, id2=NULL, area =NULL, over1 = NULL, over2 = NULL)
for (k in 1:length(id.sort))
for (op in 1:length(id.sort))
int.out = st_intersection(data.sp[data.sp$id==id.sort[k],],
data.sp[data.sp$id==id.sort[op],])
# int.out
if(nrow(int.out) != 0)
area.tmp = st_area(int.out)#/10000
over1 = area.tmp/int.out$area.ha
over2 = area.tmp/int.out$area.ha.1
else area.tmp = 0;over1=0;over2=0
df.fill.tmp = data.frame(id1 = id.sort[k], id2=id.sort[op],
area = area.tmp,
over1 = over1*100,
over2 = over2*100)
df.fill = rbind(df.fill,df.fill.tmp)
df.fill$over1 = as.numeric(df.fill$over1)
df.fill$over2 = as.numeric(df.fill$over2)
df.fill %>%
select(-c(area, over2)) %>%
pivot_wider(names_from = id2,values_from = over1,
values_fill = 0)
【问题讨论】：
嗨@M。博索莱伊。你的问题很有趣。我试图向您建议一个解决方案（请参阅下面的答案）。我希望它能满足你的需求。干杯。 【参考方案1】：简单的问题，但不是一个明显的答案！我建议的解决方案遵循与您的策略略有不同的策略，并且不涉及任何 for 循环。首先，我开发了一个函数（即area_cover_()），它生成一个交叉表，其中只有那些至少有一个交叉点的多边形。然后，在第二步中，我开发了另一个函数（即add_isolated_poly()），它在第一步生成的交叉表的末尾添加了没有交集的多边形。如果您有许多没有交集的多边形，这将更容易阅读最终表格。所以，请在下面找到代表。
注意：reprex 的输入数据对应于您的 sf 对象 s.f.2 和 area 列
Reprex
1.第一步： 创建一个交叉表，只包含至少有一个交点的多边形（不包括没有交叉点的多边形，这样读取交叉表的效率更高）。为此，我开发了函数area_cover()
area_cover() 函数的代码
library(sf)
library(dplyr)
library(tidyr)
area_cover <- function(x)
x %>%
st_intersection() %>%
filter(n.overlaps>1) %>%
mutate(area_inter = st_area(.)) %>%
unnest(., cols = c(origins, geometry)) %>%
left_join(., as.data.frame(x), by = c("origins" = "id")) %>%
mutate(cover_percent = area_inter/area.y*100) %>%
select(.,origins, area.y, area_inter, cover_percent) %>%
rename("id" = "origins", "area" = "area.y") %>%
st_drop_geometry() %>%
group_by(area_inter) %>%
mutate(poly_X_id = rev(id)) %>%
relocate(poly_X_id, .before = area_inter) %>%
xtabs(cover_percent ~ id + poly_X_id, data = .) %>%
replace(.== 0, 100) %>%
round(., digits = 1)
area_cover() 函数的输出
Results <- area_cover(s.f.2)
Results # Only polygons with at least one intersection are present in this cross table
#> poly_X_id
#> id 1 2
#> 1 100.0 31.6
#> 2 27.0 100.0
class(Results)
#> [1] "xtabs" "table" # you can convert 'Results' into matrix with 'as.matrix()' if needed.
2。第二步（可选）： 在交叉表“结果”（即上一步的结果）的末尾添加孤立的多边形（即没有交叉点）。为此，我开发了函数add_isolated_poly()，它创建了一个数据框，其中的 n 列对应于孤立多边形的 n 个 id，并用 0 填充
add_isolated_poly() 函数的代码
add_isolated_poly <- function(y, z) # 'y arg.' = s.f.2 and 'z arg.' = result of the function area_cover()
id_isolated_poly <- setdiff(y$id, colnames(z))
df_isolated_poly <- y %>%
filter(.,id %in% id_isolated_poly) %>%
st_drop_geometry() %>%
select(., id) %>%
t() %>%
as.data.frame() %>%
`colnames<-`(., id_isolated_poly) %>%
rbind(.,rep(list(rep(0, nrow(y))), length(id_isolated_poly))) %>%
slice(-c(1))
cbind.fill <- function(...)
nm <- list(...)
nm <- lapply(nm, as.matrix)
n <- max(sapply(nm, nrow))
do.call(cbind, lapply(nm, function (x)
rbind(x, matrix(0, n-nrow(x), ncol(x)))))
Results %>%
cbind.fill(., df_isolated_poly) %>%
replace(., col(.) == row(.), 100) %>%
`rownames<-`(., c(colnames(z), id_isolated_poly))
add_isolated_poly() 函数的输出
Results_2 <- add_isolated_poly(s.f.2, Results)
Results_2
#> 1 2 3
#> 1 100 31.6 0
#> 2 27 100.0 0
#> 3 0 0.0 100
class(Results_2)
#> [1] "matrix" "array
由reprex package (v2.0.1) 于 2021 年 11 月 19 日创建
重要编辑
虽然它们在建议的最小示例中产生了正确的结果，但我上面提出的两个函数不可推广并产生错误的结果。经过大量的试验和错误，这是一个简单，非常快速的功能......而且，这一次，对！所以，请在下面找到代表。
解决方案
area_cover() 函数的代码
library(sf)
library(dplyr)
area_cover <- function(x)
Results <- x %>%
st_intersection(.,.) %>%
mutate(area_inter = st_area(.),
cover = area_inter/area.1*100) %>%
st_drop_geometry() %>%
xtabs(cover ~ id.1 + id, data = ., sparse = TRUE) %>%
round(., digits = 1) %>%
as.matrix(.)
names(dimnames(Results)) <- NULL
return(Results)
area_cover() 函数的输出
area_cover(s.f.2)
#> 1 2 3
#> 1 100 31.6 0
#> 2 27 100.0 0
#> 3 0 0.0 100
由reprex package (v2.0.1) 于 2021 年 11 月 22 日创建
基准测试
我根据@M 提供的新数据集将我的函数与@wibeasley 的经过验证的解决方案进行了比较。 Beausoleil（见下面的 cmets）。
为了使比较有效，我稍微修改了@wibeasley 的函数，以便输出是一个包含百分比的矩阵（即与我的函数相同的输出）
数据set.seed(1234)
m = rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))
p = st_polygon(list(m))
n = 100
l = vector("list", n)
for (i in 1:n)
l[[i]] = p + 7 * runif(2)
s = st_sfc(l)
s.f = st_sf(s)
s.f$id = c(1,1,2,2,3,4,5,5,6,7,7,8,9,8,7,3,4,5,5,6)
s.f.2 = s.f %>% group_by(id) %>% summarise(geometry = sf::st_union(s)) # %>% summarise(area = st_area(s))
s.f.2$area = st_area(s.f.2)
library(sf)
library(dplyr)
library(tidyr)
library(bench)
bench_results <- bench::mark(
"wibeasley_validated_answer" = #!!!!! NB: lighlty modified function for the comparison to be valid!!!!!
wibeasley <- sf::st_intersection(s.f.2, s.f.2) %>%
dplyr::mutate(
area = sf::st_area(.),
proportion = area / area.1 * 100
) %>%
tibble::as_tibble() %>%
dplyr::select(
id_1 = id,
id_2 = id.1,
proportion,
) %>%
# tidyr::complete(id_1, id_2, fill = list(proportion = 0))
tidyr::pivot_wider(
names_from = id_1,
values_from = proportion,
values_fill = 0
) %>%
as.matrix(.,rownames.force = TRUE) %>%
`<-`(., .[,-c(1)]) %>%
round(.,1)
,
"lovalery_answer" =
lovalery <- area_cover(s.f.2)
,
min_iterations = 1000,
relative = TRUE,
check = TRUE)
bench_results
#> # A tibble: 2 x 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 wibeasley_validated_answer 1.13 1.10 1 1 1.16
#> 2 lovalery_answer 1 1 1.10 24.1 1
bench_results
#> # A tibble: 2 x 6
#> expression min median `itr/sec` mem_alloc `gc/sec`
#> <bch:expr> <bch:tm> <bch:tm> <dbl> <bch:byt> <dbl>
#> 1 wibeasley_validated_answer 46.4ms 49.1ms 20.2 2.54MB 1.11
#> 2 lovalery_answer 41.5ms 43.7ms 22.7 61.09MB 0.969
all.equal(wibeasley, lovalery)
#> [1] TRUE
wibeasley
#> 1 2 3 4 5 6 7 8 9
#> 1 100.0 19.0 13.3 4.6 22.6 21.2 12.8 3.7 11.6
#> 2 18.3 100.0 28.7 31.9 33.0 14.3 32.2 25.1 5.1
#> 3 14.6 32.7 100.0 23.3 35.5 7.2 28.8 7.0 13.7
#> 4 5.1 36.4 23.3 100.0 20.3 23.7 26.2 26.8 14.7
#> 5 13.0 19.8 18.7 10.7 100.0 10.3 27.0 15.6 12.4
#> 6 21.3 15.0 6.6 21.8 18.0 100.0 26.9 24.4 15.8
#> 7 9.5 24.8 19.5 17.7 34.7 19.9 100.0 11.8 17.9
#> 8 3.8 26.8 6.6 25.1 27.7 24.9 16.3 100.0 1.3
#> 9 22.1 10.0 23.7 25.5 41.0 30.0 46.0 2.3 100.0
lovalery
#> 1 2 3 4 5 6 7 8 9
#> 1 100.0 19.0 13.3 4.6 22.6 21.2 12.8 3.7 11.6
#> 2 18.3 100.0 28.7 31.9 33.0 14.3 32.2 25.1 5.1
#> 3 14.6 32.7 100.0 23.3 35.5 7.2 28.8 7.0 13.7
#> 4 5.1 36.4 23.3 100.0 20.3 23.7 26.2 26.8 14.7
#> 5 13.0 19.8 18.7 10.7 100.0 10.3 27.0 15.6 12.4
#> 6 21.3 15.0 6.6 21.8 18.0 100.0 26.9 24.4 15.8
#> 7 9.5 24.8 19.5 17.7 34.7 19.9 100.0 11.8 17.9
#> 8 3.8 26.8 6.6 25.1 27.7 24.9 16.3 100.0 1.3
#> 9 22.1 10.0 23.7 25.5 41.0 30.0 46.0 2.3 100.0
由reprex package (v2.0.1) 于 2021 年 11 月 22 日创建
【讨论】：
嗨@wibeasley，我不是OP，但是，是的，我认为你误解了目标。所需的结果是您在评论中描述的结果：它对应于相交区域的面积除以每个多边形的总面积。最后，我认为您应该将您的评论放在 OP 的问题或您的答案下，而不是我的 ;-) 干杯。 确实，@lovalery。这种计算背后的想法是重叠区域可能意味着不同的东西，具体取决于所比较的多边形。例如，如果我有 2 只动物的领土（多边形）重叠。根据每只动物的领土大小，这两个领土之间重叠的面积有不同的百分比。因此，有趣的是看看一种动物和另一种动物覆盖了多少（也许一种动物完全重叠了另一种动物的领土，但第二种动物只重叠了另一种动物的一小部分。这有意义吗？ 您好 Beausoleil 先生和@wibeasley，非常感谢您的反馈。说实话，我很难重现你得到的错误，因为在我的电脑上一切正常。也就是说，我想我已经找到了一个解决方法，让它在所有情况下都能顺利运行！我认为@wibeasley 的假设是正确的：问题在于replace_na() 函数需要一个向量或数据帧而不是矩阵，这是cbind.fill() 函数返回的。
所以我通过修改cbind.fill() 函数来编辑我上面的答案，然后不使用replace_na() 函数。请让我知道这现在是否也适合您。干杯。 （注意：我建议的解决方案应该快速运行（只要它有效！））
现在为我工作。【参考方案2】：
如果没有一个真实的示例来进行基准测试，我不确定它是否比您的解决方案更快。但它更简单、更容易理解（至少对我的大脑而言）。
sf::st_intersection() 是矢量化的。所以它会为你找到并返回第一个和第二个参数的所有交集。在这种情况下，这两个参数是同一组多边形。
sf::st_intersection(s.f.2, s.f.2) %>%
dplyr::mutate(
area = sf::st_area(.),
proportion = area / area.1
) %>%
tibble::as_tibble() %>%
dplyr::select(
id_1 = id,
id_2 = id.1,
proportion,
) %>%
# tidyr::complete(id_1, id_2, fill = list(proportion = 0))
tidyr::pivot_wider(
names_from = id_1,
values_from = proportion,
values_fill = 0
)
输出：
# A tibble: 3 x 4
id_2 `1` `2` `3`
<dbl> <dbl> <dbl> <dbl>
1 1 1 0.316 0
2 2 0.270 1 0
3 3 0 0 1
需要考虑的事项：
将面积保留为比例，而不是百分比。通常以后计算比较好。 待久一点，不要转身。通常以后计算会更好，因为您可以加入id_1 和 id_2。
如果您旋转宽度，您可能希望它作为一个矩阵，而不是一个 data.frame。
【讨论】：
我没有这方面的基准，但我可以与以前的基准进行比较，这样更快！谢谢！如果 area 列中有一些单位，请小心：由于某种原因，pivot_wider 不能使用单位...【参考方案3】：我会优化一些东西。我不能真正评估只有三个多边形的优化，我猜交叉计算是唯一真正昂贵的部分，所以我会从那里开始。
如果您不计算 (1) 多边形 A 和 B，然后 (2) 多边形 B 和 A，您将立即减少约 50%。从某种意义上说，只计算 upper triangle，并且重用/反映 lower 三角形的值。
# I think you wanted to create 5 empty columns. Use `numeric(0)` instead of `NULL`
df.fill=data.frame(id1=numeric(0), id2=numeric(0), area=numeric(0), over1=numeric(0), over2=numeric(0))
for (k in seq_along(id.sort))
for (op in seq(from = k, to = length(id.sort), by = 1)) # Avoid the lower triangle
int.out = st_intersection(
data.sp[data.sp$id==id.sort[k],],
data.sp[data.sp$id==id.sort[op],]
)
if(nrow(int.out) != 0)
area.tmp = st_area(int.out)#/10000
over1 = area.tmp/int.out$area.ha
over2 = area.tmp/int.out$area.ha.1
else area.tmp = 0;over1=0;over2=0
df.fill.tmp.upper = data.frame(id1 = id.sort[k], id2=id.sort[op],
area = area.tmp,
over1 = over1,
over2 = over2)
df.fill.tmp.lower = data.frame(id1 = id.sort[op], id2=id.sort[k],
area = area.tmp,
over1 = over2,
over2 = over1)
df.fill <-
if (k == op) rbind(df.fill, df.fill.tmp.upper)
else rbind(df.fill, df.fill.tmp.upper, df.fill.tmp.lower)
df.fill %>%
dplyr::mutate(
over1 = as.numeric(over1) * 100
over2 = as.numeric(over2) * 100
) %>%
select(-area, -over2) %>%
pivot_wider(
names_from = id2,
values_from = over1,
values_fill = 0
)
输出：
# A tibble: 3 x 4
id1 `1` `2` `3`
<dbl> <dbl> <dbl> <dbl>
1 1 100 31.6 0
2 2 27.0 100 0
3 3 0 0 100
【讨论】：
不要介意我在@lovalery 的解决方案上的cmets。我的生产相同的不对称矩阵。我把昨晚的事情搞糊涂了。谢谢你的耐心。交集只需要计算一次，除法有两次。一次由多边形 A，一次由多边形 B。【参考方案4】：如果与其余代码相比，交集计算确实很昂贵，则此解决方案可能比我的第二个解决方案更快。与@lovalery 的解决方案一样，只有一个data.frame 被传递给sf:st_intersection()。根据我的解释，主要区别在于该解决方案会保持较长时间，并使用交叉连接来枚举所有可能的组合，并在最后一步中广泛使用。当事物长而不是宽时（对我来说）更容易操作。
intersections <-
s.f.2 %>%
sf::st_intersection() %>%
dplyr::filter(2L <= n.overlaps) %>%
dplyr::mutate(
numerator = sf::st_area(.)
) %>%
tibble::as_tibble() %>%
tidyr::unnest_wider(origins) %>%
dplyr::select(
id_1 = `...1`,
id_2 = `...2`,
numerator
)
denominators <-
s.f.2 %>%
tibble::as_tibble() %>%
dplyr::select(id, area)
denominators %>%
dplyr::full_join(denominators, by = character()) %>%
dplyr::select(
id_1 = id.x,
id_2 = id.y,
area = area.y
) %>%
dplyr::left_join(intersections, by = c("id_1" = "id_1", "id_2" = "id_2")) %>%
dplyr::left_join(intersections, by = c("id_1" = "id_2", "id_2" = "id_1")) %>%
dplyr::mutate(
numerator = dplyr::coalesce(numerator.x, numerator.y, 0),
proportion = dplyr::if_else(id_1 == id_2, 1, numerator / area),
) %>%
dplyr::select(id_1, id_2, proportion) %>%
tidyr::pivot_wider(
names_from = id_1,
values_from = proportion,
values_fill = 0
)
输出：
# A tibble: 3 × 4
id_2 `1` `2` `3`
<dbl> <dbl> <dbl> <dbl>
1 1 1 0.316 0
2 2 0.270 1 0
3 3 0 0 1
【讨论】：
以上是关于如何计算多边形之间的所有成对交互以及 R 中 sf 的百分比覆盖率？的主要内容，如果未能解决你的问题，请参考以下文章
如何在 Pyspark Dataframe 中创建多列的所有成对组合？