rpurrrr-sf

`purrr` mapping of `st_intersects()` using inputs from across two lists


I have a list of geocoded point dataframes that I want to intersect with a corresponding list of polygon geometry dataframes. Below is a reprex to create the data structure:

library(dplyr)
library(purrr)
library(sf)
library(tigris)
# pull two Oregon counties and their polygons 
counties <- counties(state = "OR")[1:2,] %>% dplyr::select(NAME, GEOID, geometry)

# split the dataframe into a list to emulate my real data
or_counties <- counties %>% split(.$GEOID)
# create random point data for the two counties and add the GEOID field (a unique state + county numeric index)
wallowa_pts = as.tibble(st_jitter(st_sample(counties[1,], 10), factor=0.2)) %>% mutate(GEOID = counties[1,]$GEOID)
crook_pts = as.tibble(st_jitter(st_sample(counties[2,], 10), factor=0.2)) %>% mutate(GEOID = counties[2,]$GEOID)

# bind and convert to sf
pts <- rbind(wallowa_pts, crook_pts) %>% st_as_sf()

# split to a list
or_pts <- pts %>% split(.$GEOID)

So I now have point data in one list grouped by GEOID, and polygon data in another list, also grouped by GEOID. I need to create a new field in each list-element in or_pts that indicates the result of the evaluation of st_intersects() on the corresponding list-element of or_counties, in this case recoding T/F to colors for visualization. I can get a narrowly scoped result if I manually define the list-element in or_counties like so:

or_pts %>% 
  map(
    ~ mutate(., status =
               case_when(
                 st_intersects(., or_counties[[1]], sparse = F) == T ~ "green",
                 TRUE ~ "red"
               )
    )
  )

However, what I want is to dynamically intersect data from or_pts$41013 with or_counties$41013, or_pts$41063 with or_counties$41063, and so on for tens of thousands of points across thousands of counties to identify geocoding errors.

This seems like a nested map2/imap job but I' cant figure out where that would fit in this process.


Solution

  • Not sure where nesting comes in, so perhaps I don't quite understand the problem, but
    map2, if lists are aligned:

    map2(or_pts, or_counties, 
         \(pts, poly) mutate(pts, status = if_else(st_intersects(pts, poly, sparse = FALSE), "green", "red" )))
    

    imap, if you need to use or_pts names to access corresponding items in or_counties:

    imap(or_pts,
         \(pts, idx) mutate(pts, status = if_else(st_intersects(pts, or_counties[[idx]], sparse = FALSE), "green", "red" )))