Load in Packages

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.1.8     ✔ dplyr   1.0.9
## ✔ tidyr   1.2.0     ✔ stringr 1.4.1
## ✔ readr   2.1.2     ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
## 
## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.

Business Categories and Location of Interest

I am interested in getting a list of campground and coffee shops in White County, GA, as I want to go on a backpacking trip in the mountains of GA. I also need coffee to function so if I could find a campsite with easy coffee access that would be great.

Get County SF Data

#Load in tract geometries using the tigris package, filtered for White County, GA
whiteCountyTracts <- tracts(state = 'GA', county = 'White County')
## Retrieving data for the year 2020
#Just for fun, lets visualize the tracts
tmap_mode('plot')
## tmap mode set to plotting
tm_shape(whiteCountyTracts) + tm_borders()

Prep Search Radii for Each Tract

First we need to write a function that will calculate the radius from the centroid of each tract, to the corner of the tract’s bounding box

get_radius <- function(geo, epsg_id){
  
  #Get the bounding box of the polygon
  bb <- st_bbox(geo)
  
  #Get the coordinates of a corner on the box
  bb_corner <- st_point(c(bb[1], bb[2])) %>%
    #Adjust it to a point on the Earth's surface
    st_sfc(crs = epsg_id)
  
  #Get the centroid of the bounding box by finding the center of both lines.
  bb_center_x <- (bb[3] + bb[1]) / 2
  bb_center_y <- (bb[4] + bb[2]) / 2
  bb_centroid <- st_point(c(bb_center_x, bb_center_y)) %>%
    #Adjust the point to the earths surface
    #And convert it to an SF object
    st_sfc(crs = epsg_id) %>% st_sf()
  
  #Lastly, calculate the radius
  r <- st_distance(bb_corner, bb_centroid)
  #Multiply by 1.1 to make it slightly bigger than the bounding box
  bb_centroid$radius <- r
  return(bb_centroid)
}

Now we will calculate the radius for each tract.

#Establish the CRS
epsg_id <- 4326

whiteGeom <- whiteCountyTracts %>%
  #Get geometries
  st_geometry() %>%
  #Account for earth's surface
  st_transform(crs = epsg_id) %>%
  #Apply the get radius function for each geometry
  lapply(function(x) get_radius(x, epsg_id = epsg_id)) %>%
  #Bind the rows into a dataframe
  bind_rows() %>%
  #Get x,y coords
  mutate(x = st_coordinates(.)[,1],
         y = st_coordinates(.)[,2])

# Select the first 10 rows
whiteGeom %>% 
  # Draw a buffer centered at the centroid of Tract polygons.
  # Radius of the buffer is the radius we just calculated using loop
  st_buffer(., dist = .$radius) %>% 
  # Display this buffer in red
  tm_shape(.) + tm_polygons(alpha = 0.5, col = 'red') +
  # Display the original polygon in blue
  tm_shape(whiteCountyTracts) + tm_borders(col= 'blue')

Get Yelp Data

As noted in the practice RMD, we will need to do multiple requests for each tract and for each category. We need to write a custom function to do this.

get_yelp <- function(tracts, categories, epsg_id, key = Sys.getenv('yelp_key')){
  
  #Establish output frame
  frame = data.frame()
  
  #For each category
  for(cat in categories){
    
    #For every tract
    for(i in 1:nrow(tracts)){
      
      #Offset ticker
      n <- 1
      
      #Get the list of businesses
      businessList <- business_search(key,
                                      categories = cat,
                                      latitude = tracts[i,]$y,
                                      longitude = tracts[i,]$x,
                                      radius = round(tracts[i,]$radius),
                                      offset = (n - 1)*50,
                                      limit = 50
                                      )
      
      #Append to final dataframe
      frame <- bind_rows(frame, businessList$businesses)
      
      #If total is greater than 50
      if(length(businessList$total) > 0 & businessList$total > 50){
        
        #Iterate through the remaining businesses
        while(n < ceiling(businessList$total / 50)){
         
          #Append the business list to a final output dataframe
          
          #Adjust the offset ticker
          n <- n + 1
          
          #Get the list of businesses
          businessList <- business_search(key,
                                          categories = cat,
                                          latitude = tracts[i,]$y,
                                          longitude = tracts[i,]$x,
                                          radius = round(tracts[i,]$radius),
                                          offset = (n - 1)*50,
                                          limit = 50
                                          )
          
          #Append the list to the dataframe
          frame <- bind_rows(frame, businessList$businesses)
        
        }
      }
      
      
    }
  }
  
  #Convert x, y coords to ST points, normalized with a crs
  frameGeom <- frame %>%
    #Rename the alias field to unnest the categories
    rename('alias.business' = alias) %>%
    #Unnest and pivot into one categories column
    unnest_wider(categories, transform = ~ paste(.x, collapse = ', ')) %>%
    
    #Rename the category variables
    rename('title.category' = 'title',
           'alias.category' = 'alias'
           ) %>%
    
    #Flatten nested frame
    jsonlite::flatten() %>%
    st_as_sf(coords = c('coordinates.longitude', 'coordinates.latitude'), crs = epsg_id)
  return(frameGeom)
  
}

Now lets get the businesses

white_camping_yelp <- get_yelp(whiteGeom, c('campgrounds', 'coffee'), epsg_id) %>%
  
  #Identify if the business is a coffee shop or a campground
  mutate(coffee_or_camp = case_when(
    sum(str_detect(title.category, c('Campgrounds', 'Coffee'))) == 2 ~ 'Both',
    str_detect(title.category, 'Campgrounds') ~ 'Campground',
    str_detect(title.category, 'Coffee') ~ 'Coffee Shop'
  ))

Map Businesses on the County Map

Lets overlay the business data onto the map of white county.

#Set tmap mode to view
tmap_mode('view')
## tmap mode set to interactive viewing
#Plot the county
tm_shape(whiteCountyTracts) + tm_polygons(id = c('Tract Name' = 'NAMELSAD')) +
  tm_shape(white_camping_yelp) + tm_dots(col = c('Business Category' = 'coffee_or_camp'), jitter = .5, title = 'Business Category', id = 'name')

Final Questions

  1. What’s the county and state of your choice?
  1. How many businesses are there in total?
## [1] "67 businesses"
  1. How many businesses are there for each business category?
##   coffee_or_camp  n
## 1     Campground 24
## 2    Coffee Shop 43
  1. Upon visual inspection, can you see any noticeable spatial patterns to the way they are distributed across the county (e.g., clustering of businesses at some parts of the county)?
  1. (Optional) Are there any other interesting findings?