title: “L4 Report” author: “Jing Zhang” date: “2024-03-22” output: html_document

Lesson 4 Report

Tornado in USA

This week, we are going to explore tornado data in USA from year 2016 to 2021.

My first step was to load all the packages and data for this project.

#loading all packages
library(sf)           
library(ggplot2)      
library(dplyr)        
library(tidyr)        
library(scales)       
library(RColorBrewer) 
library(units)
library(cowplot)
#reading in tornado shapefiles
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)

Since there are many years of record in the data, we want to narrow it down to five years, which are 2016 to 2021.

Looking at data

Part one of this report will show you general tornado path and location maps.

#mapping tornado path with different color by year
pathyr<-ggplot() +
  geom_sf(data = okcounty, fill = NA) +
  geom_sf(data = tpath_16_21, 
          aes(color =as.factor(yr)), 
          size = 1) +
  scale_color_discrete(name = "Year") +
  theme_bw()
  
#mapping tornado point by year
pointyr<-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 one graph
plot_grid(pathyr, pointyr, 
          labels = c("A) Tornado Paths Map", 
                     "B) Tornado Points Map",
                     label_size = 12),
          ncol =2, 
          hjust =0, 
          label_x = 0, 
          align = "hv")
Figure 1: Left : a map of historical tornado path in USA by year. Right : a map of historical tornado events in USA by year.

Figure 1: Left : a map of historical tornado path in USA by year. Right : a map of historical tornado events in USA by year.

Density Research

Part two of this project will present tornado density in each county by year.

#preparing data for density calculation by county
countypnt <- st_join(tpoint_16_21, okcounty)
countypnt <- st_drop_geometry(countypnt)
countysum <- countypnt %>%
  group_by(GEOID,yr) %>%
  summarize(tcnt = n())
countyarea<- okcounty%>%mutate(area=st_area(okcounty))

#calculation of density by county
countymap <- countyarea %>%
  full_join(countysum, by = "GEOID") %>%
  na.omit(TRUE)%>%
  mutate(
         tdens = 10^6 * 10^3 * tcnt / area) %>%drop_units()
#mapping density of each county by year
ggplot() +geom_sf(data = okcounty, fill = NA) +
  geom_sf(data = countymap, aes(fill=tdens)) +
  scale_fill_distiller(name = expression("Tornadoes/1000 km"^2), 
                      palette = "Purples", 
                      breaks = pretty_breaks(),
                      direction = 1) +
  facet_wrap(vars(yr), ncol = 2) +
  theme_void()
Figure 2: Tornado density in each county by year.

Figure 2: Tornado density in each county by year.

Visual Explorations

Part three explores classification changes in displaying the graphs.

#setting quantile breaks for 3,4,5,6 breaks
numclas <- 3
qbrks3 <- seq(0, 1, length.out = numclas+1)

countymap3 <- countymap %>%
  mutate(tdens_c3 = cut(tdens,
                        breaks = quantile(tdens,probs = qbrks3),
                        include.lowest = T))
#4
numclas <- 4
qbrks4<- seq(0,1, length.out = numclas+1 )


countymap4 <- countymap %>%
  mutate(tdens_c4 = cut(tdens,
                        breaks = quantile(tdens,probs = qbrks4),
                        include.lowest = T))
#5
numclas <- 5
qbrks5 <- seq(0, 1, length.out = numclas +1)


countymap5 <- countymap %>%
  mutate(tdens_c5 = cut(tdens,
                        breaks = quantile(tdens,probs = qbrks5),
                        include.lowest = T))

#6
numclas <- 6
qbrks6 <- seq(0,1, length.out = numclas+1 )


countymap6 <- 
  
  countymap %>%
  mutate(tdens_c6 = cut(tdens,
                        breaks = quantile(tdens,probs = qbrks6),
                        include.lowest = T))

#mapping by different class breaks
map3<-ggplot(data = countymap3) +
  geom_sf(aes(fill = tdens_c3)) +
  scale_fill_brewer(name = expression("Tornadoes/1000 km"^2),   
                    palette = "YlOrRd") +
  theme_void() +
  theme(legend.position = "bottom") 

map4<-ggplot(data = countymap4) +
  geom_sf(aes(fill = tdens_c4)) +
  scale_fill_brewer(name = expression("Tornadoes/1000 km"^2),   
                    palette = "YlOrRd") +
  theme_void() +
  theme(legend.position = "bottom")

map5<-ggplot(data = countymap5) +
  geom_sf(aes(fill = tdens_c5)) +
  scale_fill_brewer(name = expression("Tornadoes/1000 km"^2),   
                    palette = "YlOrRd") +
  theme_void() +
  theme(legend.position = "bottom")

map6<-ggplot(data = countymap6) +
  geom_sf(aes(fill = tdens_c6)) +
  scale_fill_brewer(name = expression("Tornadoes/1000 km"^2),   
                    palette = "YlOrRd") +
  theme_void() +
  theme(legend.position = "bottom")

#ploting 4 maps on one graph
plot_grid(map3,map4,map5,map6, 
          labels = c("A) 3 Classes", 
                     "B) 4 Classes",
                     "C) 5 Classes",
                     "D) 6 Classes",
                     label_size = 12),
          ncol = 1, 
          hjust = 0, 
          label_x = 0, 
          align = "hv")
Figure 3: mapping data with different classification breaks.

Figure 3: mapping data with different classification breaks.

Setting Up A Map

Displaying tornado data on leaflet map

Finally, we have map the events on a map and you can explore the interactive map by clicking on each of the point.

#import all packages
library(leaflet)
library(sf)
library(sp)
#try to re-project data into WGS84
stpoint<- st_transform(tpoint,CRS("+proj=longlat +init=epsg:4326 +ellps=WGS84 +datum=WGS84"))
#filter data into year 2016-2021
tpoint_16_21 <- stpoint %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  select(om, yr, date, mag)
  
#setting color pal for leaflet map
pal <- colorFactor("YlOrRd", domain =unique(tpoint_16_21$mag),reverse = F)
#drawing the map with popup, legend and labels.
leaflet(data = tpoint_16_21) %>% addProviderTiles(providers$CartoDB) %>% 
  addCircleMarkers(color = ~pal(mag), fillOpacity = .5, stroke = FALSE,popup =~as.character(yr) )%>%
  addLegend(position = "bottomright", pal = pal,values = ~mag ,title = "Magnitude")

Figure 4: Map of Tornado Events from 2016 to 2021.

Conclusion

In conclusion, hope the visualizations and figures in this report can give you a informative story of tornado events happend in the USA during the year 2016 to year 2021.