Lab 4

#load the libraries
library(sf)           
library(rgdal)        
library(ggplot2)      
library(dplyr)        
library(tidyr)        
library(scales)       
library(RColorBrewer) 
library(units)
library(cowplot)
library(here)
library(leaflet)
library(tidyverse)
library(USAboundaries)
#load the data
okcounty <- st_read(here("Chapter5/Chapter5/data", "ok_counties.shp"), quiet = TRUE)
tpoint <- st_read(here("Chapter5/Chapter5/data", "ok_tornado_point.shp"), quiet = TRUE)
tpath <- st_read(here("Chapter5/Chapter5/data", "ok_tornado_path.shp"), quiet = TRUE)

Part A

Question 1

Generate a map of tornado paths where the paths from each year are displayed as a different color, similar to the map of tornado points in Figure 5.4 Links to an external site.. Create a composite figure containing the map of tornado paths and the map of tornado points using plot_grid().

#Question 1

#filter down and select relevant variables for the point and path datasets
tpoint_16_21 <- tpoint %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  select(om, yr, date)
tpath_16_21 <- tpath %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  select(om, yr, date)

#tornado path plotting
tpath_plot <- ggplot() +
  geom_sf(data = tpath_16_21, 
          aes(color = as.factor(yr))) +
  geom_sf(data = okcounty, fill = NA) +
  scale_color_discrete(name = "Year") +
  coord_sf(datum = NA) +
  theme_void()

#tornado point plotting
tpoint_plot <- ggplot() +
  geom_sf(data = tpoint_16_21, 
          aes(color = as.factor(yr))) +
  geom_sf(data = okcounty, fill = NA) +
  scale_color_discrete(name = "Year") +
  coord_sf(datum = NA) +
  theme_void()

#plotting them on a grid
plot_grid(tpath_plot, tpoint_plot, 
          labels = c("A) Tornado Paths", 
                     "B) Tornado Points",
                     label_size = 12),
          ncol = 1, 
          hjust = 0, 
          label_x = 0, 
          align = "hv")

Question 2

Summarize the density of tornado points by both county and year and generate a faceted plot that displays maps of county-level tornado density from 2016-2021.

#Question 2

#joining the county with points so the density of points per county can
#be calculated
countypnt <- st_join(tpoint_16_21, okcounty)
countypnt <- st_drop_geometry(countypnt)
countysum <- countypnt %>%
  group_by(yr, GEOID) %>%
  summarize(tcnt = n())

countymap <- okcounty %>%
  merge(countysum, by = "GEOID") %>%
  na.omit() %>%
  mutate(area = st_area(.),
         tdens = 10^6 * 10^3 * tcnt / area) %>%
  drop_units()

#plotting the density
ggplot() +
  geom_sf(data = okcounty,
          fill = NA,
          color = "gray") +
  geom_sf(data = countymap,
          aes(fill=tdens)) +
  facet_wrap(vars(yr), ncol = 2) +
  coord_sf(datum = NA) +
  scale_fill_distiller(name=expression("Tornadoes per 1000 km"^2),
                       palette = "YlOrRd", 
                       breaks = pretty_breaks(),
                       direction = 1) +
  theme_void()

Question 3

Generate four choropleth maps of tornado density based on quantile breaks with numbers of classes ranging from 3 to 6. Create a composite figure containing the four maps using plot_grid() and examine how the number of classes changes the interpretation of the map.

#Question 3

#3 quantiles
numclass3 <- 3
qbrks3 <- seq(0,1, length.out = numclass3 + 1)
countymap_quartile3 <- countymap %>%
  mutate(tdens_c3 = cut(tdens,
                        breaks = quantile(tdens, probs = qbrks3),
                        include.lowest = T))

countymap_quartile3_plot <- ggplot() +
  geom_sf(data = countymap_quartile3,
          aes(fill = tdens_c3)) +
  scale_fill_brewer(name = expression("Tornadoes per 1000 km"^2),
                    palette = "YlOrRd") +
  theme_void() +
  theme(legend.position = "right")

#4 quantiles
numclass4 <- 4
qbrks4 <- seq(0,1, length.out = numclass4 + 1)
countymap_quartile4 <- countymap %>%
  mutate(tdens_c4 = cut(tdens,
                        breaks = quantile(tdens, probs = qbrks4),
                        include.lowest = T))
countymap_quartile4_plot <- ggplot() +
  geom_sf(data = countymap_quartile4,
          aes(fill = tdens_c4)) +
  scale_fill_brewer(name = expression("Tornadoes per 1000 km"^2),
                    palette = "YlOrRd") +
  theme_void() +
  theme(legend.position = "right")


#5 quantiles
numclass5 <- 5
qbrks5 <- seq(0,1, length.out = numclass5 + 1)
countymap_quartile5 <- countymap %>%
  mutate(tdens_c5 = cut(tdens,
                        breaks = quantile(tdens, probs = qbrks5),
                        include.lowest = T))
countymap_quartile5_plot <- ggplot() +
  geom_sf(data = countymap_quartile5,
          aes(fill = tdens_c5)) +
  scale_fill_brewer(name = expression("Tornadoes per 1000 km"^2),
                    palette = "YlOrRd") +
  theme_void() +
  theme(legend.position = "right")


#6 quantiles
numclass6 <- 6
qbrks6 <- seq(0,1, length.out = numclass6 + 1)
countymap_quartile6 <- countymap %>%
  mutate(tdens_c6 = cut(tdens,
                        breaks = quantile(tdens, probs = qbrks6),
                        include.lowest = T))
countymap_quartile6_plot <- ggplot() +
  geom_sf(data = countymap_quartile6,
          aes(fill = tdens_c6)) +
  scale_fill_brewer(name = expression("Tornadoes per 1000 km"^2),
                    palette = "YlOrRd") +
  theme_void() +
  theme(legend.position = "right")

#plotting the maps on a grid
plot_grid(countymap_quartile3_plot, countymap_quartile4_plot, 
          countymap_quartile5_plot, countymap_quartile6_plot,  
          nrow = 4,
          ncol = 1, 
          label_x = 0,
          align = "hv")

Part B

Choose any dataset of your interest (point data), and create an interactive map using the leaflet package. Customize your map, add a legend, popup information, and labels, and use an adequate color palette and scheme according to the variable mapped. These are the minimum requirements for this lab. But if you want to explore other functionalities included in the library and other datasets, please do so! Include your resulting map and code in the part B section of your RMarkdown report. Then, guide your readers through the map creation process and introduce the results.

counties_va_highres <- us_counties(states = "Virginia", resolution = "high")
starbucks = read_csv('directory.csv') %>% 
    filter(City %in% c("Alexandria", "Springfield", "Lorton")) %>%
    st_as_sf(coords = c("Longitude", "Latitude"), crs = 4326) %>% 
    select(store_name = `Store Name`, phone = `Phone Number`, address = `Street Address`, city = City, brand = Brand)

pal <- colorFactor(c("lightgreen", "orange", "turquoise"), domain = starbucks$city)
pop = paste0('<a href=', starbucks$url, '>', starbucks$store_name, "</a>")

starbucks_map = function(){
leaflet() %>%
  addProviderTiles(providers$CartoDB) %>%
  setView(lng=-77.09400, lat=38.73343, zoom = 11) %>% 
  addPolygons(data = counties_va_highres, color="black", fillOpacity = .05) %>%
  addMarkers(lng=-77.09400, lat=38.73343, popup="James Jordan's Home") %>%
  addCircleMarkers(data = starbucks, color = ~pal(city), fillOpacity = .9, radius = 5, 
                   stroke = FALSE, label = ~address,
                   popup = pop) %>%
  addLegend("bottomright", pal = pal, labels = 'Cities', values = starbucks$city)
}

starbucks_map()

Due to lack of creative thinking I decided to use the Starbucks dataset from the walkthrough. I decided to display all of the Starbucks locations near where I live. I marked my home on the map (don’t make me regret this haha) and differentiated Starbucks locations by the towns they reside in. Overlaid are the boundaries of all counties in Virginia. Pop up gives the links to the location, and hovering provides the address.