使用 sf,选择包含至少一个选择的点特征的区域特征?

With sf, selecting area features that contain at least one of a selection of point features?

例如,假设我有一个 sf 对象,其中包含美国本土的 4 个城市及其坐标。然后我有一个具有 48 个特征的 sf 对象(每个可能的状态一个)。有没有办法 select 包含指定城市的州子集?类似于:

cities_sf

state_sf %>%
  filter(states s.t. there exists x in cities_sf s.t. x in states_sf) +
  ggplot() +
  ...

编辑: st_within(my_cities, my_states) 给了我

structure(list(290L, 378L, 51L, integer(0), 283L, 478L, 415L, 
380L, 489L, 64L, 189L, 184L, 311L, 488L, 66L, 73L, 49L, 1L, 
359L, 111L, 502L, 489L, 272L, 115L, 352L, 241L), predicate = "within", 
region.id = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", 
"12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", 
"24", "25", "26"), ncol = 544L, class = "sgbp")

我可以看出这 26 个索引对应于 my_states 中包含城市的多面体,但我不确定如何使用该 SGBD("sparse geometry binary predicate",根据文档) ggplot/geom_sf 条款

中的对象

编辑 2:我最终使用了 slice(states_sf, unlist(st_within(cities_sf, states_sf))),它给出了一个 sf 对象,它是我需要的子集

没有样本数据,这是我能想到的最好的。

library( sf )

#find intersecting points/polygons    
intersect <- st_intersection(x = polygons, y = points)

#and go further from there

更新

使用@Spacedman 在他的回答中提供的样本数据。

library(dplyr)
library(sf)

states %>% 
  #create ID's for the states (if they don't have one already) 
  #state ID should be equal to rownumber (fot the filter later on)
  mutate( id = row_number() ) %>%
  #filter out states that do not have any intersetcions with the points/cities
  filter( id %in% unlist( st_intersects(cities, states) ) ) %>%
  #plot
  mapview::mapview()   

使用 USAboundaries 包中的 us_states 函数,让我们制作一小组状态:

> states <- us_states(map_date = "2000-01-01", resolution = "high", states = c("CA", "OR", "WA","NV","NM","UT","CO","ID","AZ"))

这是我创建的一些点:

> pts
Simple feature collection with 4 features and 0 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -121.7663 ymin: 34.86508 xmax: -110.7263 ymax: 46.65593
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
                    geometry
1 POINT (-110.7263 34.86508)
2 POINT (-111.7345 38.64123)
3 POINT (-120.1531 46.65593)
4 POINT (-121.7663 39.37335)

要测试交集:

> st_intersects(states, pts)
although coordinates are longitude/latitude, st_intersects assumes that they are planar
Sparse geometry binary predicate list of length 9, where the predicate was `intersects'
 1: 1
 2: 4
 3: (empty)
 4: (empty)
 5: (empty)
 6: (empty)
 7: (empty)
 8: 2
 9: 3

该对象是一个列表,因此您可以获得元素的长度并找到大于零的元素 - 即其中包含一些东西:

> lengths(st_intersects(states, pts)) > 0
although coordinates are longitude/latitude, st_intersects assumes that they are planar
[1]  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE

然后以常规方式对空间多边形进行子集化:

 > plot(st_geometry(states[lengths(st_intersects(states, pts)) > 0,]))

绘制了具有四个点的四个状态。

创建子集并将其提供给 ggplot 如果这就是您绘制地图的方式。