#Laoding packages

getwd()
## [1] "/Users/devinhainje/Desktop/R"
setwd("/Users/devinhainje/Desktop/R")
library(ggplot2)
library(dbplyr)
library(sf)
## Linking to GEOS 3.8.1, GDAL 3.1.1, PROJ 6.3.1
library("dplyr")
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:dbplyr':
## 
##     ident, sql
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(RColorBrewer)
library(raster)
## Loading required package: sp
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
library(rasterVis)
## Loading required package: lattice
## Loading required package: latticeExtra
## 
## Attaching package: 'latticeExtra'
## The following object is masked from 'package:ggplot2':
## 
##     layer
library(ggmap)
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
library(ggmap)

#Read Data

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
head(arbys)
## Simple feature collection with 6 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -73.25727 ymin: 42.12308 xmax: -71.83845 ymax: 42.57244
## geographic CRS: NAD83
##   COMPANY NumRest                      Shape
## 1  ARBY'S       1 POINT (-72.57847 42.18045)
## 2  ARBY'S       1 POINT (-72.54092 42.36302)
## 3  ARBY'S       1 POINT (-72.62648 42.12308)
## 4  ARBY'S       1 POINT (-72.48863 42.14323)
## 5  ARBY'S       1 POINT (-73.25727 42.57244)
## 6  ARBY'S       1  POINT (-71.83845 42.2023)
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

#Reprojected the USA layer

usa83 <- st_transform(usa,st_crs(arbys))
st_crs(usa83)
## Coordinate Reference System:
##   User input: NAD83 
##   wkt:
## GEOGCRS["NAD83",
##     DATUM["North American Datum 1983",
##         ELLIPSOID["GRS 1980",6378137,298.257222101,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["North America - NAD83"],
##         BBOX[14.92,167.65,86.46,-47.74]],
##     ID["EPSG",4269]]
#reprojected usa to NAD83, like arbys and hardees

#Joined the Arby’s layer to the USA

z <- usa83 %>% 
  mutate(tract_area = st_area(geometry)) %>%
  st_join(arbys) %>%
  group_by(STATE_FIPS) %>% 
  summarize(n_arbys = n(),
            tract_area = unique(tract_area),
            arbys_rate = n_arbys/tract_area * 1e6)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## `summarise()` ungrouping output (override with `.groups` argument)
hist(z$arbys_rate)

#histogram of rate
plot(z["arbys_rate"], axes=TRUE, main= "Arbys Distribution", nbreaks=7,breaks="quantile", )

#ITSSSS A MAP

#Joined Hardees to USA

h <- usa83 %>% 
  mutate(tract_area = st_area(geometry)) %>%
  st_join(hardees) %>%
  group_by(STATE_FIPS) %>% 
  summarize(n_hardees = n(),
            tract_area = unique(tract_area),
            hardees_rate = n_hardees/tract_area * 1e6)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## `summarise()` ungrouping output (override with `.groups` argument)
plot(h["hardees_rate"], axes=TRUE, main= "Hardee's Distribution", nbreaks=7,breaks="quantile" )