Hey @YuJiang ,
For your question, your best friend is the st_intersects()
function which is in the {sf}
package.
I took some time to write the code to do what you asked. The codedoes some web scraping, spatial analysis and mapping. Don't hesitate to let me know if you have any questions.
# 1. LOAD PACKAGES ----
library(dplyr)
library(purrr)
library(rvest)
library(sf)
library(slugify) # remotes::install_github("hrbrmstr/slugify")
library(stringr)
library(tidygeocoder)
library(tigris)
library(tmap)
library(units)
# 2. SPATIAL ANALYSIS 1: MAIN ADDRESSES ----
# 2.1 Geocoding ----
addresses <- c("1331 Cir Park Dr, Knoxville, TN, 37916",
"1210 Varsity Dr., E. Carroll Joyner Visitor Center, Raleigh, NC, 27606",
"1600 Pennsylvania Avenue NW, Washington, DC 20500")
main_locs_geo <- tibble(
location = paste0("location", seq_along(addresses)),
address = addresses
) %>%
geocode(address = address, method = "arcgis") %>%
st_as_sf(coords = c("long", "lat"), crs = 4326)
# 2.2 Create 10-mile buffers ----
main_locs_buffers <- st_buffer(x = main_locs_geo, dist = units::set_units(x = 10, value = "mile"))
# 2.3 Download counties spatial data ----
counties_geo <- counties(cb = TRUE, keep_zipped_shapefile = TRUE) %>%
st_simplify() %>%
st_transform(crs = 4326)
# 2.4 Counties traversed by the 10-mile buffers ----
i <- st_intersects(x = main_locs_buffers, y = counties_geo)
counties2_geo_list <- map(.x = i, .f = \(x){
counties_geo %>% slice(x)
})
# 3. WEB SCRAPING OF USPS LOCATIONS ----
# 3.1 Determine counties of interest ----
counties_urls <- counties2_geo %>%
transmute(url = slugify(x = paste(NAMELSAD, STUSPS) %>% tolower())) %>%
pull(url) %>%
paste0("https://www.postlocations.com/usps-in-", .)
# 3.2 Create a function to scrape USPS locations and use it ----
scrape_usps_addresses <- function(county_url){
html <- read_html(x = county_url)
name <- html_elements(x = html, css = ".name") %>%
html_text()
addresses <- html_elements(x = html, css = ".address") %>%
html_text()
city_state <- html_elements(x = html, css = ".city-state") %>%
html_text()
zipcode <- html_elements(x = html, css = ".zipcode") %>%
html_text()
tibble(
name = name,
address = paste0(addresses, ", ", city_state, " ", zipcode)
)
}
usps_addresses <- map_dfr(.x = counties_urls, .f = scrape_usps_addresses) %>%
filter(!name %in% c("USPS Mailbox"))
usps_geo <- geocode(.tbl = usps_addresses, address = address, method = "arcgis") %>%
st_as_sf(coords = c("long", "lat"), crs = 4326)
# 5. MAPPING ----
# 5.1 Create a function to plot the maps and use it ----
tmap_mode(mode = "plot")
plot_map <- function(shp){
tm_shape(shp = shp) +
tm_polygons() +
tm_shape(shp = main_locs_buffers) +
tm_borders(col = "purple") +
tm_shape(shp = main_locs_geo) +
tm_dots(col = "red", size = 0.2) +
tm_shape(shp = usps_geo) +
tm_dots(col = "blue", size = 0.05)
}
map(.x = counties2_geo_list, .f = plot_map)
- Address in Tennessee

- Address in North Carolina

- Address in DC
