Reputation: 37
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)?
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)))
Upvotes: 1
Views: 1268
Answers (1)
Reputation: 1370
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.
Upvotes: 3
Related Questions
- How to merge multiple sf polygons?
- Create buffers in km units for global point dataset in R?
- How to unite touching polygons with equal number of users's points in sf package
- Counting points when buffers are overlapping
- Is there a function is "sf" package (R) to fully merge geometries?
- Spatial merge of two sf objects
- R: Append data when points overlap/within distance; add buffer rectangle to set1, add radius to set2
- How to combine two sf objects with the same CRS and non-overlapping area?
- How to merge dataframes of unequal length based on time with buffer intervals in R?
- Merge two polygon regions into a single polygon region in R