使用 sf 和 tidyverse 将多边形名称添加到数据框

Question

这个问题之前应该在这里提问和回答。我找不到解决方案，所以会问。如有重复问题请指出。

我正在尝试附加一个地理数据点的数据框（或者更准确地说是一个小标题），其中包含每行数据所在的多边形信息。我知道如何在 sp 和我正在尝试找到最佳的 sf/tidyverse 方式。我根据 sp 到 sf 字典 here 提出的解决方案显得不必要地复杂。

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

^{由 reprex package (v2.0.0)}

于 2021-06-18 创建

有人知道更简单的方法吗？

Answer 1

我在这里合并 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

^{由 reprex package (v2.0.0)}

于 2021-06-19 创建

使用 sf 和 tidyverse 将多边形名称添加到数据框

Add polygon name to a data frame using sf and tidyverse

gis

r

dplyr

sf