Report for Reproducible Research Course Project 2 Assignment

This report is based on the U.S. NOAA Storm Data and contains an analysis to determine which types of severe weather events are most harmful to population health and which have the greatest economic consequences. The raw data were reduced to the variables relevant for the analysis, and property and crop damage were decoded from their coded exponents (K/M/B), with the few undocumented codes impact shown to be negligible. The ~985 raw event labels were consolidated toward the 48 official NWS event categories. The timespan of the analysis was restricted to 1996-2011, the period in which NOAA recorded all event types. Additionally, a single erroneous $115 billion flood record (the 2006 Napa entry) was corrected to $115 million. With respect to population health, tornadoes cause by far the most injuries (20,670) while heat causes the most fatalities (2,034), with heat deaths highly concentrated — the top five states account for 62.1% of the national total. With respect to the greatest economic consequences, hurricanes (typhoons) cause the greatest total damage ($87.1 billion), followed by storm surge/tide and flood, with drought notable as the only leading event driven primarily by crop rather than property loss.

The report is addressing the following questions:

  1. Across the United States, which types of events (as indicated in the EVTYPE) are most harmful with respect to population health?
  2. Across the United States, which types of events have the greatest economic consequences?

Data Processing

As a first step of data processing the presence of the input data file in the working directory is checked, if not present data set is downloaded, then it is read into stormData, the code chunk is run with cache=TRUE for optimizing the performance.

# check if source data has already been downloaded
if (!file.exists("FStormData.csv.bz2")) {
    download.file("https://d396qusza40orc.cloudfront.net/repdata%2Fdata%2FStormData.csv.bz2",
        "FStormData.csv.bz2", mode = "wb")
}

# load in the data optimized - please note cache=TRUE in the header
# of this code chunk
stormData <- read.csv("FStormData.csv.bz2")

str(stormData)
## 'data.frame':    902297 obs. of  37 variables:
##  $ STATE__   : num  1 1 1 1 1 1 1 1 1 1 ...
##  $ BGN_DATE  : chr  "4/18/1950 0:00:00" "4/18/1950 0:00:00" "2/20/1951 0:00:00" "6/8/1951 0:00:00" ...
##  $ BGN_TIME  : chr  "0130" "0145" "1600" "0900" ...
##  $ TIME_ZONE : chr  "CST" "CST" "CST" "CST" ...
##  $ COUNTY    : num  97 3 57 89 43 77 9 123 125 57 ...
##  $ COUNTYNAME: chr  "MOBILE" "BALDWIN" "FAYETTE" "MADISON" ...
##  $ STATE     : chr  "AL" "AL" "AL" "AL" ...
##  $ EVTYPE    : chr  "TORNADO" "TORNADO" "TORNADO" "TORNADO" ...
##  $ BGN_RANGE : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ BGN_AZI   : chr  "" "" "" "" ...
##  $ BGN_LOCATI: chr  "" "" "" "" ...
##  $ END_DATE  : chr  "" "" "" "" ...
##  $ END_TIME  : chr  "" "" "" "" ...
##  $ COUNTY_END: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ COUNTYENDN: logi  NA NA NA NA NA NA ...
##  $ END_RANGE : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ END_AZI   : chr  "" "" "" "" ...
##  $ END_LOCATI: chr  "" "" "" "" ...
##  $ LENGTH    : num  14 2 0.1 0 0 1.5 1.5 0 3.3 2.3 ...
##  $ WIDTH     : num  100 150 123 100 150 177 33 33 100 100 ...
##  $ F         : int  3 2 2 2 2 2 2 1 3 3 ...
##  $ MAG       : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ 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  "" "" "" "" ...
##  $ WFO       : chr  "" "" "" "" ...
##  $ STATEOFFIC: chr  "" "" "" "" ...
##  $ ZONENAMES : chr  "" "" "" "" ...
##  $ LATITUDE  : num  3040 3042 3340 3458 3412 ...
##  $ LONGITUDE : num  8812 8755 8742 8626 8642 ...
##  $ LATITUDE_E: num  3051 0 0 0 0 ...
##  $ LONGITUDE_: num  8806 0 0 0 0 ...
##  $ REMARKS   : chr  "" "" "" "" ...
##  $ REFNUM    : num  1 2 3 4 5 6 7 8 9 10 ...

As a next step the number of columns is reduced to the ones relevant for this analysis, 10 of the original 37 columns are kept.

library(data.table)
# the column names that are kept for further analysis are stored in
# `keep`
keep <- c("EVTYPE", "FATALITIES", "INJURIES", "PROPDMG", "PROPDMGEXP",
    "CROPDMG", "CROPDMGEXP", "BGN_DATE", "STATE", "REFNUM")
stormData <- as.data.table(stormData)[, ..keep]

# the structure of the reduced data frame is shown
str(stormData)
## Classes 'data.table' and 'data.frame':   902297 obs. of  10 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  "" "" "" "" ...
##  $ BGN_DATE  : chr  "4/18/1950 0:00:00" "4/18/1950 0:00:00" "2/20/1951 0:00:00" "6/8/1951 0:00:00" ...
##  $ STATE     : chr  "AL" "AL" "AL" "AL" ...
##  $ REFNUM    : num  1 2 3 4 5 6 7 8 9 10 ...
##  - attr(*, ".internal.selfref")=<pointer: 0x101e90db0>

Expenses are represented with coded exponents in separate columns. In the next step these are decoded and dropped once they are not needed for further processing. The impact of the undocumented codes is reported.

In addition to these steps a single entry correction is also performed due to the error’s significant impact on the economic impact summary. One record (REFNUM 605943, a January 2006 Napa Valley flood) was coded with a B (billion) property-damage exponent, yielding an implausible $115 billion impact — more than Hurricane Katrina’s total.

The event narrative describes a localized flood, and the value is a well-known data-entry error; the exponent was corrected to M (million, ~$115M). This single fix reduces the FLOOD category total from ~$149B to ~$34B and changes the top economic event type from FLOOD to HURRICANE (TYPHOON).

# Decode, drop and then show the impact of handling unocumented
# exponents NWS 10-1605 defines only K/M/B (H is a common
# convention). Any other PROPDMGEXP/CROPDMGEXP value is undocumented.
# Report how much money is impacted by using the multiplier of 1.

documented <- c("", "H", "K", "M", "B", as.character(0:9))

# Decode the damage exponents into real dollar amounts ---
exp_to_mult <- function(e) {
    e <- toupper(trimws(as.character(e)))
    data.table::fcase(e == "H", 100, e == "K", 1000, e == "M", 1e+06, e ==
        "B", 1e+09, e %in% as.character(0:9), 10^suppressWarnings(as.numeric(e)),
        default = 1  # '', +, -, ?, NA, anything else
)
}

# --- Correct the known Napa 2006 data-entry error (REFNUM 605943)
# --- PROPDMGEXP was recorded as 'B' (billion) but the event was a
# localized Napa Valley flood; the remarks and cross-checks indicate
# ~$115M. This single record otherwise contributes $115B, ~77% of the
# FLOOD type total.
stormData[REFNUM == 605943 & PROPDMGEXP == "B", PROPDMGEXP := "M"]

stormData[, propCost := PROPDMG * exp_to_mult(PROPDMGEXP)]
stormData[, cropCost := CROPDMG * exp_to_mult(CROPDMGEXP)]

undoc_prop_cost <- stormData[!toupper(trimws(PROPDMGEXP)) %in% documented,
    sum(propCost)]
undoc_crop_cost <- stormData[!toupper(trimws(CROPDMGEXP)) %in% documented,
    sum(cropCost)]

# Drop the now-redundant raw damage columns ---
stormData[, c("PROPDMG", "PROPDMGEXP", "CROPDMG", "CROPDMGEXP") := NULL]

# VERIFY (true-dollar scale): express the undocumented-code damage as
# a fraction of the actual decoded dollar totals Digit codes 0-9 are
# treated as documented powers of ten (10^digit); only the genuinely
# undefined symbols (+, -, ?) remain 'undocumented' and are mapped to
# x1. Their combined impact ($132 property / $0 crop) is a vanishing
# fraction of the decoded totals -> no effect on rankings.

message(sprintf("Undocumented PROPDMG: $%.0f of $%.0f decoded total (%.8f%%)",
    undoc_prop_cost, stormData[, sum(propCost)], 100 * undoc_prop_cost/stormData[,
        sum(propCost)]))
## Undocumented PROPDMG: $132 of $313339868768 decoded total (0.00000004%)
message(sprintf("Undocumented CROPDMG: $%.0f of $%.0f decoded total (%.8f%%)",
    undoc_crop_cost, stormData[, sum(cropCost)], 100 * undoc_crop_cost/stormData[,
        sum(cropCost)]))
## Undocumented CROPDMG: $0 of $49104192181 decoded total (0.00000000%)

With exponent decoding done the next step is to consolidate EVTYPE values, then check its effectiveness and report residual impact.

# ============================================================
# Consolidate EVTYPE, then check the consolidation's effectiveness
# ============================================================

# Official NWS Storm Data event types (NWS Instruction 10-1605, Aug
# 2007, Table 1 'Storm Data Event Table') - the 48 permitted event
# names, uppercased to match the cleaning pipeline.
official_events <- c("ASTRONOMICAL LOW TIDE", "AVALANCHE", "BLIZZARD",
    "COASTAL FLOOD", "COLD/WIND CHILL", "DEBRIS FLOW", "DENSE FOG", "DENSE SMOKE",
    "DROUGHT", "DUST DEVIL", "DUST STORM", "EXCESSIVE HEAT", "EXTREME COLD/WIND CHILL",
    "FLASH FLOOD", "FLOOD", "FROST/FREEZE", "FUNNEL CLOUD", "FREEZING FOG",
    "HAIL", "HEAT", "HEAVY RAIN", "HEAVY SNOW", "HIGH SURF", "HIGH WIND",
    "HURRICANE (TYPHOON)", "ICE STORM", "LAKE-EFFECT SNOW", "LAKESHORE FLOOD",
    "LIGHTNING", "MARINE HAIL", "MARINE HIGH WIND", "MARINE STRONG WIND",
    "MARINE THUNDERSTORM WIND", "RIP CURRENT", "SEICHE", "SLEET", "STORM SURGE/TIDE",
    "STRONG WIND", "THUNDERSTORM WIND", "TORNADO", "TROPICAL DEPRESSION",
    "TROPICAL STORM", "TSUNAMI", "VOLCANIC ASH", "WATERSPOUT", "WILDFIRE",
    "WINTER STORM", "WINTER WEATHER")

# --- Baseline totals BEFORE consolidation (for conservation check)
# --- Consolidation only relabels EVTYPE; these sums must be
# unchanged after.
base_totals <- stormData[, .(fat = sum(FATALITIES), inj = sum(INJURIES),
    prop = sum(propCost), crop = sum(cropCost))]
base_distinct <- stormData[, uniqueN(EVTYPE)]

# --- Normalize and consolidate EVTYPE --- ~985 raw values (typos,
# casing, near-duplicates) folded toward the official NWS categories.
# RULES APPLY SEQUENTIALLY and overwrite, so specific patterns must
# precede the general ones that would swallow them (e.g. FLASH FLOOD
# before FLOOD; FREEZING FOG before FOG; RAIN after FLOOD).  NOTE:
# some rules deliberately merge DISTINCT official categories for a
# coarser harm/damage analysis - these merges are flagged inline.
stormData[, EVTYPE := toupper(trimws(EVTYPE))]
stormData[, EVTYPE := gsub("\\s+", " ", EVTYPE)]  # collapse internal spaces

# --- Thunderstorm / tornado family ---
stormData[grepl("TSTM|THUNDERSTORM", EVTYPE), EVTYPE := "THUNDERSTORM WIND"]
stormData[grepl("DRY MICROBURST|MICROBURST|DOWNBURST", EVTYPE), EVTYPE :=
    "THUNDERSTORM WIND"]
stormData[grepl("TORNADO|FUNNEL|WATERSPOUT", EVTYPE), EVTYPE := "TORNADO"]
stormData[grepl("THUNDERTORM|TUNDERSTORM|THUNERSTORM|THUNDEERSTORM", EVTYPE),
    EVTYPE := "THUNDERSTORM WIND"]

# --- Tropical ---
stormData[grepl("HURRICANE|TYPHOON", EVTYPE), EVTYPE := "HURRICANE (TYPHOON)"]
stormData[grepl("TROPICAL STORM", EVTYPE), EVTYPE := "TROPICAL STORM"]
stormData[grepl("TROPICAL DEPRESSION", EVTYPE), EVTYPE := "TROPICAL DEPRESSION"]

# --- Heat ---
stormData[grepl("EXCESSIVE HEAT|EXTREME HEAT|^HEAT|WARM|HOT", EVTYPE),
    EVTYPE := "HEAT"]


# --- Flood family (specific BEFORE general) ---
stormData[grepl("FLASH FLOOD", EVTYPE), EVTYPE := "FLASH FLOOD"]
stormData[grepl("COASTAL FLOOD|COASTAL SURGE|CSTL FLOOD|EROSION/CSTL",
    EVTYPE), EVTYPE := "COASTAL FLOOD"]
stormData[grepl("LAKESHORE FLOOD", EVTYPE), EVTYPE := "LAKESHORE FLOOD"]
stormData[grepl("^FLOOD|RIVER FLOOD|URBAN|MAJOR FLOOD|FLOODING", EVTYPE),
    EVTYPE := "FLOOD"]

# --- Surge / tide ---
stormData[grepl("SURGE|TIDE", EVTYPE), EVTYPE := "STORM SURGE/TIDE"]

# --- Cold / winter family (WINDCHILL one-word variant included) ---
stormData[grepl("COLD|WIND ?CHILL|WINDCHILL|FREEZE|FROST", EVTYPE), EVTYPE :=
    "COLD/WIND CHILL"]
stormData[grepl("BLIZZARD", EVTYPE), EVTYPE := "BLIZZARD"]
stormData[grepl("ICE STORM|^ICE|ICY|GLAZE", EVTYPE), EVTYPE := "ICE STORM"]
stormData[grepl("WINTER STORM", EVTYPE), EVTYPE := "WINTER STORM"]
stormData[grepl("WINTER WEATHER|WINTRY MIX|WINTER MIX", EVTYPE), EVTYPE :=
    "WINTER WEATHER"]
stormData[grepl("SLEET|FREEZING RAIN|FREEZING DRIZZLE", EVTYPE), EVTYPE :=
    "SLEET"]
stormData[grepl("SNOW", EVTYPE), EVTYPE := "HEAVY SNOW"]

# --- Fog (FREEZING FOG protected BEFORE generic FOG) ---
stormData[grepl("FREEZING FOG", EVTYPE), EVTYPE := "FREEZING FOG"]
stormData[grepl("FOG", EVTYPE), EVTYPE := "DENSE FOG"]

# --- Wind (non-thunderstorm) ---
stormData[grepl("HIGH WIND|STRONG WIND|GUSTY WIND|STORM FORCE WIND|^WIND",
    EVTYPE) & !grepl("THUNDERSTORM|MARINE", EVTYPE), EVTYPE := "HIGH WIND"]

# --- Surf ---
stormData[grepl("SURF", EVTYPE), EVTYPE := "HIGH SURF"]

# --- Fire, hail, lightning, drought ---
stormData[grepl("WILDFIRE|WILD.*FIRE|FOREST FIRE|GRASS FIRE|BRUSH FIRE",
    EVTYPE), EVTYPE := "WILDFIRE"]
stormData[grepl("HAIL", EVTYPE), EVTYPE := "HAIL"]
stormData[grepl("LIGHTNING|LIGHTING|LIGNTNING", EVTYPE), EVTYPE := "LIGHTNING"]
stormData[grepl("DROUGHT", EVTYPE), EVTYPE := "DROUGHT"]

# --- Landslide (officially renamed to DEBRIS FLOW in 10-1605) ---
stormData[grepl("LANDSLIDE|LANDSLUMP|MUDSLIDE|MUD SLIDE|DEBRIS FLOW", EVTYPE),
    EVTYPE := "DEBRIS FLOW"]

# --- Rip current ---
stormData[grepl("RIP CURRENT", EVTYPE), EVTYPE := "RIP CURRENT"]

# --- Rain (MUST come after flood rules so flood hybrids aren't
# stolen) ---
stormData[grepl("RAIN|WET|PRECIP|SHOWER", EVTYPE), EVTYPE := "HEAVY RAIN"]
# extend as needed for any category that still reaches the top of the
# rankings


# ============================================================ CHECK
# 1 - TOTALS CONSERVATION (relabeling must not change sums)
# ============================================================
after_totals <- stormData[, .(fat = sum(FATALITIES), inj = sum(INJURIES),
    prop = sum(propCost), crop = sum(cropCost))]
conserved <- isTRUE(all.equal(base_totals, after_totals))
message(sprintf("Totals conserved after consolidation: %s", conserved))
## Totals conserved after consolidation: TRUE
if (!conserved) {
    message("  WARNING: sums changed - a rule dropped or altered rows.")
    print(base_totals)
    print(after_totals)
}

# ============================================================ CHECK
# 2 - CARDINALITY REDUCTION (did we actually collapse labels?)
# ============================================================
message(sprintf("Distinct EVTYPE: %d -> %d", base_distinct, stormData[,
    uniqueN(EVTYPE)]))
## Distinct EVTYPE: 985 -> 227
# ============================================================ CHECK
# 3 - COVERAGE vs OFFICIAL TAXONOMY How much of the signal falls
# outside the official 48 event names?
# ============================================================
stormData[, is_official := EVTYPE %in% official_events]

tot <- stormData[, .(d = sum(FATALITIES), i = sum(INJURIES), c = sum(propCost +
    cropCost))]
uncap <- stormData[is_official == FALSE, .(d = sum(FATALITIES), i = sum(INJURIES),
    c = sum(propCost + cropCost))]
message(sprintf("Outside official taxonomy: %.2f%% deaths, %.2f%% injuries, %.2f%% dollars",
    100 * uncap$d/tot$d, 100 * uncap$i/tot$i, 100 * uncap$c/tot$c))
## Outside official taxonomy: 0.38% deaths, 0.08% injuries, 0.00% dollars
# ============================================================ CHECK
# 4 - RESIDUAL TAIL (largest non-official leftovers to fix next)
# ============================================================
print(stormData[is_official == FALSE, .(n = .N, deaths = sum(FATALITIES),
    inj = sum(INJURIES), cost = sum(propCost + cropCost)), by = EVTYPE][order(-cost)][1:15])
##                     EVTYPE     n deaths   inj    cost
##                     <char> <int>  <num> <num>   <num>
##  1:              HEAVY MIX     8      0     0 1305000
##  2: RIVER AND STREAM FLOOD     2      0     0 1200000
##  3:                  OTHER    52      0     4 1089900
##  4:              DAM BREAK     4      0     0 1002000
##  5:        COASTAL EROSION     1      0     0  766000
##  6:             HIGH WATER     6      3     0  505000
##  7:             ROCK SLIDE     2      0     0  150000
##  8:               GUSTNADO     6      0     0  103600
##  9:          BEACH EROSION     4      0     0  100000
## 10:      THUDERSTORM WINDS     2      0     0   55000
## 11:      SEVERE TURBULENCE     1      0     0   50000
## 12:        MARINE ACCIDENT     1      1     2   50000
## 13:          COASTAL STORM    10      3     2   50000
## 14:          GRADIENT WIND     9      0     0   37000
## 15:             LAKE FLOOD     1      0     0   30000
stormData[, is_official := NULL]  # tidy up the helper column

# --- Stopping rationale: consolidation is sufficient ---
message(sprintf(paste("Consolidation halted: %.2f%% of deaths, %.2f%% of injuries and %.2f%%",
    "of dollars remain \noutside the 48 official categories. The residual",
    "labels are rare, low-impact records \nthat cannot alter the top-ranked",
    "event types, so further mapping would add complexity \nwithout changing",
    "any result."), 100 * uncap$d/tot$d, 100 * uncap$i/tot$i, 100 * uncap$c/tot$c))
## Consolidation halted: 0.38% of deaths, 0.08% of injuries and 0.00% of dollars remain 
## outside the 48 official categories. The residual labels are rare, low-impact records 
## that cannot alter the top-ranked event types, so further mapping would add complexity 
## without changing any result.

As a final data-processing step the analysis is restricted to events from 1996 onward. According to NOAA’s National Centers for Environmental Information (NCEI), the Storm Events Database was created in 1996/1997; earlier records were limited to tornadoes (from 1950) and thunderstorm wind and hail (from 1955), so the pre-1996 data cannot support a fair comparison across event types. NCEI notes that “due to changes in the data collection and processing procedures over time, there are unique periods of record available depending on the event type” (NOAA NCEI Storm Events Database).

# ============================================================
# Restrict to 1996 onward
# ============================================================ NOAA
# recorded only a limited set of event types (mainly tornado, some
# thunderstorm/hail) before 1996; all 48 official event types have
# been recorded from 1996 on. Filtering to 1996+ removes the
# historical bias toward tornado and gives a fair cross-event
# comparison.
stormData[, BGN_DATE := as.Date(BGN_DATE, format = "%m/%d/%Y %H:%M:%S")]
stormData[, year := as.integer(format(BGN_DATE, "%Y"))]

before_n <- nrow(stormData)
stormData <- stormData[year >= 1996]
after_n <- nrow(stormData)

message(sprintf("Restricted to 1996+: %d -> %d rows (%.1f%% retained)",
    before_n, after_n, 100 * after_n/before_n))
## Restricted to 1996+: 902297 -> 653530 rows (72.4% retained)

Results

Storm Events with the Greatest Population Health Impact

The question to be answered by this analysis is: Across the United States, which types of events (as indicated in the EVTYPE variable) are most harmful with respect to population health?

# ============================================================
# Population health: rank event types by fatalities and injuries
# ============================================================
library(data.table)

# Aggregate both health metrics by event type
healthByType <- stormData[, .(fatalities = sum(FATALITIES), injuries = sum(INJURIES),
    casualties = sum(FATALITIES + INJURIES)  # optional combined view
),
    by = EVTYPE]

# Top 10 by fatalities
topFatal <- healthByType[order(-fatalities)][1:10]
cat("Top 10 event types by FATALITIES:\n")
## Top 10 event types by FATALITIES:
print(topFatal[, .(EVTYPE, fatalities, injuries)])
##                EVTYPE fatalities injuries
##                <char>      <num>    <num>
##  1:              HEAT       2034     7702
##  2:           TORNADO       1513    20670
##  3:       FLASH FLOOD        887     1674
##  4:         LIGHTNING        651     4141
##  5:       RIP CURRENT        542      503
##  6:             FLOOD        444     6839
##  7: THUNDERSTORM WIND        401     5189
##  8:   COLD/WIND CHILL        390      135
##  9:         HIGH WIND        368     1478
## 10:         AVALANCHE        223      156
# Top 10 by injuries
topInjury <- healthByType[order(-injuries)][1:10]
cat("\nTop 10 event types by INJURIES:\n")
## 
## Top 10 event types by INJURIES:
print(topInjury[, .(EVTYPE, injuries, fatalities)])
##                  EVTYPE injuries fatalities
##                  <char>    <num>      <num>
##  1:             TORNADO    20670       1513
##  2:                HEAT     7702       2034
##  3:               FLOOD     6839        444
##  4:   THUNDERSTORM WIND     5189        401
##  5:           LIGHTNING     4141        651
##  6:         FLASH FLOOD     1674        887
##  7:           HIGH WIND     1478        368
##  8:            WILDFIRE     1458         87
##  9: HURRICANE (TYPHOON)     1328        125
## 10:        WINTER STORM     1292        191

A figure depicting the aspects of the impact of events on population health:

    # ============================================================
    # Figure: population health impact by event type (1996+)
    # ============================================================
    library(ggplot2)
    library(data.table)
  
    # Build a long-format table of the top 10 for each metric, tagged by metric.
    # (healthByType was created in the previous chunk.)
    topFatal10  <- healthByType[order(-fatalities)][1:10,
                      .(EVTYPE, value = fatalities, metric = "Fatalities")]
    topInjury10 <- healthByType[order(-injuries)][1:10,
                      .(EVTYPE, value = injuries,   metric = "Injuries")]
  
    healthPlotData <- rbind(topFatal10, topInjury10)
  
    # Order bars by value *within* each facet using reorder_within-style trick:
    # create a per-facet ordering key so each panel sorts independently.
    healthPlotData[, ord := paste(metric, EVTYPE, sep = "___")]
    healthPlotData[, ord := reorder(ord, value)]
  
    ggplot(healthPlotData, aes(x = ord, y = value / 1000, fill = metric)) +
      geom_col(show.legend = FALSE) +
      facet_wrap(~ metric, scales = "free") +   # each panel: own events & y-scale
      coord_flip() +
      scale_x_discrete(labels = function(x) sub(".*___", "", x)) +  # strip the key
      scale_y_continuous(labels = scales::comma) +
      labs(
        title    = "Most harmful event types to population health \n(U.S., 1996 - 2011)",
        subtitle = "Top 10 event types by total fatalities and by total injuries",
        x        = NULL,
        y        = "Total persons (thousands)",
        caption  = "Source: NOAA Storm Data. Fatalities and injuries shown separately."
      ) +
      theme_minimal(base_size = 11) +
      theme(
        strip.text  = element_text(face = "bold"),
        plot.title  = element_text(face = "bold")
      )
Figure 1. Top 10 event types by total fatalities (left) and total injuries (right), United States 1996-2011. Tornadoes cause by far the most injuries while heat causes the most fatalities.

Figure 1. Top 10 event types by total fatalities (left) and total injuries (right), United States 1996-2011. Tornadoes cause by far the most injuries while heat causes the most fatalities.

The two health metrics identify different lead threats: tornadoes cause by far the most injuries (20,670), while heat causes the most fatalities (2,034, aided by consolidating “excessive heat” into “heat”). Reporting fatalities and injuries separately avoids masking this distinction behind a single combined score.

Because heat is the leading cause of weather-related fatalities, we examine its distribution across states in more detail.

Heat-related fatalities are highly concentrated geographically: the five states with the most deaths (Illinois, Pennsylvania, Texas, Missouri, and California) account for 1,263 of 2,034 heat fatalities - 62.1% of the national total - despite representing only 10% of the 50 U.S. states.

Heat related fatalities by states:

    # ============================================================
    # Figure: top 10 vs bottom 10 states by HEAT fatalities (1996+)
    # ============================================================
    library(ggplot2)
    library(data.table)
  
    # Heat fatalities by state (HEAT includes EXCESSIVE HEAT)
    heatByState <- stormData[EVTYPE == "HEAT",
                             .(fatalities = sum(FATALITIES)), by = STATE]
  
    # Bottom 10 = lowest states that still recorded >=1 heat fatality
    # (excludes the many states tied at 0, which aren't informative)
    heatByStateNZ <- heatByState[fatalities > 0][order(-fatalities)]

    # ------------------------------------------------------------
    # Concentration: share of heat fatalities in the top 5 states
    # (top 5 = top 10% of the 50 U.S. states)
    # ------------------------------------------------------------
    heatRanked     <- heatByState[order(-fatalities)]
    totalHeatFatal <- heatRanked[, sum(fatalities)]
    top5HeatFatal  <- heatRanked[1:5, sum(fatalities)]
  
    message(sprintf(
      "Top 5 states account for %d of %d heat fatalities (%.1f%%)",
      top5HeatFatal, totalHeatFatal, 100 * top5HeatFatal / totalHeatFatal))
## Top 5 states account for 1263 of 2034 heat fatalities (62.1%)
    top10    <- heatByStateNZ[1:10]
    bottom10 <- heatByStateNZ[(.N - 9):.N]
  
    top10[,    group := "Top 10 states"]
    bottom10[, group := "Bottom 10 states (with >=1 fatality)"]
  
    plotData <- rbind(top10, bottom10)
    # fix facet order so Top appears first
    plotData[, group := factor(group,
               levels = c("Top 10 states",
                          "Bottom 10 states (with >=1 fatality)"))]
  
    # per-facet ordering (reorder_within-style trick)
    plotData[, ord := paste(group, STATE, sep = "___")]
    plotData[, ord := reorder(ord, fatalities)]
  
    print(plotData[, .(group, STATE, fatalities)])
##                                    group  STATE fatalities
##                                   <fctr> <char>      <num>
##  1:                        Top 10 states     IL        357
##  2:                        Top 10 states     PA        308
##  3:                        Top 10 states     TX        285
##  4:                        Top 10 states     MO        195
##  5:                        Top 10 states     CA        118
##  6:                        Top 10 states     NY         93
##  7:                        Top 10 states     MD         89
##  8:                        Top 10 states     OK         87
##  9:                        Top 10 states     NV         62
## 10:                        Top 10 states     AZ         58
## 11: Bottom 10 states (with >=1 fatality)     WA          3
## 12: Bottom 10 states (with >=1 fatality)     IA          2
## 13: Bottom 10 states (with >=1 fatality)     RI          2
## 14: Bottom 10 states (with >=1 fatality)     CT          1
## 15: Bottom 10 states (with >=1 fatality)     IN          1
## 16: Bottom 10 states (with >=1 fatality)     NM          1
## 17: Bottom 10 states (with >=1 fatality)     OR          1
## 18: Bottom 10 states (with >=1 fatality)     AK          1
## 19: Bottom 10 states (with >=1 fatality)     UT          1
## 20: Bottom 10 states (with >=1 fatality)     SD          1
##                                    group  STATE fatalities
##                                   <fctr> <char>      <num>
    ggplot(plotData, aes(x = ord, y = fatalities, fill = group)) +
      geom_col(show.legend = FALSE) +
      geom_text(aes(label = fatalities), hjust = -0.15, size = 3) +
      facet_wrap(~ group, scales = "free") +          # own states & scale per panel
      coord_flip() +
      scale_x_discrete(labels = function(x) sub(".*___", "", x)) +
      scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
      labs(
        title    = "U.S. states most and least affected by heat-related fatalities",
        subtitle = "Top 5 states account for 62.1% of U.S. heat fatalities (1996-2011)",
        x        = "State",
        y        = "Total fatalities",
        caption  = "Source: NOAA Storm Data. Bottom 10 = lowest states with at least one fatality."
      ) +
      theme_minimal(base_size = 11) +
      theme(
        strip.text = element_text(face = "bold"),
        plot.title = element_text(face = "bold")
      )
Figure 2. States with the most (left) and fewest non-zero (right) heat-related fatalities, 1996-2011. Heat deaths are highly concentrated: the top five states account for 62.1% of the national total.

Figure 2. States with the most (left) and fewest non-zero (right) heat-related fatalities, 1996-2011. Heat deaths are highly concentrated: the top five states account for 62.1% of the national total. 

    # States/territories with ZERO heat-related fatalities (1996+)
    allStates  <- stormData[, unique(STATE)]
    heatStates <- stormData[EVTYPE == "HEAT" & FATALITIES > 0, unique(STATE)]
    zeroHeatStates <- sort(setdiff(allStates, heatStates))
    cat(zeroHeatStates, sep = ", ")
## AM, AN, AS, CO, GM, GU, HI, ID, LC, LE, LH, LM, LO, LS, MA, ME, MT, ND, NH, PH, PK, PM, PR, PZ, SL, VI, VT, WV, WY, XX

Storm Events with the Greatest Economic Impact

The question to be answered by this analysis is: Across the United States, which types of events have the greatest economic consequences?

# ============================================================
# Economic consequences: rank event types by total damage
# ============================================================
library(data.table)

# Total economic damage = property + crop cost, summed by event type
econByType <- stormData[, .(propCost = sum(propCost), cropCost = sum(cropCost),
    total = sum(propCost + cropCost)), by = EVTYPE][order(-total)]

# Top 10 economic event types
econTop10 <- econByType[1:10]
print(econTop10)
##                  EVTYPE    propCost    cropCost       total
##                  <char>       <num>       <num>       <num>
##  1: HURRICANE (TYPHOON) 81718889010  5350107800 87068996810
##  2:    STORM SURGE/TIDE 47844469000      855000 47845324000
##  3:               FLOOD 29244603200  5013161500 34257764700
##  4:             TORNADO 24622810010   283425010 24906235020
##  5:                HAIL 14595217420  2496822450 17092039870
##  6:         FLASH FLOOD 15222258910  1334901700 16557160610
##  7:             DROUGHT  1046101000 13367566000 14413667000
##  8:   THUNDERSTORM WIND  7921242880  1017007600  8938250480
##  9:      TROPICAL STORM  7642475550   677711000  8320186550
## 10:            WILDFIRE  7760449500   402255130  8162704630

Figure:

    # ============================================================
    # Figure: economic damage by event type - property + crop stacked
    # ============================================================
    library(ggplot2)
    library(data.table)
  
    # Reshape top-10 to long: one row per (event, component)
    econLong <- melt(
      econTop10[, .(EVTYPE, Property = propCost, Crop = cropCost)],
      id.vars      = "EVTYPE",
      variable.name = "component",
      value.name    = "cost"
    )
  
    # Order events by TOTAL (so the tallest full bar is on top after coord_flip)
    order_by_total <- econTop10[order(total), EVTYPE]
    econLong[, EVTYPE := factor(EVTYPE, levels = order_by_total)]
  
    ggplot(econLong, aes(x = EVTYPE, y = cost / 1e9, fill = component)) +
      geom_col() +                                   # stacked by default
      coord_flip() +
      scale_y_continuous(labels = scales::comma) +
      scale_fill_manual(values = c(Property = "steelblue", Crop = "darkgoldenrod2")) +
      labs(
        title    = "Economic consequences by event type (U.S., 1996-2011)",
        subtitle = "Property and crop damage stacked; bar height = total damage",
        x        = NULL,
        y        = "Damage ($ billions)",
        fill     = "Component",
        caption  = "Source: NOAA Storm Data"
      ) +
      theme_minimal(base_size = 11) +
      theme(
        plot.title    = element_text(face = "bold"),
        legend.position = "top"
      )
Figure 3. Total economic damage (property + crop, stacked) for the top 10 event types, 1996-2011. Hurricanes/typhoons lead; drought is the only leading event dominated by crop rather than property loss.

Figure 3. Total economic damage (property + crop, stacked) for the top 10 event types, 1996-2011. Hurricanes/typhoons lead; drought is the only leading event dominated by crop rather than property loss.

After correcting the erroneous $115B Napa flood record, hurricanes (typhoons) are the costliest event type ($87.1B), followed by storm surge/tide ($47.8B) and flood ($34.3B). Drought is distinctive as the only leading event driven by crop rather than property damage ($13.4B of $14.4B).