Simple question but not an obvious answer! The solution I suggest follows a slightly different strategy than yours and does not involve any for loops. First, I developed a function (i.e. area_cover_()) that produces a cross table with only those polygons that have at least one intersection. Then, in a second step, I developed another function (i.e. add_isolated_poly()) that adds the polygons with no intersection at the end of the cross table produced in the first step. This makes it easier to read the final table if you have many polygons without intersection. So, please find the reprex below.
NB: The input data for the reprex corresponds to your sf object s.f.2 with area column
Reprex
1. First step: create a cross table including only the polygons that have at least one intersection (not including the polygons without intersection makes the reading of the cross table more efficient). To do this, I developed the function area_cover()
- Code of the
area_cover() function
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)
}
- Output of the
area_cover() function
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. Second step (optional): add isolated polygons (i.e. with no intersection) at the end of the cross table 'Results' (i.e. result of the previous step). To do this, I developed the function add_isolated_poly() which creates a dataframe with n columns corresponding to the n id's of the isolated polygons and filled with 0s
- Code of the
add_isolated_poly() function
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))
}
- Output of the
add_isolated_poly() function
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
Created on 2021-11-19 by the reprex package (v2.0.1)
IMPORTANT EDIT
Although they produce the right result with the proposed minimal example, the two functions I proposed above are not generalizable and produce wrong results. After a lot of trial and error, here is a simple, very fast function... and, this time, right!!! So, please find the reprex below.
Solution
- Code of the
area_cover() function
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)
}
- Output of the
area_cover() function
area_cover(s.f.2)
#> 1 2 3
#> 1 100 31.6 0
#> 2 27 100.0 0
#> 3 0 0.0 100
Created on 2021-11-22 by the reprex package (v2.0.1)
Benchmarking
I compared my function with the validated solution of @wibeasley based on the new dataset provided by @M. Beausoleil (see comments below).
For the comparison to be valid, I have slightly modified the @wibeasley's function so that the output is a matrix containing percentages (i.e. same output as my function)
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)
- Code to compare the two solutions
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)
- Results of the benchmarking (relative values)
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
- Results of the benchmarking (absolute values)
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
- Final checking that the two functions produce the same result
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
Created on 2021-11-22 by the reprex package (v2.0.1)
