Add polygon name to a data frame using sf and tidyverse

user_1082004 0
Viewed 139

This question ought to be asked and answered here previously. I cannot find the solution, so will ask. Please point to the duplicate question if there is one.

I am trying to append a data frame (or a tibble to be precise) of geographic data points with information of which polygon each row of data lies in. I know how to do this in sp and am trying to find the optimal sf/tidyverse way. The solution I have come up with based on sp to sf dictionary here appears unnecessarily complicated.

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

Created on 2021-06-18 by the reprex package (v2.0.0)

A simpler way to do this, anyone?

sf dplyr gis r
1 Answers

I combine here st_join selecting just the desired variables of the polygon and renaming on tidyverse style. I think this is a bit more simpler than your proposed approach

Data

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

My approach


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

Created on 2021-06-19 by the reprex package (v2.0.0)

user_78779170
Related