найти состояние, в котором находятся точки, `sf` r - PullRequest
Новая
       41

найти состояние, в котором находятся точки, `sf` r

2 голосов
/ 05 марта 2020

У меня много точек, сгенерированных из спутниковых снимков, и я хочу найти, где находятся точки (состояния). Я искал некоторые ссылки и получил подсказку, что я должен использовать sf::st_intersects, но оказывается, что не работает. Вот минимальный пример: 

library(sf)
library(ggplot2)
library(dplyr)

# Contiguous US state boundaries
usa = st_as_sf(maps::map("state", fill = TRUE, plot = FALSE))

# Simulate some random points
pts = data.frame(
  x = c(-91.6, -74.3, -101.5),
  y = c(36.1, 42.1, 25.3)
  ) %>%
  st_as_sf(coords=c("x", "y"), crs = 4326)

ggplot() +
  geom_sf(data = usa, fill = NA) +
  geom_sf(data = pts,
          shape = 21, size = 4, fill = "red") +
  coord_sf(crs = st_crs(102003)) +
  theme_minimal()

Вот результирующий рисунок: enter image description here Я хочу получить pts data.frame с еще одной строкой, указывающей состояние Очки: 

             geometry  state
1  POINT (-91.6 36.1)  arkansas
2  POINT (-74.3 42.1)  new york
3  POINT (-101.5 25.3) NA

Я знаю, что должен показать sf::st_transform, но мне это не удалось. В идеале я хочу, чтобы пересечение было масштабируемым, поскольку у меня более 1 000 000 000 точек.

Спасибо!

Ответы [ 2 ]

2 голосов
/ 06 марта 2020

Вы можете использовать st_join. 

a <- pts %>% 
  st_join(usa)

> a 
Simple feature collection with 3 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: -101.5 ymin: 25.3 xmax: -74.3 ymax: 42.1
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
        ID            geometry
1 arkansas  POINT (-91.6 36.1)
2 new york  POINT (-74.3 42.1)
3     <NA> POINT (-101.5 25.3)
2 голосов
/ 05 марта 2020 
int = st_intersects(pts, usa)
sapply(int, function(x) if(length(x) == 0) NA else as.character(usa$ID[x]))
#[1] "arkansas" "new york" NA
...