NOAA Storm database based Health/Economic effects analysis
Elena Chernousova
Sunday, July 26, 2015
This analysis is based on data provided U.S. National Oceanic and Atmospheric Administration’s (NOAA). The file database contains characteristics of major storms and weather events in the United States, including when and where they occur, as well as estimates of any fatalities, injuries, and property damage. Expected results are set of natural events affecting population health and economic effects.
Initialsetup
Loading libriries to proper data presentation:
options("scipen"=10)
if (!require(dplyr)) {
install.packages("dplyr")
library(dplyr)
}## Loading required package: dplyr
##
## Attaching package: 'dplyr'
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionif (!require(ggplot2)) {
install.packages("ggplot2")
library(ggplot2)
}## Loading required package: ggplot2Data download
Initial NOAA data are available for public here: http://d396qusza40orc.cloudfront.net/repdata%2Fdata%2FStormData.csv.bz2
First step in the project is data download from URL above:
if (!file.exists("data")) {
dir.create("data")
}
if (!file.exists("data/noaa.csv.bz2")) {
download.file(url = "http://d396qusza40orc.cloudfront.net/repdata%2Fdata%2FStormData.csv.bz2", destfile = "data/noaa.csv.bz2", method = "auto")
}Data Processing
Second step is reading data from the file into “maindata” variable:
maindata <- read.csv("data/noaa.csv.bz2")The data have the following structure:
str(maindata)## '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 ""," Christiansburg",..: 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 ""," CANTON"," TULIA",..: 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","%SD",..: 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 "","\t","\t\t",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ REFNUM : num 1 2 3 4 5 6 7 8 9 10 ...THe data analysis must address 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?
Overall population health depends on number of injures and deaths from events. Corresponding columns in the data to represent the information are “INJURIES” and “FATALITIES”.
Total number of deaths for the period of monitoring is: 15145 Total number of injuries for the period of monitoring is: 140528
To answer the first question events should be agregated by event types. It is done by using dplyr library to total deaths and injuries. Also, data in EVTYPE sgould be refine. In the begginig, values in the column should be uppercased to eliminate different case duplicates, then trimmed spaces and group events which are started with summary in name.
maindata.health <- maindata %>%
select(FATALITIES, INJURIES, EVTYPE) %>%
mutate(EVTYPE = gsub("^SUMMARY.*", "SUMMARY", gsub("^\\s+|\\s+$", "", toupper(EVTYPE)))) %>%
group_by(EVTYPE) %>%
summarise(
fatal = sum(FATALITIES, na.rm = TRUE),
inj = sum(INJURIES, na.rm = TRUE)
) %>%
arrange(desc(fatal))
maindata.health$EVTYPE <- factor(maindata.health$EVTYPE, levels=maindata.health$EVTYPE)Economic effect of events is represented by PROPDMG which is dollar property damage and CROPDMG which crop damage. The column PROPDMGEXP and CROPDMGEXP indicate the dimension. Auxilary list would be created to determine biases by factors in PROPDMGEXP and CROPDMGEXP columns.
biases <- list("K" = 10 ^ 3, "M" = 10 ^ 6, "B" = 10 ^ 9, "0" = 1)Let’s update levels for EXP columns to match biases:
maindata$PROPDMGEXP[maindata$PROPDMGEXP=="b"] <- "B"
maindata$PROPDMGEXP[maindata$PROPDMGEXP=="m"] <- "M"
maindata$PROPDMGEXP[maindata$PROPDMGEXP=="k"] <- "K"
maindata$PROPDMGEXP[!(maindata$PROPDMGEXP %in% c("B", "M", "K"))] <- "0"
maindata$CROPDMGEXP[maindata$CROPDMGEXP=="b"] <- "B"
maindata$CROPDMGEXP[maindata$CROPDMGEXP=="m"] <- "M"
maindata$CROPDMGEXP[maindata$CROPDMGEXP=="k"] <- "K"
maindata$CROPDMGEXP[!(maindata$CROPDMGEXP %in% c("B", "M", "K"))] <- "0"After cleaning data biases and dplyr library are used to analyse economical effect and health data.
maindata.economic <- maindata %>%
select(PROPDMG, PROPDMGEXP, CROPDMG, CROPDMGEXP, EVTYPE) %>%
mutate(EVTYPE = gsub("^SUMMARY.*", "SUMMARY", gsub("^\\s+|\\s+$", "", toupper(EVTYPE)))) %>%
mutate(PROPDMG = PROPDMG * as.numeric(biases[as.character(PROPDMGEXP)])) %>%
mutate(CROPDMG = CROPDMG * as.numeric(biases[as.character(CROPDMGEXP)])) %>%
group_by(EVTYPE) %>%
summarise(
prop = sum(PROPDMG, na.rm = TRUE),
crop = sum(CROPDMG, na.rm = TRUE)
) %>%
arrange(desc(prop))
maindata.economic$EVTYPE <- factor(maindata.economic$EVTYPE, levels=maindata.economic$EVTYPE)Result
10 Most harmful events by number of death caused
maindata.health[1:10,1:2]## Source: local data frame [10 x 2]
##
## EVTYPE fatal
## 1 TORNADO 5633
## 2 EXCESSIVE HEAT 1903
## 3 FLASH FLOOD 978
## 4 HEAT 937
## 5 LIGHTNING 816
## 6 TSTM WIND 504
## 7 FLOOD 470
## 8 RIP CURRENT 368
## 9 HIGH WIND 248
## 10 AVALANCHE 224ggplot(data=maindata.health[1:10,], aes(x=EVTYPE, y=fatal)) +
geom_bar(fill="blue", color="black", stat="identity") +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
xlab("Number of deaths") + ylab("Event") +
ggtitle("Top 10 Number of deaths by event")
10 Most harmful events by number of injures caused
arrange(maindata.health, desc(inj))[1:10,c(1,3)]## Source: local data frame [10 x 2]
##
## EVTYPE inj
## 1 TORNADO 91346
## 2 TSTM WIND 6957
## 3 FLOOD 6789
## 4 EXCESSIVE HEAT 6525
## 5 LIGHTNING 5230
## 6 HEAT 2100
## 7 ICE STORM 1975
## 8 FLASH FLOOD 1777
## 9 THUNDERSTORM WIND 1488
## 10 HAIL 136110 Most harmful events by property damaged
maindata.economic[1:10,1:2]## Source: local data frame [10 x 2]
##
## EVTYPE prop
## 1 FLOOD 144657709807
## 2 HURRICANE/TYPHOON 69305840000
## 3 TORNADO 56937160779
## 4 STORM SURGE 43323536000
## 5 FLASH FLOOD 16140862067
## 6 HAIL 15732267048
## 7 HURRICANE 11868319010
## 8 TROPICAL STORM 7703890550
## 9 WINTER STORM 6688497251
## 10 HIGH WIND 5270046295ggplot(data=maindata.economic[1:10,], aes(x=EVTYPE, y=prop)) +
geom_bar(fill="blue", color="black", stat="identity") +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
xlab("Property damage") + ylab("Event") +
ggtitle("Top 10 property damage by event")
10 Most harmful events by crop damaged
arrange(maindata.economic, desc(crop))[1:10,c(1,3)]## Source: local data frame [10 x 2]
##
## EVTYPE crop
## 1 DROUGHT 13972566000
## 2 FLOOD 5661968450
## 3 RIVER FLOOD 5029459000
## 4 ICE STORM 5022113500
## 5 HAIL 3025954473
## 6 HURRICANE 2741910000
## 7 HURRICANE/TYPHOON 2607872800
## 8 FLASH FLOOD 1421317100
## 9 EXTREME COLD 1312973000
## 10 FROST/FREEZE 1094186000