Read in spatial data and assign variables

usa <- read_sf("Data/seminarClass/usa.shp")
arbys <- sf::st_read("Data/seminarClass/restaurants.gdb", layer="arbys")
## Reading layer `arbys' from data source `C:\Users\asharp5\Desktop\r-intro\Mapping\ArbysvsHardeesDistribution\Data\seminarClass\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 <- sf::st_read("Data/seminarClass/restaurants.gdb", layer="hardees")
## Reading layer `hardees' from data source `C:\Users\asharp5\Desktop\r-intro\Mapping\ArbysvsHardeesDistribution\Data\seminarClass\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
plot(usa)
## Warning: plotting the first 9 out of 54 attributes; use max.plot = 54 to plot
## all

st_crs(usa)
## Coordinate Reference System:
##   User input: USA_Contiguous_Albers_Equal_Area_Conic 
##   wkt:
## PROJCRS["USA_Contiguous_Albers_Equal_Area_Conic",
##     BASEGEOGCRS["NAD83",
##         DATUM["North American Datum 1983",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["Degree",0.0174532925199433]]],
##     CONVERSION["USA_Contiguous_Albers_Equal_Area_Conic",
##         METHOD["Albers Equal Area",
##             ID["EPSG",9822]],
##         PARAMETER["Latitude of false origin",37.5,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8821]],
##         PARAMETER["Longitude of false origin",-96,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8822]],
##         PARAMETER["Latitude of 1st standard parallel",29.5,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8823]],
##         PARAMETER["Latitude of 2nd standard parallel",45.5,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8824]],
##         PARAMETER["Easting at false origin",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8826]],
##         PARAMETER["Northing at false origin",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8827]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["USA - CONUS - onshore"],
##         BBOX[24.41,-124.79,49.38,-66.91]],
##     ID["ESRI",102003]]
st_crs(hardees)
## 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]]
st_crs(arbys)
## 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]]
#Right now the CRS are different, I want each layer to be in the same CRS.

Reprojection

#reproject usa into NAD83 CRS and assign to new object
usa_sf_nad <- usa[1:8] %>%
  st_as_sf() %>%
  st_transform(st_crs(arbys))

#check to confirm the transformation.
st_crs(usa_sf_nad)
## 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]]

Spatial Join

#join arbys to usa basemap layer
arbys_usa_join <- usa_sf_nad %>% 
  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)

Scatterplot of Arbys rate per state

plot(arbys_usa_join$arbys_rate)

## #Choropleth Map of Arbys rate per state#

plot(arbys_usa_join["arbys_rate"], axes=TRUE, main= "Arbys Distribution", nbreaks=7,breaks="quantile", )

Join hardees to usa basemap layer

Hardees_usa_join <- usa_sf_nad %>% 
  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)

Scatterplot of Hardees rate per state

plot(Hardees_usa_join$hardees_rate)

## #Choropleth Map of Hardees rate per state#

plot(Hardees_usa_join["hardees_rate"], axes=TRUE, main= "Arbys Distribution", nbreaks=7,breaks="quantile", )

### Ggplot2 map of Arbys rate per state

ggplot(arbys_usa_join)+
  geom_sf(aes(fill=arbys_rate))

### Ggplot2 map of Arbys rate per state

ggplot(Hardees_usa_join)+
  geom_sf(aes(fill=hardees_rate))