Reproducible Research: Peer Assessment 2
Dean Fantazzini
2 March 2016
Analysis of the health and economic impacts of different weather events in the USA based on the NOAA Storm Database
1. Synopsis
The U.S. National Oceanic and Atmospheric Administration’s (NOAA) storm database collects the characteristics of major storms and weather events in the United States, including the date, location, estimates of any fatalities, injuries, and property damage. This report examined which types of events are the most harmful with respect to the population health and have the greatest economic consequences. The analysis found that TORNADO is the most harmful event in terms of total number of casualties, while HEAT is the worst event in terms of fatalities. As for economic impacts, FLOOD is the worst event in terms of total economic impact, while HEAT is the worst event in terms of crop damage.
2. Data Processing
2.1 Load libraries
library(ggplot2) # for plots
library(reshape2) # for melt
library(gridExtra) # for advanced plots2.2 Download the data and load it into R
Download the file from the internet
if(!file.exists('StormData.csv.bz2')){
download.file("http://d396qusza40orc.cloudfront.net/repdata%2Fdata%2FStormData.csv.bz2",
destfile='StormData.csv.bz2')
}Uncompress the file and load the CSV file into a data file
if(file.exists('StormData.csv.bz2')){
data <- read.csv(bzfile('StormData.csv.bz2'), header = TRUE)
}2.3 Clean and prepare the dataset
The BGN_DATE column is transformed to a date format and the YEAR column is created to indicate in which year the weather event took place
data <- transform(data, YEAR = as.numeric(format(as.Date(BGN_DATE, format = "%M/%d/%Y"), "%Y")))To answer the two questions of the assignment, only the following 9 columns are needed:
- YEAR: it indicates the when the weather event took place
- REFNUM: it is the reference number of the event type
- EVTYPE: it describes the weather event (e.g. tornado, flood, etc.)
- FATALITIES: it is a measure of the harm to human health
- INJURIES: it is a a measure of harm to human health
- PROPDMG: it is a a measure of property damage in USD
- PROPDMGEXP: it is a a measure of the magnitude of property damage (e.g. thousands, millions USD, etc.)
- CROPDMG: it is a a measure of crop damage in USD
- CROPDMGEXP: it is a as a measure of the magnitude of crop damage (e.g. thousands, millions USD, etc.)
Therefore, we select only these columns to decrease the subsequent computational burden:
data <- data[ , c("REFNUM", "YEAR", "EVTYPE", "FATALITIES", "INJURIES", "PROPDMG", "PROPDMGEXP", "CROPDMG", "CROPDMGEXP")]
str(data)## 'data.frame': 902297 obs. of 9 variables:
## $ REFNUM : num 1 2 3 4 5 6 7 8 9 10 ...
## $ YEAR : num 1950 1950 1951 1951 1951 ...
## $ EVTYPE : Factor w/ 985 levels " HIGH SURF ADVISORY",..: 834 834 834 834 834 834 834 834 834 834 ...
## $ 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: Factor w/ 19 levels "","-","?","+",..: 17 17 17 17 17 17 17 17 17 17 ...
## $ CROPDMG : num 0 0 0 0 0 0 0 0 0 0 ...
## $ CROPDMGEXP: Factor w/ 9 levels "","?","0","2",..: 1 1 1 1 1 1 1 1 1 1 ...The PROPDMGEXP and CROPDMGEXP columns indicate the multiplier of the PROPDMG and CROPDMG columns, respectively:
unique(data$PROPDMGEXP)## [1] K M B m + 0 5 6 ? 4 2 3 h 7 H - 1 8
## Levels: - ? + 0 1 2 3 4 5 6 7 8 B h H K m Munique(data$CROPDMGEXP)## [1] M K m B ? 0 k 2
## Levels: ? 0 2 B k K m MMore specifically, the meaning of the previous symbols is given below:
- “B”,“b”,“9” meaning the values are in Billions -> 1 000 000 000
- “M”,“m”,“6” meaning the values are in Millions -> 1 000 000
- “K”,“k”,“3” meaning the values are in Thousands, -> 1 000
- “0”,“1”,“2”,.,“8” values are in 1, 10, 100,. ect
- “+”,“-”,“?” invalid values -> 0.
Therefore, these columns are properly transformed. First, the property damage data:
data$PROPDMGEXP <- as.character(data$PROPDMGEXP)
data$PROPDMGEXP <- gsub("\\-|\\?|\\+|h|H|0","0",data$PROPDMGEXP)
data$PROPDMGEXP <- gsub("k|K", "1000", data$PROPDMGEXP)
data$PROPDMGEXP <- gsub("m|M", "1000000", data$PROPDMGEXP)
data$PROPDMGEXP <- gsub("b|B", "1000000000", data$PROPDMGEXP)
data$PROPDMGEXP <- as.numeric(data$PROPDMGEXP)
data$PROPDMGEXP[is.na(data$PROPDMGEXP)] <- 0
#A new column is created containing the total property damage in USD
data$PROPERTIES <- data$PROPDMG * data$PROPDMGEXPThen, the same cleaning procedure is performed with the crop damage data:
data$CROPDMGEXP <- as.character(data$CROPDMGEXP)
data$CROPDMGEXP <- gsub("\\-|\\?|\\+|h|H|0","0",data$CROPDMGEXP)
data$CROPDMGEXP <- gsub("k|K", "1000", data$CROPDMGEXP)
data$CROPDMGEXP <- gsub("m|M", "1000000", data$CROPDMGEXP)
data$CROPDMGEXP <- gsub("b|B", "1000000000", data$CROPDMGEXP)
data$CROPDMGEXP <- as.numeric(data$CROPDMGEXP)
data$CROPDMGEXP[is.na(data$CROPDMGEXP)] <- 0
#A new column is created containing the total crop damage in USD
data$CROPS <- data$CROPDMG * data$CROPDMGEXP2.4 Aggregate and reshape the datasets for making plots with time series data
The total harm and damages are computed and the data are ordered:
data$HARM <- data$FATALITIES + data$INJURIES
data$DAMAGE <- data$PROPERTIES + data$CROPS
data <- data[with(data, order(YEAR, REFNUM)), ]The data are aggregated by YEAR:
Fatal<- aggregate(FATALITIES ~ YEAR, data, sum)
Injur<- aggregate(INJURIES ~ YEAR, data, sum)
Prop <- aggregate(PROPERTIES ~ YEAR, data, sum)
Crop <- aggregate(CROPS ~ YEAR, data, sum)
AnnualDamage <- cbind(Fatal, Injur, Prop, Crop)[, c(1, 2, 4, 6, 8)]2.5 Plot the time series data for population health and economic damages from 1950 till 2011
YearHRM <- melt(AnnualDamage[1:3], id.vars = "YEAR")
gyhrm <- ggplot(YearHRM, aes(x = YEAR, y = value, fill = variable))
gyhrm <- gyhrm + geom_area(stat = "identity")+theme_bw()
gyhrm <- gyhrm + labs(title = "Casualties of Adverse Weather", x = "", y = "No. of People", fill = "")
YearDMG <- melt(AnnualDamage[c(1, 4:5)], id.vars = "YEAR")
gydmg <- ggplot(YearDMG, aes(x = YEAR, y = value/10^9, fill = variable))
gydmg <- gydmg + geom_area(stat = "identity")+theme_bw()
gydmg <- gydmg + labs(title = "Economic Impact of Adverse Weather", x = "", y = "Billion USD", fill = "")
grid.arrange(gyhrm, gydmg, nrow = 2)
The previous plot shows that the number of data before 1993 is clearly limited, and particularly there is no data for crops impact before 1993. Threfore, only data from 1993 onwards will be used for the subsequent analysis.
data <- subset(data, YEAR >= 1993)2.6 Aggregate weather events
If the EVTYPE column is examined,
length(unique(data$EVTYPE))## [1] 985head(unique(data$EVTYPE), 40)## [1] WINTER STORM HAIL
## [3] THUNDERSTORM WINDS HIGH WINDS
## [5] TORNADO FUNNEL CLOUD
## [7] THUNDERSTORM WIND WALL CLOUD
## [9] FLASH FLOOD LIGHTNING
## [11] HAIL 1.75) BLIZZARD
## [13] RIVER FLOOD HEAVY SNOW
## [15] HIGH WIND WIND CHILL
## [17] HIGH WIND/ BLIZZARD ICE STORM
## [19] COASTAL FLOOD FLOOD WATCH/
## [21] RECORD HIGH TEMPERATURES HIGH WIND/SEAS
## [23] HIGH WINDS/HEAVY RAIN HIGH SEAS
## [25] RECORD COLD SEVERE TURBULENCE
## [27] RECORD RAINFALL RECORD SNOWFALL
## [29] RECORD WARMTH HEAVY SNOW/WIND
## [31] FLOOD HEAVY RAIN
## [33] GUSTY WINDS FUNNEL CLOUDS
## [35] URBAN/SMALL WILD FIRES
## [37] URBAN/SMALL FLOODING SNOW
## [39] WINTER STORMS FLOODING
## 985 Levels: HIGH SURF ADVISORY COASTAL FLOOD ... WNDit is possible to notice that there are 985 types of weather events, but several of them have a similar meaning if not completely identical, such as “THUNDER STORM WIND” and “THUNDER STORM WINDS”, etc. To aggregate these duplicated events, some keywords are selected for each weather event and replaced with 17 unique weather conditions. All remaining events are then included into “OTHERS”.
#Change to capital letters
data$EVTYPE <- toupper(data$EVTYPE)
#New general weather event types
EV <- c("FLOOD", "TORNADO", "STORM", "HEAT", "HURRICANE", "FIRE", "HIGH WIND", "LIGHTNING", "HAIL", "COLD", "SNOW", "RAIN", "HIGH SURF", "FOG", "LANDSLIDE", "TSUNAMI", "AVALANCHE")
#Substitute similar weather events with a single unique event
data$EVTYPE[grep("FLOOD|URBAN|TIDE", data$EVTYPE)] <- "FLOOD"
data$EVTYPE[grep("HURRICANE|TYPHOON", data$EVTYPE)] <- "HURRICANE"
data$EVTYPE[grep("STORM|BLIZZARD|TSTM", data$EVTYPE)] <- "STORM"
data$EVTYPE[grep("TORNADO|SPOUT|FUNNEL", data$EVTYPE)] <- "TORNADO"
data$EVTYPE[grep("HAIL", data$EVTYPE)] <- "HAIL"
data$EVTYPE[grep("HEAT|WARM|HOT|DROUGHT|DRY|HIGH TEMP", data$EVTYPE)] <- "HEAT"
data$EVTYPE[grep("FIRE", data$EVTYPE)] <- "FIRE"
data$EVTYPE[grep("WIND|BURST|DUST", data$EVTYPE)] <- "HIGH WIND"
data$EVTYPE[grep("RAIN|WET", data$EVTYPE)] <- "RAIN"
data$EVTYPE[grep("COLD|FREEZ|FROST|GLAZE|ICE|ICY|WINTER", data$EVTYPE)] <- "COLD WX"
data$EVTYPE[grep("SNOW", data$EVTYPE)] <- "SNOW"
data$EVTYPE[grep("LIGHTNING", data$EVTYPE)] <- "LIGHTNING"
data$EVTYPE[grep("SLIDE", data$EVTYPE)] <- "LANDSLIDE"
data$EVTYPE[grep("TSUNAMI", data$EVTYPE)] <- "TSUNAMI"
data$EVTYPE[grep("CURRENT|SURF|HIGH WAVES", data$EVTYPE)] <- "HIGH SURF"
data$EVTYPE[grep("FOG", data$EVTYPE)] <- "FOG"
data$EVTYPE[grep("AVALANC", data$EVTYPE)] <- "AVALANCHE"
#All remaining events are then included into "OTHERS".
others <- grep(paste(EV, collapse="|"), data$EVTYPE, invert = TRUE)
data[others, ]$EVTYPE <- "OTHERS" Aggregate the final dataset by EVTYPE and order it:
#Helth impacts:
human <- aggregate(cbind(FATALITIES, INJURIES, HARM) ~ EVTYPE, data, sum)
# sort by HARM
human <- human[with(human, order(HARM)), ]
# Save ascending order for later use
ord_h <- human$EVTYPE
# Reorder in descenting order
human <- human[with(human, order(-HARM)), ]
colnames(human)[1] <- "HEALTH"
#Economic impacts:
damage <- aggregate(cbind(PROPERTIES, CROPS, DAMAGE) ~ EVTYPE, data, sum)
# Sort by damage
damage <- damage[with(damage, order(DAMAGE)), ]
# Save ascending order for later use
ord_d <- damage$EVTYPE
# Reorder in descenting order
damage <- damage[with(damage, order(-DAMAGE)), ]
colnames(damage)[1] <- "ECONOMIC"3. Results
3.1 Which types of events are most harmful to population health?
# Melt the final data for health impacts
hmelt <- melt(human[1:3], id.vars = "HEALTH")
# Give EVTYPE a sorting order
hmelt$HEALTH <- factor(hmelt$HEALTH, levels = ord_h)
# Plot the data
hgraph <- ggplot(hmelt, aes(x = HEALTH, y = value, fill = variable))
hgraph <- hgraph + labs(title = "Casualties (Fatalities and Injuries) of Adverse Weather", x = "Event Type", y = "Number of People", fill = "")
hgraph <- hgraph + geom_bar(stat = "identity") +theme_bw()+ coord_flip()
print(hgraph)
TORNADO is the most dangerous event in terms of total number of casualties, while HEAT is the worst event in terms of fatalities.
3.2 Which types of events have the greatest economic consequences?
# Melt the final data for economic impacts
dmelt <- melt(damage[1:3], id.vars = "ECONOMIC")
# Give EVTYPE a sorting order
dmelt$ECONOMIC <- factor(dmelt$ECONOMIC, levels = ord_d)
# Plot the data
dgraph <- ggplot(dmelt, aes(x = ECONOMIC, y = value/10^9, fill = variable))
dgraph <- dgraph + labs(title = "Economic Impact of Adverse Weather", x = "Event Type", y = "Billions USD", fill = "")
dgraph <- dgraph + geom_bar(stat = "identity") +theme_bw()+ coord_flip()
print(dgraph)
FLOOD is the worst event in terms of total economic impact, while HEAT is the worst event in terms of crop damage.
4. Conclusions
The initial inspection of the NOAA dataset showed that the number of data before 1993 is very limited, without any data for crops impact before 1993. Therefore, only data from 1993 onwards were used in the final analysis. Moreover, even though the NOAA dataset contains 985 types of weather events, several of them have a similar meaning if not completely identical. The weather events with similar or identical meaning were replaced with 17 unique weather conditions, while the remaining unclassified events were assigned to the “OTHERS” event type.
The final analysis found that TORNADO is the most harmful event in terms of total number of casualties, while HEAT is the worst event in terms of fatalities. As for economic impacts, FLOOD is the worst event in terms of total economic impact, while HEAT is the worst event in terms of crop damage.