The aim of this study is to analyse Storm Data across the United States and identify the most harmful types of events with respect to population health and determine the events that had the greatest economical consequences. An end-to-end analysis has been performed from the data source via all transformation steps until the results and conclusions (with all the R codes accordingly).
The results say that in terms of health impacts (fatalities and injuries combined), the TOP 10 events accounts for the 88% (137 177 fatalities and injuries out of the total 155 673 cases) and the 3 most harmful events namely are
Regarding the economical consequences, the TOP 10 accounts for the 86% of the total caused damages (for property and crop damages aggregated, 408.6 Billion $USD out of the 477.3 Billion $USD). The 3 most harmful events namely are:
The following packages are loaded for the analysis.
# Load packages
library(readr)
library(readxl)
library(dplyr)
library(tidyverse)
library(ggplot2)Download the data to folder inst/extdata.
# Download and unzip the file ----
fileUrl <- "https://d396qusza40orc.cloudfront.net/repdata%2Fdata%2FStormData.csv.bz2"
fileName <- "storm_data.csv"
fileFolder <- "data"
dataDest <- "inst/extdata/"
if(!file.exists(paste0("inst/extdata/", fileName))){
download.file(url = fileUrl,
destfile = paste0("inst/extdata/", fileName),
method = "curl")
}Read the data set to a tibble called stormdata.
# Read the data ----
path <- paste0(getwd(), "/", dataDest)
stormdata <- read_csv( file.path(path, fileName)) %>% as_tibble()The Storm Data is an official publication of the National Oceanic and Atmospheric Administration (NOAA) which documents the occurrence of storms and other significant weather phenomena having sufficient intensity to cause loss of life, injuries, significant property damage, and/or disruption to commerce. (Source: Storm Data Documentation)
THe report is focusing on the following questions:
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?
For the analysis it’s worth to check the structure of the data set.
summary(stormdata)## STATE__ BGN_DATE BGN_TIME TIME_ZONE
## Min. : 1.0 Length:902297 Length:902297 Length:902297
## 1st Qu.:19.0 Class :character Class :character Class :character
## Median :30.0 Mode :character Mode :character Mode :character
## Mean :31.2
## 3rd Qu.:45.0
## Max. :95.0
##
## COUNTY COUNTYNAME STATE EVTYPE
## Min. : 0.0 Length:902297 Length:902297 Length:902297
## 1st Qu.: 31.0 Class :character Class :character Class :character
## Median : 75.0 Mode :character Mode :character Mode :character
## Mean :100.6
## 3rd Qu.:131.0
## Max. :873.0
##
## BGN_RANGE BGN_AZI BGN_LOCATI END_DATE
## Min. : 0.000 Length:902297 Length:902297 Length:902297
## 1st Qu.: 0.000 Class :character Class :character Class :character
## Median : 0.000 Mode :character Mode :character Mode :character
## Mean : 1.484
## 3rd Qu.: 1.000
## Max. :3749.000
##
## END_TIME COUNTY_END COUNTYENDN END_RANGE
## Length:902297 Min. :0 Mode:logical Min. : 0.0000
## Class :character 1st Qu.:0 NA's:902297 1st Qu.: 0.0000
## Mode :character Median :0 Median : 0.0000
## Mean :0 Mean : 0.9862
## 3rd Qu.:0 3rd Qu.: 0.0000
## Max. :0 Max. :925.0000
##
## END_AZI END_LOCATI LENGTH WIDTH
## Length:902297 Length:902297 Min. : 0.0000 Min. : 0.000
## Class :character Class :character 1st Qu.: 0.0000 1st Qu.: 0.000
## Mode :character Mode :character Median : 0.0000 Median : 0.000
## Mean : 0.2301 Mean : 7.503
## 3rd Qu.: 0.0000 3rd Qu.: 0.000
## Max. :2315.0000 Max. :4400.000
##
## F MAG FATALITIES INJURIES
## Min. :0.0 Min. : 0.0 Min. : 0.0000 Min. : 0.0000
## 1st Qu.:0.0 1st Qu.: 0.0 1st Qu.: 0.0000 1st Qu.: 0.0000
## Median :1.0 Median : 50.0 Median : 0.0000 Median : 0.0000
## Mean :0.9 Mean : 46.9 Mean : 0.0168 Mean : 0.1557
## 3rd Qu.:1.0 3rd Qu.: 75.0 3rd Qu.: 0.0000 3rd Qu.: 0.0000
## Max. :5.0 Max. :22000.0 Max. :583.0000 Max. :1700.0000
## NA's :843563
## PROPDMG PROPDMGEXP CROPDMG CROPDMGEXP
## Min. : 0.00 Length:902297 Min. : 0.000 Length:902297
## 1st Qu.: 0.00 Class :character 1st Qu.: 0.000 Class :character
## Median : 0.00 Mode :character Median : 0.000 Mode :character
## Mean : 12.06 Mean : 1.527
## 3rd Qu.: 0.50 3rd Qu.: 0.000
## Max. :5000.00 Max. :990.000
##
## WFO STATEOFFIC ZONENAMES LATITUDE
## Length:902297 Length:902297 Length:902297 Min. : 0
## Class :character Class :character Class :character 1st Qu.:2802
## Mode :character Mode :character Mode :character Median :3540
## Mean :2875
## 3rd Qu.:4019
## Max. :9706
## NA's :47
## LONGITUDE LATITUDE_E LONGITUDE_ REMARKS
## Min. :-14451 Min. : 0 Min. :-14455 Length:902297
## 1st Qu.: 7247 1st Qu.: 0 1st Qu.: 0 Class :character
## Median : 8707 Median : 0 Median : 0 Mode :character
## Mean : 6940 Mean :1452 Mean : 3509
## 3rd Qu.: 9605 3rd Qu.:3549 3rd Qu.: 8735
## Max. : 17124 Max. :9706 Max. :106220
## NA's :40
## REFNUM
## Min. : 1
## 1st Qu.:225575
## Median :451149
## Mean :451149
## 3rd Qu.:676723
## Max. :902297
## After some checks in the aforementioned Storm Data Documentation, it can be concluded that at first glance for the purpose of the report, the most important data fields are:
EVTYPE: this field contains the different storm event
types. In the data set there are 977 different events.FATALITIES: number of the weather-related fatalities
(direct and indirect)INJURIES: number of the weather-related injuries
(direct and indirect)PROPDMG: nominal value of the property damagePROPDMGEXP: category of the property damage. There are
19 different categories.CROPDMG: nominal value of the crop damageCROPDMGEXP: category of the crop damage. There are 9
different categories.Based on Section 2.1 in the Storm Data Documentation and based on the
instruction of the we can assume that the main categories are correctly
entered apart from immaterial discrepancies, hence using the data field
without any modification doesn’t lead to significant differences. Even
though I reckon further cleansing would be a nice to have exercise, it
is out of the scope of this analysis. Variable EVTYPE will
be used for the categorization.
Based on the Storm Data Documentation (Section 2.6), the
FATALITIES and INJURIES (direct and indirect)
are stored in the data fields and actually no extra manipulation is
needed to get the exact amounts from them, they are measured on the same
scale. To get all the harmful health impacts on the population, we can
sum up these variables and check them in total by categories.
Regarding the damages, let’s refer again to the documentation to get
some explanation about variables PROPDMG,
PROPDMGEXP, CROPDMG and
CROPDMGEXP. Estimates should be in the form of US
Dollar values and rounded to three significant digits, followed by the
magnitude of the value (i.e., 1.55 Billion $USD for $1,550,000,000).
Values used to signify magnitude include: Thousand $USD, Million $USD,
and Billion $USD. If additional precision is available, it may be
provided in the narrative part of the entry. When damage is due to more
than one element of the storm, indicate, when possible, the amount of
damage caused by each element. If the dollar amount of damage is
unknown, or not available, check the “no information available”
box. (Section 2.7)
Let’s have a closer look at on the different measures that are
applied in the “...EXP” variables.
# List of different damage identifiers
unique(stormdata$PROPDMGEXP) %>% sort()## [1] "-" "?" "+" "0" "1" "2" "3" "4" "5" "6" "7" "8" "B" "h" "H" "K" "m" "M"unique(stormdata$CROPDMGEXP) %>% sort()## [1] "?" "0" "2" "B" "k" "K" "m" "M"We see that there are expected magnitude identifiers. In
PROPDMGEXP these are: “B” “h” “H” “K” “m” “M”, meanwhile in
CROPDMGEXP the expected ones are: “B” “k” “K” “m” “M”. For
these we can use the multipliers as indicated by the letters. There are
number indicators (“0” “1” “2” “3” “4” “5” “6” “7” “8”) as well. For
these based on the documentation we can suppose that these are
multipliers of the power of 10 (“followed by the magnitude of the
value”). Regarding the other values (“-” “?” “+”) I will suppose
that the damages are entered as $USD, so the multiplier is 1. Lastly I
consider NA as 0.
Hence, the following conversions will be used for the magnitude identifiers:
| Value | Multiplier |
|---|---|
| NA | \(0\) |
| “-” “?” “+” | \(1\) |
| \(n=\) “0” “1” “2” “3” “4” “5” “6” “7” “8” | \(10^n\) |
| “h” and “H” | \(100\) |
| “k” and “K” | \(1000\) |
| “m” and “M” | \(1000000\) |
| “B” | \(1000000000\) |
stormdata_mod <- stormdata %>%
select(EVTYPE, FATALITIES, INJURIES, contains("DMG")) %>%
mutate(
total_health = FATALITIES + INJURIES,
pdmg_mult = case_when(
PROPDMGEXP %in% c("-", "?", "+") ~ 1,
PROPDMGEXP %in% c("0", "1", "2", "3", "4", "5", "6", "7", "8") ~ 10^as.integer(PROPDMGEXP),
PROPDMGEXP %in% c("h", "H") ~ 100,
PROPDMGEXP %in% c("k", "K") ~ 1000,
PROPDMGEXP %in% c("m", "M") ~ 1000000,
PROPDMGEXP %in% c("b", "B") ~ 1000000000,
TRUE ~ 0
),
cdmg_mult = case_when(
CROPDMGEXP %in% c("-", "?", "+") ~ 1,
CROPDMGEXP %in% c("0", "1", "2", "3", "4", "5", "6", "7", "8") ~ 10^as.integer(CROPDMGEXP),
CROPDMGEXP %in% c("h", "H") ~ 100,
CROPDMGEXP %in% c("k", "K") ~ 1000,
CROPDMGEXP %in% c("m", "M") ~ 1000000,
CROPDMGEXP %in% c("b", "B") ~ 1000000000,
TRUE ~ 0
),
pdmg_val = PROPDMG * pdmg_mult,
cdmg_val = CROPDMG * cdmg_mult,
totdmg_val = pdmg_val + cdmg_val
)
# Health related data
stormdata_mod %>% select(EVTYPE, FATALITIES, INJURIES, total_health) %>% print(n = 10)## # A tibble: 902,297 × 4
## EVTYPE FATALITIES INJURIES total_health
## <chr> <dbl> <dbl> <dbl>
## 1 TORNADO 0 15 15
## 2 TORNADO 0 0 0
## 3 TORNADO 0 2 2
## 4 TORNADO 0 2 2
## 5 TORNADO 0 2 2
## 6 TORNADO 0 6 6
## 7 TORNADO 0 1 1
## 8 TORNADO 0 0 0
## 9 TORNADO 1 14 15
## 10 TORNADO 0 0 0
## # ℹ 902,287 more rows# Damage related data
stormdata_mod %>% select(EVTYPE, contains("DMG"), contains("dmg") ) %>% print(n = 10)## # A tibble: 902,297 × 10
## EVTYPE PROPDMG PROPDMGEXP CROPDMG CROPDMGEXP pdmg_mult cdmg_mult pdmg_val
## <chr> <dbl> <chr> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 TORNADO 25 K 0 <NA> 1000 0 25000
## 2 TORNADO 2.5 K 0 <NA> 1000 0 2500
## 3 TORNADO 25 K 0 <NA> 1000 0 25000
## 4 TORNADO 2.5 K 0 <NA> 1000 0 2500
## 5 TORNADO 2.5 K 0 <NA> 1000 0 2500
## 6 TORNADO 2.5 K 0 <NA> 1000 0 2500
## 7 TORNADO 2.5 K 0 <NA> 1000 0 2500
## 8 TORNADO 2.5 K 0 <NA> 1000 0 2500
## 9 TORNADO 25 K 0 <NA> 1000 0 25000
## 10 TORNADO 25 K 0 <NA> 1000 0 25000
## # ℹ 902,287 more rows
## # ℹ 2 more variables: cdmg_val <dbl>, totdmg_val <dbl>With this data table we are ready to produce the results.
Let’s aggregating the result and check the most harmful events. First let’s create a table with the TOP 10 events. Then plot it to a bar chart where the categories are the events and the order is descending by the health impact.
# Get the health data needed for the summary
stormdata_mod_health <- stormdata_mod %>%
select(EVTYPE, total_health) %>%
group_by(EVTYPE) %>%
summarize(health_impact = sum(total_health)) %>%
arrange(desc(health_impact)) %>%
mutate(category = "Total Health") %>%
head(10)
stormdata_mod_health## # A tibble: 10 × 3
## EVTYPE health_impact category
## <chr> <dbl> <chr>
## 1 TORNADO 96979 Total Health
## 2 EXCESSIVE HEAT 8428 Total Health
## 3 TSTM WIND 7461 Total Health
## 4 FLOOD 7259 Total Health
## 5 LIGHTNING 6046 Total Health
## 6 HEAT 3037 Total Health
## 7 FLASH FLOOD 2755 Total Health
## 8 ICE STORM 2064 Total Health
## 9 THUNDERSTORM WIND 1621 Total Health
## 10 WINTER STORM 1527 Total Health# Reorder the EVTYPE based on health_impact
stormdata_mod_health$EVTYPE <- factor(stormdata_mod_health$EVTYPE, levels = unique(stormdata_mod_health$EVTYPE[order(stormdata_mod_health$health_impact, decreasing = TRUE)]))the TOP 10 most harmful events with respect to population health across the US are:
# Create a bar chart to represent the results
g1 <- ggplot(stormdata_mod_health, aes(x = EVTYPE, y = health_impact)) +
geom_bar(stat = "identity", fill = "#0099f9") +
geom_text(aes(label = health_impact), size = 4, vjust = -0.3, colour = "blue") +
ggtitle("Most harmful events for the population health") +
xlab("Storm event type") +
ylab("Harmful cases (fatalities + injuries)") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
g1
From the bar chart above, we can see that
# Get the damages for the summary
stormdata_mod_econ <- stormdata_mod %>%
select(EVTYPE, contains("DMG"), contains("dmg") ) %>%
group_by(EVTYPE) %>%
summarize(economical_impact = round(sum(totdmg_val) / 10^9, 1)) %>%
arrange(desc(economical_impact)) %>%
head(10)
stormdata_mod_econ## # A tibble: 10 × 2
## EVTYPE economical_impact
## <chr> <dbl>
## 1 FLOOD 150.
## 2 HURRICANE/TYPHOON 71.9
## 3 TORNADO 57.4
## 4 STORM SURGE 43.3
## 5 HAIL 18.8
## 6 FLASH FLOOD 18.2
## 7 DROUGHT 15
## 8 HURRICANE 14.6
## 9 RIVER FLOOD 10.1
## 10 ICE STORM 9# Reorder the EVTYPE based on economical_impact
stormdata_mod_econ$EVTYPE <- factor(stormdata_mod_econ$EVTYPE, levels = unique(stormdata_mod_econ$EVTYPE[order(stormdata_mod_econ$economical_impact, decreasing = TRUE)]))the TOP 10 most events that had the greatest economic consequenses in the US are:
# Create a chart that represents the results
g2 <- ggplot(stormdata_mod_econ, aes(x = EVTYPE, y = economical_impact)) +
geom_bar(stat = "identity", fill = "#CC6666") +
geom_text(aes(label = economical_impact), size = 4, vjust = -0.3, colour = "red") +
ggtitle("Most harmful events for the economy") +
xlab("Storm event type") +
ylab("Damage caused (property + crop, Billion $USD") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
g2
This part contains summary statistic calculations for the Synopsis.
# TOP10 health proportion and nominal values
round(sum(stormdata_mod_health$health_impact) / sum(stormdata_mod$total_health), 2)## [1] 0.88sum(stormdata_mod_health$health_impact)## [1] 137177sum(stormdata_mod$total_health)## [1] 155673round(as.numeric(stormdata_mod_health[1,2]) / sum(stormdata_mod$total_health), 2)## [1] 0.62round(as.numeric(stormdata_mod_health[2,2]) / sum(stormdata_mod$total_health), 2)## [1] 0.05round(as.numeric(stormdata_mod_health[3,2]) / sum(stormdata_mod$total_health), 2)## [1] 0.05# TOP10 economical
round(sum(stormdata_mod_econ$economical_impact)*10^9/sum(stormdata_mod$totdmg_val),2)## [1] 0.86stormdata_mod_econ$economical_impact## [1] 150.3 71.9 57.4 43.3 18.8 18.2 15.0 14.6 10.1 9.0round(stormdata_mod_econ$economical_impact[1]*10^9 / sum(stormdata_mod$totdmg_val), 2)## [1] 0.31round(stormdata_mod_econ$economical_impact[2]*10^9 / sum(stormdata_mod$totdmg_val), 2)## [1] 0.15round(stormdata_mod_econ$economical_impact[3]*10^9 / sum(stormdata_mod$totdmg_val), 2)## [1] 0.12