Search code examples
rpolygonintersectionr-sfgtfs

R intersecting multiple (2+) overlapping polygons


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.

enter image description here

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


Solution

  • 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)