This is an analysis of the National Oceanic and Atmospheric Administration (NOAA) Storm data to investigate the impact of weather events on public health (as measured by average number of injuries and/or fatalaties) and on the economy (as measured by the average amount of crop damage and/or property damage). This allows the goverment to allocate resources and effort effectively.
The aim of the analysis is to investigate the following questions.
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?
The was downloaded as a zip file from the NOAA Storm Data. The documentation for the data can be found in the Storm Data Documentation. The data was read into R from the zipepd comma separated file.
library(tidyverse)
if (!file.exists('repdata%2Fdata%2FStormData.csv.bz2')) {
## Download the zip file
download.file("https://d396qusza40orc.cloudfront.net/repdata%2Fdata%2FStormData.csv.bz2",
'repdata%2Fdata%2FStormData.csv.bz2')
}
input_data <- read.csv('repdata%2Fdata%2FStormData.csv.bz2')
After reading in the data we investigate the first few files and attributes of the dataset.
dim(input_data)
## [1] 902297 37
str(input_data)
## '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 ...
head(input_data)
## STATE__ BGN_DATE BGN_TIME TIME_ZONE COUNTY COUNTYNAME STATE
## 1 1 4/18/1950 0:00:00 0130 CST 97 MOBILE AL
## 2 1 4/18/1950 0:00:00 0145 CST 3 BALDWIN AL
## 3 1 2/20/1951 0:00:00 1600 CST 57 FAYETTE AL
## 4 1 6/8/1951 0:00:00 0900 CST 89 MADISON AL
## 5 1 11/15/1951 0:00:00 1500 CST 43 CULLMAN AL
## 6 1 11/15/1951 0:00:00 2000 CST 77 LAUDERDALE AL
## EVTYPE BGN_RANGE BGN_AZI BGN_LOCATI END_DATE END_TIME COUNTY_END
## 1 TORNADO 0 0
## 2 TORNADO 0 0
## 3 TORNADO 0 0
## 4 TORNADO 0 0
## 5 TORNADO 0 0
## 6 TORNADO 0 0
## COUNTYENDN END_RANGE END_AZI END_LOCATI LENGTH WIDTH F MAG FATALITIES
## 1 NA 0 14.0 100 3 0 0
## 2 NA 0 2.0 150 2 0 0
## 3 NA 0 0.1 123 2 0 0
## 4 NA 0 0.0 100 2 0 0
## 5 NA 0 0.0 150 2 0 0
## 6 NA 0 1.5 177 2 0 0
## INJURIES PROPDMG PROPDMGEXP CROPDMG CROPDMGEXP WFO STATEOFFIC ZONENAMES
## 1 15 25.0 K 0
## 2 0 2.5 K 0
## 3 2 25.0 K 0
## 4 2 2.5 K 0
## 5 2 2.5 K 0
## 6 6 2.5 K 0
## LATITUDE LONGITUDE LATITUDE_E LONGITUDE_ REMARKS REFNUM
## 1 3040 8812 3051 8806 1
## 2 3042 8755 0 0 2
## 3 3340 8742 0 0 3
## 4 3458 8626 0 0 4
## 5 3412 8642 0 0 5
## 6 3450 8748 0 0 6
Next we will check if there any missing values in the data especially in the INJURIES and FATALITIES fields
sum(is.na(input_data$INJURIES))
## [1] 0
sum(is.na(input_data$INJURIES))
## [1] 0
There are no missing values in the data for injuries and fatalities.
The range of the data captured is from following dates:
library(lubridate)
range(mdy_hms(as.character(input_data$BGN_DATE)))
## [1] "1950-01-03 UTC" "2011-11-30 UTC"
The following questions were investigated and attempted to be addressed.
In order to answer this question we will look at the total number of injuries and fatalities caused by each event. For clarity purposes we will only count event types which have more than 10 average injuries and more than 1 average fatalities.
We will group the data by event type (variable = EVTYPE) and calculate the average number of injuries and pick only events with avg. injuries more than 10 and sort the data by the avg number of injuries and plot the data in a bar plot.
g_data <- input_data %>% group_by(EVTYPE) %>% summarise(inj = mean(INJURIES, na.rm = TRUE)) %>% filter(inj > 10) %>% arrange(inj) %>% mutate(EVTYPE = factor(EVTYPE, levels = unique(EVTYPE)))
g <- ggplot(data = g_data)
g <- g + geom_bar(mapping = aes(x = EVTYPE, y = inj), stat = "identity") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + ylab("Average number of injuries") + xlab("Event Type") + labs(title = "Average number of injuries")
print(g)
Heatwaves caused the most injuries on average followed by Tropical storm and Wildfires.
Next we will do a similar investigation of the number of fatalities on average by event type filtering only events with more than 1 fatality on average.
g_data2 <- input_data %>% group_by(EVTYPE) %>% summarise(fatal = mean(FATALITIES, na.rm = TRUE)) %>% filter(fatal > 1) %>% arrange(fatal) %>% mutate(EVTYPE = factor(EVTYPE, levels = unique(EVTYPE)))
g2 <- ggplot(data = g_data2)
g2 <- g2 + geom_bar(mapping = aes(x = EVTYPE, y = fatal), stat = "identity") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + ylab("Average number of Fatalities") + xlab("Event Type") + labs(title = "More an 1 avg. num of fatalities")
print(g2)
Tornadoes/T storms/Hail caused the most number of fatalities on average followed by Cold and snow.
In order to estimate the economic consequences of a weather event we will look at 2 variables, CROP damage and Property damage. First we will look at crop damage, showing only the events where the average damage was more than 10 million dollars. Transformation of the data is done for the dollar amount, replacing the levels in CROPDMGEXP with appropriate multiplication factor (e.g. 1000000 for “M”)
levels.CROPDMG <- levels(input_data$CROPDMGEXP)
levels.CROPDMG.rep <- c(0,0,0,0,1000000000,1000,1000,1000000,1000000)
crop_data <- input_data %>% select(EVTYPE, CROPDMG, CROPDMGEXP) %>%
filter(CROPDMG > 0) %>% mutate(CROPMDG_rep = plyr::mapvalues(CROPDMGEXP, from = levels.CROPDMG, to = levels.CROPDMG.rep)) %>% mutate(CROPDMG_USD = CROPDMG * as.numeric(paste(CROPMDG_rep)))
crop_damage <- crop_data %>% group_by(EVTYPE) %>% summarise(CROPDMG_mean = mean(CROPDMG_USD, na.rm = TRUE)) %>% arrange(CROPDMG_mean) %>% mutate(EVTYPE = factor(EVTYPE, levels = unique(EVTYPE))) %>% filter(CROPDMG_mean > 10000000)
g3 <- ggplot(data = crop_damage)
g3 <- g3 + geom_bar(mapping = aes(x = EVTYPE, y = CROPDMG_mean), stat = "identity") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + ylab("Average Crop damage") + xlab("Event Type") + labs(title = "Average Crop damage (more than 10 million USD)")
print(g3)
River flood seems to have caused the most crop damage on average.
Next we will investigate the Property damage focussing only on damage more than 10 million dollars on average. Transformation of the data is done for the dollar amount, replacing the levels in PROPDMGEXP with appropriate multiplication factor (e.g. 1000000 for “M”)
prop_data <- input_data %>% select(EVTYPE, PROPDMG, PROPDMGEXP) %>%
filter(PROPDMG > 0) %>% filter(!(PROPDMGEXP %in% c("","-","?","+","0","1","2","3","4", "5","6","7","8")))
prop_levels <- c("B","h","H","K","m","M")
prop_levels_rep <- c(1000000000,100,100,1000,1000000,1000000)
prop_data <- prop_data %>% mutate(PROPDMG_REP = plyr::mapvalues(PROPDMGEXP, from = prop_levels, to = prop_levels_rep)) %>% mutate(PROPDMG_USD = PROPDMG * as.numeric(paste(PROPDMG_REP)))
prop_damage <- prop_data %>% group_by(EVTYPE) %>% summarise(PROPDMG_mean = mean(PROPDMG_USD, na.rm = TRUE)) %>% arrange(PROPDMG_mean) %>% mutate(EVTYPE = factor(EVTYPE, levels = unique(EVTYPE))) %>% filter(PROPDMG_mean > 10000000)
g4 <- ggplot(data = prop_damage)
g4 <- g4 + geom_bar(mapping = aes(x = EVTYPE, y = PROPDMG_mean), stat = "identity") + theme(axis.text.x = element_text(angle = 90, hjust = 1)) + ylab("Average Property damage") + xlab("Event Type") + labs(title = "Average Property damage (more than 10 million USD)")
print(g4)
Heavy rains/Severe weather seems to have caused the most property damage on average followed by Tornadoes/T Storms/Hail.