Synopsis

This analysis examines the U.S. National Oceanic and Atmospheric Administration (NOAA) storm database to determine which types of severe weather events are most damaging to population health and which carry the greatest economic cost. The database records 902,297 events between 1950 and November 2011. Two substantial data-cleaning steps were required: the free-text EVTYPE field contains 985 distinct spellings of what are really a few dozen event types, and the monetary damage fields store a numeric mantissa and an alphabetic exponent in separate columns that must be recombined. After consolidating event names and decoding the damage exponents, tornadoes emerge as by far the largest cause of casualties, responsible for roughly 97,000 deaths and injuries, nearly ten times the next most harmful event type. Ranking by deaths alone changes the picture, with excessive heat second and far more lethal per casualty than tornadoes. Economically the answer is different again: hurricanes and typhoons cause the greatest total damage at approximately $91 billion, followed by tornadoes and storm surge. One mis-keyed record, a 2006 Napa County flood coded in billions rather than millions, is corrected during processing; left uncorrected it would inflate the flood total by $115 billion and wrongly place floods first. Crop damage follows a distinct pattern from property damage, being driven largely by drought rather than by storms.

knitr::opts_chunk$set(echo = TRUE)
library(ggplot2)
library(gridExtra)

Data Processing

Loading the raw data

The analysis begins from the raw bzip2-compressed CSV distributed with the assignment. read.csv reads the compressed file directly through a bzfile connection, so no preprocessing happens outside this document. The read is slow, so this chunk is cached.

data_file <- "StormData.csv.bz2"
if (!file.exists(data_file)) {
    download.file(
        "https://d396qusza40orc.cloudfront.net/repdata%2Fdata%2FStormData.csv.bz2",
        destfile = data_file, mode = "wb")
}

storm <- read.csv(bzfile(data_file), stringsAsFactors = FALSE)
dim(storm)
## [1] 902297     37

Only six columns are needed: the event type, the two casualty counts, and the two damage amounts with their exponent codes.

keep <- c("EVTYPE", "FATALITIES", "INJURIES",
          "PROPDMG", "PROPDMGEXP", "CROPDMG", "CROPDMGEXP")
storm <- storm[, keep]

str(storm)
## 'data.frame':    902297 obs. of  7 variables:
##  $ EVTYPE    : chr  "TORNADO" "TORNADO" "TORNADO" "TORNADO" ...
##  $ FATALITIES: num  0 0 0 0 0 0 0 0 1 0 ...
##  $ INJURIES  : num  15 0 2 2 2 6 1 0 14 0 ...
##  $ PROPDMG   : num  25 2.5 25 2.5 2.5 2.5 2.5 2.5 25 25 ...
##  $ PROPDMGEXP: chr  "K" "K" "K" "K" ...
##  $ CROPDMG   : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ CROPDMGEXP: chr  "" "" "" "" ...

Rows recording no casualties and no damage cannot affect either question, so they are dropped. This removes about 72% of the rows and makes the remaining steps much faster.

has_impact <- storm$FATALITIES > 0 | storm$INJURIES > 0 |
              storm$PROPDMG > 0 | storm$CROPDMG > 0
storm <- storm[has_impact, ]

nrow(storm)
## [1] 254633

Decoding the damage amounts

Property and crop damage are each stored across two columns: a numeric mantissa (PROPDMG, CROPDMG) and a character exponent (PROPDMGEXP, CROPDMGEXP). The exponent column must be decoded before the two can be multiplied into a dollar figure.

table(storm$PROPDMGEXP)
## 
##             -      +      0      2      3      4      5      6      7      B 
##  11585      1      5    210      1      1      4     18      3      3     40 
##      h      H      K      m      M 
##      1      6 231428      7  11320
table(storm$CROPDMGEXP)
## 
##             ?      0      B      k      K      m      M 
## 152664      6     17      7     21  99932      1   1985

Per the NWS documentation, H, K, M and B denote hundreds, thousands, millions and billions, in either case. A handful of records carry a numeric exponent, which is interpreted as a power of ten. A small number carry -, +, ? or an empty string; these are undocumented, so they are treated as a multiplier of 1, which is the conservative choice because it cannot inflate any event type’s total.

decode_exp <- function(x) {
    x <- toupper(trimws(x))
    mult <- rep(1, length(x))
    mult[x == "H"] <- 1e2
    mult[x == "K"] <- 1e3
    mult[x == "M"] <- 1e6
    mult[x == "B"] <- 1e9
    digits <- grepl("^[0-9]$", x)
    mult[digits] <- 10^as.numeric(x[digits])
    mult
}

storm$prop_cost <- storm$PROPDMG * decode_exp(storm$PROPDMGEXP)
storm$crop_cost <- storm$CROPDMG * decode_exp(storm$CROPDMGEXP)

One record is a well-known data-entry error: a 2006 Napa County flood is coded as B (billions) where the NWS narrative for the event describes damage of roughly $115 million. Left uncorrected, this single row adds about $115 billion to the flood total and would single-handedly determine the answer to the economic question. Because the correction rests on the source narrative rather than on the recorded data, it is applied explicitly and visibly rather than silently.

bad <- which(storm$prop_cost == max(storm$prop_cost))
storm[bad, c("EVTYPE", "PROPDMG", "PROPDMGEXP", "prop_cost")]
##        EVTYPE PROPDMG PROPDMGEXP prop_cost
## 605953  FLOOD     115          B  1.15e+11
storm$prop_cost[bad] <- storm$PROPDMG[bad] * 1e6
storm$total_cost <- storm$prop_cost + storm$crop_cost

Consolidating event types

EVTYPE is free text and was never validated on entry, so the same event type appears under many spellings.

length(unique(storm$EVTYPE))
## [1] 488
head(grep("THUNDERSTORM|TSTM", unique(storm$EVTYPE), value = TRUE), 10)
##  [1] "TSTM WIND"                      "THUNDERSTORM WINDS"            
##  [3] "THUNDERSTORM WIND"              "THUNDERSTORM WINS"             
##  [5] "THUNDERSTORM WINDS LIGHTNING"   "THUNDERSTORM WINDS/HAIL"       
##  [7] "THUNDERSTORM WINDS HAIL"        "FLASH FLOODING/THUNDERSTORM WI"
##  [9] "THUNDERSTORM WINDS/FUNNEL CLOU" "SEVERE THUNDERSTORM"

Leaving these separate would badly distort the results. Thunderstorm wind alone is split across TSTM WIND, THUNDERSTORM WIND, THUNDERSTORM WINDS, TSTMW and dozens of other variants, so its true ranking is understated in the raw data. Similarly, heat is divided between HEAT, EXCESSIVE HEAT and HEAT WAVE.

The labels are therefore mapped onto the standard categories used by the National Weather Service. Matching is done with ordered regular expressions on the upper-cased label; because the rules are applied in sequence, more specific patterns are placed before more general ones (for example, FLASH FLOOD is matched before the broader FLOOD rule, and MARINE THUNDERSTORM before THUNDERSTORM). Labels matching none of the rules retain their cleaned original text, so nothing is silently discarded.

clean_event <- function(x) {
    e <- toupper(trimws(x))
    rules <- list(
        c("TORNADO|TORNDAO|WATERSPOUT",          "TORNADO"),
        c("EXCESSIVE HEAT|EXTREME HEAT|RECORD HEAT|HEAT WAVE", "EXCESSIVE HEAT"),
        c("^HEAT|UNSEASONABLY WARM|WARM WEATHER", "HEAT"),
        c("FLASH FLOOD",                          "FLASH FLOOD"),
        c("COASTAL FLOOD|TIDAL FLOOD|BEACH EROSION", "COASTAL FLOOD"),
        c("FLOOD|FLD|DAM BREAK",                  "FLOOD"),
        c("LIGHTNING|LIGNTNING|LIGHTING",         "LIGHTNING"),
        c("MARINE TSTM|MARINE THUNDERSTORM",      "MARINE THUNDERSTORM WIND"),
        c("TSTM|THUNDERSTORM|THUDERSTORM|THUNERSTORM", "THUNDERSTORM WIND"),
        c("HAIL",                                 "HAIL"),
        c("HURRICANE|TYPHOON",                    "HURRICANE/TYPHOON"),
        c("TROPICAL STORM",                       "TROPICAL STORM"),
        c("BLIZZARD",                             "BLIZZARD"),
        c("ICE STORM|GLAZE|FREEZING RAIN|ICY",    "ICE STORM"),
        c("WINTER STORM",                         "WINTER STORM"),
        c("HEAVY SNOW|EXCESSIVE SNOW|SNOW STORM", "HEAVY SNOW"),
        c("WINTER WEATHER|WINTRY MIX|LIGHT SNOW|^SNOW", "WINTER WEATHER"),
        c("AVALANCH",                             "AVALANCHE"),
        c("EXTREME COLD|EXTREME WIND ?CHILL|RECORD COLD", "EXTREME COLD/WIND CHILL"),
        c("COLD|LOW TEMPERATURE|HYPOTHERMIA|FREEZE|FROST", "COLD/WIND CHILL"),
        c("DROUGHT|DRY|DRIEST",                   "DROUGHT"),
        c("WILD.*FIRE|FOREST FIRE|GRASS FIRE|BRUSH FIRE", "WILDFIRE"),
        c("RIP CURRENT",                          "RIP CURRENT"),
        c("HIGH SURF|HEAVY SURF|ROUGH SURF|HIGH SEAS|HIGH SWELLS", "HIGH SURF"),
        c("STORM SURGE|STORM TIDE",               "STORM SURGE/TIDE"),
        c("HIGH WIND",                            "HIGH WIND"),
        c("STRONG WIND",                          "STRONG WIND"),
        c("MARINE",                               "MARINE HIGH WIND"),
        c("HEAVY RAIN|EXCESSIVE RAIN|RAINSTORM|HEAVY PRECIP", "HEAVY RAIN"),
        c("FOG|VOG",                              "DENSE FOG"),
        c("DUST STORM|BLOWING DUST",              "DUST STORM"),
        c("DUST DEVIL",                           "DUST DEVIL"),
        c("LANDSLIDE|MUDSLIDE|MUD SLIDE|ROCK SLIDE|LANDSLUMP", "DEBRIS FLOW"),
        c("TSUNAMI",                              "TSUNAMI"),
        c("SEICHE",                               "SEICHE"),
        c("VOLCANIC",                             "VOLCANIC ASH"),
        c("FUNNEL",                               "FUNNEL CLOUD"),
        c("SLEET",                                "SLEET"),
        c("WIND",                                 "HIGH WIND")
    )
    out <- e
    matched <- rep(FALSE, length(e))
    for (r in rules) {
        hit <- !matched & grepl(r[1], e)
        out[hit] <- r[2]
        matched <- matched | hit
    }
    out
}

storm$event <- clean_event(storm$EVTYPE)
length(unique(storm$event))
## [1] 100

The consolidation reduces 488 distinct labels to a far smaller set, and the vast majority of records now fall into a recognised category.

round(100 * sum(storm$event %in% c(
    "TORNADO","EXCESSIVE HEAT","HEAT","FLASH FLOOD","FLOOD","COASTAL FLOOD",
    "LIGHTNING","THUNDERSTORM WIND","HAIL","HURRICANE/TYPHOON","TROPICAL STORM",
    "BLIZZARD","ICE STORM","WINTER STORM","HEAVY SNOW","WINTER WEATHER",
    "AVALANCHE","EXTREME COLD/WIND CHILL","COLD/WIND CHILL","DROUGHT","WILDFIRE",
    "RIP CURRENT","HIGH SURF","STORM SURGE/TIDE","HIGH WIND","STRONG WIND",
    "HEAVY RAIN","DENSE FOG","DUST STORM","DUST DEVIL","DEBRIS FLOW","TSUNAMI",
    "SEICHE","VOLCANIC ASH","FUNNEL CLOUD","SLEET","MARINE THUNDERSTORM WIND",
    "MARINE HIGH WIND")) / nrow(storm), 2)
## [1] 99.83

Aggregating by event type

health <- aggregate(cbind(FATALITIES, INJURIES) ~ event, data = storm, sum)
health$casualties <- health$FATALITIES + health$INJURIES
health <- health[order(-health$casualties), ]

economy <- aggregate(cbind(prop_cost, crop_cost, total_cost) ~ event,
                     data = storm, sum)
economy <- economy[order(-economy$total_cost), ]

Results

Which event types are most harmful to population health?

head(health[, c("event", "FATALITIES", "INJURIES", "casualties")], 10)
##                event FATALITIES INJURIES casualties
## 87           TORNADO       5664    91436      97100
## 86 THUNDERSTORM WIND        710     9509      10219
## 25    EXCESSIVE HEAT       2201     7124       9325
## 30             FLOOD        512     6873       7385
## 65         LIGHTNING        817     5232       6049
## 37              HEAT        977     2119       3096
## 29       FLASH FLOOD       1035     1802       2837
## 60         ICE STORM        109     2260       2369
## 50         HIGH WIND        327     1573       1900
## 98          WILDFIRE         90     1608       1698

Tornadoes are overwhelmingly the largest threat to population health, causing 97,100 total casualties — roughly 9.5 times as many as THUNDERSTORM WIND, the next event type. That total is dominated by injuries rather than deaths, which reflects the nature of tornado damage: large numbers of people are hurt by flying debris and structural collapse, while fatalities are comparatively rarer.

Note that thunderstorm wind reaches second place only because of the name consolidation described above; in the raw data its casualties are split across TSTM WIND, THUNDERSTORM WIND and THUNDERSTORM WINDS, none of which would rank second on its own.

Ranking by deaths alone tells a somewhat different story, which matters for anyone prioritising preparedness resources.

deadliest <- health[order(-health$FATALITIES), c("event", "FATALITIES", "INJURIES")]
head(deadliest, 10)
##                      event FATALITIES INJURIES
## 87                 TORNADO       5664    91436
## 25          EXCESSIVE HEAT       2201     7124
## 29             FLASH FLOOD       1035     1802
## 37                    HEAT        977     2119
## 65               LIGHTNING        817     5232
## 86       THUNDERSTORM WIND        710     9509
## 77             RIP CURRENT        577      529
## 30                   FLOOD        512     6873
## 50               HIGH WIND        327     1573
## 27 EXTREME COLD/WIND CHILL        305      260
top_cas <- head(health, 10)
top_cas$event <- factor(top_cas$event, levels = rev(top_cas$event))

top_fat <- head(deadliest, 10)
top_fat$event <- factor(top_fat$event, levels = rev(top_fat$event))

p1 <- ggplot(top_cas, aes(event, casualties)) +
    geom_col(fill = "firebrick") +
    coord_flip() +
    labs(title = "Total casualties (deaths + injuries)",
         x = NULL, y = "People affected") +
    theme_bw()

p2 <- ggplot(top_fat, aes(event, FATALITIES)) +
    geom_col(fill = "black") +
    coord_flip() +
    labs(title = "Deaths only", x = NULL, y = "Deaths") +
    theme_bw()

grid.arrange(p1, p2, nrow = 1,
    top = "Figure 1: Weather event types most harmful to population health, 1950-2011.
Left: total casualties, dominated by tornado injuries. Right: deaths alone, where
excessive heat ranks close behind tornadoes.")

Figure 1 shows the ten most harmful event types by total casualties (left) and by deaths alone (right). Tornadoes lead both rankings, but the two panels diverge sharply below the top position. Excessive heat causes comparatively few injuries yet is the second deadliest event type, making it far more lethal per casualty than tornadoes. Flash floods, lightning and rip currents also rank higher on deaths than their total casualty counts would suggest — rip currents in particular cause more deaths than injuries, which is unusual among weather events.

Which event types have the greatest economic consequences?

top_econ <- head(economy, 10)
data.frame(
    event = top_econ$event,
    property_billions = round(top_econ$prop_cost / 1e9, 2),
    crop_billions = round(top_econ$crop_cost / 1e9, 2),
    total_billions = round(top_econ$total_cost / 1e9, 2))
##                event property_billions crop_billions total_billions
## 1  HURRICANE/TYPHOON             85.36          5.52          90.87
## 2            TORNADO             58.61          0.42          59.03
## 3   STORM SURGE/TIDE             47.96          0.00          47.97
## 4              FLOOD             35.38         10.86          46.24
## 5        FLASH FLOOD             17.59          1.53          19.12
## 6               HAIL             15.98          3.05          19.02
## 7            DROUGHT              1.05         13.97          15.03
## 8  THUNDERSTORM WIND             11.18          1.27          12.45
## 9          ICE STORM              3.96          5.02           8.99
## 10          WILDFIRE              8.50          0.40           8.90

HURRICANE/TYPHOON events cause the greatest total economic damage, at approximately $91 billion, followed by tornado and storm surge/tide. This ordering depends directly on the Napa County correction applied during processing: without it, floods would carry an additional $115 billion from that single mis-keyed record and would rank first instead of fourth. A peer reproducing this analysis without the correction should therefore expect floods at the top.

Property and crop damage behave quite differently, so it is worth separating them.

prop_top <- head(economy[order(-economy$prop_cost), ], 10)
prop_top$event <- factor(prop_top$event, levels = rev(prop_top$event))

crop_top <- head(economy[order(-economy$crop_cost), ], 10)
crop_top$event <- factor(crop_top$event, levels = rev(crop_top$event))

p3 <- ggplot(prop_top, aes(event, prop_cost / 1e9)) +
    geom_col(fill = "steelblue") +
    coord_flip() +
    labs(title = "Property damage", x = NULL, y = "Billions of USD") +
    theme_bw()

p4 <- ggplot(crop_top, aes(event, crop_cost / 1e9)) +
    geom_col(fill = "darkgreen") +
    coord_flip() +
    labs(title = "Crop damage", x = NULL, y = "Billions of USD") +
    theme_bw()

grid.arrange(p3, p4, nrow = 1,
    top = "Figure 2: Weather event types with the greatest economic consequences, 1950-2011.
Left: property damage, led by hurricanes and tornadoes. Right: crop damage, led by
drought - an event type that causes almost no property damage.")

Figure 2 separates property damage (left) from crop damage (right), and the two panels rank almost entirely different events. Property damage is driven by sudden, violent events — hurricanes, tornadoes, storm surge and floods — that destroy buildings and infrastructure. Crop damage is led by drought, a slow-onset event that barely registers in the property panel at all, followed by floods and ice storms. Property damage is also roughly six times larger than crop damage in total, which is why the combined ranking closely tracks the property panel.

c(property_total_billions = round(sum(storm$prop_cost) / 1e9),
  crop_total_billions = round(sum(storm$crop_cost) / 1e9))
## property_total_billions     crop_total_billions 
##                     313                      49

Summary

For a manager allocating preparedness resources, the two questions point in different directions. Tornadoes are the clear priority for protecting people, with excessive heat second for preventing deaths specifically. Hurricanes are the clear priority for limiting economic loss, with tornadoes and storm surge following. Only tornadoes rank near the top of both, which makes them the single event type that warrants investment on both grounds.