## ── 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.
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.
#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()
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')
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'
))
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')
## [1] "67 businesses"
## coffee_or_camp n
## 1 Campground 24
## 2 Coffee Shop 43