GEO 588 Lab 4 Part A
Richard Wetzel
2024-03-22
#Load Libraries-----
library(sf)## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUElibrary(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)
library(scales)
library(RColorBrewer)
library(units)## udunits database from C:/Users/rwetzel/AppData/Local/Programs/R/R-4.3.2/library/units/share/udunits/udunits2.xmllibrary(cowplot)
library(here)## here() starts at C:/Users/rwetzel/OneDrive - Granger Waste Services/Documents/Richard Wetzel/Administrative/MGIS/588/Lab 4library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.5
## ✔ ggplot2 3.5.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ readr::col_factor() masks scales::col_factor()
## ✖ purrr::discard() masks scales::discard()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ lubridate::stamp() masks cowplot::stamp()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsIntroduction
This exercise focuses on conducting basic spatial analysis using geographic data and producing some basic maps of the results of that analysis. Specifically, the focus is on analyzing and mapping tornadoes occurring in Oklahoma between 2016 and 2021.
#Read in 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 datasets to only include the desired years (2016-2021)
#Select relevant fields
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)
#Perform a spatial join of the tornado points layer with the county layer
#Summarize the number of tornadoes in each county and calculate the density
#of tornadoes per square kilometer
countypnt <- st_join(tpoint_16_21, okcounty)
countypnt <- st_drop_geometry(countypnt)
countysum <- countypnt %>%
group_by(GEOID) %>%
summarize(tcnt =n())
countymap <- okcounty %>%
left_join(countysum, by = "GEOID") %>%
replace(is.na(.),0) %>% #Replace null values with 0 for counties with no tornadoes
mutate(area=st_area(okcounty),
tdens=10^6 * 10^3 * tcnt/area) %>%
drop_units() #Removes the units designator from the fieldExercise 1
The first exercise involves creating a composite figure containing two separate maps: one containing the touchdown points for each tornado and another containing the tornado’s path, with data symbolized by the year in which they occurred.
tpointcolored <-ggplot() +
geom_sf(data = okcounty,
fill = NA,
color = "gray") +
geom_sf(data = tpoint_16_21, size = 0.75,
aes(color=as.factor(yr))) +
scale_color_discrete(name="Year")+
coord_sf(datum=NA)+
theme_void()
tpathcolored <-ggplot() +
geom_sf(data = okcounty,
fill = NA,
color = "gray") +
geom_sf(data = tpath_16_21, size = 0.75,
aes(color=as.factor(yr))) +
scale_color_discrete(name="Year")+
coord_sf(datum=NA)+
theme_void()
plot_grid (tpathcolored, tpointcolored,
labels=c("Tornado Paths",
"Tornado Points",
label_size=10),
ncol=1,
hjust=0,
label_x=0,
align="hv")
Figure 1. Tornados appear to occur more frequently in a band of Oklahoma running from the southwest corner to the northeast. Their paths also seem to align with this.
Exercise 2
The second exercise consists of summarizing the density of tornado points by both county and year and creating a faceted map in which each year’s values are plotted separately in one image.
#Perform spatial join of tornado point and county layer
#Summarize the number of tornadoes by county and by year to calculate density of
#tornadoes by square kilometer per year
countypnt <- st_join(tpoint_16_21, okcounty)
countypnt <- st_drop_geometry(countypnt)
countysum <- countypnt %>%
group_by(GEOID,yr) %>%
summarize(tcnt = n())## `summarise()` has grouped output by 'GEOID'. You can override using the
## `.groups` argument.countymap2 <- okcounty %>%
left_join(countysum, by = "GEOID") %>%
replace(is.na(.), 0) %>%
filter(yr!=0)%>%
mutate(area = st_area(.),
tdens=10^6*10^3*tcnt/area)%>%
drop_units()
#Create a faceted map which displays the density of tornadoes by county, by year.
ggplot(data=countymap2) +
ggtitle("Tornados per 1000 SqKm by Year")+
geom_sf(data=okcounty, fill=NA)+
geom_sf(aes(fill=tdens))+
scale_fill_distiller(name=expression("Tornadoes per 1000 SqKm"),
palette = 'YlOrRd',
breaks=pretty_breaks(),
direction=1)+
theme_void()+
theme(legend.position="bottom")+
facet_wrap(~yr)
Figure 2. 2019 appears to have areas with higher densities of tornados by area, although this may be a function of aggregating the data to irregularly sized counties. Performing a spatial analysis in which tornado density is calculated as part of a moving kernel may produce different results
Exercise 3
The third exercise consists of creating four choropleth maps of tornado density based on quantile breaks with varying numbers of classes, ranging from 3 to 6 classes. Combining these maps into one image allows for a visual comparison of how the change in number of classes affects how the data is interpreted.
#Create 4 sets of break values for chloropleth mapping with the number of classes #between 3 and 6.
numclass_3 <- 3
qbreaks_3 <- seq(0,1,length.out = numclass_3 +1)
qbreaks_3## [1] 0.0000000 0.3333333 0.6666667 1.0000000numclass_4 <- 4
qbreaks_4 <- seq(0,1, length.out= numclass_4 +1)
qbreaks_4## [1] 0.00 0.25 0.50 0.75 1.00numclass_5 <- 5
qbreaks_5 <-seq(0,1, length.out=numclass_5 +1)
qbreaks_5## [1] 0.0 0.2 0.4 0.6 0.8 1.0numclass_6 <-6
qbreaks_6 <-seq(0,1, length.out=numclass_6 +1)
qbreaks_6## [1] 0.0000000 0.1666667 0.3333333 0.5000000 0.6666667 0.8333333 1.0000000#Create 4 different maps, one for each of the four quantile breaks created above
countymap_3 <-countymap %>%
mutate(tdens_c1=cut(tdens,
breaks=quantile(tdens, probs=qbreaks_3),
include.lowest=T))
q3_map <-ggplot(data=countymap_3)+
geom_sf(aes(fill=tdens_c1))+
scale_fill_brewer(name='', palette= "YlOrRd")+
theme_void()+
theme(legend.position="bottom",
text=element_text(color="525252"))
countymap_4 <-countymap %>%
mutate(tdens_c1=cut(tdens,
breaks=quantile(tdens, probs=qbreaks_4),
include.lowest=T))
q4_map <-ggplot(data=countymap_4)+
geom_sf(aes(fill=tdens_c1))+
scale_fill_brewer(name='', palette= "YlOrRd")+
theme_void()+
theme(legend.position="bottom",
text=element_text(color="525252"))
countymap_5 <-countymap %>%
mutate(tdens_c1=cut(tdens,
breaks=quantile(tdens, probs=qbreaks_5),
include.lowest=T))
q5_map <-ggplot(data=countymap_5)+
geom_sf(aes(fill=tdens_c1))+
scale_fill_brewer(name='', palette= "YlOrRd")+
theme_void()+
theme(legend.position="bottom",
text=element_text(color="525252"))
countymap_6 <-countymap %>%
mutate(tdens_c1=cut(tdens,
breaks=quantile(tdens, probs=qbreaks_6),
include.lowest=T))
q6_map <-ggplot(data=countymap_6)+
geom_sf(aes(fill=tdens_c1))+
scale_fill_brewer(name='', palette= "YlOrRd")+
theme_void()+
theme(legend.position="bottom",
text=element_text(color="525252"))
#Plot the four maps in a 2x2 grid to compare how using each set of quantile breaks
#affects the interpretation of the map
plot_grid(q3_map, q4_map, q5_map, q6_map,
labels=c ("Three Classes",
"Four Classes",
"Five Classes",
"Six Classes",
label_size=10),
ncol = 2,
label_x = 0,
align = "hv")
Figure 3. Quantile classification in choropleth mappping result in features being distributed evenly across each available class. By adding more breaks, this can result in data appearing to be more variable between counties than is actually true.