# Mapping Tornado Alley with R

boB Rudis Originally published at rud.is on ・6 min read

I caught a re-tweet of this tweet by @harry_stevens:

THREAD: I wrote a post on @observablehq about a map I made today. It shows a typical day in the life of a graphics journalist: You never know what problems you'll have to solve on deadline! https://t.co/yRhW1wbLxN #d3js #dataviz 1/7 pic.twitter.com/7N6mmK0nz3

— Harry Stevens (@Harry_Stevens) May 16, 2019

Harry’s thread and Observable post are great on their own and both show power and utility of Observable javascript notebooks.

However, the re-tweet (which I’m not posting because it’s daft) took a swipe at both Python & R. Now, I’m all for a good swipe at Python (mostly to ensure we never forget all those broken spacebars and tab keys that language has caused) but I’ll gladly defend it and R together when it comes to Getting Things Done, even on deadline.

Let’s walk through what one of us might have done had we been in the same scenario as Harry.

So, we have to create a map of historical tornado frequency trends on deadline.

We emailed researchers and received three txt files. One is a set of latitudes, another longitudes, and the final one is the trend value. It’s gridded data.

Download that ZIP and pretend you got three files in email vs a nice ZIP and make a new RStudio project called “tornado” and put those three files in a local-to-the-project-root data/ directory. Let’s read them in and look at them:

library(hrbrthemes) # not 100% necessary but i like my ggplot2 theme(s) :-)
library(tidyverse) # data wrangling & ggplot2

tibble(
lat = scan(here::here("data/lats.txt")),
lon = scan(here::here("data/lons.txt")),
trend = scan(here::here("data/trends.txt"))


You very likely never directly use the base::scan() function, but it's handy here since we just have files of doubles with each value separated by whitespace. Now, let's see what we have:

tornado

geom_point(aes(color = trend))



ggplot(tornado, aes(trend)) +
geom_histogram() +
scale_x_continuous(breaks = seq(-0.5, 0.5, 0.05))


Since we’re looking for trends (either direction) in just the United States the latitude and longitude ranges will need to be shrunk down a bit (it does indeed look like globally gridded data) and we’ll be able to shrink the data set a bit more since we only want to look at large or small tends.

We don’t really need modern R/ggplot2 mapping idioms for this project (i.e. the new {sf} ecosystem), so we’ll keep it “simple” (scare quotes since that’s a loaded term) and just use the built in maps and geom_map(). First, let’s get the U.S. states and extract their bounding boxes/limits:

maps::map("state", ".", exact = FALSE, plot = FALSE, fill = TRUE) %>%
fortify(map_obj) %>%
as_tibble() -> state_map

xlim <- range(state_map$long) ylim <- range(state_map$lat)



NOTE: I tend not to use the handy ggplot::map_data() function since it ends up clobbering purrr::map() which I use heavily (though not in this post). I also try to use {sf} these days so this tends to not be an issue anymore anyway.

Now, let’s focus in on the target area in the original paper and the Axios article:

filter(
between(lon, -107, xlim[2]), between(lat, ylim[1], ylim[2]), # -107 gets us ~left-edge of TX
((trend < -0.07) | (trend > 0.07)) # approximates notebook selection range

geom_point(aes(color = trend))



Now we’re getting close to our final solution.

As stated in the Observable notebook and implied by the word “grid” these dots are centroids of grid rectangles. This means we really want boxes, not points. The article got all fancy but it’s not really necessary since we can use ggplot2::geom_tile() to get us said boxes:

ggplot(tornado, aes(lon, lat)) +
geom_tile(aes(fill = trend, color = trend))



Now, we just need to add in map layers, and tweak some aesthetics to make it look like a map. We’ll start naively:

ggplot() +
geom_tile(
aes(lon, lat, fill = trend, color = trend)
) +
geom_map(
data = state_map, map = state_map,
aes(long, lat, map_id = region),
color = "black", size = 0.125, fill = NA
)



Our gridded data is definitely covering the right/same areas so we just need to make this more suitable for an article. We’ll use Harry’s palette and layer in U.S. state borders, an overall country border, and approximate the title and legend aesthetics:

c(
"#023858", "#045a8d", "#0570b0", "#3690c0", "#74a9cf",
"#a6bddb", "#d0d1e6", "#ece7f2", "#fff7fb", "#ffffff",
"#ffffcc", "#ffeda0", "#fed976", "#feb24c", "#fd8d3c",
"#fc4e2a", "#e31a1c", "#bd0026", "#800026"
) -> grad_cols # colors from article

ggplot() +

# tile layer

geom_tile(
aes(lon, lat, fill = trend, color = trend)
) +

# state borders

geom_map(
data = state_map, map = state_map,
aes(long, lat, map_id = region),
color = ft_cols$slate, size = 0.125, fill = NA ) + # usa border borders("usa", colour = "black", size = 0.5) + # color scales scale_colour_gradientn( colours = grad_cols, labels = c("Fewer", rep("", 4), "More"), name = "Change in tornado frequency, 1979-2017" ) + scale_fill_gradientn( colours = grad_cols, labels = c("Fewer", rep("", 4), "More"), name = "Change in tornado frequency, 1979-2017" ) + # make it Albers-ish and ensure we can fit the borders in coord_map( projection = "polyconic", xlim = scales::expand_range(range(tornado$lon), add = 2),
) +

# tweak legend aesthetics

guides(
colour = guide_colourbar(
title.position = "top", title.hjust = 0.5
),
fill = guide_colourbar(
title.position = "top", title.hjust = 0.5
)
) +
labs(
x = NULL, y = NULL
) +
theme_ipsum_rc(grid="") +
theme(axis.text = element_blank()) +
theme(legend.position = "top") +
theme(legend.title = element_text(size = 16, hjust = 0.5)) +
theme(legend.key.width = unit(4, "lines")) +
theme(legend.key.height = unit(0.5, "lines"))



### FIN

I went through some extra steps for folks new to R but the overall approach was at the very least equally as expedient as the Observable one and — despite the claims by the quite daft retweet — this is no less “shareable” or “reusable” than the Observable notebook. You can clone the repo (https://git.sr.ht/~hrbrmstr/tornado) and reuse this work immediately.

If you take a stab at an alternate approach — especially if you do use {sf} — definitely blog about it and drop a link here or on Twitter.

Posted on by:

### boB Rudis

Don't look at me…I do what he does—just slower. #rstats avuncular • Chef • Cyclist • Christian • [Master] Chef des Donné