Brief exploratory analysis of the impact of storms and weather events in The United States

Synopsis

This analysis uses data from The US National Oceanic and Atmospheric Administration (NOAA) to briefly explore the impact of storms and weather events on the health of the US population and economy.
The data set that will be loaded into R to perform the analysis contains information on meteorological events from the year 1950 to November 2011. For the first years of the database there are usually fewer events recorded, most likely due to the lack of good records. More recent years can be considered more complete.
The following shows how the data loading and preprocessing was carried out, and subsequently the results of the exploratory analysis.

Data loading and preprocessing

This first section of the analysis shows how to load the data into R and how the processing of said data is carried out prior to analyzing the results in the following section.
First, the tidyverse library is loaded into R, then the data is downloaded to the project directory and then read and loaded into R:

library(tidyverse)
url <- "https://d396qusza40orc.cloudfront.net/repdata%2Fdata%2FStormData.csv.bz2"

if(file.exists("data_NOAA.csv.bz2")){
    raw_data <- read_csv(file = "data_NOAA.csv.bz2")
} else {
    download.file(url, destfile = "data_NOAA.csv.bz2")
    raw_data <- read_csv(file = "data_NOAA.csv.bz2")
}

The data structure is as follows:

glimpse(raw_data)

## Rows: 902,297
## Columns: 37
## $ STATE__    <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 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:0…
## $ BGN_TIME   <chr> "0130", "0145", "1600", "0900", "1500", "2000", "0100", "09…
## $ TIME_ZONE  <chr> "CST", "CST", "CST", "CST", "CST", "CST", "CST", "CST", "CS…
## $ COUNTY     <dbl> 97, 3, 57, 89, 43, 77, 9, 123, 125, 57, 43, 9, 73, 49, 107,…
## $ COUNTYNAME <chr> "MOBILE", "BALDWIN", "FAYETTE", "MADISON", "CULLMAN", "LAUD…
## $ STATE      <chr> "AL", "AL", "AL", "AL", "AL", "AL", "AL", "AL", "AL", "AL",…
## $ EVTYPE     <chr> "TORNADO", "TORNADO", "TORNADO", "TORNADO", "TORNADO", "TOR…
## $ BGN_RANGE  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ BGN_AZI    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ BGN_LOCATI <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ END_DATE   <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ END_TIME   <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ COUNTY_END <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ COUNTYENDN <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ END_RANGE  <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ END_AZI    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ END_LOCATI <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ LENGTH     <dbl> 14.0, 2.0, 0.1, 0.0, 0.0, 1.5, 1.5, 0.0, 3.3, 2.3, 1.3, 4.7…
## $ WIDTH      <dbl> 100, 150, 123, 100, 150, 177, 33, 33, 100, 100, 400, 400, 2…
## $ F          <dbl> 3, 2, 2, 2, 2, 2, 2, 1, 3, 3, 1, 1, 3, 3, 3, 4, 1, 1, 1, 1,…
## $ MAG        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ FATALITIES <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 4, 0, 0, 0, 0,…
## $ INJURIES   <dbl> 15, 0, 2, 2, 2, 6, 1, 0, 14, 0, 3, 3, 26, 12, 6, 50, 2, 0, …
## $ PROPDMG    <dbl> 25.0, 2.5, 25.0, 2.5, 2.5, 2.5, 2.5, 2.5, 25.0, 25.0, 2.5, …
## $ PROPDMGEXP <chr> "K", "K", "K", "K", "K", "K", "K", "K", "K", "K", "M", "M",…
## $ CROPDMG    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ CROPDMGEXP <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ WFO        <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ STATEOFFIC <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ ZONENAMES  <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ LATITUDE   <dbl> 3040, 3042, 3340, 3458, 3412, 3450, 3405, 3255, 3334, 3336,…
## $ LONGITUDE  <dbl> 8812, 8755, 8742, 8626, 8642, 8748, 8631, 8558, 8740, 8738,…
## $ LATITUDE_E <dbl> 3051, 0, 0, 0, 0, 0, 0, 0, 3336, 3337, 3402, 3404, 0, 3432,…
## $ LONGITUDE_ <dbl> 8806, 0, 0, 0, 0, 0, 0, 0, 8738, 8737, 8644, 8640, 0, 8540,…
## $ REMARKS    <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
## $ REFNUM     <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, …

The columns of interest for the analysis are the following:

EVTYPE, of type character. This variable comprises all weather events and storm types for which NOAA collects data.
FATALITIES, of double type. Records the number of deaths due to weather events and storms.
INJURIES, of double type. Records the number of people injured due to weather events and storms.
PROPDMG, of type double. Records the USD amount of property damage caused by weather events and storms.
CROPDMG, of type double. Records the amount in USD of crop damage caused by weather events and storms.
PROPDMGEXP, of type character. This variable reports whether property damage in USD recorded in PROPDMG is in hundreds (H), thousands (K), millions (M), or billions (B).
CROPDMGEXP, of type character. This variable reports whether damage in USD to crops recorded in CROPDMG is in hundreds (H), thousands (K), millions (M) or billions (B).

Actually, in the FATALITIES and INJURIES columns, all values are integers:

all(round(raw_data$FATALITIES, 0) == raw_data$FATALITIES)

## [1] TRUE

all(round(raw_data$INJURIES, 0) == raw_data$INJURIES)

## [1] TRUE

The values stored in PROPDMGEXP and CROPDMGEXP are as follows:

table(raw_data$PROPDMGEXP, useNA = "always")

## 
##      -      ?      +      0      1      2      3      4      5      6      7 
##      1      8      5    216     25     13      4      4     28      4      5 
##      8      B      h      H      K      m      M   <NA> 
##      1     40      1      6 424665      7  11330 465934

table(raw_data$CROPDMGEXP, useNA = "always")

## 
##      ?      0      2      B      k      K      m      M   <NA> 
##      7     19      1      9     21 281832      1   1994 618413

At this point all the letters in the two previous columns will be written in capital letters, the rest of the values that appear (numbers and symbols) will be considered missing values:

raw_data$PROPDMGEXP <- if_else(str_detect(raw_data$PROPDMGEXP, regex("[a-z]", ignore_case = T)),
                           str_to_upper(raw_data$PROPDMGEXP),
                           NA)
raw_data$CROPDMGEXP <- if_else(str_detect(raw_data$CROPDMGEXP, regex("[a-z]", ignore_case = T)),
                           str_to_upper(raw_data$CROPDMGEXP),
                           NA)

Now the values stored in the PROPDMGEXP and CROPDMGEXP columns are:

table(raw_data$PROPDMGEXP, useNA = "always")

## 
##      B      H      K      M   <NA> 
##     40      7 424665  11337 466248

table(raw_data$CROPDMGEXP, useNA = "always")

## 
##      B      K      M   <NA> 
##      9 281853   1995 618440

Let’s now look at the proportions of missing data for each of the variables of interest:

mean(is.na(raw_data$EVTYPE))

## [1] 0

mean(is.na(raw_data$FATALITIES))

## [1] 0

mean(is.na(raw_data$INJURIES))

## [1] 0

mean(is.na(raw_data$PROPDMG))

## [1] 0

mean(is.na(raw_data$CROPDMG))

## [1] 0

mean(is.na(raw_data$PROPDMGEXP))

## [1] 0.5167345

mean(is.na(raw_data$CROPDMGEXP))

## [1] 0.6854062

There is a lot of missing data in PROPDMGEXP and CROPDMGEXP!

All PROPDMGEXP and CROPDMGEXP values will now be expressed in USD units:

analytical_data <- raw_data %>%
  select(EVTYPE, FATALITIES, INJURIES, PROPDMG, PROPDMGEXP, CROPDMG, CROPDMGEXP) %>%
  mutate(
    PROPDMG = case_when(
      PROPDMGEXP == "H" ~ PROPDMG*10^2,
      PROPDMGEXP == "K" ~ PROPDMG*10^3,
      PROPDMGEXP == "M" ~ PROPDMG*10^6,
      PROPDMGEXP == "B" ~ PROPDMG*10^9
      ),
    CROPDMG = case_when(
      CROPDMGEXP == "K" ~ CROPDMG*10^3,
      CROPDMGEXP == "M" ~ CROPDMG*10^6,
      CROPDMGEXP == "B" ~ CROPDMG*10^9
      )
  )

Finally, a variable called ECONDMG will be created to obtain the total economic damages and will be equal to the sum of PROPDMG and CROPDMG:

analytical_data$ECONDMG <- analytical_data$PROPDMG + analytical_data$CROPDMG

This concludes the preprocessing of the data.

Results

1. Damage to people’s health

First of all, the events that cause the most fatalities are obtained with the following code:

fatalities <- analytical_data %>%
  group_by(EVTYPE) %>%
  summarise(FATALITIES = sum(FATALITIES)) %>%
  top_n(n = 10, wt = FATALITIES) %>%
  arrange(desc(FATALITIES))

Below are the above results in a bar graph:

library(scales)
fatalities %>%
  mutate(
    EVTYPE = case_when(
      EVTYPE == "RIP CURRENT" ~ "RIP current",
      EVTYPE == "TSTM WIND" ~ "TSTM wind",
      .default = str_to_sentence(EVTYPE))
  ) %>%
  ggplot(aes(x = fct_inorder(EVTYPE) %>% fct_rev(), y = FATALITIES)) +
    geom_col(fill = 18) +
    labs(title = "Fatalities by type of meteorological event in the USA",
         subtitle = "Values according to NOAA records, period 1950 - 2011",
         x = "Event type",
         y = "Fatalities") +
    geom_text(aes(y = FATALITIES + 200), 
              label = format(fatalities$FATALITIES, big.mark = ",")) +
    scale_y_continuous(labels = label_number(big.mark = ",")) +
    coord_flip() +
    theme_light()

Secondly, the events that leave the most people injured are obtained with the following code:

injuries <- analytical_data %>%
  group_by(EVTYPE) %>%
  summarise(INJURIES = sum(INJURIES)) %>%
  top_n(n = 10, wt = INJURIES) %>%
  arrange(desc(INJURIES))

Below are the above results in a bar graph:

injuries %>%
  mutate(
    EVTYPE = case_when(
      EVTYPE == "TSTM WIND" ~ "TSTM wind",
      .default = str_to_sentence(EVTYPE))
  ) %>%
  ggplot(aes(x = fct_inorder(EVTYPE) %>% fct_rev(), y = INJURIES)) +
    geom_col(fill = 29) +
    labs(title = "Injuries by type of meteorological event in the USA",
         subtitle = "Values according to NOAA records, period 1950 - 2011",
         x = "Event type",
         y = "Injuries") +
    geom_text(aes(y = INJURIES + 4000), 
              label = format(injuries$INJURIES, big.mark = ",")) +
    scale_y_continuous(limits = c(0, 100000), labels = label_number(big.mark = ",")) +
    coord_flip() +
    theme_light()

2. Economic damages

The events that produce the most economic damage are obtained with the following code:

econ_dmg <- analytical_data %>%
  group_by(EVTYPE) %>%
  summarise(ECONDMG = sum(ECONDMG)) %>%
  top_n(n = 10, wt = ECONDMG) %>%
  arrange(desc(ECONDMG))

From the previous results, the following graph is obtained:

econ_dmg %>%
  mutate(
    EVTYPE = case_when(
      EVTYPE == "TORNADOES, TSTM WIND, HAIL" ~ "Tornadoes, TSTM wind, hail",
      .default = str_to_sentence(EVTYPE)
    )
  ) %>%
  ggplot(aes(x = fct_inorder(EVTYPE) %>% fct_rev(), y = ECONDMG)) +
    geom_col(fill = 71) +
    labs(title = "Economic damages by type of meteorological event in the USA",
         subtitle = "Values according to NOAA records, period 1950 - 2011",
         x = "Event type",
         y = "Damages in USD") +
    geom_text(aes(y = ECONDMG + 170000000), 
              label = dollar(econ_dmg$ECONDMG, scale = 1/10^6, suffix = "M")) +
    scale_y_continuous(limits = c(0, 2000000000),
                       labels = label_dollar(scale = 1/10^6, suffix = "M")) +
    coord_flip() +
    theme_light()

Conclusions

The event that causes the most fatalities is the tornado type event.
The event that leaves the most people injured is the tornado.
The events that cause the most economic damage are tornado-type events, storm winds, and hail.