How to get the number of USPS stores and their distances to a specific address within a certain mile

Hello I have a question about how to get the number of USPS stores and their distances to a specific address with a certain mile. For example, I have three addresses, "1331 Cir Park Dr, Knoxville, TN, 37916", "1210 Varsity Dr., E. Carroll Joyner Visitor Center, Raleigh, NC, 27606" and "1600 Pennsylvania Avenue NW, Washington, DC 20500". I want to know if I use these three addresses as three circle centers, and then I use 10 miles as the radius to draw a circle, how many USPS stores can be included inside this circle and for each USPS store, what's the distance between the USPS store and the address (the circle center)? Really appreciate your help.

The diagram is shown below. The blue soli circle is an address and the organ circles are USPS stores.
Screenshot 2023-08-22 at 17.31.37

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

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)
  1. Address in Tennessee

tennessee

  1. Address in North Carolina

north_carolina

  1. Address in DC

dc

Thanks so much for your help. It does help me out of my problem!!! REALLY APPRECIATE IT.

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.

If you have a query related to it or one of the replies, start a new topic and refer back with a link.