Introduction

This document presents an analysis of tornado data in Oklahoma, including the generation of maps depicting tornado paths and points, summarizing tornado density by county and year, and creating choropleth maps to visualize tornado density variation.

The code begins by loading necessary libraries such as sf for spatial data handling, here for file path management, units for unit conversions, dplyr for data manipulation, ggplot2 for plotting, and cowplot for arranging plots.

# Load required libraries
library(sf) # For working with spatial data
library(here) # For easy file path management
library(units) # For unit conversions
library(dplyr) # For data manipulation
library(ggplot2) # For creating plots
library(cowplot) # For arranging plots

This section imports shapefile data containing Oklahoma county boundaries (ok_counties.shp), tornado paths (ok_tornado_path.shp), and tornado points (ok_tornado_point.shp).

# Import Oklahoma county, tornado paths, and points data
oklahoma_counties <- st_read(here("r_data/chapter5", "ok_counties.shp"), quiet = TRUE)
tornado_paths <- st_read(here("r_data/chapter5", "ok_tornado_path.shp"), quiet = TRUE)
tornado_points <- st_read(here("r_data/chapter5", "ok_tornado_point.shp"), quiet = TRUE)

Data Preparation

The tornado paths data for recent years (2016-2021) were extracted and filtered from the dataset. The data include the tornado identification number (om), year (yr), and date (date) of occurrence.

# Filter tornado path data for recent years (2016-2021)
tornado_paths_recent <- tornado_paths %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  select(om, yr, date)

head(tornado_paths_recent)
## Simple feature collection with 6 features and 3 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -101.988 ymin: 35.6916 xmax: -94.4104 ymax: 36.8317
## Geodetic CRS:  NAD83
##       om   yr       date                       geometry
## 1 613662 2016 2016-03-23 LINESTRING (-94.5042 35.691...
## 2 613675 2016 2016-03-30 LINESTRING (-96.0151 36.215...
## 3 613676 2016 2016-03-30 LINESTRING (-95.5523 36.729...
## 4 613677 2016 2016-03-30 LINESTRING (-95.6384 36.284...
## 5 613678 2016 2016-03-30 LINESTRING (-95.3178 36.815...
## 6 613727 2016 2016-04-15 LINESTRING (-101.988 36.74,...

The tornado point data for recent years (2016-2021) were processed to calculate tornado density by county and year. The data were joined with county boundary data to facilitate spatial analysis.

# Filter tornado point data for recent years (2016-2021)
tornado_points_recent <- tornado_points %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  select(om, yr, date)

head(tornado_points_recent)
## Simple feature collection with 6 features and 3 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -101.988 ymin: 35.6916 xmax: -94.5042 ymax: 36.8151
## Geodetic CRS:  NAD83
##       om   yr       date                 geometry
## 1 613662 2016 2016-03-23 POINT (-94.5042 35.6916)
## 2 613675 2016 2016-03-30 POINT (-96.0151 36.2151)
## 3 613676 2016 2016-03-30 POINT (-95.5523 36.7291)
## 4 613677 2016 2016-03-30  POINT (-95.6384 36.284)
## 5 613678 2016 2016-03-30 POINT (-95.3178 36.8151)
## 6 613727 2016 2016-04-15   POINT (-101.988 36.74)

1. Generate a map of tornado paths by year

This section visualizes tornado paths in Oklahoma from 2016 to 2021. The main objective is to understand the spatial distribution and trends of tornado occurrences over recent years.

# Plot tornado paths with different colors for each year
# Plot tornado paths with different colors for each year
path_plot <- ggplot() +
  geom_sf(data = tornado_paths_recent, aes(color = as.factor(yr))) +
  geom_sf(data = oklahoma_counties, fill = NA) +
  scale_color_discrete(name = "Year") +
  coord_sf(datum = NA) +
  theme(plot.caption = element_text(hjust = 0.5)) +
  labs(title = "Tornado Paths by Year",
       x = "Longitude",
       y = "Latitude",
       caption = "The generated map illustrates tornado paths for each year within the specified timeframe (2016-2021). Each year \n is represented by a distinct color, allowing for easy visual differentiation of tornado occurrences. The map provides insights \ninto the spatial distribution and clustering of tornadoes over recent years, aiding in the identification of high-risk areas.") +
  theme_void()

# Display tornado paths by year
path_plot

2. Summarize tornado density by county and year

This section analyzes tornado density by county and year in Oklahoma from 2016 to 2021. The objective is to identify regions with higher tornado occurrence rates and explore temporal variations in tornado density.

# Join tornado point data with county boundary data to calculate tornado density
tornado_density <- st_join(oklahoma_counties, tornado_points_recent) %>%
  filter(!is.na(yr)) %>%  # Remove NA from yr
  group_by(COUNTYFP, yr) %>%
  summarise(total_tornado_points = n()) %>%
  mutate(area = st_area(geometry),
         tornado_density = 10^6 * 10^3 * total_tornado_points / area) %>%
  drop_units()

# Display tornado density data summary
glimpse(tornado_density)
## Rows: 217
## Columns: 6
## Groups: COUNTYFP [75]
## $ COUNTYFP             <chr> "001", "001", "001", "001", "001", "005", "005", …
## $ yr                   <dbl> 2016, 2017, 2019, 2020, 2021, 2016, 2019, 2017, 2…
## $ total_tornado_points <int> 1, 1, 1, 1, 2, 2, 1, 1, 1, 2, 3, 1, 3, 1, 1, 1, 1…
## $ geometry             <POLYGON [°]> POLYGON ((-94.8073 35.6386,..., POLYGON (…
## $ area                 <dbl> 1493619771, 1493619771, 1493619771, 1493619771, 1…
## $ tornado_density      <dbl> 0.6695144, 0.6695144, 0.6695144, 0.6695144, 1.339…
# Generate a faceted plot to display maps of county-level tornado density from 2016-2021
density_plot <- ggplot() +
  geom_sf(data = oklahoma_counties, fill = "white", color = "black") +  # Add county polygons
  geom_sf(data = tornado_density, aes(fill = tornado_density)) +  # Add tornado density data
  scale_fill_gradient(name = "Tornado Density (points/km²)",
                      low = "lightblue",
                      high = "darkred") + # Gradient scale
  facet_wrap(~ yr, labeller = labeller(yr = as.character)) + # Facet by year
  theme_bw() + # Theme settings
  theme(legend.position = "bottom", plot.caption = element_text(hjust = 0)) + 
  labs(title = "County-Level Tornado Density from 2016-2021",
       x = "Longitude",
       y = "Latitude",
       caption = "The faceted plot displays county-level tornado density from 2016 to 2021.Each facet represents a year, \nwith tornado density depicted using a color gradient. Regions with higher tornado density are indicated \nby warmer colors (e.g., red), allowing for easy identification of areas prone to tornado occurrences.")

# Display tornado density plots
density_plot

3. Generate choropleth maps of tornado density

Choropleth Maps iterates over different numbers of classes (3 to 6) to create choropleth maps of tornado density. It defines breaks for tornado density classes based on quantiles and classifies tornado density accordingly. Furthermore, in this section choropleth maps are generated with varying numbers of classes, each representing tornado density using distinct colors.

# Create choropleth maps with different numbers of classes
maps <- list()
for (num_classes in 3:6) {
  # Define breaks for tornado density classes
  breaks <- quantile(unique(tornado_density$tornado_density), probs = seq(0, 1, length.out = num_classes + 1))
  
  # Classify tornado density
  tornado_density <- tornado_density %>%
    mutate(tornado_density_class = cut(tornado_density, breaks = breaks, include.lowest = TRUE))
  
  # Generate choropleth map
  map <- ggplot(tornado_density) +
    geom_sf(aes(fill = tornado_density_class)) +
    scale_fill_brewer(name = "Tornado Density (classes)",
                      palette = "YlOrRd") +
    theme_void() +
    theme(legend.position = "right", plot.title = element_text(hjust = 0.5), plot.caption = element_text(hjust = 0)) +
    labs(title = paste("Choropleth Map with", num_classes, "Classes"),
         x = "Longitude",
         y = "Latitude",
         caption = paste("Tornado density is categorized into", num_classes, "classes based on quantile breaks. Choropleth map were generated \n to visualize tornado density in Oklahoma. The maps categorize tornado density into distinct classes \n based on quantile breaks, allowing for the identification of regions with similar tornado occurrence rates. \n The number of classes affects the granularity of tornado density representation, providing flexibility in data interpretation."))
  
  maps[[num_classes]] <- map
  
  # Display the map
  print(map)
}

Choropleth maps with varying numbers of classes were generated to visualize tornado density in Oklahoma. The maps categorize tornado density into distinct classes based on quantile breaks, allowing for the identification of regions with similar tornado occurrence rates. The number of classes affects the granularity of tornado density representation, providing flexibility in data interpretation.

Conclusion

The analysis of tornado data in Oklahoma provides valuable insights into the spatial and temporal patterns of tornado occurrences. By visualizing tornado paths, density by county and year, and creating choropleth maps, we have gained a better understanding of tornado activity in the region. The findings can inform disaster preparedness and mitigation strategies, helping to minimize the impact of tornadoes on communities. Further research could explore additional factors influencing tornado occurrence patterns and refine predictive models for enhanced risk management.