右 | sf:我在每个点周围都有点和 2 个缓冲区。如果较大的缓冲区重叠(但较小的不重叠),如何将这些点组合成单个单元?
Posted
技术标签:
【中文标题】右 | sf:我在每个点周围都有点和 2 个缓冲区。如果较大的缓冲区重叠(但较小的不重叠),如何将这些点组合成单个单元?【英文标题】: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)? 【发布时间】:2021-12-18 06:10:53 【问题描述】:我目前在地图(学校)上有点。每个在点(学校)周围都有两个缓冲区。一个是450m,一个是250m。如果点重叠,我想将它们视为一个单元(因为否则事情会变得复杂),但我希望它们保持它们所覆盖的几何形状/区域。
因此,在此处给出的示例地图上,我希望将前三所学校/分数合并为一个“单元”。我希望他们保留他们所覆盖的区域,但只需将 R 视为一个单位。如果我使用“st_union”函数,我必须对较小和较大的缓冲区都这样做。但这会导致单元总数的差异,因为更多的较大缓冲区重叠。所以基本上我希望较小的缓冲区与其他较小的缓冲区合并(成一个单元),如果较大的缓冲区重叠并因此合并。
最终目标是,我想计算一家商店在距离学校 250m 范围内的次数与 450m 的距离。我认为,如上所述,由于学校之间存在大量重叠,它最容易合并。
我有sample data here.
我的代码如下,但请注意,由于上述原因,它会导致学校/单位的数量不同。
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)))
【问题讨论】:
我建议在重叠缓冲区多边形上使用st_intersection()
来查找重叠区域(以及哪些原始缓冲区/学校)对重叠区域有贡献。这个答案***.com/a/69293465/13478749 描述了这个过程(针对不同的问题),这将是这个问题的一个好的开始。
我知道该怎么做,但我不太清楚的是,当且仅当较大的缓冲区相交时,如何将较小的缓冲区合并为一个(即使它们不相交)。我想这是某种 if/then 语句,但我无法在 sf 中弄清楚。
【参考方案1】:
最初我建议使用st_intersection()
,但开始变得复杂,所以我采用了这种方法,而不是依赖于st_join()
。
首先,加载库和数据。
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
然后,在学校周围建立缓冲区。
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
使用st_union()
加入所有450 m缓冲区,以确定“单元”的身份。
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
请注意,在 156 所学校中,我们只有 39 个单位。当然,当多个 450 m 缓冲区链接在一起时,这些单元可能会影响深远 - 请参阅最后的地图。
然后,我们使用st_join()
来确定每个单元内有哪些学校(或其缓冲区)。指定join = st_within
是这里的关键。
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
然后将同一单元的 250 m 缓冲区(即使它们不接触)组合成一个 MULTIPOLYGON
,我们使用 group_by(unit)
并使用 summarise(do_union = TRUE)
(默认值)。这类似于st_union()
,但尊重group_by()
(实际上使用st_union()
可能得到完全相同的结果,但这也有效)。
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
希望这已经回答了您的问题,就像我们现在一样:
每个单元属于哪个学校schools_units
联合 250 m 缓冲区schools_250_units
联合 450 m 缓冲区schools_450_units
就寻找靠近单位/学校的商店而言,例如在您的other question 中,您可能会考虑这些单位的分布范围——当 450 m 缓冲区链接在一起时。以下面的地图为例。
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')
看看第 3 单元 - 它非常大,几乎吞没了第 22 单元。由您决定要做什么以及是否合适,但这是我基于快速地图的初步想法之一。
【讨论】:
您可能是对的,另一个问题可能更适用,但这个问题仍然有用,我可能仍会以某种方式使用它,即使只是作为敏感性分析。请在另一个问题中看到我的评论,你很有帮助!以上是关于右 | sf:我在每个点周围都有点和 2 个缓冲区。如果较大的缓冲区重叠(但较小的不重叠),如何将这些点组合成单个单元?的主要内容,如果未能解决你的问题,请参考以下文章