右 | 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 个缓冲区。如果较大的缓冲区重叠(但较小的不重叠),如何将这些点组合成单个单元?的主要内容,如果未能解决你的问题,请参考以下文章

sf中标准的分页功能介绍

使用 tmap 和 sf 循环记录一系列地图

poj 1189 钉子和小球

(转)Eclipse的jar里的class文件打断点

(转)Eclipse的jar里的class文件打断点

钉子与小球