##load packages and data

library(data.table)   # fast read + group-by
library(dplyr)        # tidy transforms
library(ggplot2)      # plotting
library(stringr)      # text cleaning
library(scales)       # axis labels
setwd("C:/Users/golnaz/OneDrive/Documents/Agro-ecosystem/repdata_data_StormData.csv")
file <- "repdata_data_StormData.csv"
dt <- data.table::fread(file)

# Peek at key fields
str(dt[, .(BGN_DATE, EVTYPE, FATALITIES, INJURIES,
           PROPDMG, PROPDMGEXP, CROPDMG, CROPDMGEXP)])
## Classes 'data.table' and 'data.frame':   902297 obs. of  8 variables:
##  $ 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" ...
##  $ 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  "" "" "" "" ...
##  - attr(*, ".internal.selfref")=<externalptr>

##transformations with justification

# Normalize EVTYPE
norm_evtype <- function(x){
  x <- toupper(str_trim(x))
  x <- str_replace_all(x, "\\s+", " ")
  x <- str_replace_all(x, "TSTM", "THUNDERSTORM")
  x <- str_replace_all(x, "THUNDERSTORMS", "THUNDERSTORM")
  x <- str_replace_all(x, "HVY", "HEAVY")
  x <- str_replace_all(x, "FLD", "FLOOD")
  x
}

# Convert exponent to multiplier
exp_to_mult <- function(e){
  e <- toupper(trimws(as.character(e)))
  dplyr::case_when(
    e == "K" ~ 1e3,
    e == "M" ~ 1e6,
    e == "B" ~ 1e9,
    e %in% as.character(0:9) ~ 10^as.numeric(e),
    e %in% c("", "+", "-", "?", NA) ~ 1,
    TRUE ~ 1
  )
}

dt <- dt %>%
  mutate(
    EVTYPE_CLEAN = norm_evtype(EVTYPE),
    PROP_MULT = exp_to_mult(PROPDMGEXP),
    CROP_MULT = exp_to_mult(CROPDMGEXP),
    PROP_DMG_USD = PROPDMG * PROP_MULT,
    CROP_DMG_USD = CROPDMG * CROP_MULT,
    ECON_DMG_USD = PROP_DMG_USD + CROP_DMG_USD,
    YEAR = as.integer(format(as.Date(BGN_DATE, "%m/%d/%Y %H:%M:%S"), "%Y"))
  )

##Aggregation

# Health impacts
health_by_event <- dt %>%
  group_by(EVTYPE_CLEAN) %>%
  summarise(
    fatalities = sum(FATALITIES, na.rm = TRUE),
    injuries   = sum(INJURIES,   na.rm = TRUE),
    health_total = fatalities + injuries
  ) %>%
  arrange(desc(health_total))

# Economic impacts
econ_by_event <- dt %>%
  group_by(EVTYPE_CLEAN) %>%
  summarise(
    prop = sum(PROP_DMG_USD, na.rm = TRUE),
    crop = sum(CROP_DMG_USD, na.rm = TRUE),
    econ_total = prop + crop
  ) %>%
  arrange(desc(econ_total))

# Sensitivity: more complete reporting period
health_1995 <- dt %>%
  filter(YEAR >= 1995) %>%
  group_by(EVTYPE_CLEAN) %>%
  summarise(health_total = sum(FATALITIES + INJURIES, na.rm = TRUE)) %>%
  arrange(desc(health_total))

econ_1995 <- dt %>%
  filter(YEAR >= 1995) %>%
  group_by(EVTYPE_CLEAN) %>%
  summarise(econ_total = sum(ECON_DMG_USD, na.rm = TRUE)) %>%
  arrange(desc(econ_total))

##Which events are most harmful to population health?

top_n <- 10
health_by_event %>%
  slice_head(n = top_n) %>%
  ggplot(aes(x = reorder(EVTYPE_CLEAN, health_total), y = health_total)) +
  geom_col() +
  coord_flip() +
  labs(
    title = paste("Top", top_n, "Event Types by Fatalities + Injuries (1950–2011)"),
    x = "Event Type", y = "Total (Fatalities + Injuries)"
  ) +
  scale_y_continuous(labels = comma)

##Which events have the greatest economic consequences?

econ_by_event %>%
  slice_head(n = top_n) %>%
  ggplot(aes(x = reorder(EVTYPE_CLEAN, econ_total), y = econ_total)) +
  geom_col() +
  coord_flip() +
  labs(
    title = paste("Top", top_n, "Event Types by Economic Damage (1950–2011)"),
    x = "Event Type", y = "Economic Damage (USD)"
  ) +
  scale_y_continuous(labels = dollar_format())

##Sensitivity (more complete years: 1995–2011)

econ_1995 %>%
  slice_head(n = 5) %>%
  ggplot(aes(x = reorder(EVTYPE_CLEAN, econ_total), y = econ_total)) +
  geom_col() +
  coord_flip() +
  labs(
    title = "Economic Damage (Top 5), 1995–2011",
    x = "Event Type", y = "Economic Damage (USD)"
  ) +
  scale_y_continuous(labels = dollar_format())

sessionInfo()
## R version 4.3.2 (2023-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 26100)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_Europe.utf8  LC_CTYPE=English_Europe.utf8   
## [3] LC_MONETARY=English_Europe.utf8 LC_NUMERIC=C                   
## [5] LC_TIME=English_Europe.utf8    
## 
## time zone: Europe/Paris
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] scales_1.3.0      stringr_1.5.1     ggplot2_3.4.4     dplyr_1.1.4      
## [5] data.table_1.14.8
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.4      jsonlite_1.8.7    highr_0.10        compiler_4.3.2   
##  [5] tidyselect_1.2.0  jquerylib_0.1.4   yaml_2.3.7        fastmap_1.1.1    
##  [9] R6_2.5.1          labeling_0.4.3    generics_0.1.3    knitr_1.45       
## [13] tibble_3.2.1      munsell_0.5.0     bslib_0.6.1       pillar_1.9.0     
## [17] rlang_1.1.2       utf8_1.2.4        cachem_1.0.8      stringi_1.8.2    
## [21] xfun_0.41         sass_0.4.7        cli_3.6.1         withr_3.0.2      
## [25] magrittr_2.0.3    digest_0.6.33     grid_4.3.2        rstudioapi_0.15.0
## [29] lifecycle_1.0.4   vctrs_0.6.5       evaluate_1.0.4    glue_1.6.2       
## [33] farver_2.1.1      codetools_0.2-19  fansi_1.0.5       colorspace_2.1-0 
## [37] rmarkdown_2.25    tools_4.3.2       pkgconfig_2.0.3   htmltools_0.5.7