This is an R Markdown document to explore the NOAA Storm Database and answer some basic questions about severe weather events. There will be a data extraction from database and a cleansing. Then with some analysis we will try to find the most harmful events for health population. All the code and information about each code chunk is available for future research reproduction. Requirements :
First of all, data must be loaded from the database. Link below shows the address of storm database. We load Storm Data into variable “data” through the code below,
url<-"http://d396qusza40orc.cloudfront.net/repdata%2Fdata%2FStormData.csv.bz2"
download.file(url=url, destfile="stormdata.csv.bz2")
bz2stream<-bzfile("stormdata.csv.bz2", "r")
data<-read.table(bz2stream, sep = ",", header=TRUE, row.names = NULL)
close(bz2stream)
If date is needed anytime, it's better for it to be of date class
data$BGN_DATE<-as.Date(as.character(data$BGN_DATE), "%m/%d/%Y")
data$END_DATE<-as.Date(as.character(data$END_DATE), "%m/%d/%Y")
In documents explaining data we have just 48 event types. But in the data the event type field has about 1000 distinct values. First of all we must decrease these 1000 values as much as possible without losing any specific event. I had done the following steps :
1. Put all the 48 events in a vector called validEVTYPEList
This is a matrix which contains valid event types in its first column du to “National Weather Service Storm Data Documentation” and their corresponding words that will be seen in real data in its second column
## Event-Official-Name Event-Seen-Words
## [1,] "Astronomical Low Tide" "Astronomical Low Tide"
## [2,] "Avalanche" "Avalanche"
## [3,] "Blizzard" "Blizzard"
## [4,] "Coastal Flood" "Coastal Flood"
## [5,] "Cold/Wind Chill" "Cold/Wind Chill"
## [6,] "Debris Flow" "Debris Flow"
## [7,] "Dense Fog" "Dense Fog"
## [8,] "Dense Smoke" "Dense Smoke"
## [9,] "Drought" "Drought"
## [10,] "Dust Devil" "Dust Devil"
## [11,] "Dust Storm" "Dust Storm"
## [12,] "Excessive Heat" "Excessive Heat"
## [13,] "Extreme Cold/Wind Chill" "Extreme Cold/Wind Chill"
## [14,] "Flash Flood" "Flash Flood"
## [15,] "Flood" "Flood"
## [16,] "Frost/Freeze" "Frost/Freeze"
## [17,] "Funnel Cloud" "Funnel Cloud"
## [18,] "Freezing Fog" "Freezing Fog"
## [19,] "Hail" "Hail"
## [20,] "Heat" "Heat"
## [21,] "Heavy Rain" "Heavy Rain"
## [22,] "Heavy Snow" "Heavy Snow"
## [23,] "High Surf" "High Surf"
## [24,] "High Wind" "High Wind"
## [25,] "Hurricane" "Hurricane"
## [26,] "Ice Storm" "Ice Storm"
## [27,] "Lake-Effect Snow" "Lake-Effect Snow"
## [28,] "Lakeshore Flood" "Lakeshore Flood"
## [29,] "Lightning" "Lightning"
## [30,] "Rip Current" "Rip Current"
## [31,] "Seiche" "Seiche"
## [32,] "Sleet" "Sleet"
## [33,] "Storm Surge/Tide" "Storm Surge/Tide"
## [34,] "Strong Wind" "Strong Wind"
## [35,] "Thunderstorm Wind" "Thunderstorm Wind"
## [36,] "Tornado" "Tornado"
## [37,] "Tropical Depression" "Tropical Depression"
## [38,] "Tropical Storm" "Tropical Storm"
## [39,] "Tsunami" "Tsunami"
## [40,] "Volcanic Ash" "Volcanic Ash"
## [41,] "Waterspout" "Waterspout"
## [42,] "wild/forest fire" "wild/forest fire"
## [43,] "Winter Storm" "Winter Storm|winter-storm|winterstorm"
## [44,] "Winter Weather" "Winter Weather|wintry mix"
## [45,] "Marine Hail" "Marine Hail"
## [46,] "Marine High Wind" "Marine High Wind"
## [47,] "Marine Strong Wind" "Marine Strong Wind"
## [48,] "Marine Thunderstorm Wind" "Marine Thunderstorm Wind"
as it seems,
2. A new attribute for data named “validEVTYPE”“ will be defined and filled with the corresponding validEVTYPE due to each EVTYPE in the data
Using the lapply I loop over all validEVTYPEList values and if find any similar words in EVTYPE using grep function, put the validEVTYPEList value as the EVTYPE's corresponding validEVTYPE. Cause, "data” is a non-environmental variable for lapply, to assign “data” a value we must use <<- sign
tmp<-apply(validEVTYPEList, 1, function(x) {
data[grep(x[2], data$EVTYPE, ignore.case=TRUE),
"validEVTYPE"]<<-x[1]})
rm(tmp)
3. Check the result of validEVTYPE using the following codes
length(data[is.na(data$validEVTYPE), "EVTYPE"])
to see how changes to the validEVTYPE vector causes NAs to vanish and lessen and also check
notConverted<-tapply(data[is.na(data$validEVTYPE), "EVTYPE"],
data[is.na(data$validEVTYPE), "EVTYPE"], length)
head(sort(notConverted, decreasing=TRUE))
to see the most dense event types.
4. do some changes to the validEVTYPEList and run sections 2 and 3 again These changes include :
For example I
added the events :
edited some valid events with the supposition of :
edited some valid events for misspellings :
After doing steps 2,3, and 4 for several times, the final valid event type list will be as follows:
## Event-Official-Name
## [1,] "WARM"
## [2,] "WIND"
## [3,] "DRY"
## [4,] "Astronomical Low Tide"
## [5,] "Avalanche"
## [6,] "Blizzard"
## [7,] "Coastal Flood"
## [8,] "Cold/Wind Chill"
## [9,] "Debris Flow"
## [10,] "Dense Fog"
## [11,] "Dense Smoke"
## [12,] "Drought"
## [13,] "Dust Devil"
## [14,] "Dust Storm"
## [15,] "Excessive Heat"
## [16,] "Extreme Cold/Wind Chill"
## [17,] "Flash Flood"
## [18,] "Flood"
## [19,] "Frost/Freeze"
## [20,] "Funnel Cloud"
## [21,] "Freezing Fog"
## [22,] "Hail"
## [23,] "Heat"
## [24,] "Heavy Rain"
## [25,] "Heavy Snow"
## [26,] "High Surf"
## [27,] "High Wind"
## [28,] "Hurricane"
## [29,] "Ice Storm"
## [30,] "Lake-Effect Snow"
## [31,] "Lakeshore Flood"
## [32,] "Lightning"
## [33,] "Rip Current"
## [34,] "Seiche"
## [35,] "Sleet"
## [36,] "Storm Surge/Tide"
## [37,] "Strong Wind"
## [38,] "Thunderstorm Wind"
## [39,] "Tornado"
## [40,] "Tropical Depression"
## [41,] "Tropical Storm"
## [42,] "Tsunami"
## [43,] "Volcanic Ash"
## [44,] "Waterspout"
## [45,] "wild/forest fire"
## [46,] "Winter Storm"
## [47,] "Winter Weather"
## [48,] "URBAN/SML STREAM FLD"
## [49,] "Landslide"
## [50,] "Summary"
## [51,] "Marine Hail"
## [52,] "Marine High Wind"
## [53,] "Marine Strong Wind"
## [54,] "Marine Thunderstorm Wind"
## Event-Seen-Words
## [1,] "warm"
## [2,] "wind"
## [3,] "dry"
## [4,] "Astronomical Low Tide"
## [5,] "Avalanche"
## [6,] "Blizzard"
## [7,] "Coastal Flood"
## [8,] "Cold/Wind Chill|cold|Wind Chill"
## [9,] "Debris Flow"
## [10,] "Dense Fog|fog"
## [11,] "Dense Smoke|smoke"
## [12,] "Drought"
## [13,] "Dust Devil|DustDevil|Dust-Devil|Dust devel"
## [14,] "Dust Storm|DustStorm|Dust-Storm"
## [15,] "Excessive Heat|ExcessiveHeat|Excessive-Heat"
## [16,] "Extreme Cold|Extreme Cold/Wind Chill|extreme windchill|extreme wind chill"
## [17,] "Flash Flood"
## [18,] "Flood"
## [19,] "Freezing|Freeze|Frost|Frost/Freeze"
## [20,] "Funnel Cloud|FunnelCloud|Funnel-Cloud"
## [21,] "Freezing Fog|ice fog"
## [22,] "Hail"
## [23,] "Heat|Hot|very Warm"
## [24,] "Heavy Rain|rain"
## [25,] "Heavy Snow|snow"
## [26,] "High Surf"
## [27,] "High Wind"
## [28,] "Hurricane"
## [29,] "Ice Storm"
## [30,] "Lake-Effect Snow|Lake Effect Snow"
## [31,] "Lakeshore Flood"
## [32,] "Lightning"
## [33,] "Rip Current"
## [34,] "Seiche"
## [35,] "Sleet"
## [36,] "Tide|Storm Surge|Storm Surge/Tide"
## [37,] "Strong Wind|gusty wind"
## [38,] "TSTM|TSTM Wind|Thunderstorm|Thunderstorm Wind|TSTMW|TSTM WND|THUNDEERSTORM WINDS|THUNDERESTORM WINDS|THUNERSTORM WINDS|THUNDERSTORM W INDS|thunderstorm wins|thundestorm|thuderstorm|thunderstorms"
## [39,] "Tornado"
## [40,] "Tropical Depression"
## [41,] "Tropical Storm"
## [42,] "Tsunami"
## [43,] "Volcanic Ash"
## [44,] "Waterspout|WAYTERSPOUT|WATER SPOUT"
## [45,] "wild/forest fire|Wildfire|forest fire"
## [46,] "Winter Storm|winter-storm|winterstorm"
## [47,] "Winter Weather|wintry mix"
## [48,] "URBAN/SML STREAM FLD"
## [49,] "landslide"
## [50,] "Summary"
## [51,] "Marine Hail"
## [52,] "Marine High Wind|marine high-wind"
## [53,] "Marine Strong Wind|marine strong-wind"
## [54,] "Marine Thunderstorm Wind|Marine Thunderstorm"
This matrix can easily be compared with the basic one for changes.
5. Put the remaining EVTYPEs as “OTHERS” in data$validEVTYPE
data[is.na(data$validEVTYPE), "validEVTYPE"]<-"Others"
So if we repeat the replacement of valid event types for the last time we will have a data of 50 event types. Some of the official event types doesn't seem to appear in real data.
## [1] "Tornado" "Thunderstorm Wind"
## [3] "Hail" "Heavy Rain"
## [5] "Heavy Snow" "Ice Storm"
## [7] "Winter Storm" "Hurricane"
## [9] "Cold/Wind Chill" "Lightning"
## [11] "Dense Fog" "Rip Current"
## [13] "Flood" "High Wind"
## [15] "Funnel Cloud" "Heat"
## [17] "WIND" "Others"
## [19] "Waterspout" "Extreme Cold/Wind Chill"
## [21] "Blizzard" "Frost/Freeze"
## [23] "Storm Surge/Tide" "Avalanche"
## [25] "WARM" "Dust Storm"
## [27] "Sleet" "Dust Devil"
## [29] "Strong Wind" "High Surf"
## [31] "DRY" "Winter Weather"
## [33] "Drought" "Tropical Storm"
## [35] "wild/forest fire" "Lake-Effect Snow"
## [37] "Landslide" "URBAN/SML STREAM FLD"
## [39] "Freezing Fog" "Summary"
## [41] "Volcanic Ash" "Seiche"
## [43] "Tropical Depression" "Dense Smoke"
## [45] "Marine Hail" "Marine High Wind"
## [47] "Tsunami" "Lakeshore Flood"
## [49] "Marine Thunderstorm Wind" "Marine Strong Wind"
We must omit Summary from the list of events cause as it seems, it is just a monthly summary of injuries and damages. In this stage, data is ready for injuries computations.
To compute the real amounts of damages in USD (United states Dollar) we must use the exponential attribute of them.
The official document divides the exponentials to three groups of
But getting the distinct value of exponential for property damage we have :
unique(data$PROPDMGEXP)
## [1] K M B m + 0 5 6 ? 4 2 3 h 7 H - 1 8
## Levels: - ? + 0 1 2 3 4 5 6 7 8 B h H K m M
and getting the distinct value of exponential for cropping damage we have :
unique(data$CROPDMGEXP)
## [1] M K m B ? 0 k 2
## Levels: ? 0 2 B k K m M
For clarifying the amount of damage, I do the following assumptions :
Big and Small Case letters are equal so for example “M” and “m” both mean Million
If value of exponential attribute is a number, it shows the power of 10 to be crossed with the main value.
Anything except digits and the special characters are considered as NA
Codes below add two new attributes to data which contain the amount of dollars wasted for damages of both types “Crop” and “Property” damages.
The code checks in “PROPDMGEXP” and “CROPDMGEXP” for
Any numbers to use the number for the money zeros
If no, letters “k”, “K” to put 3 zeros infront of the number
If no, letters “m”, “M” to put 6 zeros infront of the number
If no, letters “b”, “B” to put 9 zeros infront of the number
Else to set the record as Not Available one (NA)
data$PROPDMGUSD<-data$PROPDMG*ifelse(grepl("[[:digit:]]",data$PROPDMGEXP) & !grepl("[[:alpha:]]",data$PROPDMGEXP),
10^as.numeric(as.character(data$PROPDMGEXP)),
ifelse(data$PROPDMGEXP %in% c("k", "K"), 1000,
ifelse(data$PROPDMGEXP %in% c("M", "m"), 1000000,
ifelse(data$PROPDMGEXP %in% c("B", "b"), 1000000000, NA))))
data$CROPDMGUSD<-data$CROPDMG*ifelse(grepl("[[:digit:]]",data$CROPDMGEXP) & !grepl("[[:alpha:]]",data$CROPDMGEXP),
10^as.numeric(as.character(data$CROPDMGEXP)),
ifelse(data$CROPDMGEXP %in% c("k", "K"), 1000,
ifelse(data$CROPDMGEXP %in% c("M", "m"), 1000000,
ifelse(data$CROPDMGEXP %in% c("B", "b"), 1000000000, NA))))
After reducing Number of events near their officially reported ones, it is now possible to compare various events to each other due to
Across the United States, which types of events are most harmful with respect to population health? Code below computes sum of injuries/fatalities each event has caused
injuriesSum<-aggregate(cbind(data$INJURIES+data$FATALITIES,
data$INJURIES,
data$FATALITIES),
by=list(data$validEVTYPE), sum, na.rm=TRUE)
names(injuriesSum)<-c("event-name", "total-damage", "property-damage", "crop-damage")
If we sort the injuries in a descending order so that the most causable event is the first, we can see that Tornado is the worst Event causing greatest amount of injury. After that, come Heat and Thunderstorm Wind as the enemies of human health. Below you can see the ten most devil events in United States in order.
injuriesSumSorted<-injuriesSum[order(injuriesSum[2], decreasing=TRUE), ]
head(injuriesSumSorted, 10)
## event-name total-damage property-damage crop-damage
## 39 Tornado 97022 91364 5658
## 16 Heat 12362 9224 3138
## 38 Thunderstorm Wind 10237 9518 719
## 11 Flood 10126 8602 1524
## 26 Lightning 6048 5231 817
## 22 Ice Storm 2081 1992 89
## 20 High Wind 1802 1507 295
## 49 Winter Storm 1570 1353 217
## 47 wild/forest fire 1543 1456 87
## 21 Hurricane 1463 1328 135
But if we want to see what is the average of injuries/fatalities per occurance, we must divide amount of “injuries+fatalities” by number of recorded occurance of each event. code below computes this ratio.
injuriesOccurance<-tapply(data$INJURIES+data$FATALITIES, data$validEVTYPE, length)
injuriesRatioSum<-cbind(injuriesSum[1], injuriesSum[2:4]/injuriesOccurance)
injuriesRatioSorted<-injuriesRatioSum[order(injuriesRatioSum[2], decreasing=TRUE), ]
head(injuriesRatioSorted, 10)
## event-name total-damage property-damage crop-damage
## 42 Tsunami 8.1000 6.4500 1.65000
## 21 Hurricane 5.0799 4.6111 0.46875
## 16 Heat 4.6352 3.4586 1.17660
## 39 Tornado 1.5988 1.5055 0.09323
## 32 Rip Current 1.4234 0.6808 0.74260
## 9 Dust Storm 1.0794 1.0280 0.05140
## 22 Ice Storm 1.0246 0.9808 0.04382
## 1 Avalanche 1.0207 0.4404 0.58031
## 31 Others 0.9068 0.7945 0.11231
## 29 Marine Strong Wind 0.7500 0.4583 0.29167
As you can see, there are big differences between these two result. Tsunami and Hurricane cause more injuries considering their frequency of occurance.
It means that in average, Tsunami and Hurricane cause more injuries than for example Tornado each time they occure.
Something like “Dust Storm” isn't even in the first ten most injury-causings. But the result is actually because “Dust Storms” are rare in occurance. Due to both these investigates Heat is always a great danger to Population Health.
Here we can see the five worst events in both “ratio”“ and "raw sum” format
par(mfrow=c(2,1))
toDraw<-matrix(injuriesSumSorted[1:5, 2], nrow=1, ncol=5)
colnames(toDraw)<-injuriesSumSorted[1:5, 1]
barplot(toDraw,
xlab="Event", ylab="Number of Injury/Fatility", main="Five Events causing the Worst Injuries across U.S.")
toDraw<-matrix(injuriesRatioSorted[1:5, 2], nrow=1, ncol=5)
colnames(toDraw)<-injuriesRatioSorted[1:5, 1]
barplot(toDraw,
xlab="Event", ylab="Ratio of Injury/Fatility", main="Five Events causing Most Injuries per Occurance across U.S.")
Across the United States, which types of events have the greatest economic consequences?
Code below computes sum of damage amount that each event has caused. Then the events are sorted due to the amount of their overal damages
damageSum<-aggregate(cbind(data$PROPDMGUSD+data$CROPDMGUSD,
data$PROPDMGUSD,
data$CROPDMGUSD),
by=list(data$validEVTYPE), sum, na.rm=TRUE)
names(damageSum)<-c("event-name", "total-damage", "property-damage", "crop-damage")
damageSumSorted<-damageSum[order(damageSum[2], decreasing=TRUE), ]
head(damageSumSorted, 10)
## event-name total-damage property-damage crop-damage
## 11 Flood 1.578e+11 1.682e+11 1.227e+10
## 21 Hurricane 4.431e+10 8.476e+10 5.515e+09
## 39 Tornado 1.817e+10 5.855e+10 4.175e+08
## 15 Hail 1.005e+10 1.598e+10 3.047e+09
## 22 Ice Storm 5.925e+09 3.946e+09 5.022e+09
## 38 Thunderstorm Wind 5.513e+09 1.118e+10 1.272e+09
## 35 Storm Surge/Tide 4.645e+09 4.797e+10 8.550e+05
## 47 wild/forest fire 3.839e+09 7.872e+09 4.033e+08
## 20 High Wind 3.296e+09 6.005e+09 6.868e+08
## 6 Drought 1.886e+09 1.046e+09 1.397e+10
As you can see Flood has the most total economic cost on US Government which is more than 1000 billion in all days from 1950 till 2011. Hurricane and Tornado are in the next levels. As the data shows Hurricane and Tornado are overally one of the worst evil events, Humanly and Economically.
The ratio part can also be computed exactly like the injuries ratio. But economics is something different from injuries. Computing economic cost per each occurance does not to be enough to compare Events. But if it can be useful code below computes the ratio of economic cost per each occurance :
damageOccurance<-tapply(data$INJURIES+data$FATALITIES, data$validEVTYPE, length)
damageRatioSum<-cbind(damageSum[1], damageSum[2:4]/damageOccurance)
damageRatioSorted<-damageRatioSum[order(damageRatioSum[2], decreasing=TRUE), ]
head(damageRatioSorted, 10)
## event-name total-damage property-damage crop-damage
## 21 Hurricane 153866131 294292292 19150322
## 42 Tsunami 7204100 7203100 1000
## 35 Storm Surge/Tide 6712049 69327268 1236
## 22 Ice Storm 2917356 1942653 2472729
## 41 Tropical Storm 2195627 11067992 996981
## 11 Flood 1908449 2034589 148404
## 47 wild/forest fire 907888 1861983 95384
## 6 Drought 758206 420461 5615983
## 13 Frost/Freeze 613170 12578 1298479
## 24 Lakeshore Flood 327826 327826 0