1 Introduction

This report completes the three practice exercises from Wimberly’s chapter on vector geospatial data. The analyses use tornado locations and paths in Oklahoma from 2016 through 2021. The report first compares the spatial distribution of tornado points and paths, then maps annual county-level tornado density, and finally examines how quantile classification with different numbers of classes affects a choropleth map.

1.1 Data and required objects

2 Part A

2.1 Exercise 1: Map tornado points and paths by year

The first exercise compares the locations of individual tornado observations with the paths associated with those tornadoes. A separate color is assigned to each year so that the spatial distribution can be compared across the 2016–2021 study period. County boundaries are added as outlines to provide geographic context without obscuring the tornado features.

2.1.1 Map tornado points

The point map uses the tornado point layer and shows each year as a different color.

# create map of tornado points
pointmap <- 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() + 
  theme(palette.colour.discrete = scales::pal_brewer(palette = "Set1"),
        plot.margin = margin(t = 20, r = 15, b = 20, l = 15, unit = "mm")) +
  labs(title = "Oklahoma tornado locations",
       subtitle = "2016-2021")

2.1.2 Map tornado paths

The tornado path map shows trajectories of tornadoes colorized by year.

# plot of tornado paths with years as different colors

pathmap <- ggplot() +
  geom_sf(data = tpath_16_21, 
          aes(color = as.factor(yr)), size = 2) +
  geom_sf(data = okcounty, fill = NA) +
  scale_color_discrete(name = "Year") +
  coord_sf(datum = NA) +
  theme_void() + 
  theme(palette.colour.discrete = scales::pal_brewer(palette = "Set1"),
        plot.margin = margin(t = 20, r = 15, b = 20, l = 15, unit = "mm")) +
  labs (title = "Oklahoma tornado paths",
        subtitle = "2016-2021")

2.1.3 Combine the maps

The two maps are arranged vertically with plot_grid() so that the maps are easily compared.

plot_grid(
  pointmap,
  pathmap,
  labels = c("A) Tornado points", "B) Tornado paths"),
  label_size = 12,
  ncol = 1,
  hjust = 0,
  label_x = 0,
  align = "v")
Figure 1. Locations of tornadoes in Oklahoma from 2016–2021 shown as (A) individual tornado points and (B) tornado paths. Colors represent the year of occurrence.

Figure 1. Locations of tornadoes in Oklahoma from 2016–2021 shown as (A) individual tornado points and (B) tornado paths. Colors represent the year of occurrence.

Results Both maps show the locations of tornadoes colorized by year, but the path map shows the start and end points of each tornado connected with a line. Some tornadoes did not move a long distance and are difficult to see on the path map whereas all tornado locations are equally visible on the point map.

2.2 Exercise 2: County-level tornado density by year

The second exercise summarizes tornado point counts by county and year. County area is then used to calculate the tornado density per county. Counties with no recorded tornado points in a given year are retained and assigned a count of zero. This is important because otherwise those counties would appear as missing rather than as having no recorded tornado points.

2.2.1 Count tornado points by county and year

The data are grouped by GEOID and yr, summarized by number of tornadoes. complete() adds any missing county-year combinations and assigns them a count of zero.

tornadoyrs <- countypnt %>%
  st_drop_geometry() %>%
  group_by(GEOID, yr) %>%
  summarise(tcnt = n(), .groups = "drop") %>%
  complete(GEOID = okcounty$GEOID, yr, fill = list(tcnt = 0))

2.2.2 Calculate county area and tornado density

The area of each county is calculated from the polygon geometry using st_area(). The area multiplied by the number of tornadoes per county is used to calculate the density of tornadoes per 1,000 square kilometers.

# Calculate county area in square meters.
okcounty <- okcounty %>% 
  mutate(area = st_area(geometry)) %>% 
  drop_units()

# Join the counts and county polygons, then calculate density.
countyyrs <- tornadoyrs %>% 
  replace(is.na(.), 0) %>%
  left_join(okcounty, 
            by = "GEOID") %>% 
  mutate(tdensyr = 10^6 * 10^3 * tcnt / area) %>%
  st_as_sf()

2.2.3 Faceted maps of annual density

The maps below use the same continuous color scale across all six years, which allows the density values to be compared.

# density of tornado points by county and year as faceted maps

ggplot() +
  geom_sf(data = countyyrs,
          aes(fill = tdensyr)) +
  scale_fill_distiller(name = expression("Tornadoes/1000 km"^2), 
                       palette = "Spectral", 
                       breaks = pretty_breaks(),
                       direction = -1) +
  theme_void() +
  theme(legend.position = "bottom") + 
  facet_wrap(vars(yr), ncol = 3) + 
  labs(title = "Maps of tornado density by county in Oklahoma",
       subtitle = "(Tornadoes/km2) for each year from 2016 to 2021")
Figure 2. County-level tornado point density in Oklahoma for each year from 2016–2021. Density is the number of tornado points per 1,000 km² of county area. Counties without recorded tornado points are represented by zero.

Figure 2. County-level tornado point density in Oklahoma for each year from 2016–2021. Density is the number of tornado points per 1,000 km² of county area. Counties without recorded tornado points are represented by zero.

Results Because the maps use the same color scale, counties can be compared across years. Density shows the standardized tornado impact per county, regardless of county size. It appears that 2019 had higher tornado activity across the north and eastern portions of the state and other years had moderate activity in some counties.

2.3 Exercise 3: Compare quantile classifications with 3–6 classes

The third exercise maps tornado density for the full 2016–2021 period using different quantile classification. The four maps use three, four, five, and six classes, respectively.

For this exercise, density is calculated across the entire study period rather than separately for each year. County counts are summed across all years, joined to county area, and converted to tornado points per 1,000 km².

2.3.1 Calculate overall county tornado density

# Count all tornado point records in each county across the study period.
county_counts <- countypnt %>%
  st_drop_geometry() %>%
  group_by(GEOID) %>%
  summarise(tcnt = n(), .groups = "drop")

# calculate tornado density per county
countymap <- okcounty %>%
  left_join(county_counts, by = "GEOID") %>%
  replace(is.na(.), 0) %>%
  mutate(area = st_area(okcounty),
         tdens = 10^6 * 10^3 * tcnt / area) %>%
  drop_units()

2.3.2 Calculate quantile breaks

The probability values below define the requested number of quantile classes. For example, three classes use the 0th, 33rd, 67th, and 100th percentiles. The same approach is repeated for four, five, and six classes.

Tied density values can produce repeated quantile break values. The code therefore removes duplicate breaks before applying cut(). If many counties have the same density, particularly zero, the resulting map may contain fewer distinct classes than requested. The actual number of classes should therefore be checked when interpreting the maps.

# quantile breaks for 3 classes
qbrks3 <- seq(0, 1, length.out = 3 + 1)

# quantile breaks for 4 classes
qbrks4 <- seq(0, 1, length.out = 4 + 1)

# quantile breaks for 5 classes
qbrks5 <- seq(0, 1, length.out = 5 + 1)

# quantile breaks for 6 classes
qbrks6 <- seq(0, 1, length.out = 6 + 1)

# Calculate unique quantile break values for each classification.
breaks3 <- unique(quantile(countymap$tdens, probs = qbrks3, na.rm = TRUE))
breaks4 <- unique(quantile(countymap$tdens, probs = qbrks4, na.rm = TRUE))
breaks5 <- unique(quantile(countymap$tdens, probs = qbrks5, na.rm = TRUE))
breaks6 <- unique(quantile(countymap$tdens, probs = qbrks6, na.rm = TRUE))

# Check the number of distinct intervals available.
c(
  classes_3_requested = length(breaks3) - 1,
  classes_4_requested = length(breaks4) - 1,
  classes_5_requested = length(breaks5) - 1,
  classes_6_requested = length(breaks6) - 1
)
## classes_3_requested classes_4_requested classes_5_requested classes_6_requested 
##                   3                   4                   5                   6

2.3.3 Assign counties to quantile classes

Each county is assigned to a class based on its overall tornado density. include.lowest = TRUE ensures that the lowest observed value is included in the first interval.

# assign breaks to data
countymap3 <- countymap %>%
  mutate(tdens_c1 = cut(tdens,
    breaks = quantile(tdens,
    probs = qbrks3,
    na.rm = TRUE),
    include.lowest = TRUE))

countymap4 <- countymap %>%
  mutate(tdens_c1 = cut(
    tdens,
    breaks = quantile(tdens,
    probs = qbrks4,
    na.rm = TRUE),
    include.lowest = TRUE))

countymap5 <- countymap %>%
  mutate(tdens_c1 = cut(
    tdens,
    breaks = quantile(tdens,
    probs = qbrks5,
    na.rm = TRUE),
    include.lowest = TRUE))

countymap6 <- countymap %>%
  mutate(tdens_c1 = cut(tdens, breaks = quantile(tdens, 
    probs = qbrks6, na.rm = TRUE),
    include.lowest = TRUE))

2.3.4 Create the four choropleth maps

The four maps use the same color palette but a different number of quantile classes. The legend intervals show the density values assigned to each class.

map3 <- ggplot(data = countymap3) +
  geom_sf(aes(fill = tdens_c1)) +
  scale_fill_brewer(
    name = expression("Tornadoes/1000 km"^2),
    palette = "Spectral", direction = -1) +
  theme_void() +
  theme(legend.position = "bottom")
  

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


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


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

2.3.5 Combine the four maps

The maps are arranged in a two-by-two layout using plot_grid(). This makes it easier to see how the number of classes changes the pattern and the amount of detail shown.

library(patchwork)

plot_grid(map3, map4, map5, map6,
          labels = c("A) 3 Breaks", 
                     "B) 4 Breaks",
                     "C) 5 Breaks",
                     "D) 6 Breaks",
                     label_size = 12),
          ncol = 1, 
          hjust = 0, 
          label_x = 0, 
          align = "hv") +
  plot_annotation(
    title = "Tornado Density by County",
    subtitle = "Comparison of Quantile Classification Methods")
Figure 3. County-level tornado point density across 2016–2021 classified using quantile breaks with (A) three, (B) four, (C) five, and (D) six classes.

Figure 3. County-level tornado point density across 2016–2021 classified using quantile breaks with (A) three, (B) four, (C) five, and (D) six classes.

Results Quantile maps have equal numbers of counties in each class. Increasing the number of classes from three to six reveals distinct differences between methods, but each shows the same overall patterns. Fewer classes results in less distinction between high and low values, while more classes can be difficult to interpret.

3 Overall conclusion

There are many methods of mapping the same data, each with unique results and different implications for interpretation of the data. Even the same style of map can appear vastly differnt depending on the classification groups used, as in example three. I would say that some methods are more appropriate for this type of data than others, and even the color scheme chosen can have a great impact on the results. Becoming familiar with different types of mapping techniques and data classification or normalization methods is the best way to know which mapping technique to choose for your type of data, audience, and research question.