Assignment 4
Madison Quill
2026-03-27
Part A
The first part of this assignment involved using ggplot to create static maps of tornado locations and paths from 2016 to 2021. Three different figures were created in this section.
Load Libraries
Before the analysis can be done it is necessary to load all relevant packages for this analysis.
library(sf)
library(ggplot2)
library(dplyr)
library(tidyr)
library(scales)
library(RColorBrewer)
library(units)
library(cowplot)Tornado Paths and Points by Year
This part of the assignment involves creating a map of tornado paths where the paths from each year are displayed as a different color and creating a composite figure containing the map of tornado paths and the map of tornado points by year.
To accomplish this, path and point tornado data first needs to be loaded and filtered to the appropriate years. County polygons will also be loaded for background reference.
#Load Data
okcounty <- st_read("ok_counties.shp", quiet = TRUE)
tpoint <- st_read("ok_tornado_point.shp", quiet = TRUE)
tpath <- st_read("ok_tornado_path.shp", quiet = TRUE)
#Filter to relevant dates for points and paths
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)Next, maps of tornado path and points need to be created using ggplot and the sf package. These maps will be saved for use in the composite figure later.
#Create map of tornado points by year
ptgraph <- 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()
#Create map of tornado paths by year
linegraph <- ggplot() +
geom_sf(data = tpath_16_21, size = 1,
aes(color = as.factor(yr))) +
geom_sf(data = okcounty, fill = NA) +
scale_color_discrete(name = "Year") +
coord_sf(datum = NA) +
theme_void()Finally, a composite figure containing both of the created maps is generated using plot_grid().
#Create a composite figure containing the map of tornado paths and the map of tornado points using plot_grid().
plot_grid(ptgraph, linegraph,
labels = c("A) Tornado Points by Year",
"B) Tornado Paths by Year",
label_size = 12),
ncol = 1,
hjust = 0,
label_x = 0,
align = "hv")
Tornado Path Densities by Year
The second part of part A involves summarizing the density of tornado points by both county and year and using these densities to generate a faceted plot that displays maps of county-level tornado density from 2016-2021. We can use the previously created and filtered tornado points dataset and the previously loaded in county dataset in order to create a choropleth map of tornado density for each year.
First, the area for each county needs to be calculated. Then the tornado point data needs to be joined to the county data to calculate the county which each point is in. This joined dateset will be used to calculate a summary of tornados by year and county.
#Add Area to each county
okcounty <- okcounty %>%
mutate(area = st_area(okcounty))%>%
drop_units()
#Generate counts of tornado's by county and year
countypnt <- st_join(tpoint_16_21, okcounty)#join counties to tornados
countypnt <- st_drop_geometry(countypnt)#drop geometry so summary data can be created
countysum <- countypnt %>% #create a summary of tornados by county and year
group_by(GEOID,yr) %>%
summarize(tcnt = n())Next, the summary table needs to be joined to the county geometry and then the density for each county and year can be calculated using the counts of tornados by year and county and the area for each county.
#Join county counts to geometry and calculate density
countymap <- okcounty %>%
inner_join(countysum, by = "GEOID") %>%
replace(is.na(.), 0) %>%
ungroup()%>%
mutate(tdens = 10^6 * 10^3 * tcnt / area) %>%
drop_units()Finally, a plot can be created showing the density of tornados by county for each year.
#Create map
ggplot(data = countymap) +
geom_sf(data = okcounty, fill = NA) +#add background county data
geom_sf(aes(fill = tdens)) +
facet_wrap(vars(yr), ncol = 2)+#separate maps by year
scale_fill_distiller(name = expression("Tornadoes/1000 km"^2),
palette = "YlOrRd",
breaks = pretty_breaks(),
direction = 1) +
theme_void()+
labs(title = "Density of tornado points by both county and year 2016 to 2021")#add titles
Tornado density using different class breaks.
This part of the analysis involves generating different choropleth maps using different class breaks and adding them to the same plot for comparison of how class breaks change how a map can be interpreted.
First, total county densities needed to be generated in a similar way as in the previous section except only by county and not year. The previously calculated county areas were used for this calculation as well.
#Generate Density Numbers
countypnttot <- st_join(tpoint_16_21, okcounty)
countypnttot <- st_drop_geometry(countypnt)
countysumtot <- countypnttot %>%
group_by(GEOID) %>%
summarize(tcnt = n())
countybreaks <- okcounty %>%
left_join(countysumtot, by = "GEOID") %>%
replace(is.na(.), 0) %>%
mutate(area = st_area(okcounty),
tdens = 10^6 * 10^3 * tcnt / area) %>%
drop_units()Next class breaks can be created for numbers from 3-6. A sequence for the class breaks can be created and then used to cut the density numbers at specific quantile intervals and store these classes in a new column of the data. This new column can be used to create a map plot with class breaks. This plot is saved for later use in a composite plot.
# Create Three Breaks
numclas3 <- 3
qbrks3 <- seq(0, 1, length.out = numclas3 + 1)
#Add class breaks field to dataset
countybreaks <- countybreaks %>%
mutate(tdens_c3 = cut(tdens,
breaks = quantile(tdens, probs = qbrks3),
include.lowest = T))
#Create map using 3 class breaks
breaks3 <- ggplot(data = countybreaks) +
geom_sf(aes(fill = tdens_c3)) +
scale_fill_brewer(name="Density",
palette = "YlOrRd") +
theme_void()
#Four Breaks
numclas4 <- 4
qbrks4 <- seq(0, 1, length.out = numclas4 + 1)
countybreaks <- countybreaks %>%
mutate(tdens_c4 = cut(tdens,
breaks = quantile(tdens, probs = qbrks4),
include.lowest = T))
breaks4 <- ggplot(data = countybreaks) +
geom_sf(aes(fill = tdens_c4)) +
scale_fill_brewer(name="Density",
palette = "YlOrRd") +
theme_void()
#Five Breaks
numclas5 <- 5
qbrks5 <- seq(0, 1, length.out = numclas5 + 1)
countybreaks <- countybreaks %>%
mutate(tdens_c5 = cut(tdens,
breaks = quantile(tdens, probs=qbrks5),
include.lowest = T))
breaks5 <- ggplot(data = countybreaks) +
geom_sf(aes(fill = tdens_c5)) +
scale_fill_brewer(name="Density",
palette = "YlOrRd") +
theme_void()
#Six Breaks
numclas6 <- 6
qbrks6 <- seq(0, 1, length.out = numclas6 + 1)
countybreaks <- countybreaks %>%
mutate(tdens_c6 = cut(tdens,
breaks = quantile(tdens, probs = qbrks6),
include.lowest = T))
breaks6 <- ggplot(data = countybreaks) +
geom_sf(aes(fill = tdens_c6)) +
scale_fill_brewer(name="Density",
palette = "YlOrRd") +
theme_void()Finally a composite plot showing density maps with 3, 4, 5, and 6 class breaks is created as you can see the distribution of the data can look different based of where class breaks are made. More class breaks allows for more distinguishment between areas.
plot_grid(breaks3, breaks4, breaks5, breaks6,
labels = c("A) Three Breaks",
"B) Four Breaks",
"C) Five Breaks",
"D) Six Breaks",
label_size = 12),
ncol = 2,
hjust = 0,
label_x = 0,
align = "hv")
Part B
For the second part of this assignment I used a csv of point data of flash flooding events in the state of Florida to create an interactive map using R’s leaflet package.
Loading Libraries
A few libraries will be used in part B that were not used in part A. These are loaded below.
library(leaflet)
library(tidyverse)
library(leafem)
library(leafpop)Creating an Interactive Leaflet Map
To create this map, the flooding data needs to be read in as spatial data using the appropriate latitude and longitude fields. Additionally, only the necessary fields will be selected for this map.
#Read in point data
(floods = read_csv('floods.csv') |>
st_as_sf(coords = c("BEGIN_LON", "BEGIN_LAT"), crs = 4326) |>
select(location = `BEGIN_LOCATION`, bdate = `BEGIN_DATE`, damage = `DAMAGE_PROPERTY_NUM`))## Simple feature collection with 500 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -87.49 ymin: 25.5687 xmax: -80.0501 ymax: 30.9876
## Geodetic CRS: WGS 84
## # A tibble: 500 × 4
## location bdate damage geometry
## <chr> <chr> <dbl> <POINT [°]>
## 1 JENSEN BEACH 2/28/2015 105000 (-80.2 27.26)
## 2 LEMON CITY 2/28/2015 3000 (-80.1954 25.8111)
## 3 CHIPLEY 4/12/2015 0 (-85.537 30.7881)
## 4 STEELE CITY 4/12/2015 2000 (-85.4158 30.7244)
## 5 COTTONDALE 4/12/2015 2000 (-85.3613 30.8221)
## 6 COTTONDALE 4/12/2015 2000 (-85.4149 30.797)
## 7 COTTONDALE 4/12/2015 2000 (-85.3467 30.8385)
## 8 SIMSVILLE 4/12/2015 5000 (-85.2401 30.6681)
## 9 KYNESVILLE 4/12/2015 5000 (-85.3111 30.6771)
## 10 ROUND LAKE 4/12/2015 5000 (-85.3954 30.6597)
## # ℹ 490 more rowsWith this data, we can create an interactive map showing all the locations of flash floods. Symbol color is used to represent the level of damage caused by each event. Hovering over an event shows the beginning date of the event. Clicking on a point shows a popup containing information on the pertinent fields.
#Color map creation
pal <- colorNumeric(palette = "YlOrRd",domain = floods$damage)
#Create Map
leaflet(data=floods) |>
addTiles() |> #Add open street map basemap
addCircleMarkers(
color = ~pal(damage), fillOpacity = .50, #add color pallete based on damage value for symbols
stroke = FALSE,
label = ~bdate, #add label
popup = leafpop::popupTable(st_drop_geometry(floods[,1:4]), feature.id = FALSE, row.numbers = FALSE) #add popup table
) |>
addLegend("bottomleft",pal=pal, values = ~damage)#Add legend