有没有一种有效的方法来计算 R 中 sf 中多边形的所有成对交集(不使用 for 循环)?
Is there an efficient way to calculate all pairwise intersection of polygons in sf (without using for loops) in R?
我想根据数据框中的变量计算多边形的成对交集。
这是我要找的东西:
set.seed(131)
library(sf)
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)
plot(s.f.2, col = s.f.2$id)
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)
all.ids.pol = unique(s.f.2$id)
df = data.frame(NULL)
for (i in 1:length(all.ids.pol)) {
for (j in 1:length(all.ids.pol)) {
ol1 = st_intersection(s.f.2[c(all.ids.pol[i],all.ids.pol[j]),])
if (dim(ol1)[1]>2 | i == j) {
#add in an area count column (area of each arable poly in intersect layer)
ol1$areaover <- st_area(ol1$geometry)
} else {ol1$areaover = 0}
#run the intersect function, converting the output to a tibble in the process
if (i == j) {
ol1.1 <- as_tibble(ol1)[1,]
} else {ol1.1 <- as_tibble(ol1)[2,] }
id.names = c(all.ids.pol[i],all.ids.pol[j])
my.df =data.frame(area = ol1.1$areaover,id1 = id.names[1],id2 =id.names[2])
df = rbind(df,my.df)
}
}
intersected.areas = df %>%
arrange(id1) %>%
mutate(area.ha = units::set_units(area, ha)) %>%
select(-area) %>%
pivot_wider(names_from = id1, values_from = area.ha) %>%
arrange(id2) %>%
rename(ID = id2)
最后我想要一个矩阵,在行和列中显示“ID”,矩阵中的每个单元格显示多边形之间的交叉区域。
# A tibble: 3 × 4
ID `1` `2` `3`
<dbl> [ha] [ha] [ha]
1 1 1.59 0.501 0
2 2 0.501 1.86 0
3 3 0 0 1
所以代码有效,但我想知道 sf 中是否存在一个函数可以比使用 for 循环更有效地执行此操作(这对于大型数据集来说变得非常不切实际)
library(sf)
library(tidyverse)
s.f.2 %>%
st_intersection(s.f.2) %>%
mutate(intersect_area = st_area(.)) %>%
st_drop_geometry() %>%
pivot_wider(id_cols = id, names_from = id.1,
values_from = intersect_area,
values_fill = 0)
# # A tibble: 3 x 4
# id `1` `2` `3`
# <dbl> <dbl> <dbl> <dbl>
# 1 1 1.59 0.501 0
# 2 2 0.501 1.86 0
# 3 3 0 0 1
我想根据数据框中的变量计算多边形的成对交集。
这是我要找的东西:
set.seed(131)
library(sf)
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)
plot(s.f.2, col = s.f.2$id)
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)
all.ids.pol = unique(s.f.2$id)
df = data.frame(NULL)
for (i in 1:length(all.ids.pol)) {
for (j in 1:length(all.ids.pol)) {
ol1 = st_intersection(s.f.2[c(all.ids.pol[i],all.ids.pol[j]),])
if (dim(ol1)[1]>2 | i == j) {
#add in an area count column (area of each arable poly in intersect layer)
ol1$areaover <- st_area(ol1$geometry)
} else {ol1$areaover = 0}
#run the intersect function, converting the output to a tibble in the process
if (i == j) {
ol1.1 <- as_tibble(ol1)[1,]
} else {ol1.1 <- as_tibble(ol1)[2,] }
id.names = c(all.ids.pol[i],all.ids.pol[j])
my.df =data.frame(area = ol1.1$areaover,id1 = id.names[1],id2 =id.names[2])
df = rbind(df,my.df)
}
}
intersected.areas = df %>%
arrange(id1) %>%
mutate(area.ha = units::set_units(area, ha)) %>%
select(-area) %>%
pivot_wider(names_from = id1, values_from = area.ha) %>%
arrange(id2) %>%
rename(ID = id2)
最后我想要一个矩阵,在行和列中显示“ID”,矩阵中的每个单元格显示多边形之间的交叉区域。
# A tibble: 3 × 4
ID `1` `2` `3`
<dbl> [ha] [ha] [ha]
1 1 1.59 0.501 0
2 2 0.501 1.86 0
3 3 0 0 1
所以代码有效,但我想知道 sf 中是否存在一个函数可以比使用 for 循环更有效地执行此操作(这对于大型数据集来说变得非常不切实际)
library(sf)
library(tidyverse)
s.f.2 %>%
st_intersection(s.f.2) %>%
mutate(intersect_area = st_area(.)) %>%
st_drop_geometry() %>%
pivot_wider(id_cols = id, names_from = id.1,
values_from = intersect_area,
values_fill = 0)
# # A tibble: 3 x 4
# id `1` `2` `3`
# <dbl> <dbl> <dbl> <dbl>
# 1 1 1.59 0.501 0
# 2 2 0.501 1.86 0
# 3 3 0 0 1