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