0

I would like to get the maximum percentage polygons of a sf object are overlapped by other polygons of the same sf object. The idea is to then subset the sf object for polygons which are overlapped by e.g. less than 30%.

Example data:

set.seed(5)

poly <- st_sfc(st_polygon(list( cbind(c(0,100,100,0,0),c(0,0,100,100,0))
                              )))

points <- st_sample(poly, size = 20, type = "random" )

buffpoints <- st_buffer(points, dist = 10)

plot(poly)
plot(points, add = T)
plot(buffpoints, add = T)

I found this post. However, I am currently not able to use it for my purpose.

How to compute all pairwise interaction between polygons and the the percentage coverage in R with sf?

It feels like there should be a relatively simple solution, however, I can´t get my brain wrapped around it...

RasK
  • 11
  • 6
  • Hi. Just to be clear, you want to retain only those polygons that share at least a certain percentage of their area with any other polygon in a list of polygons, right? – agila Apr 06 '22 at 17:53
  • @agila: Less than a certain percentage was the idea. But yes. See the solution based on the answer by mrhellmann – RasK Apr 07 '22 at 06:50
  • @RasK Poor form to post my working solution in your question. – mrhellmann Apr 07 '22 at 10:37
  • @mrhellmann: Sorry. Thought it would be a bad way to edit your solution in your answer as it might be working for others... Therefore I also mentioned you and your answer. I deleted the edit in my question... – RasK Apr 07 '22 at 12:20
  • @RasK Are you having a problem with the solution below? It works as posted for me. From you previous edit, it looks like you missed the step of creating `buff_sf` from `buffpoints`. – mrhellmann Apr 07 '22 at 12:32

1 Answers1

0

Below is (an admittedly ugly) solution using a for loop. It checks each circle from 'buffpoints' for the total amount of area overlapped by other circles in 'buffpoints'.

There's likely a way to do this with purrr, similar to this solution to a similar problem, but I couldn't quite get there. Hopefully someone else can.

library(sf)
library(tidyverse)

##### example data
set.seed(5)
poly <- st_sfc(st_polygon(list( cbind(c(0,100,100,0,0),c(0,0,100,100,0))
)))
points <- st_sample(poly, size = 20, type = "random" )
buffpoints <- st_buffer(points, dist = 10)
#####

# make buffpoints an sf object & add an id column
buff_sf <- buffpoints %>% st_as_sf() %>% mutate(id = row_number())

# for each row(polygon), find the area covered by the union of buff_sf, except for that row(polygon).  
for(i in 1:nrow(buff_sf)){
 covered_area <- buff_sf[i,] %>% 
                   st_intersection(., st_union(buff_sf[-i,])) %>% 
                   st_area()
  #fix problem with numeric(0)
 buff_sf[i, 'covered'] <- ifelse(length(covered_area) == 0, 0, covered_area)
}

#Add percent covered
buff_sf <- buff_sf %>% 
            mutate(area = st_area(.), coverage_pct = covered / area)

#Plotting to eyeball coverage, making break points at the coverage %
#  you're looking for will help
ggplot(buff_sf) + 
  geom_sf(aes(fill = coverage_pct), alpha = .2) + 
  scale_fill_viridis_c()

Looks about right. Highly covered polygons are yellow, uncovered polygons are purple. The most covered polygon is at ~74%, 4 are covered 0% and by eye that looks to be the case.

Created on 2022-04-06 by the reprex package (v2.0.1)

>head(buff_sf)
#Simple feature collection with 6 features and 4 fields
#Geometry type: POLYGON
#Dimension:     XY
#Bounding box:  xmin: 0.4650128 ymin: 3.99837 xmax: 101.6876 ymax: 99.02071
#CRS:           NA
#                               x id     covered     area coverage_pct
#1 POLYGON ((30.02145 89.02071...  1 104.5753253 314.0157 0.3330257463
#2 POLYGON ((78.52186 72.0701,...  2   0.0000000 314.0157 0.0000000000
#3 POLYGON ((101.6876 21.13403...  3 158.5424259 314.0157 0.5048868799
#4 POLYGON ((38.43995 22.57173...  4 195.9680109 314.0157 0.6240706676
#5 POLYGON ((20.46501 13.99837...  5   0.0968619 314.0157 0.0003084619
#6 POLYGON ((80.10575 47.99139...  6  54.2027628 314.0157 0.1726116124

And filtering based on % coverage:

buff_sf %>% filter(coverage_pct < .3)
mrhellmann
  • 5,069
  • 11
  • 38
  • Your solution almost did the job. However, due to an error: "Error in UseMethod("st_intersection"): not applicable method for 'st_intersection' applied to object of class "c('tbl_df', 'tbl', 'data.frame')" I had to change the code slightly. See above – RasK Apr 07 '22 at 06:43
  • @RasK Looks like you applied my solution to a data.frame, not an sf object. The solution I posted works. – mrhellmann Apr 07 '22 at 10:36