Introduction

The story this map tells is where tornadoes touched down in Oklahoma between 2016 and 2021 and how strong they were. Most tornadoes are weak, and the strongest ones are rare, so a map that colors each tornado by its magnitude shows a state-wide field of pale, low-intensity points with a handful of darker exceptions. A second part of the story is what we do not know: about a quarter of the tornadoes have no recorded magnitude, and the map keeps them visible in gray instead of hiding them. I use the same NOAA Storm Prediction Center data as in Part A and build an interactive leaflet map, so that readers can zoom into dense areas, click a tornado to read its details, and hover over a county to see how many tornadoes it had.

Setup

# Spatial data
library(sf)

# Data wrangling
library(dplyr)

# Interactive maps
library(leaflet)

I read the tornado touchdown points and the county boundaries with st_read().

# Folder that holds the unzipped Chapter 5 shapefiles
data_dir <- "C:/Users/Wiest/Documents/chapter5_data/Chapter5/"

okcounty <- st_read(paste0(data_dir, "ok_counties.shp"), quiet = TRUE)
tpoint <- st_read(paste0(data_dir, "ok_tornado_point.shp"), quiet = TRUE)

Data preparation

Tornado points

The data cover 1950-2021, so I filter to 2016-2021 and keep the columns needed for the map: date, magnitude (mag), injuries (inj), and fatalities (fat).

In the NOAA data, a magnitude of -9 means that the magnitude is unknown. If left as -9, leaflet would treat it as the lowest value and stretch the color scale; I recode it to NA instead. About a quarter of the tornadoes (105 of 434) have an unknown magnitude, so I keep them on the map and show them in gray rather than dropping them.

tornado_pts <- tpoint |>
  filter(yr >= 2016, yr <= 2021) |>
  select(om, yr, date, mag, inj, fat) |>
  # -9 codes an unknown magnitude
  mutate(mag = na_if(mag, -9))

# Check: tornadoes per magnitude, with the unknown ones counted as NA
table(tornado_pts$mag, useNA = "ifany")
## 
##    0    1    2    3    4 <NA> 
##  118  178   27    5    1  105

County counts

For context, the map includes county outlines that also report how many tornadoes each county had. A spatial join (st_join()) finds the county of each tornado, and then tornadoes are counted per county.

county_counts <- tornado_pts |>
  # Attach county attributes to each tornado point (join by location)
  st_join(okcounty) |>
  st_drop_geometry() |>
  count(GEOID, name = "tcnt")

okcounty_cnt <- okcounty |>
  left_join(county_counts, by = "GEOID") |>
  # Counties missing from the counts had no tornadoes, so set them to zero
  mutate(tcnt = coalesce(tcnt, 0L))

Coordinate reference system

Leaflet expects longitude and latitude in WGS84 (EPSG 4326), but the shapefiles use the NAD83 datum. I transform both layers with st_transform().

tornado_map <- st_transform(tornado_pts, 4326)
okcounty_ll <- st_transform(okcounty_cnt, 4326)

# Check: should print 4326
st_crs(tornado_map)$epsg
## [1] 4326

Building the map

Palette

Magnitude is an ordered variable, so I use a sequential palette. colorFactor() gives each of the five magnitudes (0-4) its own color from the YlOrRd ramp: darker and redder means a stronger tornado. A diverging or qualitative palette would be inappropriate because the categories are ordered and have no meaningful midpoint. Unknown magnitude is shown in gray, which is neutral and clearly not part of the ramp.

pal <- colorFactor(
  palette = "YlOrRd",
  domain = 0:4,
  na.color = "gray"
)

Map layers

The map is built in layers, added in this order:

  1. A light gray basemap (Esri.WorldGrayCanvas) keeps the point colors easy to see.
  2. County polygons with a transparent fill are drawn first so they sit under the points. Hovering over a county shows its name and tornado count.
  3. Tornado points are drawn as circle markers. The fill is colored by magnitude and a thin black outline keeps pale points visible. Hovering shows the magnitude, and clicking opens a popup with the date, magnitude, injuries, and fatalities.
  4. A legend (with an “Unknown” entry), a title, and a starting view that fits the whole state complete the map.
  5. A layers control lets readers switch the county outlines and the tornado points on and off. Each layer has a group name, and addLayersControl() lists those groups as checkboxes.
leaflet(tornado_map) |>
  # Basemap
  addProviderTiles(providers$Esri.WorldGrayCanvas) |>
  # County outlines (transparent fill so that hovering inside a county works)
  addPolygons(
    data = okcounty_ll,
    group = "County outlines",
    fillOpacity = 0,
    color = "gray20",
    weight = 1.5,
    opacity = 0.8,
    label = ~paste0(
      NAME, " County: ", tcnt,
      ifelse(tcnt == 1, " tornado", " tornadoes")
    ),
    highlightOptions = highlightOptions(weight = 3, color = "black")
  ) |>
  # Tornado points colored by magnitude
  addCircleMarkers(
    group = "Tornadoes",
    fillColor = ~pal(mag),
    color = "black",
    weight = 0.5,
    radius = 4,
    fillOpacity = 0.9,
    label = ~paste0("Magnitude: ", ifelse(is.na(mag), "Unknown", mag)),
    popup = ~paste0(
      "Date: ", date, "<br>",
      "Magnitude: ", ifelse(is.na(mag), "Unknown", mag), "<br>",
      "Injuries: ", inj, "<br>",
      "Fatalities: ", fat
    )
  ) |>
  # Legend: values include NA so the "Unknown" entry is shown
  addLegend(
    position = "bottomleft",
    pal = pal,
    values = tornado_map$mag,
    title = "Tornado magnitude (EF)",
    na.label = "Unknown"
  ) |>
  # Checkboxes to toggle layers
  addLayersControl(
    overlayGroups = c("County outlines", "Tornadoes"),
    options = layersControlOptions(collapsed = FALSE)
  ) |>
  # Title
  addControl("Oklahoma tornadoes, 2016-2021", position = "bottomright") |>
  # Starting view that fits all of Oklahoma
  fitBounds(lng1 = -103.1, lat1 = 33.6, lng2 = -94.4, lat2 = 37.1)

Discussion of results

The map shows tornadoes concentrated in central and eastern Oklahoma, with fewer in the western Panhandle, which matches the static maps in Part A. Among the 329 tornadoes with a known magnitude, most are weak: 118 are EF0 and 178 are EF1, which appear as pale yellow and light orange. Only 27 are EF2, 5 are EF3, and 1 is EF4, so the darker red points that signal strong tornadoes are rare and easy to pick out against the pale background. This is why a sequential palette works well here: the eye is drawn to the few strong events.

Two cautions apply when reading the map. First, the gray points (105 tornadoes, about 24% of the total) have an unknown magnitude and should not be read as weak; the legend identifies them separately. Because they are scattered across the state, the pattern among the known magnitudes should not be assumed to hold for them. Second, in the densest areas, points overlap, so some are partly hidden; zooming in or switching the county layer off helps. The county counts in the hover labels cover 2016-2021 and include tornadoes of unknown magnitude, so they will not match the colored points exactly. Each tornado is counted in the county where it touched down, even if its path crossed into other counties.