Severe Weather Events in US and their Consequences
Thymios C
2 December 2016
Synopsis
In this report I will explore the U.S. National Oceanic and Atmospheric Administration’s (NOAA) storm database and try to extract useful information from it. The questions that I’ll try to answer are the following:
Across the United States, which types of events are most harmful with respect to population health?
Across the United States, which types of events have the greatest economic consequences?
Moreover I will show which States are more prone to severe weather events.
Data Processing
For this report I have used R version 3.3.2
The data that we are going to analyse were downloaded from here(47Mb)
Load libraries.
library(ggplot2)
library(plyr)
library(dplyr)
library(RecordLinkage)
library(lubridate)
library(datasets)
library(scales)
library(gridExtra)
library(stringr)Import the Data set.
#download.file("https://d396qusza40orc.cloudfront.net/repdata%2Fdata%2FStormData.csv.bz2",
# "StormData.csv.bz2")
rData <- read.csv("StormData.csv.bz2")Explore the structure of the data set.
dim(rData)## [1] 902297 37str(rData)## 'data.frame': 902297 obs. of 37 variables:
## $ STATE__ : num 1 1 1 1 1 1 1 1 1 1 ...
## $ BGN_DATE : Factor w/ 16335 levels "1/1/1966 0:00:00",..: 6523 6523 4242 11116 2224 2224 2260 383 3980 3980 ...
## $ BGN_TIME : Factor w/ 3608 levels "00:00:00 AM",..: 272 287 2705 1683 2584 3186 242 1683 3186 3186 ...
## $ TIME_ZONE : Factor w/ 22 levels "ADT","AKS","AST",..: 7 7 7 7 7 7 7 7 7 7 ...
## $ COUNTY : num 97 3 57 89 43 77 9 123 125 57 ...
## $ COUNTYNAME: Factor w/ 29601 levels "","5NM E OF MACKINAC BRIDGE TO PRESQUE ISLE LT MI",..: 13513 1873 4598 10592 4372 10094 1973 23873 24418 4598 ...
## $ STATE : Factor w/ 72 levels "AK","AL","AM",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ EVTYPE : Factor w/ 985 levels " HIGH SURF ADVISORY",..: 834 834 834 834 834 834 834 834 834 834 ...
## $ BGN_RANGE : num 0 0 0 0 0 0 0 0 0 0 ...
## $ BGN_AZI : Factor w/ 35 levels ""," N"," NW",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ BGN_LOCATI: Factor w/ 54429 levels "","- 1 N Albion",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ END_DATE : Factor w/ 6663 levels "","1/1/1993 0:00:00",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ END_TIME : Factor w/ 3647 levels ""," 0900CST",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ 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 : Factor w/ 24 levels "","E","ENE","ESE",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ END_LOCATI: Factor w/ 34506 levels "","- .5 NNW",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ 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: 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 ...
## $ WFO : Factor w/ 542 levels ""," CI","$AC",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ STATEOFFIC: Factor w/ 250 levels "","ALABAMA, Central",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ ZONENAMES : Factor w/ 25112 levels ""," "| __truncated__,..: 1 1 1 1 1 1 1 1 1 1 ...
## $ 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 : Factor w/ 436781 levels "","-2 at Deer Park\n",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ REFNUM : num 1 2 3 4 5 6 7 8 9 10 ...For this report I’m only going to use the following columns:
- STATE
- EVTYPE
- FATALITIES
- INJURIES
- PROPDMG
- CROPDMG
Remove the rest of the columns
tData <- rData[,c(7,8,23,24,25,27)]Look for missing values
table(is.na(tData))##
## FALSE
## 5413782There are no missing values.
In this report I’m only interested for events that were harmful with respect to population health or to economy so I will exclude all the rest rows.
tData <- with(tData, tData[FATALITIES >0 | INJURIES >0 | PROPDMG >0 | CROPDMG >0,])Explore columns
STATE
STATE is a factor variable with 72 levels that is obviously wrong. Since state is not an important element of this report I will just map the misspelled states to a “similar one”.
# Using mutate (dplyr) to replace states with the correct ones from state.abb (datasets)
# using levenshteinDist (RecordLinkage) to count distance between strings.
tData$STATE <- sapply(tData$STATE,
function(x) state.abb[which.min(levenshteinDist(x,state.abb))])
#Convert to factor & rename
tData$STATE <- as.factor(tData$STATE)
names(tData)[1] <- "State"EVTYPE
EVTYPE is a factor variable with 985 levels (488 in the tidy data set) while according to the Storm Data Documentation there are only 48 Storm Data Events.
Before I map the “wrong” event types to the correct ones based on text similarity I will find the most harmful events per category and see if they have valid descriptions.
sumD <- ddply(.data = tData, .(EVTYPE), .fun = summarise, FAT = sum(FATALITIES),
INJ = sum(INJURIES), PROP = sum(PROPDMG), CROP = sum(CROPDMG))
head(sumD[order(sumD$FAT, decreasing = T),],10)## EVTYPE FAT INJ PROP CROP
## 407 TORNADO 5633 91346 3212258.2 100018.52
## 61 EXCESSIVE HEAT 1903 6525 1460.0 494.40
## 73 FLASH FLOOD 978 1777 1420124.6 179200.46
## 151 HEAT 937 2100 298.5 662.70
## 258 LIGHTNING 816 5230 603351.8 3580.61
## 423 TSTM WIND 504 6957 1335965.6 109202.60
## 86 FLOOD 470 6789 899938.5 168037.88
## 306 RIP CURRENT 368 232 1.0 0.00
## 200 HIGH WIND 248 1137 324731.6 17283.21
## 11 AVALANCHE 224 170 1623.9 0.00head(sumD[order(sumD$INJ, decreasing = T),],10)## EVTYPE FAT INJ PROP CROP
## 407 TORNADO 5633 91346 3212258.16 100018.52
## 423 TSTM WIND 504 6957 1335965.61 109202.60
## 86 FLOOD 470 6789 899938.48 168037.88
## 61 EXCESSIVE HEAT 1903 6525 1460.00 494.40
## 258 LIGHTNING 816 5230 603351.78 3580.61
## 151 HEAT 937 2100 298.50 662.70
## 238 ICE STORM 89 1975 66000.67 1688.95
## 73 FLASH FLOOD 978 1777 1420124.59 179200.46
## 364 THUNDERSTORM WIND 133 1488 876844.17 66791.45
## 134 HAIL 15 1361 688693.38 579596.28head(sumD[order(sumD$PROP, decreasing = T),],10)## EVTYPE FAT INJ PROP CROP
## 407 TORNADO 5633 91346 3212258.2 100018.52
## 73 FLASH FLOOD 978 1777 1420124.6 179200.46
## 423 TSTM WIND 504 6957 1335965.6 109202.60
## 86 FLOOD 470 6789 899938.5 168037.88
## 364 THUNDERSTORM WIND 133 1488 876844.2 66791.45
## 134 HAIL 15 1361 688693.4 579596.28
## 258 LIGHTNING 816 5230 603351.8 3580.61
## 381 THUNDERSTORM WINDS 64 908 446293.2 18684.93
## 200 HIGH WIND 248 1137 324731.6 17283.21
## 481 WINTER STORM 206 1321 132720.6 1978.99head(sumD[order(sumD$CROP, decreasing = T),],10)## EVTYPE FAT INJ PROP CROP
## 134 HAIL 15 1361 688693.38 579596.28
## 73 FLASH FLOOD 978 1777 1420124.59 179200.46
## 86 FLOOD 470 6789 899938.48 168037.88
## 423 TSTM WIND 504 6957 1335965.61 109202.60
## 407 TORNADO 5633 91346 3212258.16 100018.52
## 364 THUNDERSTORM WIND 133 1488 876844.17 66791.45
## 49 DROUGHT 0 4 4099.05 33898.62
## 381 THUNDERSTORM WINDS 64 908 446293.18 18684.93
## 200 HIGH WIND 248 1137 324731.56 17283.21
## 159 HEAVY RAIN 98 251 50842.14 11122.80All event types are valid except “TSTM WIND”.
Replace TSTM with Thunderstorm. Then following the technique that I used with the States map the rest of event types to the correct ones. To do that I first create a vector with the valid events as they appear in the Storm Data Documentation.
tData$EVTYPE <- str_replace_all(tData$EVTYPE, "TSTM", "Thunderstorm")
# Creat a vector with the correct names
stormD <- 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")
#Replacing evtypes with their correct similar ones
tData$EVTYPE <- sapply(tData$EVTYPE,
function(x) stormD[
which.min(levenshteinDist(tolower(x),
tolower(stormD)))])
# Convert to factor and Rename
tData$EVTYPE <- as.factor(tData$EVTYPE)
names(tData)[2] <- "EvType"FATALITIES, INJURIES, PROPDMG, CROPDMG
Rename.
names(tData)[3:6] <- c("Fatalities", "Injuries", "PropDmg", "CropDmg")Reevaluate totals.
sumD2 <- ddply(.data = tData, .(EvType), .fun = summarise,
tFat = sum(Fatalities),
tInj = sum(Injuries),
tPropD = sum(PropDmg),
tCropD = sum(CropDmg))Importance of first ten event types per category.
topFat <- head(sumD2[order(sumD2$tFat, decreasing = T),],10)
percent(sum(topFat$tFat)/sum(tData$Fatalities))## [1] "86.2%"topInj <- head(sumD2[order(sumD2$tInj, decreasing = T),],10)
percent(sum(topInj$tInj)/sum(tData$Injuries))## [1] "92.7%"topProp <- head(sumD2[order(sumD2$tPropD, decreasing = T),],10)
percent(sum(topProp$tPropD)/sum(tData$PropDmg))## [1] "95.6%"topCrop <- head(sumD2[order(sumD2$tCropD, decreasing = T),],10)
percent(sum(topCrop$tCropD)/sum(tData$CropDmg))## [1] "96.6%"It’s obvious that there is no reason to care about all the rest event types.
Plots
plot1 <- ggplot(data=topFat,
aes(x=EvType, y=tFat)) +
geom_bar(stat="identity", fill = "darkblue", colour = "white") +
ylab("#Fatalities") +xlab("Events") +
scale_y_continuous(labels=comma) +
theme(legend.position="none",axis.text.x = element_text(angle = 45, hjust = 1))
plot2 <- ggplot(data=topInj,
aes(x=EvType, y=tInj/1000)) +
geom_bar(stat="identity", fill = "darkgreen", colour = "white")+
ylab("#Injuries(in K$)") +xlab("Events")+
theme(legend.position="none",axis.text.x = element_text(angle = 45, hjust = 1))
plot3 <- ggplot(data=topProp,
aes(x=EvType, y=tPropD/1e+06)) +
geom_bar(stat="identity", fill = "maroon", colour = "white")+
ylab("Property Damage(in M$)") +xlab("Events")+
theme(legend.position="none",axis.text.x = element_text(angle = 45, hjust = 1))
plot4 <- ggplot(data=topCrop,
aes(x=EvType, y=tCropD/1000)) +
geom_bar(stat="identity", fill = "orange", colour = "white")+
ylab("Crop Damage(in K$)") +xlab("Events")+
theme(legend.position="none",axis.text.x = element_text(angle = 45, hjust = 1))
grid.arrange(plot1, plot2, plot3, plot4, nrow = 2, ncol = 2, top = "Most Harmful Events in US (1950 - 2011)")
From the above plots it is obvious that torando is the event type that causes most fatalities, injuries and property damages whereas hail is responsible for most crop damages.
Find the 5 states with the biggest problems caused by the worst event types per category (tornado for fatalities, injuries, property damages, hail for crop damages) and plot.
tornD <- tData[tData$EvType == "Tornado",]
tornFA <- with(tornD, tapply(Fatalities, State, sum))
tornFS <- head(sort(tornFA, decreasing = T), 5)
tornFDF <- data.frame(State = names(tornFS), Fatalities = tornFS[1:5])
pl1 <- ggplot(data = tornFDF, aes(x = State, y = Fatalities)) +
geom_bar(stat="identity", fill = "lightblue", colour = "darkgrey") +
ylab("Fatalities \n from Tornadoes") +xlab("States") +
theme(legend.position="none")
tornIA <- with(tornD , tapply(Injuries, State, sum))
tornIS <- head(sort(tornIA, decreasing = T), 5)
tornIDF <- data.frame(State = names(tornIS), Injuries = tornIS[1:5])
pl2 <- ggplot(data = tornIDF, aes(x = State, y = Injuries)) +
geom_bar(stat="identity", fill = "lightgreen", colour = "darkgrey") +
ylab("Injuries \n from Tornadoes") +xlab("States") +
theme(legend.position="none")
tornPA <- with(tornD , tapply(PropDmg, State, sum))
tornPS <- head(sort(tornPA, decreasing = T), 5)
tornPDF <- data.frame(State = names(tornPS), PropDmg = tornPS[1:5])
pl3 <- ggplot(data = tornPDF, aes(x = State, y = PropDmg/1000)) +
geom_bar(stat="identity", fill = "#F8766D", colour = "darkgrey") +
ylab("Prop. Damage from \n Tornadoes(in thousands $)") +xlab("States") +
theme(legend.position="none")
hailD <- tData[tData$EvType == "Hail",]
tornCA <- with(hailD, tapply(CropDmg, State, sum))
tornCS <- head(sort(tornCA, decreasing = T), 5)
tornCDF <- data.frame(State = names(tornCS), CorpDmg = tornCS[1:5])
pl4 <- ggplot(data = tornCDF, aes(x = State, y = CorpDmg/1000)) +
geom_bar(stat="identity", fill = "khaki", colour = "darkgrey") +
ylab("Corp Dmg from \n Hail (in thousands $)") +xlab("States") +
theme(legend.position="none")
grid.arrange(pl1, pl2, pl3, pl4, nrow = 2, ncol = 2, top = "States with most Public Health/Economic Problems \n Caused from Severe Weather Events")
The above plot shows the 5 states with most fatalities / injuries / property damages and crop damages caused by the worst weather event by category.
Results
The above analysis gave the following results:
The weather event that causes most fatalities, injuries and property damage is by far the tornado.
The weather event that causes most crop damages is the hail.
Alabama and Texas are the States with most fatalities and injuries, caused by severe weather events.
Texas and Nebraska are the States with most property and crop damages, caused by severe weather events.