2 Data Processing

2.1 Loading the data

library(dplyr)

cfile <- "StormData.csv.bz2"
url <- paste("https://d396qusza40orc.cloudfront.net/repdata/data/", cfile, sep="")

(!file.exists(cfile)) && { download.file(url, destfile=cfile, method="curl", quiet=TRUE) }

## [1] FALSE

df <- tibble::as.tibble(read.csv(cfile, header=TRUE, na.strings="NA", stringsAsFactors=FALSE))

dim(df)

## [1] 902297     37

head(df, 5)

2.2 Tidy up

Looking at the dataset and the Storm Data Event Table, pd01016005curr.pdf, section 2.1.1, it looks like some values in the event type column are spelled differently (sometimes, even the case is different) but identical in meaning and should be put together.

Anticipating on the required columns to address the two main questions, I will only work on a subset of the original dataset, namely the following columns:

EVTYPE, FATALITIES, INJURIES, PROPDMG, PROPDMGEXP, CROPDMG, CROPDMGEXP

df <- subset(df, select=c("EVTYPE", "FATALITIES", "INJURIES", "PROPDMG", "PROPDMGEXP", "CROPDMG", "CROPDMGEXP"))

Before applying transformations on the column EVTYPE, we can see that we have 985 different event types:

length(unique(df$EVTYPE))

## [1] 985

For example, the values beginning with FLOOD are the following (notice the spelling, case, plural/singular, repetition… Typical of how, we human interprets, reports things slightly differenlty):

unique(df[grep("^FLOOD", df$EVTYPE, ignore.case=TRUE, perl=TRUE),]$EVTYPE)

##  [1] "FLOODING"             "FLOOD"                "FLOOD WATCH/"        
##  [4] "FLOODING/HEAVY RAIN"  "FLOOD/FLASH FLOOD"    "FLOOD/RAIN/WINDS"    
##  [7] "FLOOD/RAIN/WIND"      "FLOODS"               "FLOOD/RIVER FLOOD"   
## [10] "FLOOD/FLASH FLOODING" "FLOOD FLASH"          "FLOOD FLOOD/FLASH"   
## [13] "FLOOD/FLASH"          "FLOOD & HEAVY RAIN"   "FLOOD/FLASHFLOOD"    
## [16] "Flood/Flash Flood"    "Flood"                "Flood/Strong Wind"   
## [19] "FLOOD/FLASH/FLOOD"

The following functions work on the EVTYPE column of our dataset, trying to clean up that column.
Several pipelines are defined on EVTYPE to:

remove all extra spaces (whether at the edges or betwwen words),
replace / with -,
remove some special characters (double back-quote, ampersand…)
capitalize the words,
suppress pluralization,
replace abbreviation (when known),
swap order,
trying to consolidate synonym expressions

library(stringr)

capitalize <- function(str) {
  s <- strsplit(str, " ")[[1]]
  paste(toupper(substring(s, 1, 1)), tolower(substring(s, 2)), sep="", collapse=" ")
}

# character replacement
s_normalize <- function(v) {
  v %>%
    str_trim() %>%
    str_replace("AND|&", "/") %>%
    str_replace("\\\\", "/") %>% 
    str_replace_all("/", " / ") %>%
    str_replace_all(" / | /|/ ", "-") %>%
    sapply(capitalize, simplify=TRUE)
}

# applying regexp to replace some values
s_replace <- function(d_f, colname, pattern, repl, icase=TRUE, p=TRUE) { 
  ixs <- grep(pattern, d_f[[colname]], ignore.case=icase, perl=p)
  if (length(ixs) != 0) {
    d_f[ixs,][[colname]] = repl
  }
  d_f
} 

# define a sequence of changes - this can be extended....
wr_replace <- function(d_f) {
  d_f %>%
    s_replace("EVTYPE", "[ ]+/$", "") %>%
    s_replace("EVTYPE", "^River Flood", "Flood") %>%  # 2.1.1 refer only to Flood
    s_replace("EVTYPE", "Flooding|Floods", "Flood") %>%
    s_replace("EVTYPE", "Flood Flood", "Flood") %>%
    s_replace("EVTYPE", "Flood.Flash.*$", "Flash-Flood") %>%
    s_replace("EVTYPE", "Flood-Rain/wind.?", "Flood-Rain-Wind") %>%
    s_replace("EVTYPE", "^Heavy Rains$", "Heavy Rain") %>%
    s_replace("EVTYPE", "^Tstm Win|^Thunderstorm Win|^Thunderstormw", "Thunderstorm Wind") %>%
    s_replace("EVTYPE", "Typhoon|Hurricane", "Typhoon-Hurricane") %>%
    s_replace("EVTYPE", "^Hail.?(?:\\d+|\\d+\\.\\d+)\\)?$", "Hail") %>%
    s_replace("EVTYPE", "^High Winds? \\d+$", "High Wind")
}

# call previous defnied functions
combine <- function(df) {
  df$EVTYPE <- s_normalize(df$EVTYPE)
  wr_replace(df)
}

ndf <- combine(df)

Now the column EVTYPE, has 726 consolidated event types:

length(unique(ndf$EVTYPE))

## [1] 726

For example, the values begining with FLOOD are the following:

unique(ndf[grep("^Flood", ndf$EVTYPE, perl=TRUE),]$EVTYPE)

## [1] "Flood"              "Flood Watch-"       "Flood-rain-winds"  
## [4] "Flood-rain-wind"    "Flood-river Flood"  "Flood - Heavy Rain"
## [7] "Flood-strong Wind"

Note: Some more work can be done on the remaining values on the EVTYPE column, however this does not have impact on the top n results (n being 10) of my answers to the two main questions.

2.3 Most harmful meteorological events for the population health

According to datasset (Storm Data Event Table, pd01016005curr.pdf, section 2.1.1) the relevant fields (columns) to address this question are the following:

EVTYPE which we will rename as event_type,
FATALITIES,
INJURIES.

No missing values, as we can see from:

sapply(subset(df, select=c('EVTYPE', 'FATALITIES', 'INJURIES')), function(x) {sum(is.na(x))})

##     EVTYPE FATALITIES   INJURIES 
##          0          0          0

In the following my interpretation is to consider fatalities, then injuries in this order to determine the most harmful events.

To summarize the data we use a pipeline (using R package dplyr) to

only retain the columns EVTYPE, FATALITIES and INJURIES (a projection onto these columns)
a grouping by EVTYPE
a summation over fatalities and injuries columns (the na.rm=TRUE can be omitted),
a sort to reorder the dataset per fatalities, then injuries in descending order,
a final renaming

The code is as follows (and less verbose):

ndf_harm <- ndf %>%
  select(EVTYPE, FATALITIES, INJURIES) %>%
  group_by(EVTYPE) %>%
  summarize(fatalities=sum(FATALITIES, na.rm=TRUE), injuries=sum(INJURIES, na.rm=TRUE)) %>%
  arrange(desc(fatalities), desc(injuries)) %>%
  transmute(event_type=EVTYPE, fatalities, injuries)

The top 10 worst meteorological events in term of fatalities and injuries are the following:

head(ndf_harm, 10)

We can now build a plot with two graphs representing fatalities and injuries (in this order and next to each other) for the top 10 worst events.

library(ggplot2)
library(ggpubr)  # for ggarrange

N <- 10

# helper functions

# ggplot, using bar geometry
gp <- function(ndf, yaxe, yaxe_lab, n=N) {
  ggplot(ndf, aes(y=yaxe, x=event_type)) +
    geom_bar(stat="identity", aes(fill=event_type), alpha=0.7, position="dodge") +
    geom_label(aes(label=yaxe), group=ndf$event_type, size=2.2, colour="black", show.legend=FALSE) + 
    ylab(yaxe_lab) + xlab("Event Type") +
    coord_flip() +
    ggtitle(paste0("Top", n, " most harmful events"))
}

# customize the theme a bit
gtheme <- function(yaxe_title=TRUE, title=TRUE, hjust=0.9, 
                   axis.text.y=element_blank()) {
  thm <- theme_classic() +
    theme(
      text = element_text(size=9),
      panel.background = element_rect(fill="azure1", colour="black",
                                      size=0.5, linetype="solid"),
      strip.background = element_rect(fill="lightskyblue3", colour="black",
                                      size=0.5),
      axis.text.y = axis.text.y
    )
  if (!yaxe_title) {
    thm <- thm + 
      theme(axis.title.y = element_blank())
  }
  if (title) {
    thm <- thm + 
      theme(plot.title = element_text(color="black", size=12, face="bold", hjust=hjust))
  }
  else {
    thm <- thm + 
      theme(plot.title = element_text(color="transparent"))
  }
  return(thm)
}

# re-organizing the data (for the 10 top, so that order in plot reflects order in data frame)
ndf_top_n <- head(ndf_harm, N) %>%
  arrange(fatalities, injuries) %>%
  mutate(event_type=as.character(event_type)) %>%
  mutate(event_type=factor(event_type, levels=unique(event_type), ordered=TRUE))

# arrange the plots in one grid witjh a common legend
gp11 <- gp(ndf_top_n, ndf_top_n$fatalities, '#Fatalities')
gp12 <- gp(ndf_top_n, ndf_top_n$injuries, '#Injuries')
ggarrange(gp11 + gtheme(), 
          gp12 + gtheme(yaxe_title=FALSE, title=FALSE), 
          ncol=2, nrow=1,
          common.legend = TRUE, legend = "right")

2.4 Event type causing the greatest economic consequences

According to the datasset (Storm Data Event Table, pd01016005curr.pdf, section 2.1.1) the relevant fields (columns) to address this question are the following:

EVTYPE renamed as event_type,
PROPDMG (property damage), evaluated in $US,
PROPDMGEXP (property damage exponent), a scale factor
CROPDMG (crop damage), evaluated in $US
CROPDMGEXP (crop damage exponent), a scale factor

the PROPDMGEXP and CROPDMGEXP columns have the following values:

unique(ndf$PROPDMGEXP)

##  [1] "K" "M" ""  "B" "m" "+" "0" "5" "6" "?" "4" "2" "3" "h" "7" "H" "-"
## [18] "1" "8"

unique(ndf$CROPDMGEXP)

## [1] ""  "M" "K" "m" "B" "?" "0" "k" "2"

We need to transform the 2 columns PROPDMGEXP and CROPDMGEXP into numbers thant we can use to scale (multiply) the value of the counterpart columns, respectively PROPDMG and CROPDMG.

Section 2.7 of pd01016005curr.pdf, states that “B” stands for billions, “M” for millions… We use this rule regardless of the case.
We transform the integer into power of 10, so 1 will become 10^1, 2 will become 10^2…
"" and "?" and 0 will be turn into 0 value (Note: 0 could have been turned into 10^0 == 1)
Less obvious is what to do with + and - symbols which should be explained in the narrative (this is alluded to in section 2.7 of pd01016005curr.pdf), however I could not find any clues. I decided to ignore and it does not affect my results for the top 10. However, beyond that I would have to reconsider this choice.

Here is my proposal for the conversion:

p_convert <- function(sym) {
  if (sym %in% c("", "+", "?", "-", "0")) { # WARN: +, - may need to be reconsidered
    0
  }
  else if (sym %in% c("B", "b")) {
    1e9
  }
  else if (sym %in% c("M", "m")) {
    1e6
  }
  else if (sym %in% c("K", "k")) {
    1e3
  }
  else if (sym %in% c('H', 'h' )) {
    1e2
  }
  else if (sym %in% c("1", "2", "3", "4", "5", "6", "7", "8")) { # also "0" ?
    10^as.numeric(sym)
  }
}

Next, we are now ready to extract the relevant data to answer the question about the greatest economic consequences.
We define a new pipeline that:

transform the existing columns by multiplying together the two relevant columns for property damage and crop damage respectively,
group by event type then
summarize by summing over property damage and crop damage respectively,
rename the columns

As the ranking and finanacial estimate expressed in $US for property damages and crop damages are different, I decided to make separate plots.

Here is the snippet for the property damage:

ndf_eco <- ndf %>%
  select(EVTYPE, PROPDMG, PROPDMGEXP, CROPDMG, CROPDMGEXP) %>%
  transmute(event_type=EVTYPE,
            property_damage = PROPDMG * sapply(PROPDMGEXP, p_convert),
            crop_damage = CROPDMG * sapply(CROPDMGEXP, p_convert)) %>%
  group_by(event_type) %>%
  summarize(property_damage=sum(property_damage, na.rm=TRUE),
            crop_damage=sum(crop_damage, na.rm=TRUE))

We can then look at the 10 greatest economic consequences for the properties (expressed in billions of $US, hence the division by 1e+09 and rounding):

N <- 10

ndf_top_n_eco_prop <- ndf_eco %>%
  select(event_type, property_damage) %>%
  arrange(desc(property_damage)) %>%
  head(N) %>%
  mutate(property_damage = round(property_damage / 1e+09, 3))

ndf_top_n_eco_prop

Which can be visualized with the following plot (using the same “reoganization” for the plot as for the first plot which shows a decreading order this tile from left to right):

gp2 <- function(ndf, yaxe, xaxe, yaxe_lab, title) {
  ggplot(ndf, aes(y=yaxe, x=xaxe)) +
    geom_bar(stat="identity", aes(fill=xaxe), alpha=0.7, position="dodge") +
    geom_label(aes(label=yaxe), group=ndf$event_type, size=2.2, colour="black", show.legend=FALSE) + 
    ylab(yaxe_lab) + xlab("Event Type") +
    ggtitle(title)
}

reorg <- function(ndf) {
  ndf %>%
    mutate(event_type=as.character(event_type)) %>%
    mutate(event_type=factor(event_type, levels=unique(event_type), ordered=TRUE))
}

gp_ndf_top_n_eco_prop <- reorg(ndf_top_n_eco_prop)

ggarrange(gp2(gp_ndf_top_n_eco_prop, gp_ndf_top_n_eco_prop$property_damage,
              gp_ndf_top_n_eco_prop$event_type, 'Property Damage (Billions $US)',
              title=paste0("Top", N, " greatest economic consequence for properties")) +
          gtheme(yaxe_title=TRUE, hjust=0.5, axis.text.y=element_text()),
          ncol=1, nrow=1, common.legend = TRUE, legend = "none")

And then look at the 10 greatest economic consequences for the crops (expressed in billions of $US, hence the division by 1e+09 and rounding):

ndf_top_n_eco_crop <- ndf_eco %>%
  select(event_type, crop_damage) %>%
  arrange(desc(crop_damage)) %>%
  head(N) %>%
  mutate(crop_damage = round(crop_damage / 1e+09, 3))

ndf_top_n_eco_crop

This can be visualized with the following plot (same remarks as for property damages apply here):

gp_ndf_top_n_eco_crop <- reorg(ndf_top_n_eco_crop)

ggarrange(gp2(gp_ndf_top_n_eco_crop, gp_ndf_top_n_eco_crop$crop_damage,
              gp_ndf_top_n_eco_crop$event_type, 'Crop Damage (Billions $US)',
              title=paste0("Top", N, " greatest economic consequence for crops")) +
          gtheme(yaxe_title=TRUE, hjust=0.5, axis.text.y=element_text()),
          ncol=1, nrow=1, common.legend = TRUE, legend = "none")

Storm Data Case Study

Pascal P

30 July 2018

1 Synopsis

2 Data Processing

2.1 Loading the data

2.2 Tidy up

2.3 Most harmful meteorological events for the population health

2.4 Event type causing the greatest economic consequences

3 Results

3.1 Most harmful meteorological events for the population health

3.2 Event type causing the greatest economic consequences