library(ggplot2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(sf)
## Linking to GEOS 3.8.1, GDAL 3.1.1, PROJ 6.3.1
library(RColorBrewer)
library(classInt)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
restaurants <- st_read("restaurants.gdb")
## Multiple layers are present in data source /Users/devinhainje/Desktop/R/restaurants.gdb, reading layer `arbys'.
## Use `st_layers' to list all layer names and their type in a data source.
## Set the `layer' argument in `st_read' to read a particular layer.
## Warning in evalq((function (..., call. = TRUE, immediate. = FALSE, noBreaks. =
## FALSE, : automatically selected the first layer in a data source containing more
## than one.
## Reading layer `arbys' from data source `/Users/devinhainje/Desktop/R/restaurants.gdb' using driver `OpenFileGDB'
## Simple feature collection with 2590 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -158.0057 ymin: 19.69826 xmax: -68.01207 ymax: 61.21599
## geographic CRS: NAD83
restaurants_layers <- st_layers(dsn = "restaurants.gdb") 
restaurants_layers                              
## Driver: OpenFileGDB 
## Available layers:
##   layer_name geometry_type features fields
## 1      arbys         Point     2590      2
## 2    hardees         Point     3262      2
arbys <- st_read("restaurants.gdb", layer = "arbys")
## Reading layer `arbys' from data source `/Users/devinhainje/Desktop/R/restaurants.gdb' using driver `OpenFileGDB'
## Simple feature collection with 2590 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -158.0057 ymin: 19.69826 xmax: -68.01207 ymax: 61.21599
## geographic CRS: NAD83
hardees <- st_read("restaurants.gdb", layer = "hardees")
## Reading layer `hardees' from data source `/Users/devinhainje/Desktop/R/restaurants.gdb' using driver `OpenFileGDB'
## Simple feature collection with 3262 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -117.0076 ymin: 25.78868 xmax: -73.34611 ymax: 48.80142
## geographic CRS: NAD83
usa <- st_read("usa.shp")
## Reading layer `usa' from data source `/Users/devinhainje/Desktop/R/usa.shp' using driver `ESRI Shapefile'
## Simple feature collection with 51 features and 54 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -6293474 ymin: -1294964 xmax: 2256319 ymax: 4592025
## projected CRS:  USA_Contiguous_Albers_Equal_Area_Conic
usa83 <- st_transform(usa,st_crs(arbys))

a <- usa83 %>% 
  select(STATE_NAME, STATE_ABBR) %>%
  mutate(stateArea = st_area(geometry)) %>%
  st_join(arbys) %>%
  group_by(STATE_ABBR)%>%
  summarise(arbys_rate=n()) %>%
  st_join(hardees) %>%
  group_by(STATE_ABBR, arbys_rate)%>%
  summarise(hardees_rate=n())
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## `summarise()` ungrouping output (override with `.groups` argument)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## `summarise()` regrouping output by 'STATE_ABBR' (override with `.groups` argument)
breaks_a <- classIntervals(c(min(a$arbys_rate)-.00001,
                             a$arbys_rate), n=7,
                           style = "quantile")
breaks_h <- classIntervals(c(min(a$hardees_rate)-.00001,
                             a$hardees_rate), n=7,
                           style = "quantile")

a_cat <- a %>%
  mutate(arbys_cat=cut(arbys_rate, breaks_a$brks),
         hardees_cat=cut(hardees_rate, breaks_h$brks))

x <- ggplot(a_cat) + geom_sf(aes(fill = arbys_cat)) +
  scale_fill_brewer(palette = "PuBu")
y <- ggplot(a_cat) + geom_sf(aes(fill = hardees_cat)) +
  scale_fill_brewer(palette = "PuBu")
x

y

grid.arrange(x,y, nrow=2)