Storm Data Analysis
Antonio Castro
2025-12-07
Storm Data Report
knitr::opts_chunk$set(echo = TRUE)
## to unpack bzip2 algorithm
# install.packages("R.utils")
library(R.utils)
## graphics
# install.packages("ggplot2")
library(ggplot2)Synopsis
The present work is a Coursera course - Reproducible Research project. The main objective is to analise Storm Data, from the US National Weather Service, to answer two specific questions:
Which types of events are most harmful with respect to population health?
Which types of events have the greatest economic consequences?
Data Processing
The data is available here: Storm Data. It comes as a comma-separated-value file compressed via the bzip2 algorithm to reduce its size.
Unpack and read the file:
bunzip2(filename="repdata_data_StormData.csv.bz2", destname="repdata_data_StormData.csv", remove=FALSE, overwrite=TRUE)
dataset <- read.csv("repdata_data_StormData.csv", header=TRUE)As we have just readen the file, we have to try to understand its contents. First of all we need to know each attributes it has, and the type of each one:
str(dataset)## '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 ...Wow! There are a lot of observations and attributes. There is some
documentation available at:
Storm
Data Documentation, and at the National
Climatic Data Center Storm Events - FAQ.
From the documentation and the dataset attribute names it is possible to understand which attributes are rellevant for our analysis.
To answer both questions we need to find out which attribute defines the type of the event. The EVTYPE should be the one.
For the first question, it’s necessary to understand which attributes are related to population health. These are: INJURIES and FATALITIES.
And for the second question PROPDMG, PROPDMGEXP, CROPDMG, CROPDMGEXP are the attributes related to economics. PROPDMG, and CROPDMG are the values rounded to three significant digits, while PROPDMGEXP, and CROPDMGEXP use characters to signify magnitude include “K” for thousands, “M” for millions, and “B” for billions.
Goind direct to the point, as we will not need the other attributes, we can keep only the needed attributes:
dataset <- dataset[, c("EVTYPE", "INJURIES", "FATALITIES", "PROPDMG", "PROPDMGEXP", "CROPDMG", "CROPDMGEXP")]Let’s tke a look of the first ten lines of the dataset, to get understand better the attributes.
head(dataset, n=10)## EVTYPE INJURIES FATALITIES PROPDMG PROPDMGEXP CROPDMG CROPDMGEXP
## 1 TORNADO 15 0 25.0 K 0
## 2 TORNADO 0 0 2.5 K 0
## 3 TORNADO 2 0 25.0 K 0
## 4 TORNADO 2 0 2.5 K 0
## 5 TORNADO 2 0 2.5 K 0
## 6 TORNADO 6 0 2.5 K 0
## 7 TORNADO 1 0 2.5 K 0
## 8 TORNADO 0 0 2.5 K 0
## 9 TORNADO 14 1 25.0 K 0
## 10 TORNADO 0 0 25.0 K 0We can go further, and create a specific attribute for the information about first question that sums INJURIES and FATALITIES, and remove these two attributes:
dataset$health_impact <- dataset$INJURIES + dataset$FATALITIES
dataset$INJURIES <- NULL
dataset$FATALITIES <- NULLSo we have a new design:
head(dataset, n=10)## EVTYPE PROPDMG PROPDMGEXP CROPDMG CROPDMGEXP health_impact
## 1 TORNADO 25.0 K 0 15
## 2 TORNADO 2.5 K 0 0
## 3 TORNADO 25.0 K 0 2
## 4 TORNADO 2.5 K 0 2
## 5 TORNADO 2.5 K 0 2
## 6 TORNADO 2.5 K 0 6
## 7 TORNADO 2.5 K 0 1
## 8 TORNADO 2.5 K 0 0
## 9 TORNADO 25.0 K 0 15
## 10 TORNADO 25.0 K 0 0For the economic values we can do this in two steps: transfor PROPDMG, and CROPDMG in the full value, depending ont its magnitude, that we can take from PROPDMGEXP, and CROPDMGEXP. Before we do this, we have to check if PROPDMGEXP, and CROPDMGEXP include only the characters described on the documentation: “K” for thousands, “M” for millions, and “B” for billions.
print("Characters for *PROPDMGEXP* :")## [1] "Characters for *PROPDMGEXP* :"print(unique(dataset$PROPDMGEXP))## [1] "K" "M" "" "B" "m" "+" "0" "5" "6" "?" "4" "2" "3" "h" "7" "H" "-" "1" "8"print("Characters for *CROPDMGEXP* :")## [1] "Characters for *CROPDMGEXP* :"print(unique(dataset$CROPDMGEXP))## [1] "" "M" "K" "m" "B" "?" "0" "k" "2"Oh, there are some values that are not exactly what we needed. It is necessary to take a decision here. We will take the upper case of both this attributes, so “M” and “m” will be the same. Than take its magnitude to multiply for the value. Finally, we will sum these values in a single attribute, and remove the original values
dataset$PROPDMGEXP <- toupper(dataset$PROPDMGEXP)
dataset[dataset$PROPDMGEXP=="K", "PROPDMG"] <- dataset[dataset$PROPDMGEXP=="K", "PROPDMG"] * 100
dataset[dataset$PROPDMGEXP=="M", "PROPDMG"] <- dataset[dataset$PROPDMGEXP=="M", "PROPDMG"] * 1000
dataset[dataset$PROPDMGEXP=="B", "PROPDMG"] <- dataset[dataset$PROPDMGEXP=="B", "PROPDMG"] * 1000000
dataset$CROPDMGEXP <- toupper(dataset$CROPDMGEXP)
dataset[dataset$CROPDMGEXP=="K", "CROPDMG"] <- dataset[dataset$CROPDMGEXP=="K", "CROPDMG"] * 100
dataset[dataset$CROPDMGEXP=="M", "CROPDMG"] <- dataset[dataset$CROPDMGEXP=="M", "CROPDMG"] * 1000
dataset[dataset$CROPDMGEXP=="B", "CROPDMG"] <- dataset[dataset$CROPDMGEXP=="B", "CROPDMG"] * 1000000
dataset$economic_impact <- dataset$PROPDMG + dataset$CROPDMG
dataset$PROPDMG <- NULL
dataset$PROPDMGEXP <- NULL
dataset$CROPDMG <- NULL
dataset$CROPDMGEXP <- NULLLet’s see the changes:
So we have a new design:
head(dataset, n=10)## EVTYPE health_impact economic_impact
## 1 TORNADO 15 2500
## 2 TORNADO 0 250
## 3 TORNADO 2 2500
## 4 TORNADO 2 250
## 5 TORNADO 2 250
## 6 TORNADO 6 250
## 7 TORNADO 1 250
## 8 TORNADO 0 250
## 9 TORNADO 15 2500
## 10 TORNADO 0 2500To answer the questions we have to aggregate the data by EVTYPE, for each context, health_impact, and economic_impact:
evtype_health <- aggregate(dataset$health_impact, list(dataset$EVTYPE), FUN=sum)
colnames(evtype_health) <- c("EVTYPE", "health_impact")
evtype_economic <- aggregate(dataset$economic_impact, list(dataset$EVTYPE), FUN=sum)
colnames(evtype_economic) <- c("EVTYPE", "economic_impact")Let’s see the results for health:
head(evtype_health, n=10)## EVTYPE health_impact
## 1 HIGH SURF ADVISORY 0
## 2 COASTAL FLOOD 0
## 3 FLASH FLOOD 0
## 4 LIGHTNING 0
## 5 TSTM WIND 0
## 6 TSTM WIND (G45) 0
## 7 WATERSPOUT 0
## 8 WIND 0
## 9 ? 0
## 10 ABNORMAL WARMTH 0And for economics:
head(evtype_economic, n=10)## EVTYPE economic_impact
## 1 HIGH SURF ADVISORY 20000
## 2 COASTAL FLOOD 0
## 3 FLASH FLOOD 5000
## 4 LIGHTNING 0
## 5 TSTM WIND 18000
## 6 TSTM WIND (G45) 800
## 7 WATERSPOUT 0
## 8 WIND 0
## 9 ? 500
## 10 ABNORMAL WARMTH 0Oops, we have a question mark and some values with zero total. Let’s remove it:
evtype_health <- evtype_health[!(evtype_health$health_impact==0),]
evtype_health <- evtype_health[!(evtype_health$EVTYPE=="?"),]
evtype_economic <- evtype_economic[!(evtype_economic$economic_impact==0),]
evtype_economic <- evtype_economic[!(evtype_economic$EVTYPE=="?"),]So, all the necessary chages have been made.
Results
The answer to the questions:
- Which types of events are most harmful to population health?
Let’s see the 5 most impacting event types for people health:
data <- head(evtype_health[order(evtype_health[, "health_impact"], decreasing=TRUE), ], n=5)
ggplot(data, aes(EVTYPE, health_impact)) +
geom_bar(stat = "identity") +
labs(title = "Health Impact by Event Type",
x = "Event Type",
y = "People") +
theme_minimal()
Which one is the worst?
maximun <- max(evtype_health$health_impact)
evtype_health[evtype_health$health_impact==maximun, "EVTYPE"]## [1] "TORNADO"- Which types of events have the greatest economic consequences?
Let’s see the 5 most impacting event types for economics:
data <- head(evtype_economic[order(evtype_economic[, "economic_impact"], decreasing=TRUE), ], n=5)
ggplot(data, aes(EVTYPE, economic_impact)) +
geom_bar(stat = "identity") +
labs(title = "Health Impact by Event Type",
x = "Event Type",
y = "US $") +
theme_minimal()
Which one is the worst?
maximun <- max(evtype_economic$economic_impact)
evtype_economic[evtype_economic$economic_impact==maximun, "EVTYPE"]## [1] "TORNADO"