# load packages
library(sf)
## Warning: package 'sf' was built under R version 4.0.2
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.0.2
library(sp)
## Warning: package 'sp' was built under R version 4.0.2
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.0.2
##
## 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(RColorBrewer)
# bring everything in
arbys <- sf::st_read("seminarClass/restaurants.gdb", layer="arbys")
## Reading layer `arbys' from data source `C:\Users\mwtro\OneDrive\Documents\r spatial\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("seminarClass/restaurants.gdb", layer="hardees")
## Reading layer `hardees' from data source `C:\Users\mwtro\OneDrive\Documents\r spatial\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
usa <- st_read("seminarClass/usa.shp")
## Reading layer `usa' from data source `C:\Users\mwtro\OneDrive\Documents\r spatial\data\seminarClass\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
# check classes, are they sf?
class(usa)
## [1] "sf" "data.frame"
class(arbys)
## [1] "sf" "data.frame"
class(hardees) #they are all sf, have to combine using lapply
## [1] "sf" "data.frame"
# check projections
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]]
# need to reproject...
usa <- st_transform(usa, 4269)
st_crs(usa) #they have the same projection now, both NAD83
## Coordinate Reference System:
## User input: EPSG:4269
## 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]]
# combine
usaarbys <- st_join(usa, arbys)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
usahardees <- st_join(usa, hardees)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
# ARBYS
arbysgroups <- group_by(usaarbys,STATE_NAME)
arbysdensity <- summarize(arbysgroups, NumRest = n(), STATE_NAME = unique(STATE_NAME))
## `summarise()` ungrouping output (override with `.groups` argument)
# HARDEES
hardeesgroups <- group_by(usahardees,STATE_NAME)
hardeesdensity <- summarize(hardeesgroups, NumRest = n(), STATE_NAME = unique(STATE_NAME))
## `summarise()` ungrouping output (override with `.groups` argument)
# plotting
plot(arbysdensity["NumRest"])
plot(hardeesdensity["NumRest"])
# Tweak Plots using PLOT
pal <- brewer.pal(7, "Reds")
plot(arbysdensity["NumRest"],
main = "Distribution of Arby's Throughout the USA",
breaks = "quantile", nbreaks = 7,
pal = pal)
# Plot with GGPLOT
ggplot(arbysdensity) +
geom_sf(aes(fill=NumRest))
ggplot(hardeesdensity) +
geom_sf(aes(fill=NumRest))
# Manipulate? will show if there are more arbys than hardees or more hardees than arbys
arbysdensity$difference <- (arbysdensity$NumRest-hardeesdensity$NumRest)
pal2 <- brewer.pal(3, "RdBu")
plot(arbysdensity["difference"],
main = "Are There More Arby's (blue) or Hardee's (red) per State?",
breaks = "quantile", nbreaks = 3,
pal = pal2)
