R intersecting multiple (2+) overlapping polygons

Question

I am having some issues with st_intersection. I am trying to intersect polygons from bus stop buffers in preparation for areal interpolation. Here are the data: Here is the data: https://realtime.commuterpage.com/rtt/public/utility/gtfs.aspx

Here is my code:

ART2019Path <- file.path(GTFS_path, "2019-10_Arlington.zip")
ART2019GTFS <- read_gtfs(ART2019Path)
ART2019StopLoc <- stops_as_sf(ART2019GTFS$stops) ### Make a spatial file for stops
ART2019Buffer <- st_buffer(ART2019StopLoc, dist = 121.92) ### Make buffer with 400ft (121.92m) radius

It creates something that looks like the image below (created using mapview); as you can see, there are multiple overlapping buffers.

I tried intersecting the polygons using the following:

BufferIntersect <- st_intersection(ART2019Buffer, ART2019Buffer)
BufferIntersect <- st_make_valid(BufferIntersect) ### Fix some of the polygons that didn't quite work

But it only intersects two layers of polygons, meaning there is still overlap. How do I make all buffers intersect?

I have looked at similar questions on here like this: Loop to check multiple polygons overlap in r SF package

But there is no answer.

One of the comments suggested the following links:

https://r-spatial.org/r/2017/12/21/geoms.html

https://r-spatial.github.io/sf/reference/geos_binary_ops.html#details-1

But I can't get either to work. Any help would be greatly appreciated.

---

Edit

Couple of clarifying points in response to some comments.

1. I am interested in the area of each unique polygon within bus stop buffers as I will be using these polygons in an areal interpolation with census data to estimate population with access to bus stops

2. 400ft walking distance is standard practice for bus stop accessibility

Answer 1

It sounds like you just want the buffer(s), but without having to deal with all of the overlapping sections. It doesn't matter if a person is within 400ft of one bus-stop or three, right?

If so, you can use the st_union function to "blend" the buffers together.

library(tidytransit)
library(sf)
library(mapview)
library(ggplot2)

# s2 true allows buffering in meters, s2 off later speeds things up
sf::sf_use_s2(TRUE)
ART2019Path <- file.path("/your/file/path/")
ART2019GTFS <- read_gtfs(ART2019Path)
ART2019StopLoc <- stops_as_sf(ART2019GTFS$stops) ### Make a spatial file for stops
ART2019Buffer <- st_buffer(ART2019StopLoc, dist = 121.92) ### Make buffer with 400ft (121.92m) radius

# might be needed due to some strange geometries in buffer, and increase speed
sf::sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off

# MULTIPOLYGON  sfc object covering only the buffered areas,
#  there are no 'overlaps'.
buff_union <- st_union(st_geometry(ART2019Buffer))
#> although coordinates are longitude/latitude, st_union assumes that they are planar

buff_union
#> Geometry set for 1 feature 
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -77.16368 ymin: 38.83828 xmax: -77.04768 ymax: 38.9263
#> Geodetic CRS:  WGS 84
#> MULTIPOLYGON (((-77.08604 38.83897, -77.08604 3...

# Non-overlapping buffer & stops
ggplot() + 
  geom_sf(data = buff_union, fill = 'blue', alpha = .4) +
  geom_sf(data = ART2019StopLoc, color = 'black') + 
  coord_sf(xlim = c(-77.09, -77.07),
           ylim = c(38.885, 38.9)) 

# Overlapping buffer & stops
ggplot() +
  geom_sf(data = ART2019Buffer, fill = 'blue', alpha = .4) +
  geom_sf(data = ART2019StopLoc, color = 'black') +
  coord_sf(xlim = c(-77.09, -77.07),
           ylim = c(38.885, 38.9)) 

# Back to original settings
sf::sf_use_s2(TRUE)
#> Spherical geometry (s2) switched on

^{Created on 2022-04-18 by the reprex package (v2.0.1)}

Answer 2

Something like this works for me, albeit a bit slowly. Here I loop through each stop buffer and run the intersection process on an object containing all other stop buffers excluding that stop buffer.

library(sf)
library(tidyverse)

df<-read.csv("YOUR_PATH/google_transit/stops.txt")

# Read data
ART2019StopLoc <- st_as_sf(df, coords=c('stop_lon', 'stop_lat'))
ART2019StopLoc <- st_set_crs(ART2019StopLoc, value=4326)

# Make buffer
ART2019Buffer <- st_buffer(ART2019StopLoc, dist=121.92)

# Create empty data frame to store results
results <- data.frame()

# Loop through each stop and intersect with other stops
for(i in 1:nrow(ART2019Buffer)) {
  
  # Subset to stop of interest
  stop <- ART2019Buffer[i,]

  # Subset to all other stops excl. stop of interest
  stop_check <- ART2019Buffer[-i,]
  
  # Intersect and make valid
  stop_intersect <- st_intersection(stop, stop_check) %>% 
    st_make_valid()
  
  # Create one intersected polygon
  stop_intersect <- st_combine(stop_intersect) %>% 
    st_as_sf() %>% 
    mutate(stop_name=stop$stop_name)
  
  # Combine into one results object
  results <- rbind(results, stop_intersect)
  print(i)
  
}

ggplot() +
  geom_sf(data=ART2019Buffer %>% filter(stop_name %in% results$stop_name),
          fill='gray70') +
  geom_sf(data=results, aes(fill=stop_name), alpha=0.5)

The plot below shows the results for the first 8 stops. The gray circles are the original stop buffers and the colored buffers show the intersection with adjacent buffers.