# Load necessary libraries
library(knitr)
library(markdown)
library(rmarkdown)
library(plyr)
library(stats)
# Load the data
storm <- read.csv(file = "repdata_data_StormData.csv", header = TRUE, sep = ",")
dim(storm)
## [1] 902297 37
names(storm)
## [1] "STATE__" "BGN_DATE" "BGN_TIME" "TIME_ZONE" "COUNTY"
## [6] "COUNTYNAME" "STATE" "EVTYPE" "BGN_RANGE" "BGN_AZI"
## [11] "BGN_LOCATI" "END_DATE" "END_TIME" "COUNTY_END" "COUNTYENDN"
## [16] "END_RANGE" "END_AZI" "END_LOCATI" "LENGTH" "WIDTH"
## [21] "F" "MAG" "FATALITIES" "INJURIES" "PROPDMG"
## [26] "PROPDMGEXP" "CROPDMG" "CROPDMGEXP" "WFO" "STATEOFFIC"
## [31] "ZONENAMES" "LATITUDE" "LONGITUDE" "LATITUDE_E" "LONGITUDE_"
## [36] "REMARKS" "REFNUM"
str(storm)
## '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 ...
# Select relevant variables
varsNedeed <- c("EVTYPE", "FATALITIES", "INJURIES", "PROPDMG", "PROPDMGEXP", "CROPDMG", "CROPDMGEXP")
storm <- storm[varsNedeed]
dim(storm)
## [1] 902297 7
names(storm)
## [1] "EVTYPE" "FATALITIES" "INJURIES" "PROPDMG" "PROPDMGEXP"
## [6] "CROPDMG" "CROPDMGEXP"
str(storm)
## 'data.frame': 902297 obs. of 7 variables:
## $ EVTYPE : chr "TORNADO" "TORNADO" "TORNADO" "TORNADO" ...
## $ 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 "" "" "" "" ...
# Refactor PROPDMGEXP variable
unique(storm$PROPDMGEXP)
## [1] "K" "M" "" "B" "m" "+" "0" "5" "6" "?" "4" "2" "3" "h" "7" "H" "-" "1" "8"
storm$PROPDMGEXP <- mapvalues(storm$PROPDMGEXP, from = c("K", "M", "", "B", "m", "+", "0", "5", "6", "?", "4", "2", "3", "h", "7", "H", "-", "1", "8"), to = c(10^3, 10^6, 1, 10^9, 10^6, 0, 1, 10^5, 10^6, 0, 10^4, 10^2, 10^3, 10^2, 10^7, 10^2, 0, 10, 10^8))
storm$PROPDMGEXP <- as.numeric(as.character(storm$PROPDMGEXP))
storm$PROPDMGTOTAL <- (storm$PROPDMG * storm$PROPDMGEXP) / 1000000000
# Refactor CROPDMGEXP variable
unique(storm$CROPDMGEXP)
## [1] "" "M" "K" "m" "B" "?" "0" "k" "2"
storm$CROPDMGEXP <- mapvalues(storm$CROPDMGEXP, from = c("", "M", "K", "m", "B", "?", "0", "k", "2"), to = c(1, 10^6, 10^3, 10^6, 10^9, 0, 1, 10^3, 10^2))
storm$CROPDMGEXP <- as.numeric(as.character(storm$CROPDMGEXP))
storm$CROPDMGTOTAL <- (storm$CROPDMG * storm$CROPDMGEXP) / 1000000000
# Calculate total fatalities per event
sumFatalities <- aggregate(FATALITIES ~ EVTYPE, data = storm, FUN = "sum")
dim(sumFatalities)
## [1] 985 2
# Top 10 events with the highest fatalities
fatalities10events <- sumFatalities[order(-sumFatalities$FATALITIES), ][1:10, ]
dim(fatalities10events)
## [1] 10 2
fatalities10events
## EVTYPE FATALITIES
## 834 TORNADO 5633
## 130 EXCESSIVE HEAT 1903
## 153 FLASH FLOOD 978
## 275 HEAT 937
## 464 LIGHTNING 816
## 856 TSTM WIND 504
## 170 FLOOD 470
## 585 RIP CURRENT 368
## 359 HIGH WIND 248
## 19 AVALANCHE 224
# Barplot for top 10 events with the highest fatalities
par(mfrow = c(1, 1), mar = c(12, 4, 3, 2), mgp = c(3, 1, 0), cex = 0.8)
barplot(fatalities10events$FATALITIES, names.arg = fatalities10events$EVTYPE, las = 3, main = "10 Fatalities Highest Events", ylab = "Number of Fatalities")
dev.copy(png, "fatalities-events.png", width = 480, height = 480)
## png
## 3
dev.off()
## png
## 2
# Calculate total injuries per event
sumInjuries <- aggregate(INJURIES ~ EVTYPE, data = storm, FUN = "sum")
dim(sumInjuries)
## [1] 985 2
# Top 10 events with the highest injuries
injuries10events <- sumInjuries[order(-sumInjuries$INJURIES), ][1:10, ]
dim(injuries10events)
## [1] 10 2
injuries10events
## EVTYPE INJURIES
## 834 TORNADO 91346
## 856 TSTM WIND 6957
## 170 FLOOD 6789
## 130 EXCESSIVE HEAT 6525
## 464 LIGHTNING 5230
## 275 HEAT 2100
## 427 ICE STORM 1975
## 153 FLASH FLOOD 1777
## 760 THUNDERSTORM WIND 1488
## 244 HAIL 1361
# Barplot for top 10 events with the highest injuries
par(mfrow = c(1, 1), mar = c(12, 6, 3, 2), mgp = c(4, 1, 0), cex = 0.8)
barplot(injuries10events$INJURIES, names.arg = injuries10events$EVTYPE, las = 3, main = "10 Injuries Highest Events", ylab = "Number of Injuries")
dev.copy(png, "injuries-events.png", width = 480, height = 480)
## png
## 3
dev.off()
## png
## 2
# Calculate total property damage per event
sumPropertyDamage <- aggregate(PROPDMGTOTAL ~ EVTYPE, data = storm, FUN = "sum")
dim(sumPropertyDamage)
## [1] 985 2
# Top 10 events with the highest property damage
propdmg10Total <- sumPropertyDamage[order(-sumPropertyDamage$PROPDMGTOTAL), ][1:10, ]
propdmg10Total
## EVTYPE PROPDMGTOTAL
## 170 FLOOD 144.657710
## 411 HURRICANE/TYPHOON 69.305840
## 834 TORNADO 56.947381
## 670 STORM SURGE 43.323536
## 153 FLASH FLOOD 16.822674
## 244 HAIL 15.735268
## 402 HURRICANE 11.868319
## 848 TROPICAL STORM 7.703891
## 972 WINTER STORM 6.688497
## 359 HIGH WIND 5.270046
# Barplot for top 10 events with the highest property damage
par(mfrow = c(1, 1), mar = c(12, 6, 3, 2), mgp = c(3, 1, 0), cex = 0.8)
barplot(propdmg10Total$PROPDMGTOTAL, names.arg = propdmg10Total$EVTYPE, las = 3, main = "10 Property Damages Highest Events", ylab = "Damage Property Values (in Billions)")
dev.copy(png, "propdmg-total.png", width = 480, height = 480)
## png
## 3
dev.off()
## png
## 2
# Calculate total crop damage per event
sumCropDamage <- aggregate(CROPDMGTOTAL ~ EVTYPE, data = storm, FUN = "sum")
dim(sumCropDamage)
## [1] 985 2
# Top 10 events with the highest crop damage
cropdmg10Total <- sumCropDamage[order(-sumCropDamage$CROPDMGTOTAL), ][1:10, ]
cropdmg10Total
## EVTYPE CROPDMGTOTAL
## 95 DROUGHT 13.972566
## 170 FLOOD 5.661968
## 590 RIVER FLOOD 5.029459
## 427 ICE STORM 5.022113
## 244 HAIL 3.025954
## 402 HURRICANE 2.741910
## 411 HURRICANE/TYPHOON 2.607873
## 153 FLASH FLOOD 1.421317
## 140 EXTREME COLD 1.292973
## 212 FROST/FREEZE 1.094086
# Barplot for top 10 events with the highest crop damage
par(mfrow = c(1, 1), mar = c(10, 6, 3, 2), mgp = c(3, 1, 0), cex = 0.6)
barplot(cropdmg10Total$CROPDMGTOTAL, names.arg = cropdmg10Total$EVTYPE, las = 2, main = "10 Crop Damages Highest Events", ylab = "Damage Crop Values (in Billions)")
dev.copy(png, "cropdmg-total.png", width = 480, height = 480)
## png
## 3
dev.off()
## png
## 2