С помощью sf, выбирать объекты области, которые содержат хотя бы один из набора точечных объектов?

Question

Например, допустим, у меня есть объект sf, содержащий 4 города в смежных Соединенных Штатах и ​​их координаты. И тогда у меня есть sf объект с 48 функциями (по одному на каждое возможное состояние). Есть ли способ выбрать подмножество штатов, которые содержат назначенные города? Что-то вроде:

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 («двоичным предикатом разреженной геометрии», в соответствии с документацией) в ggplot / geom_sf условия

edit 2 : я закончил с использованием slice(states_sf, unlist(st_within(cities_sf, states_sf))), который дает sf объект, который является подмножеством, которое мне нужно

Answer 1

Используя функцию us_states из пакетов USAboundaries, давайте создадим небольшой набор состояний:

> 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, если вы так рисуете карту.

Answer 2

без образцов данных, это лучшее, что я могу придумать.

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()