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)
