R | sf: I have points and 2 buffers around each point. How to combine the points into single units if larger buffer overlaps (but smaller does not)?

Question

I currently have points on a map (schools). Each has two buffers around the point (school). One that is 450m and one that is 250m. If points are overlapping, I want to think of them as a single unit (because otherwise things get complicated), but I want them to keep their geometries/area that they cover.

So on the sample map given here, I would want the top three schools/points to be combined into one "unit." I want them to keep the area that they cover, but just count in R as a single unit. If I use the "st_union" function, I have to do it for both the smaller and the larger buffers. But this leads to a difference in the total number of units, because more of the larger buffers overlap. So basically I want the smaller buffer to combine with the other smaller buffers (into one unit) if the larger buffers overlap and therefore combine.

The ultimate goal is that I want to count how many times a store falls within 250m from a school compared to 450m. Its easiest, I think, to combine as described above, because of the vast overlap between schools.

I have sample data here.

The code I have is below, but note that it results in a differential number of schools/units due to the above explained reason.

county.sf <- get_acs(state = "MO",
                     county = c( "St. Louis City"),
                     geography = "tract",
                     variables = "B03002_001", 
                     output="wide", 
                     geometry = TRUE) %>%
  sf::st_transform(crs = "ESRI:102003")
  
class(county.sf)

# School data
school <- read.csv("C:\\myfile1.csv")
school.sf <- st_as_sf(school, coords = c("long", "lat"), crs = "epsg:4326") 
school.sf.utm <- st_transform(school.sf, crs = "ESRI:102003")


# Store data
store <- import("C:\\myfile2.csv")
store.sf <- st_as_sf(store, coords = c("XCoord", "YCoord"), crs = "ESRI:102696") 
store.sf.utm <- st_transform(store.sf, crs = "ESRI:102003")


elem.buff <-st_buffer(school.sf.utm, 250)     
elem.buff2 <-st_buffer(school.sf.utm, 450) 

pts_com<-st_union(elem.buff)
pts_pol<-st_cast(pts_com, "POLYGON")

pts_com2<-st_union(elem.buff2)
pts_pol2<-st_cast(pts_com2, "POLYGON")


#unmerged zone map
ex.map<- tm_shape(county.sf) +
  tm_polygons() + 
  
  tm_shape(elem.buff) +
  tm_borders(col="red") +  
  
  tm_shape(school.sf.utm) +
  tm_dots(col = "red") +
  
  tm_shape(elem.buff2) +
  tm_borders(col="blue") + 
    
  tm_shape(pts_pol) +
  tm_borders(col="black") +
  
  tm_shape(store.sf.utm) +
  tm_dots() 
ex.map




#merged zones map

ex.map<- tm_shape(county.sf) +
  tm_polygons() + 
  
  #(elem.buff) +
  #tm_borders(col="red") +  
  
  tm_shape(school.sf.utm) +
  tm_dots(col = "red") +
  
  #tm_shape(elem.buff2) +
  #tm_borders(col="blue") + 
  
  tm_shape(pts_pol) +
  tm_borders(col="red") +
  
  tm_shape(store.sf.utm) +
  tm_dots() +

  tm_shape(pts_pol2) +
  tm_borders(col="blue")
ex.map



lengths(st_intersects(elem.buff, store.sf.utm)) #gives per school but ignores overlapping
lengths(st_intersects(pts_pol, store.sf.utm))

lengths(st_intersects(elem.buff2, store.sf.utm)))
lengths(st_intersects(pts_pol2, store.sf.utm)))

I would suggest using `st_intersection()` on the overlapping buffer polygons to find where you have overlaps (and which original buffers/schools) contribute to the overlapping zones. This answer https://stackoverflow.com/a/69293465/13478749 describes the process (for a different problem), and would be a good start for this question. — hugh-allan, Nov 04 '21 at 04:32
I understand how to do that, but what is less clear to me is how to combine the smaller buffers into one (even though they don't intersect) if and only if the larger buffers do intersect. I imagine it is some kind of if/then statement, but I can't figure it out in sf. — revere2323, Nov 04 '21 at 13:41

score 2 · Accepted Answer · answered Nov 05 '21 at 05:41

Initially I suggested using st_intersection(), but that started getting complicated so I went with this method instead which hinges on st_join().

First, load libraries and data.

library(readxl)
library(tidyverse)
library(sf)
library(tmap)
tmap_mode('view')

read_xlsx('Schools and Stores_all.xlsx', sheet = 1) %>% 
  st_as_sf(., coords = c("long", "lat"), crs = "epsg:4326") %>% 
  st_transform(crs = "ESRI:102003") %>% 
  {. ->> schools}

schools

# Simple feature collection with 156 features and 1 field
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 208163.7 ymin: -81868.36 xmax: 565039.3 ymax: 151922.7
# Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic
# # A tibble: 156 x 2
#    School                                               geometry
# *  <chr>                                             <POINT [m]>
# 1  AcademyOf Envt Sci/math Elementary School   (497569.1 143116)
# 2  AcademyOf Envt Sci/math Middle School     (498610.1 140067.7)
# 3  Adams Elementary School                   (495000.6 141023.2)
# 4  Ames Visual/perf. Arts                    (500120.2 144483.3)
# 5  Ashland Elementary And Br.                (496514.9 146392.7)
# 6  Aspire Academy                              (495458 149014.6)
# 7  Beaumont Cte High School                  (497748.6 145417.7)
# 8  Bryan Hill Elementary School              (495926.2 144709.6)
# 9  Buder Elementary School                   (492994.6 136888.8)
# 10 Busch Ms Character Athletics              (491906.9 135569.1)
# # ... with 146 more rows

Then, make the buffer zones around the schools.

schools %>% 
  st_buffer(250) %>% 
  {. ->> schools_250}

schools_250

# Simple feature collection with 156 features and 1 field
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: 207913.7 ymin: -82118.36 xmax: 565289.3 ymax: 152172.7
# Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic
# # A tibble: 156 x 2
#    School                                                                                           geometry
# *  <chr>                                                                                       <POLYGON [m]>
# 1  AcademyOf Envt Sci/math Elementary School ((497819.1 143116, 497818.8 143103, 497817.7 143089.9, 497816 ~
# 2  AcademyOf Envt Sci/math Middle School     ((498860.1 140067.7, 498859.7 140054.7, 498858.7 140041.6, 498~
# 3  Adams Elementary School                   ((495250.6 141023.2, 495250.3 141010.1, 495249.3 140997, 49524~
# 4  Ames Visual/perf. Arts                    ((500370.2 144483.3, 500369.9 144470.2, 500368.9 144457.2, 500~
# 5  Ashland Elementary And Br.                ((496764.9 146392.7, 496764.6 146379.6, 496763.6 146366.5, 496~
# 6  Aspire Academy                            ((495708 149014.6, 495707.6 149001.5, 495706.6 148988.4, 49570~
# 7  Beaumont Cte High School                  ((497998.6 145417.7, 497998.3 145404.7, 497997.3 145391.6, 497~
# 8  Bryan Hill Elementary School              ((496176.2 144709.6, 496175.9 144696.5, 496174.9 144683.5, 496~
# 9  Buder Elementary School                   ((493244.6 136888.8, 493244.2 136875.7, 493243.2 136862.7, 493~
# 10 Busch Ms Character Athletics              ((492156.9 135569.1, 492156.6 135556, 492155.5 135543, 492153.~
# # ... with 146 more rows

schools %>% 
  st_buffer(450) %>% 
  {. ->> schools_450}

schools_450

# Simple feature collection with 156 features and 1 field
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: 207713.7 ymin: -82318.36 xmax: 565489.3 ymax: 152372.7
# Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic
# # A tibble: 156 x 2
#    School                                                                                           geometry
# *  <chr>                                                                                       <POLYGON [m]>
# 1  AcademyOf Envt Sci/math Elementary School ((498019.1 143116, 498018.5 143092.5, 498016.7 143069, 498013.~
# 2  AcademyOf Envt Sci/math Middle School     ((499060.1 140067.7, 499059.5 140044.2, 499057.6 140020.7, 499~
# 3  Adams Elementary School                   ((495450.6 141023.2, 495450 140999.6, 495448.2 140976.1, 49544~
# 4  Ames Visual/perf. Arts                    ((500570.2 144483.3, 500569.6 144459.8, 500567.8 144436.3, 500~
# 5  Ashland Elementary And Br.                ((496964.9 146392.7, 496964.3 146369.1, 496962.5 146345.6, 496~
# 6  Aspire Academy                            ((495908 149014.6, 495907.4 148991, 495905.5 148967.5, 495902.~
# 7  Beaumont Cte High School                  ((498198.6 145417.7, 498198 145394.2, 498196.2 145370.7, 49819~
# 8  Bryan Hill Elementary School              ((496376.2 144709.6, 496375.6 144686.1, 496373.8 144662.6, 496~
# 9  Buder Elementary School                   ((493444.6 136888.8, 493444 136865.2, 493442.1 136841.8, 49343~
# 10 Busch Ms Character Athletics              ((492356.9 135569.1, 492356.3 135545.5, 492354.4 135522, 49235~
# # ... with 146 more rows

Join all the 450 m buffer zones using st_union(), to determine the identities of the "units".

schools_450 %>% 
  st_union %>% 
  st_cast('POLYGON') %>% 
  st_sf %>% 
  mutate(
    unit = row_number()
  ) %>% 
  {. ->> schools_450_units}

schools_450_units

# Simple feature collection with 39 features and 1 field
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: 207713.7 ymin: -82318.36 xmax: 565489.3 ymax: 152372.7
# Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic
# First 10 features:
#                          geometry unit
# 1  POLYGON ((565450.4 -82051.3...    1
# 2  POLYGON ((499821.5 138524.2...    2
# 3  POLYGON ((498299.1 138974.4...    3
# 4  POLYGON ((500735.6 142090.4...    4
# 5  POLYGON ((499872.4 142927.1...    5
# 6  POLYGON ((496870.2 135641.5...    6
# 7  POLYGON ((495561.7 135210, ...    7
# 8  POLYGON ((498291.9 141786.4...    8
# 9  POLYGON ((496337.3 144526.6...    9
# 10 POLYGON ((208574.8 -53382.3...   10

Notice that from 156 schools, we have only 39 units. And of course the units can be quite far reaching when multiple 450 m buffers are chained together - see map at the end.

Then, we use st_join() to determine which schools (or their buffer zones) are inside each unit. Specifying join = st_within is the key here.

schools %>% 
  st_join(schools_450_units, join = st_within) %>% 
  {. ->> schools_units}

schools_units

# Simple feature collection with 156 features and 2 fields
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 208163.7 ymin: -81868.36 xmax: 565039.3 ymax: 151922.7
# Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic
# # A tibble: 156 x 3
#    School                                               geometry  unit
# *  <chr>                                             <POINT [m]> <int>
# 1  AcademyOf Envt Sci/math Elementary School   (497569.1 143116)     5
# 2  AcademyOf Envt Sci/math Middle School     (498610.1 140067.7)     3
# 3  Adams Elementary School                   (495000.6 141023.2)     3
# 4  Ames Visual/perf. Arts                    (500120.2 144483.3)     5
# 5  Ashland Elementary And Br.                (496514.9 146392.7)    32
# 6  Aspire Academy                              (495458 149014.6)    37
# 7  Beaumont Cte High School                  (497748.6 145417.7)    33
# 8  Bryan Hill Elementary School              (495926.2 144709.6)     9
# 9  Buder Elementary School                   (492994.6 136888.8)    18
# 10 Busch Ms Character Athletics              (491906.9 135569.1)    13
# # ... with 146 more rows

Then to combine the 250 m buffers of the same unit (even if they aren't touching) into a single MULTIPOLYGON, we use group_by(unit) and use summarise(do_union = TRUE) (the default). This is like st_union(), but respects group_by() (it may in fact be possible to get exactly the same result using st_union(), but this works too).

schools_units %>% 
  `st_geometry<-`(schools_250 %>% st_geometry) %>% 
  group_by(unit) %>% 
  summarise(do_union = TRUE) %>% 
  {. ->> schools_250_units}

schools_250_units

# Simple feature collection with 39 features and 1 field
# Geometry type: GEOMETRY
# Dimension:     XY
# Bounding box:  xmin: 207913.7 ymin: -82118.36 xmax: 565289.3 ymax: 152172.7
# Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic
# # A tibble: 39 x 2
#    unit                                                                                             geometry
#   <int>                                                                                       <GEOMETRY [m]>
# 1     1 POLYGON ((565289.3 -81868.36, 565289 -81881.45, 565288 -81894.49, 565286.2 -81907.47, 565283.9 -8...
# 2     2 MULTIPOLYGON (((499659 138681.1, 499657.3 138668.1, 499654.9 138655.3, 499651.9 138642.5, 499648....
# 3     3 MULTIPOLYGON (((496446.7 137764.8, 496442.9 137752.3, 496438.6 137739.9, 496433.6 137727.9, 49642...
# 4     4 MULTIPOLYGON (((500574.2 142260.4, 500573.1 142247.3, 500571.4 142234.4, 500569 142221.5, 500566 ...
# 5     5 MULTIPOLYGON (((497408.5 142703.5, 497404.7 142690.9, 497400.4 142678.6, 497395.4 142666.5, 49738...
# 6     6 MULTIPOLYGON (((496993.6 135888.6, 496982.8 135881.2, 496971.6 135874.4, 496960.1 135868.1, 49694...
# 7     7 MULTIPOLYGON (((496252.1 134426.9, 496250.4 134413.9, 496248 134401.1, 496244.9 134388.3, 496241....
# 8     8 POLYGON ((498130.8 141969.4, 498130.4 141956.3, 498129.4 141943.3, 498127.7 141930.3, 498125.3 14...
# 9     9 MULTIPOLYGON (((496175.9 144696.5, 496174.9 144683.5, 496173.1 144670.5, 496170.8 144657.6, 49616...
# 10    10 POLYGON ((208413.7 -53199.34, 208413.3 -53212.42, 208412.3 -53225.47, 208410.6 -53238.45, 208408....
# # ... with 29 more rows

Hopefully this has answered your question, as now we have:

Which school belongs to each unit schools_units
Unioned 250 m buffers schools_250_units
Unioned 450 m buffers schools_450_units

In terms of finding stores which are close to units/schools, like in your other question, you might consider just how widespread the units can be - when 450 m buffers are chained together. Take the map below for example.

tm_shape(schools_450_units, bbox = schools_450_units %>% filter(unit == 19) %>% st_buffer(2000) %>% st_bbox)+
  tm_polygons(col = 'unit', border.col = 'blue', legend.show = FALSE)+
  tm_shape(schools_250_units)+
  tm_polygons(col = 'unit', border.col = 'red', legend.show = FALSE)+
  tm_shape(schools_units)+
  tm_text(text = 'unit')

Look at unit 3 - it's pretty big and almost engulfs unit 22. Up to you to decide what you want to do and if this is suitable, but this was one of my initial thoughts based on a quick map.

You are probably right, the other question is probably more applicable, but this one is still helpful and I might still use it in some way, even if just as a sensitivity analysis. Please seem my comment in the other question, you have been so helpful! — revere2323, Nov 05 '21 at 17:14

R | sf: I have points and 2 buffers around each point. How to combine the points into single units if larger buffer overlaps (but smaller does not)?

1 Answers1

Linked