Блог

Этот вопрос следует задать и ответить на него здесь ранее. Я не могу найти решение, поэтому спрошу. Пожалуйста, укажите на дублирующий вопрос, если он есть.

Я пытаюсь добавить фрейм данных (или, точнее, фрагмент) географических точек данных с информацией о том, в каком полигоне находится каждая строка данных. Я знаю, как это сделать sp , и пытаюсь найти оптимальный sf / tidyverse способ. Решение, которое я придумал на основе sp sf словаря to здесь, кажется излишне сложным.

 library(sf)
#> Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1
library(tidyverse)

polys <- sf::st_read(system.file("shape/nc.shp", package ="sf")) %>% 
  sf::st_transform(4326)
#> Reading layer `nc' from data source 
#>   `/Library/Frameworks/R.framework/Versions/4.1/Resources/library/sf/shape/nc.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 100 features and 14 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27

set.seed(1)
points <- tibble::tibble(
  lon = sample(seq(st_bbox(polys)[c("xmin")], st_bbox(polys)[c("xmax")], 0.01), 
               size = 1e3,
               replace = TRUE
  ),
  lat = sample(seq(st_bbox(polys)[c("ymin")], st_bbox(polys)[c("ymax")], 0.01), 
               size = 1e3,
               replace = TRUE
  ),
  somedata = sample(letters, size = 1e3, replace = TRUE)
)

# Here is a way to do what I want. The code seems unnecessarily complicated, however.

points %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>% 
  mutate(county = 
           polys[
             sapply(st_intersects(.,polys), 
                        function(z) if (length(z)==0) NA_integer_ else z[1]),
             ]$NAME
  ) 
#> Simple feature collection with 1000 features and 2 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -84.32377 ymin: 33.88212 xmax: -75.47377 ymax: 36.58212
#> Geodetic CRS:  WGS 84
#> # A tibble: 1,000 x 3
#>    somedata             geometry county     
#>  * <chr>             <POINT [°]> <chr>      
#>  1 m        (-75.97377 34.19212) <NA>       
#>  2 z        (-77.54377 36.38212) Northampton
#>  3 m        (-83.04377 33.94212) <NA>       
#>  4 h        (-79.24377 36.04212) Orange     
#>  5 x        (-79.62377 36.21212) Guilford   
#>  6 m        (-81.34377 36.37212) Ashe       
#>  7 w        (-81.63377 34.00212) <NA>       
#>  8 g        (-82.46377 34.35212) <NA>       
#>  9 f        (-81.26377 34.38212) <NA>       
#> 10 x        (-78.36377 35.02212) Sampson    
#> # … with 990 more rows
 

^{Создано 2021-06-18 пакетом reprex (v2.0.0)}

Кто-нибудь знает более простой способ сделать это?

Я комбинирую здесь st_join выбор только нужных переменных полигона и переименование в стиле tidyverse. Я думаю, что это немного проще, чем предложенный вами подход

Данные

 library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(tidyverse)

polys <- sf::st_read(system.file("shape/nc.shp", package = "sf")) %>%
  sf::st_transform(4326)
#> Simple feature collection with 100 features and 14 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27


set.seed(1)
points <- tibble::tibble(
  lon = sample(seq(st_bbox(polys)[c("xmin")], st_bbox(polys)[c("xmax")], 0.01),
    size = 1e3,
    replace = TRUE
  ),
  lat = sample(seq(st_bbox(polys)[c("ymin")], st_bbox(polys)[c("ymax")], 0.01),
    size = 1e3,
    replace = TRUE
  ),
  somedata = sample(letters, size = 1e3, replace = TRUE)
)

 

Мой подход

 
points %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  # Spatial join selecting just NAME variable
  st_join(polys[, "NAME"]) %>%
  # Rename tidyverse style
  rename(county = NAME)
#> Simple feature collection with 1000 features and 2 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -84.32377 ymin: 33.88212 xmax: -75.47377 ymax: 36.58212
#> Geodetic CRS:  WGS 84
#> # A tibble: 1,000 x 3
#>    somedata             geometry county     
#>    <chr>             <POINT [°]> <chr>      
#>  1 m        (-75.97377 34.19212) <NA>       
#>  2 z        (-77.54377 36.38212) Northampton
#>  3 m        (-83.04377 33.94212) <NA>       
#>  4 h        (-79.24377 36.04212) Orange     
#>  5 x        (-79.62377 36.21212) Guilford   
#>  6 m        (-81.34377 36.37212) Ashe       
#>  7 w        (-81.63377 34.00212) <NA>       
#>  8 g        (-82.46377 34.35212) <NA>       
#>  9 f        (-81.26377 34.38212) <NA>       
#> 10 x        (-78.36377 35.02212) Sampson    
#> # ... with 990 more rows
 

^{Создано 2021-06-19 пакетом reprex (v2.0.0)}