Reproducible Research Course Project 2
Mark Chang
2026-10-19
Introduction / Synopsis
This project involves exploring the U.S. National Oceanic and Atmospheric Administration’s (NOAA) storm database. We will answer several questions related to a subset of the data in this database.
Data
The data for this project comes from NOAA; it is downloadable as a bzip2 compressed file.
Information on how the variables are constructed can be found from the following links :-
- National Weather Service Storm Data Documentation
- National Climatic Data Center Storm Events
The events in the database start in the year 1950 and end in November 2011. Less data is available in the earlier years of the database.
Assignment
The data analysis must address the following qwuestions:-
Across the United States, which types of events (as indicated in the EVTYPE variable) are most harmful with respect to population health?
Across the United States, which types of events have the greatest economic consequences?
Software Environment
sessionInfo()## R version 4.6.1 (2026-06-24)
## Platform: aarch64-apple-darwin23
## Running under: macOS Sequoia 15.7.9
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Pacific/Auckland
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.39 R6_2.6.1 fastmap_1.2.0 xfun_0.61
## [5] cachem_1.1.0 knitr_1.52 htmltools_0.5.9 rmarkdown_2.32
## [9] lifecycle_1.0.5 cli_3.6.6 sass_0.4.10 jquerylib_0.1.4
## [13] compiler_4.6.1 rstudioapi_0.19.0 tools_4.6.1 evaluate_1.0.5
## [17] bslib_0.12.0 yaml_2.3.12 otel_0.2.0 jsonlite_2.0.0
## [21] rlang_1.3.0Initialisation
Set to show R code
library(data.table)
library(ggplot2)
library(knitr)
library(dplyr)
library(tidyr)
opts_chunk$set(echo=TRUE)Data Processing
We will download and read the data into a data table
stormdata
# Windows
#setwd("C:/Users/mchang/Projects/ProgrammingAssignment7")
# Mac
setwd("/Users/mchang/Projects/ProgrammingAssignment7")
path <- getwd()
datasetpath <- file.path(path,"dataset")
datafilename <- "dataFiles.bz2"
# Download file
url <- "https://d396qusza40orc.cloudfront.net/repdata%2Fdata%2FStormData.csv.bz2"
#download.file(url, file.path(path, datafilename))
#unzip(zipfile = datafilename, exdir = datasetpath)
stormdata <- data.table::fread("repdata_data_StormData.csv", fill=TRUE,header=TRUE)We will take a look at the data
head(stormdata)## STATE__ BGN_DATE BGN_TIME TIME_ZONE COUNTY COUNTYNAME STATE
## <char> <char> <char> <char> <char> <char> <char>
## 1: 1.00 4/18/1950 0:00:00 0130 CST 97.00 MOBILE AL
## 2: 1.00 4/18/1950 0:00:00 0145 CST 3.00 BALDWIN AL
## 3: 1.00 2/20/1951 0:00:00 1600 CST 57.00 FAYETTE AL
## 4: 1.00 6/8/1951 0:00:00 0900 CST 89.00 MADISON AL
## 5: 1.00 11/15/1951 0:00:00 1500 CST 43.00 CULLMAN AL
## 6: 1.00 11/15/1951 0:00:00 2000 CST 77.00 LAUDERDALE AL
## EVTYPE BGN_RANGE BGN_AZI BGN_LOCATI END_DATE END_TIME COUNTY_END COUNTYENDN
## <char> <char> <char> <char> <char> <char> <char> <char>
## 1: TORNADO 0.00 0.00
## 2: TORNADO 0.00 0.00
## 3: TORNADO 0.00 0.00
## 4: TORNADO 0.00 0.00
## 5: TORNADO 0.00 0.00
## 6: TORNADO 0.00 0.00
## END_RANGE END_AZI END_LOCATI LENGTH WIDTH F MAG FATALITIES INJURIES
## <char> <char> <char> <char> <char> <char> <char> <char> <char>
## 1: 0.00 14.00 100.00 3 0.00 0.00 15.00
## 2: 0.00 2.00 150.00 2 0.00 0.00 0.00
## 3: 0.00 0.10 123.00 2 0.00 0.00 2.00
## 4: 0.00 0.00 100.00 2 0.00 0.00 2.00
## 5: 0.00 0.00 150.00 2 0.00 0.00 2.00
## 6: 0.00 1.50 177.00 2 0.00 0.00 6.00
## PROPDMG PROPDMGEXP CROPDMG CROPDMGEXP WFO STATEOFFIC ZONENAMES LATITUDE
## <num> <char> <num> <char> <char> <char> <char> <num>
## 1: 25.0 K 0 3040
## 2: 2.5 K 0 3042
## 3: 25.0 K 0 3340
## 4: 2.5 K 0 3458
## 5: 2.5 K 0 3412
## 6: 2.5 K 0 3450
## LONGITUDE LATITUDE_E LONGITUDE_ REMARKS REFNUM
## <num> <num> <num> <char> <num>
## 1: 8812 3051 8806 1
## 2: 8755 0 0 2
## 3: 8742 0 0 3
## 4: 8626 0 0 4
## 5: 8642 0 0 5
## 6: 8748 0 0 6str(stormdata)## Classes 'data.table' and 'data.frame': 902297 obs. of 37 variables:
## $ STATE__ : chr "1.00" "1.00" "1.00" "1.00" ...
## $ 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 : chr "97.00" "3.00" "57.00" "89.00" ...
## $ COUNTYNAME: chr "MOBILE" "BALDWIN" "FAYETTE" "MADISON" ...
## $ STATE : chr "AL" "AL" "AL" "AL" ...
## $ EVTYPE : chr "TORNADO" "TORNADO" "TORNADO" "TORNADO" ...
## $ BGN_RANGE : chr "0.00" "0.00" "0.00" "0.00" ...
## $ BGN_AZI : chr "" "" "" "" ...
## $ BGN_LOCATI: chr "" "" "" "" ...
## $ END_DATE : chr "" "" "" "" ...
## $ END_TIME : chr "" "" "" "" ...
## $ COUNTY_END: chr "0.00" "0.00" "0.00" "0.00" ...
## $ COUNTYENDN: chr "" "" "" "" ...
## $ END_RANGE : chr "0.00" "0.00" "0.00" "0.00" ...
## $ END_AZI : chr "" "" "" "" ...
## $ END_LOCATI: chr "" "" "" "" ...
## $ LENGTH : chr "14.00" "2.00" "0.10" "0.00" ...
## $ WIDTH : chr "100.00" "150.00" "123.00" "100.00" ...
## $ F : chr "3" "2" "2" "2" ...
## $ MAG : chr "0.00" "0.00" "0.00" "0.00" ...
## $ FATALITIES: chr "0.00" "0.00" "0.00" "0.00" ...
## $ INJURIES : chr "15.00" "0.00" "2.00" "2.00" ...
## $ 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 ...
## - attr(*, ".internal.selfref")=<pointer: 0x12a0114e0>Based on the information above, we have 902,297 rows of observations and 37 columns of variables names. Since the questions are only interested in harmful events with respect to population health and economic consequences, we will focus on only the relevant variables.
names(stormdata)## [1] "STATE__" "BGN_DATE" "BGN_TIME" "TIME_ZONE" "COUNTY"
## [6] "COUNTYNAME" "STATE" "EVTYPE" "BGN_RANGE" "BGN_AZI"
## [11] "BGN_LOCATI" "END_DATE" "END_TIME" "COUNTY_END" "COUNTYENDN"
## [16] "END_RANGE" "END_AZI" "END_LOCATI" "LENGTH" "WIDTH"
## [21] "F" "MAG" "FATALITIES" "INJURIES" "PROPDMG"
## [26] "PROPDMGEXP" "CROPDMG" "CROPDMGEXP" "WFO" "STATEOFFIC"
## [31] "ZONENAMES" "LATITUDE" "LONGITUDE" "LATITUDE_E" "LONGITUDE_"
## [36] "REMARKS" "REFNUM"General
From the links regarding the variable construction and naming, we choose only the following variables.
- evtype : storm event type
- fatalities: amount of fatalities per event
- injuries : amount of injuries per event
- propdmg : property damage amount
- propdmgexp: property damage in exponents
- cropdmg : crop damage amount
- cropdmgexp: crop damage in exponents
Only these variables are chosen; everything else is discarded to make
analysis simpler. For readability, all variables are transformed to
lowercase and the fatalities and injuries columns made numeric. We will
read this into a new data table stormdata2
stormdata2 <- stormdata %>% select(c("EVTYPE","FATALITIES","INJURIES","PROPDMG","PROPDMGEXP","CROPDMG","CROPDMGEXP")) %>% rename_all(tolower)
stormdata2$fatalities <- as.numeric(stormdata2$fatalities)
stormdata2$injuries <- as.numeric(stormdata2$injuries)
str(stormdata2)## Classes 'data.table' and '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 "" "" "" "" ...
## - attr(*, ".internal.selfref")=<pointer: 0x12a0114e0>stormevents <- length(unique(stormdata2$evtype))We have 985 unique types of storm events.
Injuries and Fatalities Involved
We can compute the total number of fatalities and injuries and group
them by storm event type into a new data frame
stormhealth1. We transform the data by using the
gather function to create a single column called type
(which contains the type “property” or “crop”) and a single value column
called value which contains the counts. This format is required for
ggplot2 later on.
stormhealth1 <- stormdata2 %>%
select("evtype","fatalities","injuries") %>%
group_by(evtype) %>%
summarise(totalfatalities = sum(fatalities), totalinjuries = sum(injuries), .groups='drop') %>%
arrange(desc(totalfatalities),desc(totalinjuries)) %>%
slice(1:10) %>%
gather(key=type,value=value,totalfatalities,totalinjuries)Economic Costs involved
propexp <- unique(stormdata2$propdmgexp)
cropexp <- unique(stormdata2$cropdmgexp)Looking at the *exp fields, these seem to carry either the exponents
or the place value indicators. propexp carries these values
K, M, , B, m, +, 0, 5, 6, ?, 4, 2, 3, h, 7, H, -, 1, 8 and
cropexp carries these values , M, K, m, B, ?, 0, k, 2.
We will transform them into values that can be used for actual calculations. We set all NA values to 0.
dmgexpkey <- c("-" = 10^0,
"+" = 10^0,
"0" = 10^0,
"1" = 10^1,
"2" = 10^2,
"3" = 10^3,
"4" = 10^4,
"5" = 10^5,
"6" = 10^6,
"7" = 10^7,
"8" = 10^8,
"9" = 10^9,
"H" = 10^2,
"h" = 10^2,
"K" = 10^3,
"k" = 10^3,
"M" = 10^6,
"m" = 10^6,
"B" = 10^9,
"b" = 10^9
)
stormdata2[, propdmgexp := dmgexpkey[as.character(stormdata2[,propdmgexp])]]
stormdata2[, cropdmgexp := dmgexpkey[as.character(stormdata2[,cropdmgexp])]]
stormdata2[is.na(propdmgexp), propdmgexp := 10^0]
stormdata2[is.na(cropdmgexp), cropdmgexp := 10^0]We can now compute the actual property damage and crop damage costs
by adding calculated columns propdmgvalue and
cropdmgvalue into a new data frame
stormdamage1. We then proceed to group all similar evtypes
together into stormdamage2. We transform the data by using
the gather function to create a single column called type
(which contains the type evtype) and a single value column called value
which contains the dollar values. This format is required for ggplot2
later on.
stormdamage1 <- stormdata2 %>%
select("evtype", "propdmg", "propdmgexp", "cropdmg", "cropdmgexp") %>%
mutate(propdmgvalue = propdmg * propdmgexp) %>%
mutate(cropdmgvalue = cropdmg * cropdmgexp)
stormdamage2 <- stormdamage1 %>%
group_by(evtype) %>%
summarise(totalprop = sum(propdmgvalue), totalcrop = sum(cropdmgvalue), .groups='drop') %>%
arrange(desc(totalprop),desc(totalcrop)) %>%
slice(1:10) %>%
gather(key=type,value=value,totalprop,totalcrop)Results
We are now ready to answer the questions posed in the assignment.
Population Health Analysis
Question 1 : Across the United States, which types of events (as indicated in the EVTYPE variable) are most harmful with respect to population health?
ggplot(stormhealth1,aes(reorder(evtype,-value), value, fill=type)) + geom_bar(stat="identity", position="dodge") + labs(x="Event Type", y="Count",title="Storm Damage") + theme_bw() + theme(axis.text.x = element_text(angle=20,vjust=1,hjust=1))
knitr::kable(stormhealth1)| evtype | type | value |
|---|---|---|
| TORNADO | totalfatalities | 5633 |
| EXCESSIVE HEAT | totalfatalities | 1903 |
| FLASH FLOOD | totalfatalities | 978 |
| HEAT | totalfatalities | 937 |
| LIGHTNING | totalfatalities | 816 |
| TSTM WIND | totalfatalities | 504 |
| FLOOD | totalfatalities | 470 |
| RIP CURRENT | totalfatalities | 368 |
| HIGH WIND | totalfatalities | 248 |
| AVALANCHE | totalfatalities | 224 |
| TORNADO | totalinjuries | 91346 |
| EXCESSIVE HEAT | totalinjuries | 6525 |
| FLASH FLOOD | totalinjuries | 1777 |
| HEAT | totalinjuries | 2100 |
| LIGHTNING | totalinjuries | 5230 |
| TSTM WIND | totalinjuries | 6957 |
| FLOOD | totalinjuries | 6789 |
| RIP CURRENT | totalinjuries | 232 |
| HIGH WIND | totalinjuries | 1137 |
| AVALANCHE | totalinjuries | 170 |
Answer : Based on the bar plot, tornadoes have the highest impact on population health since it causes the most fatalities and injuries. It is followed by excessive heat and floods.
Economic Consequences Analysis
ggplot(stormdamage2,aes(reorder(evtype,-value), value / 10^6, fill=type)) + geom_bar(stat="identity", position="dodge") + labs(x="Event Type", y="Damage (in millions)",title="Storm Damage") + theme_bw() + theme(axis.text.x = element_text(angle=20,vjust=1,hjust=1))
knitr::kable(stormdamage2)| evtype | type | value |
|---|---|---|
| FLOOD | totalprop | 144657709807 |
| HURRICANE/TYPHOON | totalprop | 69305840000 |
| TORNADO | totalprop | 56947380676 |
| STORM SURGE | totalprop | 43323536000 |
| FLASH FLOOD | totalprop | 16822673978 |
| HAIL | totalprop | 15735267513 |
| HURRICANE | totalprop | 11868319010 |
| TROPICAL STORM | totalprop | 7703890550 |
| WINTER STORM | totalprop | 6688497251 |
| HIGH WIND | totalprop | 5270046295 |
| FLOOD | totalcrop | 5661968450 |
| HURRICANE/TYPHOON | totalcrop | 2607872800 |
| TORNADO | totalcrop | 414953270 |
| STORM SURGE | totalcrop | 5000 |
| FLASH FLOOD | totalcrop | 1421317100 |
| HAIL | totalcrop | 3025954473 |
| HURRICANE | totalcrop | 2741910000 |
| TROPICAL STORM | totalcrop | 678346000 |
| WINTER STORM | totalcrop | 26944000 |
| HIGH WIND | totalcrop | 638571300 |
Answwer : Based on the bar plot, floods have the highest economic impact on both property and crop damages. It is followed by hurricanes/typhoons and tornadoes.