---
title: "Details of `noaastormevents`"
author: "Brooke Anderson and Ziyu Chen"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Details of noaastormevents}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r echo = FALSE}
not_cran <- Sys.getenv("NOT_CRAN") == "true"
hasData <- requireNamespace("hurricaneexposuredata", quietly=TRUE)
if (!not_cran | !hasData) {
    knitr::opts_chunk$set(eval = FALSE)
    msg <- paste("Note: Examples in this vignette are set to not run on CRAN. If you would",
                 "like to build this vignette locally, you can do so by first setting the",
                 "environmental variable 'NOT_CRAN' to 'true' on your computer and then ",
                 "rebuilding the vignette.")
    msg <- paste(strwrap(msg), collapse="\n")
    message(msg)
}
```

```{r echo = FALSE, message = FALSE, warning = FALSE}
library(ggplot2)
library(stringr)
library(lubridate)
library(forcats)
library(dplyr)
library(tidyr)
library(viridis)
library(magrittr)

library(noaastormevents)
```


This vignette provides more details on how the `noaastormevents` package interacts with the online NOAA Storm Events database to pull storm event listings based on user queries. 

## Structure of NOAA Storm Events data

The NOAA Storm Events data is available online at https://www.ncdc.noaa.gov/stormevents/. That website includes documentation on the data, as well as a page that allows bulk data download of yearly csv files either through ftp or http (https://www.ncdc.noaa.gov/stormevents/ftp.jsp). Data is available from January 1950 and tends to be updated to within a few months of present. 

Data is stored in bulk by year in compressed comma-separated files (`.csv.gz` files). Each year has three compressed files available: 

- `StormEvents_details` file 
- `StormEvents_fatalities` file 
- `StormEvents_loations` file 

File names for each file include both the year of the data (e.g., "1950") and the date the file was last modified (e.g., "20170120"). Files are given regular names other than these two specifications. This regular naming scheme allows us to use regular expressions in code within the `noaastormevents` package on all listed file names to identify the exact name of a file for a specific year, as explained in the next section. 

The size of all three file types has increased with time (see figure below; note that the y-axis is log 10). The largest file for any given year is the "Details" file. Most file sizes increased substantially in 1996 (dotted vertical line), when the database dramatically expanded the types of events it included. Before 1996, the database covered tornadoes and, for some years, a few other types of events. From 1996, the database expanded to include events like floods, tropical storms, snow storms, etc. While "Locations" files exist in the database for early years, they contain no information until 1996. See the documentation at the NOAA Storm Events database website for more information on the coverage of the database at different times across its history.

```{r echo = FALSE, fig.width = 6, fig.height = 4, fig.align = "center"}
url <- paste0("http://www1.ncdc.noaa.gov/pub/data/swdi/",
                "stormevents/csvfiles/")
times_translation <- data_frame(times = c("K", "M", "G"),
                                times_num = c(10^3, 10^6, 10^9))
htmltab::htmltab(doc = url, which = 1, rm_nodata_cols = FALSE) %>%
  select(Name, Size) %>%
  filter(!(Size %in% "-")) %>%
  mutate(file_type = str_extract(Name, ".*?-ftp"),
         file_type = str_replace(file_type, "StormEvents_", ""),
         file_type = str_replace(file_type, "-ftp", ""),
         base = str_extract(Size, "[0-9.]+"),
         times = str_extract(Size, "[A-Z]")) %>%
  filter(file_type %in% c("details", "fatalities", "locations")) %>%
  mutate(year = str_extract(Name, "_d[0-9]{1,4}_"),
         year = str_extract(year, "[0-9]+"),
         year = as.numeric(year)) %>%
  select(file_type, base, times, year) %>%
  left_join(times_translation, by = "times") %>%
  mutate(base = as.numeric(base),
         times_num = ifelse(is.na(times_num), 1, times_num),
         file_size = base * times_num,
         file_type = str_to_title(file_type)) %>%
  select(-base, -times, -times_num) %>%
  ggplot(aes(x = year, y = file_size, color = file_type)) + 
  geom_hline(aes(yintercept = 10^3), lty = 2, col = "lightgray") + 
  annotate("text", x = 1952, y = 10^3 + 500, label = "1 KB") + 
  geom_hline(aes(yintercept = 10^6), lty = 2, col = "lightgray") + 
  annotate("text", x = 1952, y = 10^6 + 500000, label = "1 MB") + 
  geom_vline(aes(xintercept = 1996), lty = 3) + 
  geom_line() + 
  theme_classic() + 
  labs(x = "Year", y = "File size", color = "File type") +
  scale_y_log10()
```


## Downloading NOAA Storm Events data for a year

The database data is stored in files separated by year, so the file for an entire year is identified and downloaded when a user asks for event listings from any time or any type of event that year. For example, if a user wants to list flood events from the week of Hurricane Floyd in 1999, functions in the `noaastormevents` package would first identify and download the full "Details" data file for 1999 and then filter down to flood events starting in the correct week. 

To identify the online file path for a specific year, the `find_file_name` function in the `noaastormevents` package uses the `htmltab` function (from the package of the same name) to create a dataframe listing all files available for download from the NOAA Storm Events database. The function then uses regular expressions to identify the file name in that listing for the requested year. For example, the name of the file with "Details" information for 1999 can be determined with:

```{r}
find_file_name(year = "1999", file_type = "detail")
```

Here is the full definition of the `find_file_name` function:

```{r}
find_file_name
```

Typically, this function will only be used internally rather than called directly by a user.

Once the file name has been determined, a function in the package then downloads that file to the user's computer. For some years, files are very large, so this download can take a little while. To avoid downloading data from the same year more than once within an R session, the downloading function stores the downloaded data for that year in a temporary environment in the R user's session. In later requests for the same year, the function will first check for data from this year in the temporary environment and only download the data from the online database if it is not already available on the user's computer. 

This environment is created to be temporary, which means that it is deleted at the end of the current R session. While some packages that access online databases cache any downloaded data in a way that persists between R sessions, we chose not to do that and instead only cache within an R session, but delete all data at the close of the R session. This is because some of the Storm Event files are very large, and most users will likely only want to keep a small subset of the data for a given year (e.g., only flood events during the week of Hurricane Floyd). It would be wasteful of memory to cache all the 1999 data indefinitely on the user's computer in this case; instead, the user should use our package to create the desired subset of the data and then explicitly store that subset locally to use in future analysis. 

The function for downloading the file for a year is called `download_storm_data`. Here is it's full definition: 

```{r}
noaastormevents:::download_storm_data
```

Finally, the `noaastormevents` package allows a user to query storm events either by a date range or by a named historical tropical storm, rather than a year. The `create_storm_data` function inputs either a date range or a storm name, as well as the requested file type, and downloads data for the appropriate year or years. If the user requests a date range, the function will download yearly data files for all years included in that range. If the user requests a tropical storm, the function will pull the data for that particular year. Here is the full definition of `create_storm_data`:

```{r}
create_storm_data
```


As a note, many of the functions in the `noaastormevents` package that allow linking events with tropical storms rely on historical data for the storms, including storm tracks, estimated distances to eastern U.S. counties, and dates when the storm was closest to each county. The package pulls this historical data from the `hurricaneexposuredata` package, through the interfacing package `hurricaneexposure`. The hurricane data goes from 1988 to (currently) 2015 and includes all Atlantic basin tropical storms that came within 250 km of at least one U.S. county. The following storms are included in that package and so available to be used for functions in `noaastormevents`: 

```{r echo = FALSE}
library(hurricaneexposuredata)
library(hurricaneexposure)

hurr_tracks %>%
        select(storm_id) %>%
        distinct() %>%
        separate(storm_id, into = c("storm_name", "Year"), sep = "-") %>%
        group_by(Year) %>%
        summarize(Storms = paste(storm_name, collapse = ", ")) %>% 
        as.data.frame() %>% 
        pander::pander(justify = "cl", split.cells = c("20%","80%"))
```

## Structure of "Details" data files

While the `noaastormevent` package focuses on higher-level functions, which result in a simplified and cleaned version of this storm events data, a user can use the `create_storm_data` function to pull the full dataset for a year into R and work with the raw, uncleaned version. For example, here is a call that pulls the raw data for 2015 into an R object called `events_2015`. 

```{r message = FALSE}
events_2015 <- create_storm_data(date_range = c("2015-01-01", "2015-12-31"))
slice(events_2015, 1:3)
```

This raw data has 51 columns. This includes:

- the date and time of the beginning and end of the event (`BEGIN_YEARMONTH`, `BEGIN_DAY`, `BEGIN_TIME`, `BEGIN_DATE_TIME`, `END_YEARMONTH`, `END_DAY`, `END_TIME`, `END_DATE_TIME`)
- unique IDs for the event and the episode, which can be used to link listings for events and episodes across the three types of files ("details", "locations", and "fatalities") (`EPISODE_ID`, `EVENT_ID`). Note that there are more unique event IDs (57,779 for the 2015 events data) than unique episode IDs (9,511 for the 2015 events data)
- the type of event (`EVENT_TYPE`)
- the location of the event (`CZ_TYPE`, `CZ_FIPS`, `CZ_NAME`, `BEGIN_AZIMUTH`, `BEGIN_LOCATION`, `END_LOCATION`, `END_AZIMUTH`, `BEGIN_LAT`, `BEGIN_LON`, `END_LAT`, `END_LON`)
- description and impacts of the event (`EPISODE_NARRATIVE`, `EVENT_NARRATIVE`, `INJURIES_DIRECT`, `INJURIES_INDIRECT`, `DEATHS_DIRECT`, `DEATHS_INDIRECT`, `DAMAGE_PROPERTY`, `DAMAGE_CROPS`)
- information on how the event was reported (`SOURCE`)
- some metrics that are specific to event types (e.g., `FLOOD_CAUSE`, `CATEGORY`, `TOR_F_SCALE`, `TOR_LENGTH`, `TOR_WIDTH`)

### Event types

The following sections provide some summary statistics for data from this database for a single year (2015), to help users better understand the available data. Users may want to conduct similar data analysis themselves with the set of data they pull from the NOAA Storm Events database relevant to a particular research project. The code from this vignette (available at the package's GitHub repository) can serve as a starting point for that. 

In the 2015 event listings, here are the types of events and the number of reported events for each: 


```{r echo = FALSE}
events_2015 %>%
  group_by(EVENT_TYPE) %>%
  summarize(N = n()) %>%
  arrange(desc(N)) %>%
  mutate(N = prettyNum(N, big.mark = ",")) %>%
  knitr::kable(col.names = c("Event type", "Number of events in 2015"))
```

Here are how the start dates for listings for each event type are distributed over the year (event types are ordered by decreasing total count during the year; note that the y-axes vary depending on the range of events by date for each event type): 

```{r fig.width = 8, fig.height = 8, echo = FALSE, fig.align = "center"}
events_2015 %>%
  select(BEGIN_DATE_TIME, EVENT_TYPE) %>%
  mutate(date = str_extract(BEGIN_DATE_TIME, ".+\\ "),
         date = dmy(date)) %>%
  group_by(date, EVENT_TYPE) %>%
  summarize(n = n()) %>%
  ungroup() %>%
  group_by(EVENT_TYPE) %>%
  mutate(tot_n = sum(n)) %>%
  ungroup() %>%
  arrange(desc(tot_n)) %>%
  mutate(EVENT_TYPE = factor(EVENT_TYPE, levels = unique(EVENT_TYPE))) %>%
  ggplot(aes(x = date, y = n)) + 
  geom_point(alpha = 0.3, size = 1) + 
  facet_wrap(~ EVENT_TYPE, scales = "free_y", ncol = 5) + 
  scale_x_date(date_labels = "%m/%d", date_breaks = "4 months") + 
  theme_classic()
```

Many event types are clearly seasonal (e.g., winter weather, winter storms, heavy snow, cold, extreme cold, blizzards, ice storms, lake-effect snow, and avalanches are all much more common during winter months, while tropical depressions and tropical storms are all limited to the hurricane season). However, for some events, reporting seasonal patterns might be based not just on the true pattern of events but also on the timing of important exposures and impacts of the events. For example, rip currents have many more listings during the spring and summer, which may be related to events being more likely to be listed when more people are swimming. Frost event listings are particularly high at the start and end of the frost season, rather than in the middle of winter, which may be related to the impacts of frost on crops being higher in spring and fall than during the winter. If working with this data, it important to keep in mind that the data are based on reporting, and there may be related influences on the probability of an event being reported and included in the data that differ from using data from something like a weather station. 

### Episodes versus events

"Episodes" seem to collect related "events", where events can vary in the type or location of the event, while an "episode" collects events that belong to the same large system. The following graph shows, for each episode listed in 2015, the number of events listed for the episode (x-axis) and the size (in days) of the range of begin dates across events in the episode.

```{r fig.height = 3.5, fig.width = 5, echo = FALSE, fig.align = "center"}
events_2015 %>%
  mutate(date = substring(BEGIN_DATE_TIME, 1, 9),
         date = dmy(date),
         yday = yday(date)) %>%
  group_by(EPISODE_ID) %>%
  summarize(n = n(),
            day_diff = max(yday) - min(yday)) %>%
  ggplot(aes(x = n, y = day_diff)) +
  geom_jitter(height = 0.1, alpha = 0.25) + 
  theme_classic() + 
  labs(x = "Number of events per episode", 
       y = "Largest difference in begin dates\nacross events in the episode")
```

An episode will never include events in more than one state, so a large weather system could potentially be described by multiple episodes in different states: 

```{r}
events_2015 %>% 
  select(EPISODE_ID, STATE) %>% 
  group_by(EPISODE_ID) %>% 
  summarize(n_states = length(unique(STATE))) %>% 
  ungroup() %>% 
  summarize(max_n_states = max(n_states))
```

Here are maps with the beginning locations of events in the episodes with the most events in 2015. Note that the beginning latitude and longitude are not listed for every event, resulting in one of the episodes not having any points on the map. From the other maps, it is clear that events within the episode were fairly close together. 

```{r message = FALSE, warning = FALSE, echo = FALSE, fig.align = "center", fig.width = 6, fig.height = 10}
top_episodes <- events_2015 %>%
  group_by(EPISODE_ID) %>%
  summarize(n = n()) %>%
  arrange(desc(n)) %>%
  slice(1:10)

events_2015 %>%
  filter(EPISODE_ID %in% top_episodes$EPISODE_ID) %>%
  mutate(EPISODE_ID = factor(EPISODE_ID, levels = unique(top_episodes$EPISODE_ID))) %>%
  ggplot(aes(x = BEGIN_LON, y = BEGIN_LAT)) +  
  geom_polygon(data= map_data("state"),
               aes(x = long, y = lat, group = group),
               colour = "black", fill = NA) + 
  geom_point(color = "red", alpha = 0.3, size = 0.5) + 
  facet_wrap(~ EPISODE_ID, ncol = 2) +
  coord_map() + 
  theme_void()
```

For these episodes with the most events in 2015, the following graph shows the number of events reported for the episode. One of the episodes was a winter storm, another was heavy rains and floods, while the rest of the episodes included high winds, hail, tornadoes, rain, and / or flooding. 

```{r echo = FALSE, fig.align = "center", fig.width = 6, fig.height = 4}
top_episodes <- events_2015 %>%
  group_by(EPISODE_ID) %>%
  summarize(n = n()) %>%
  arrange(desc(n)) %>%
  slice(1:10)

a <- events_2015 %>%
  filter(EPISODE_ID %in% top_episodes$EPISODE_ID) %>%
  group_by(EPISODE_ID, EVENT_TYPE) %>%
  summarize(n = n()) %>%
  ungroup() %>%
  spread(key = EPISODE_ID, value = n, fill = 0) %>%
  tibble::remove_rownames() %>%
  as.data.frame() %>%
  tibble::column_to_rownames(var = "EVENT_TYPE") %>%
  dist("canberra") %>%
  hclust() 
```

```{r echo = FALSE, fig.align = "center", fig.width = 6, fig.height = 3.5}
events_2015 %>%
  filter(EPISODE_ID %in% top_episodes$EPISODE_ID) %>%
  group_by(EPISODE_ID, EVENT_TYPE) %>%
  summarize(n = n()) %>%
  ungroup() %>%
  mutate(EPISODE_ID = factor(EPISODE_ID, levels = unique(top_episodes$EPISODE_ID)),
         EVENT_TYPE = factor(EVENT_TYPE, levels = unique(a$labels[a$order]))) %>%
  ggplot(aes(x = EPISODE_ID, y = EVENT_TYPE)) + 
  geom_tile(aes(fill = n)) + 
  theme_classic() + 
  labs(x = "", y = "", fill = "Number of events reported") + 
  theme(axis.text.x  = element_text(angle=45, vjust = 1, hjust = 0),
        legend.position = "top",
        plot.margin = unit(c(0, 1, 0, 0), "cm")) +
  scale_x_discrete(position = "top") +
  scale_fill_viridis(option = "A", direction = -1)
```

Once we removed event types with less that 50 listings in 2015, we did a cluster analysis of event types, to group events that are more likely to occur together within an episode. The following plot shows the resulting cluster structure of these event types.

```{r echo = FALSE, fig.align = "center", fig.width = 6, fig.height = 10}
b <- events_2015 %>%
  group_by(EVENT_TYPE) %>%
  mutate(n_events = n()) %>%
  ungroup() %>%
  filter(n_events >= 50) %>%
  group_by(EPISODE_ID, EVENT_TYPE) %>%
  summarize(n = n()) %>%
  ungroup() %>%
  spread(key = EPISODE_ID, value = n, fill = 0) %>%
  tibble::remove_rownames() %>%
  as.data.frame() %>%
  tibble::column_to_rownames(var = "EVENT_TYPE") %>%
  dist("canberra") %>%
  hclust() 
plot(as.dendrogram(b, hang = 0.25), horiz = TRUE, axes = FALSE)
```

The next graph shows the number of events of each event type (excluding event types with less than 50 total listings in 2015). Each row represents an episode.

```{r echo = FALSE, fig.height = 6, fig.width = 8, fig.align = "center"}
events_2015 %>%
  group_by(EVENT_TYPE) %>%
  mutate(n_events = n()) %>%
  ungroup() %>%
  filter(n_events >= 100) %>%
  group_by(EPISODE_ID, EVENT_TYPE) %>%
  summarize(n = n()) %>%
  ungroup() %>%
  mutate(EPISODE_ID = factor(EPISODE_ID),
         EVENT_TYPE = factor(EVENT_TYPE, levels = unique(b$labels[b$order]))) %>%
  ggplot(aes(x = EVENT_TYPE, y = EPISODE_ID)) + 
  geom_tile(aes(fill = n)) + 
  theme_classic() + 
  theme(axis.text.x  = element_text(angle = 45, vjust = 1, hjust = 0),
        legend.position = "top",
        plot.margin = unit(c(0, 1, 0, 0), "cm"),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank()) +
  scale_x_discrete(position = "top") +
  scale_fill_viridis(option = "A", direction = -1, begin = 0.1, end = 0.9) + 
  labs(x = "", y = "", fill = "Number of reported events in the episode")
```

### How events are reported

The `SOURCE` column of the raw data gives information on how each event was reported. 

```{r echo = FALSE}
events_2015 %>%
  group_by(SOURCE) %>%
  summarize(N = n()) %>%
  arrange(desc(N)) %>%
  mutate(N = prettyNum(N, big.mark = ",")) %>%
  knitr::kable(col.names = c("Source of event report", "Number of events reported in 2015"))
```

The majority of events in this database, at least for 2015, were reported by either a trained spotter or the public.

The following graph shows, for each type of event in 2015, the percent reported by each source. For some types of events, reporting is dominated by a specific source. For example, most high surf reports come from trained spotters, while most drought reports come from drought monitors and most tornado reports come from the NWS Storm Survey. For other types of events, reporting sources are more diversified. Both axes of the plot are ordered by overall frequency (i.e., overall number of each type of event and overall number of reports from each source).

```{r fig.height = 7, fig.width = 10, echo = FALSE, message = FALSE, warning = FALSE, fig.align = "center"}
tot_events <- events_2015 %>%
  group_by(EVENT_TYPE) %>%
  summarize(tot_events = n()) %>%
  arrange(desc(tot_events))
tot_source <- events_2015 %>%
  group_by(SOURCE) %>%
  summarize(tot_source = n()) %>%
  arrange(tot_source)

events_2015 %>%
  group_by(EVENT_TYPE, SOURCE) %>%
  summarize(N = n()) %>%
  ungroup() %>%
  right_join(tot_events, by = "EVENT_TYPE") %>%
  mutate(perc_events = 100 * N / tot_events) %>%
  mutate(EVENT_TYPE = factor(EVENT_TYPE, levels = unique(EVENT_TYPE)),
         SOURCE = factor(SOURCE, levels = unique(tot_source$SOURCE))) %>%
  ggplot(aes(y = SOURCE, x = EVENT_TYPE)) + 
  geom_tile(aes(fill = perc_events)) + 
  theme_classic() + 
  labs(x = "", y = "", fill = "% events reported by a given source") + 
  theme(axis.text.x  = element_text(angle=45, vjust = 1, hjust = 0),
        legend.position = "top",
        plot.margin = unit(c(0, 1, 0, 0), "cm")) +
  scale_x_discrete(position = "top") +
  scale_fill_viridis(option = "C", direction = -1)
```

### Event locations

Each event has a state listed for the event (`STATE`). The following graph gives the number of reported events in each state for 2015:

```{r echo = FALSE, fig.height = 8, fig.width = 6, fig.align = "centers"}
events_2015 %>%
  group_by(STATE) %>%
  summarize(n = n()) %>%
  mutate(STATE = str_to_title(STATE)) %>%
  ggplot(aes(x = n, y = fct_reorder(STATE, n))) + 
  geom_point() + 
  labs(x = "Number of listed events in 2015", y = "") + 
  theme_classic()
```

Note that "states" include bodies of water (e.g., specific Great Lakes, the Hawaii waters, the Gulf of Mexico) and territories (American Samoa, Guam, Puerto Rico, Virgin Islands).

For some event types, the latitude and longitude of the beginning of the event is included with the event listing. 

```{r echo = FALSE}
events_2015 %>%
  group_by(EVENT_TYPE) %>%
  summarize(n = n(),
            non_miss = sum(!is.na(BEGIN_LAT)),
            perc_non_miss = non_miss / n,
            perc_non_miss = factor(perc_non_miss, levels = c(1, 0),
                                   labels = c("With latitude / longitude",
                                              "Without latitude / longitude"))) %>%
  ungroup() %>%
  group_by(perc_non_miss) %>%
  summarize(with_loc = paste(EVENT_TYPE, collapse = ", ")) %>%
  rename(` ` = perc_non_miss, 
         `Event types` = with_loc) %>%
  as.data.frame() %>% 
  pander::pander(split.cells = c(20, 60),
                 justify = "rl")
```

Of the 2015 events with a latitude and longitude listed for the beginning of the event and that are in the continental U.S., here are those locations by month: 

```{r echo = FALSE, fig.width = 9, fig.height = 6, fig.align = "center"}
events_2015 %>%
  filter(!is.na(BEGIN_LAT) & !(is.na(BEGIN_LON))) %>%
  filter(!(STATE %in% c("GULF OF MEXICO", "ALASKA", "GUAM", "ATLANTIC NORTH",
                        "LAKE HURON", "LAKE ST CLAIR", "HAWAII", "AMERICAN SAMOA",
                        "LAKE SUPERIOR", "ATLANTIC SOUTH", "LAKE MICHIGAN", 
                        "HAWAII WATERS", "PUERTO RICO", "E PACIFIC", "LAKE ERIE",
                        "LAKE ONTARIO", "VIRGIN ISLANDS"))) %>%
  mutate(MONTH_NAME = factor(MONTH_NAME, levels = month.name)) %>%
  ggplot(aes(x = BEGIN_LON, y = BEGIN_LAT)) + 
  geom_polygon(data = map_data("usa"), aes(x = long, y = lat, group = group),
               color = "black", fill = NA) + 
  geom_point(size = 1, alpha = 0.15, color = "darkred") + 
  coord_map() + 
  theme_void() + 
  facet_wrap(~ MONTH_NAME)
```

Here are those locations by event type: 

```{r echo = FALSE, fig.width = 9, fig.height = 6, fig.align = "center"}
events_2015 %>%
  filter(!is.na(BEGIN_LAT) & !(is.na(BEGIN_LON))) %>%
  filter(!(STATE %in% c("GULF OF MEXICO", "ALASKA", "GUAM", "ATLANTIC NORTH",
                        "LAKE HURON", "LAKE ST CLAIR", "HAWAII", "AMERICAN SAMOA",
                        "LAKE SUPERIOR", "ATLANTIC SOUTH", "LAKE MICHIGAN", 
                        "HAWAII WATERS", "PUERTO RICO", "E PACIFIC", "LAKE ERIE",
                        "LAKE ONTARIO", "VIRGIN ISLANDS"))) %>%
  group_by(EVENT_TYPE) %>%
  mutate(n = n()) %>%
  ungroup() %>%
  mutate(EVENT_TYPE = fct_reorder(EVENT_TYPE, n, .desc = TRUE)) %>%
  ggplot(aes(x = BEGIN_LON, y = BEGIN_LAT)) + 
  geom_polygon(data = map_data("usa"), aes(x = long, y = lat, group = group),
               color = "black", fill = NA) + 
  geom_point(size = 0.8, alpha = 0.15, color = "darkred") + 
  coord_map() + 
  theme_void() + 
  facet_wrap(~ EVENT_TYPE)
```

Some events have different latitudes and longitudes for the beginning and ending locations. For example, here are maps for one state (Arkansas) of events with different starting and ending locations:

```{r echo = FALSE, fig.width = 7, fig.height = 5, fig.align = "center"}
events_2015 %>%
  filter(!is.na(BEGIN_LAT) & !(is.na(BEGIN_LON)) & 
        !is.na(END_LAT) & !(is.na(END_LON))) %>%
  mutate(diff_lon = BEGIN_LON - END_LON,
         diff_lat = BEGIN_LAT - END_LAT) %>%
  filter(abs(diff_lon) > 0.1 | abs(diff_lat) > 0.1) %>%
  filter(STATE == "ARKANSAS") %>%
  ggplot(aes(x = BEGIN_LON, y = BEGIN_LAT)) + 
  geom_polygon(data = map_data("county", region = "arkansas"), 
               aes(x = long, y = lat, group = group),
               color = "lightgray", fill = NA) + 
  geom_polygon(data = map_data("state", region = "arkansas"),
               aes(x = long, y = lat, group = group),
                   color = "black", fill = NA) + 
  geom_segment(aes(xend = END_LON, yend = END_LAT), color = "red", size = 1) + 
  theme_void() + 
  coord_map() + 
  facet_wrap(~ EVENT_TYPE, ncol = 3)
```

Some events are reported by forecast zone (`CZTYPE` of "Z") rather than county (`CZTYPE` of "C"). Specific types of events are typically either always reported for a county or always reported for a forecast zone (see table below). Events typically reported by county include floods ("Flash Flood", "Flood", "Debris Flow"), tornado-like events ("Tornado", "Funnel Cloud", "Dust Devil"), and a few other events often related to thunderstorms ("Thunderstorm Wind", "Hail", "Heavy Rain", "Lightning"). Events typically reported by forecast zone include severe winter weather ("Winter Weather", "Winter Storm", "Heavy Snow", "Cold/Wind Chill", "Extreme Cold/Wind Chill", "Blizzard", "Frost/Freeze", "Ice Storm", "Sleet", "Lake-Effect Snow", "Avalanche", "Freezing Fog"), extreme heat ("Heat", "Excessive Heat", "Drought"), events related to the water or coast ("Marine Thunderstorm Wind", "High Surf", "Coastal Flood", "Waterspout", "Astronomical Low Tide", "Rip Current", "Tropical Storm", "Marine High Wind", "Marine Hail", "Marine Strong Wind", "Hurricane", "Seiche", "Storm Surge/Tide", "Tropical Depression", "Marine Dense Fog", "Sneakerwave", "Tsunami"), and a few others ("High Wind", "Dense Fog", "Strong Wind", "Wildfire", "Dust Storm", "Dense Smoke").

```{r echo = FALSE}
events_2015 %>%
  group_by(CZ_TYPE, EVENT_TYPE) %>%
  summarize(n_events = n()) %>%
  ungroup() %>%
  spread(key = CZ_TYPE, value = n_events, fill = 0) %>%
  mutate(total = C + Z,
         perc_county = 100 * C / total) %>%
  arrange(desc(perc_county), desc(total)) %>%
  mutate(C = prettyNum(C, big.mark = ","),
         Z = prettyNum(Z, big.mark = ","),
         total = prettyNum(total, big.mark = ","),
         perc_county = paste0(round(perc_county), "%")) %>%
  knitr::kable(col.names = c("Event type", "County", "Forecast Zone", "Total", "% county"),
               align = "lcccc")
```

For events reported by county, here are maps showing distributions in the number of events reported in 2015:

```{r message = FALSE, warning = FALSE, fig.width = 9, fig.height = 6, fig.align = "center", echo = FALSE}
events_2015 %>%
  filter(CZ_TYPE == "C") %>%
  filter(EVENT_TYPE != "Heat") %>%
  filter(!(STATE %in% c("GULF OF MEXICO", "ALASKA", "GUAM", "ATLANTIC NORTH",
                        "LAKE HURON", "LAKE ST CLAIR", "HAWAII", "AMERICAN SAMOA",
                        "LAKE SUPERIOR", "ATLANTIC SOUTH", "LAKE MICHIGAN", 
                        "HAWAII WATERS", "PUERTO RICO", "E PACIFIC", "LAKE ERIE",
                        "LAKE ONTARIO", "VIRGIN ISLANDS"))) %>%
  mutate(CZ_FIPS = str_pad(CZ_FIPS, width = 3, pad = "0"),
         STATE_FIPS = str_pad(STATE_FIPS, width = 2, pad = "0")) %>%
  unite(fips, STATE_FIPS, CZ_FIPS, sep = "") %>%
  mutate(fips = as.numeric(fips)) %>%
  group_by(fips, EVENT_TYPE) %>%
  dplyr::summarize(value = n()) %>%
  ungroup() %>%
  left_join(maps::county.fips, by = "fips") %>%
  left_join(unite(map_data("county"), polyname, region, subregion, sep = ","),
            by = "polyname") %>%
  ggplot(aes(x = long, y = lat, group = group, fill = value)) +
  geom_polygon(color = NA) + 
  geom_polygon(data = map_data("usa"), 
               aes(x = long, y = lat, group = group),
               color = "black", fill = NA) + 
  theme_void() + 
  coord_map() + 
  scale_fill_viridis(option = "A", direction = -1) + 
  facet_wrap(~ EVENT_TYPE) + 
  theme(legend.position = "top") + 
  labs(fill = "Number of events")
```

Here is a sample of events that are instead reported by forecast zone, with the state, `CZ_NAME`, and event narrative included. Note that the county name is often provided by the `CZ_NAME` column, although the `CZ_FIPS` value is the forecast zone for any event listed by forecast zone. We use code to try to match `STATE` and `CZ_NAME` listings to a table of U.S. county names and associated county FIPS for each event listed by forecast zone, to allow these events to be included in event listings and maps created by functions in `noaastormevents`.

```{r echo = FALSE}
events_2015 %>%
  filter(CZ_TYPE == "Z" & EVENT_NARRATIVE != "") %>%
    filter(!(STATE %in% c("GULF OF MEXICO", "ALASKA", "GUAM", "ATLANTIC NORTH",
                        "LAKE HURON", "LAKE ST CLAIR", "HAWAII", "AMERICAN SAMOA",
                        "LAKE SUPERIOR", "ATLANTIC SOUTH", "LAKE MICHIGAN", 
                        "HAWAII WATERS", "PUERTO RICO", "E PACIFIC", "LAKE ERIE",
                        "LAKE ONTARIO", "VIRGIN ISLANDS"))) %>%
  select(STATE, CZ_NAME, EVENT_NARRATIVE) %>%
  sample_n(8) %>%
  mutate(STATE = str_to_title(STATE),
         CZ_NAME = str_to_title(CZ_NAME)) %>%
  as.data.frame() %>% 
  pander::pander(split.cells = c(15, 15, 50), justify = "ccl")
```

Within the code in `noaastormevents`, the `match_forecast_county` is used to try to match a county FIPS to each of the events listed by forecast zone. To get the full code for that function, you can run `match_forecast_county` (i.e., the function name, without parentheses after). To match an events listed by forecast zone to a county, this function tries the following to try to match all or part of the `state_name` and `cz_name` columns in the storm events data to the state and county names in the `county.fips` dataframe that comes with the `maps` package: 

1. Tries to match `cz_name` to the county name in `county.fips` after removing any periods or apostrophes in `cz_name`.
2. Next, for county names with 'county' in them, try to match the word before 'county'
   to county name in `county.fips`. Then check the two words before 'county', then the
   one and two words before 'counties'.
3. Next, pull out the last word in `cz_name` and try to match it to the county name
  in `county.fips`. Then check the last two words in `cz_name`, then check the last three words in `cz_name`.
4. Next, pull any words right before a slash and check that against the county name.
5. Finally, try removing anything in parentheses in `cz_name` before matching.

In addition, there are a few final steps in cleaning the data. First, all listings with "Utah" in the `cz_name` are set to missing-- while there is a Utah County, Utah, from inspection of event listings in 2015, events in Utah with "Utah" in `cz_name` often referred to parts of the state, rather than to the county. Further, any event with "National Park" listed in `cz_name` is set to not match with a county FIPS. In Wyoming, Park County was being matched to a `cz_name` for Yellowstone National Park, and this could be a problem in other states, so this extra check was included. 

```{r echo = FALSE}
z_events_2015 <- events_2015 %>%
  dplyr::select_(~ BEGIN_YEARMONTH, ~ BEGIN_DAY, ~ END_YEARMONTH, ~ END_DAY,
                   ~ EPISODE_ID, ~EVENT_ID, ~ STATE, ~ CZ_TYPE, ~ CZ_NAME,
                    ~ DEATHS_DIRECT, ~ DEATHS_INDIRECT, ~ INJURIES_DIRECT,
                  ~ INJURIES_INDIRECT, ~ DAMAGE_PROPERTY, ~ DAMAGE_CROPS,
                   ~ EVENT_TYPE, ~ STATE_FIPS, ~ CZ_FIPS, ~ SOURCE,
                   ~ EPISODE_NARRATIVE, ~ EVENT_NARRATIVE) %>%
  setNames(tolower(names(.))) %>%
  dplyr::filter_(~ cz_type == "Z") %>%
  match_forecast_county()
```

In the 2015 events data, there are 22,664 events listed by forecast zone rather than county. Once the `match_forecast_county` function is applied to these events, 16,405 events were linked to a county (72%) while 6,259 events (28%) could not be matched to a county.

```{r echo = FALSE}
outside_events <- z_events_2015 %>%
  filter(state %in% c("GULF OF MEXICO", "GUAM", "ATLANTIC NORTH",
                                    "LAKE HURON", "LAKE ST CLAIR", "AMERICAN SAMOA",
                                    "LAKE SUPERIOR", "ATLANTIC SOUTH", "LAKE MICHIGAN",
                                    "HAWAII WATERS", "PUERTO RICO", "E PACIFIC", "LAKE ERIE",
                                    "LAKE ONTARIO", "VIRGIN ISLANDS", "HAWAII", "ALASKA")) %>%
  filter(is.na(fips)) 
```

Of the events not matched to a county, 2,997 events were outside the continental U.S. (i.e., in Hawaii, Alaska, U.S. territories, or waters): 

```{r echo = FALSE}
outside_events %>%
  mutate(state = stringr::str_to_title(state)) %>%
  group_by(state) %>%
  dplyr::summarize(n = n()) %>%
  arrange(desc(n)) %>%
  knitr::kable(col.names = c("State", "Number of events listed by forecast zone"))
```

```{r echo = FALSE}
continental_events <- z_events_2015 %>%
  filter(!(state %in% c("GULF OF MEXICO", "GUAM", "ATLANTIC NORTH",
                                    "LAKE HURON", "LAKE ST CLAIR", "AMERICAN SAMOA",
                                    "LAKE SUPERIOR", "ATLANTIC SOUTH", "LAKE MICHIGAN",
                                    "HAWAII WATERS", "PUERTO RICO", "E PACIFIC", "LAKE ERIE",
                                    "LAKE ONTARIO", "VIRGIN ISLANDS", "HAWAII", "ALASKA"))) %>%
  filter(is.na(fips)) 
```

This left 3,262 events in the continental U.S. that were listed by forecast zone but could not be linked to a county by the `match_forecast_county` function. Most of these events had a value for `cz_name` with words related to mountains, water, adjacency (e.g., "area of", "vicinity of"), or a few other word types (e.g., "desert", "hwy"). The following table gives the number of these remaining unmatched events with words from each of these categories (note: the `cz_name` for an event may have words from more than one of these categories, in which case it would be counted in this table under both categories). 

```{r echo = FALSE}
mountain_words <- c("mountains", "mountain", "valley", "valleys", "divide", "range",
                    "ridge", "ridges", "gap", "front", "cascades", "foothills", "plateau",
                    "highlands", "mtns", "hills", "mesa", "mesas", "vly", "peak",
                    "gorge", "rim", "arches", "hole", "canyon")
water_words <- c("beaches", "beach", "river", "basin", "coast", "delta", "lake",
                 "shoreline", "bay")
adjacent_words <- c("area", "region", "vicinity", "country")
other_words <- c("desert", "hwy", "panhandle", "corridor", "nation", "plains", 
                 "lowlands", "prarie")

continental_events <- continental_events %>%
  mutate(cz_name = str_to_lower(cz_name),
         has_mountain_words = str_detect(cz_name, paste0(mountain_words, collapse = "|")),
         has_water_words = str_detect(cz_name, paste0(water_words, collapse = "|")),
         has_adjacent_words = str_detect(cz_name, paste0(adjacent_words, collapse = "|")),
         has_other_words = str_detect(cz_name, paste0(other_words, collapse = "|")))

continental_events %>% 
  summarize(mountain_words = sum(has_mountain_words),
            water_words = sum(has_water_words),
            adjacent_words = sum(has_adjacent_words),
            other_words = sum(has_other_words)) %>% 
  gather(key = word_type, value = n_events) %>%
  mutate(words = c(paste(mountain_words, collapse = ", "),
                   paste(water_words, collapse = ", "),
                   paste(adjacent_words, collapse = ", "),
                   paste(other_words, collapse = ", "))) %>%
  select(word_type, words, n_events) %>%
  mutate(word_type = str_to_title(str_replace(word_type, "_words", ""))) %>%
  dplyr::rename(`Word type` = word_type,
                Words = words,
                `# events` = n_events) %>%
  as.data.frame() %>% 
  pander::pander(split.cells = c(10, 50, 10), big.mark = ",", justify = "llc")
```

The following table summarizes the number of events that could not be linked to a county that contained at least one of each of these types of words:

```{r echo = FALSE}
continental_events %>% 
  dplyr::group_by(event_type) %>%
  dplyr::summarize(mountain_words = sum(has_mountain_words),
                   water_words = sum(has_water_words),
                   adjacent_words = sum(has_adjacent_words),
                   other_words = sum(has_other_words)) %>% 
  dplyr::arrange(desc(mountain_words)) %>%
  knitr::kable(col.names = c("Event type", "Mountain words", "Water words", "Adjacent words",
                             "Other words"))
```

For 2015, here are the `cz_name` values for the events that could not be matched to a county and did not include the words listed above:

```{r echo = FALSE}
continental_events %>% 
  filter(!has_mountain_words & !has_water_words & !has_adjacent_words & !has_other_words) %>%
  select(cz_name, state) %>% 
  mutate(cz_name = str_to_title(cz_name),
         state = str_to_title(state)) %>%
  distinct() %>%
  knitr::kable()
```

For the events listed by forecast zone that could be successfully matched to a county, here are the geographic distributions in event counts in 2015:

```{r message = FALSE, warning = FALSE, fig.width = 9, fig.height = 6, fig.align = "center", echo = FALSE}
z_events_2015 %>%
  filter(!is.na(fips)) %>%
  filter(!(state %in% c("GULF OF MEXICO", "ALASKA", "GUAM", "ATLANTIC NORTH",
                        "LAKE HURON", "LAKE ST CLAIR", "HAWAII", "AMERICAN SAMOA",
                        "LAKE SUPERIOR", "ATLANTIC SOUTH", "LAKE MICHIGAN", 
                        "HAWAII WATERS", "PUERTO RICO", "E PACIFIC", "LAKE ERIE",
                        "LAKE ONTARIO", "VIRGIN ISLANDS"))) %>%
  group_by(event_type) %>%
  dplyr::mutate(n_events = n()) %>%
  ungroup() %>%
  filter(n_events >= 100) %>%
  select(-n_events) %>%
  group_by(fips, event_type) %>%
  dplyr::summarize(value = n()) %>%
  ungroup() %>%
  left_join(maps::county.fips, by = "fips") %>%
  left_join(unite(map_data("county"), polyname, region, subregion, sep = ","),
            by = "polyname") %>%
  ggplot(aes(x = long, y = lat, group = group, fill = value)) +
  geom_polygon(color = NA) + 
  geom_polygon(data = map_data("usa"), 
               aes(x = long, y = lat, group = group),
               color = "black", fill = NA) + 
  theme_void() + 
  coord_map() + 
  scale_fill_viridis(option = "A", direction = -1) + 
  facet_wrap(~ event_type) + 
  theme(legend.position = "top") + 
  labs(fill = "Number of events")
```

When using this function, and for event listings generated using this function, the user may want to hand-check that event listings with names like "Lake" and "Mountain" in the `cz_name` column are not erroneously matched to counties with names like "Lake County" and "Mountain County". Code like the following can be used for these checks (in this case, checking a dataframe of event listings named `z_events_2015` that is output from `match_forecast_county` and so has `fips` added for each event listing, if a match could be found): 

```{r eval = FALSE}
z_events_2015 %>% 
  filter(cz_type == "Z") %>% 
  select(cz_name, state, fips) %>% 
  mutate(cz_name = str_to_title(cz_name)) %>% 
  filter(str_detect(cz_name, "Lake") & !is.na(fips))  %>% 
  distinct()
```

### Impact measurements

The "Details" datasets for each year include six measurements of the impacts of each event: 

- `DEATHS_DIRECT`
- `DEATHS_INDIRECT`
- `INJURIES_DIRECT`
- `INJURIES_INDIRECT`
- `DAMAGE_PROPERTY`
- `DAMAGE_CROPS`

```{r echo = FALSE}
c_events_2015 <- events_2015 %>%
  dplyr::select_(~ BEGIN_YEARMONTH, ~ BEGIN_DAY, ~ END_YEARMONTH, ~ END_DAY,
                   ~ EPISODE_ID, ~EVENT_ID, ~ STATE, ~ CZ_TYPE, ~ CZ_NAME,
                  ~ DEATHS_DIRECT, ~ DEATHS_INDIRECT, ~ INJURIES_DIRECT,
                  ~ INJURIES_INDIRECT, ~ DAMAGE_PROPERTY, ~ DAMAGE_CROPS,
                   ~ EVENT_TYPE, ~ STATE_FIPS, ~ CZ_FIPS, ~ SOURCE,
                   ~ EPISODE_NARRATIVE, ~ EVENT_NARRATIVE) %>%
  setNames(tolower(names(.))) %>%
  dplyr::filter_(~ cz_type == "C") %>%
  dplyr::mutate_(fips = ~ as.numeric(paste0(state_fips, str_pad(cz_fips, 3, pad = "0"))))

county_events_2015 <- bind_rows(c_events_2015, filter(z_events_2015, !is.na(fips)))
```

Many of the impact values are given using abbreviations for amounts. For example, this listing for a tornado gives the property damage as "5.00K":

```{r}
county_events_2015 %>% filter(event_id == "582970")
```

Therefore, the `noaastormevents` package uses a function to pull out these abbreviations and convert associated impact values to numeric values (e.g., 5000 for "5.00K"). The conversions conducted are: 

- "K" = 1,000
- "M" = 1,000,000
- "B" = 1,000,000,000
- "T" = 1,000,000,000,000

This is done by the function `parse_damage`:

```{r}
parse_damage
```


```{r echo = FALSE}
county_damage_2015 <- county_events_2015 %>%
  dplyr::select(episode_id , event_id, fips, event_type, contains("damage")) %>% 
  tidyr::gather(key = impact, value = amount, -event_type, -event_id, -episode_id, -fips) %>%
  dplyr::mutate(amount = parse_damage(amount)) 
```

Usually, it seems that the cost per events within an episode do not overlap, so that the costs from all events in an episode can be summed to generate a total damage cost. While in some cases, different event listings within the same county in the same episode have the same damage cost, it often seems that the total cost was divided in these cases across events, as in this case, where the total estimated cost of $6,186,909 in the county (see the narratives) was spread evenly over the two event listings:

```{r echo = FALSE}
county_events_2015 %>%
  filter(episode_id == "100823" & fips == 45035) %>%
  select(event_type, contains("damage"), cz_name, state, event_narrative, episode_narrative) %>% 
  as.data.frame() %>% 
  pander::pander()
```

However, there are a few cases where it appears that costs might be duplicated over different events within the same episode in a county. For example, in the following listing, it looks like the total estimated damage cost of $1.7 million (see that narrative) is repeated over the two listed events, which would mean that adding damages across events in the episode would lead to a county damage cost of twice the true amount. 

```{r echo = FALSE}
county_events_2015 %>%
  filter(episode_id == "99142" & fips == 2023) %>%
  select(event_type, contains("damage"), cz_name, state, event_narrative, episode_narrative) %>% 
  as.data.frame() %>% 
  pander::pander()
```

```{r echo = FALSE}
a <- county_damage_2015 %>%
  group_by(event_id) %>%
  summarize(episode_id = first(episode_id),
            fips = first(fips),
            total_amount = sum(amount)) %>%
  ungroup() %>% 
  group_by(episode_id, fips) %>% 
  summarize(n_events = n(),
            average_amount = mean(total_amount),
            max_amount = max(total_amount),
            range_amount = diff(range(total_amount))) %>%
  filter(!is.na(max_amount) & max_amount > 0 & n_events > 1) %>% 
  ungroup() %>%
  mutate(zero_range = range_amount == 0) %>%
  group_by(zero_range) %>%
  count()
```

Of the times when a county had two or more events as part of the same episode in 2015, in 3,136 (87%) cases the damage costs were not identical across the event listings for the county and episode, while in 452 cases (13%) the cost listings were identical across event listings (as shown in the cases just given). If using damage cost estimates for research, it may be wise to hand-check for cases where damage cost estimates are duplicated across different events in a way that prevents summing across events to get the cumulative episode total cost. 

When costs were summed over all event and county listings in an episode, here are the ten episodes in the 2015 dataset with the highest total damage costs (combining property and crop damage):

```{r echo = FALSE, warning = FALSE, message = FALSE}
county_damage_2015 %>%
  left_join(county_events_2015 %>% select(fips, cz_name, state, episode_narrative, event_id)) %>%
  unite(location, cz_name, state, sep = ", ") %>% 
  group_by(episode_id) %>%
  summarize(location = paste(unique(stringr::str_to_title(location)), collapse = "; "),
            episode_damage = sum(amount, na.rm = TRUE),
            event_types = paste(unique(event_type), collapse = ", "),
            episode_narrative = first(episode_narrative)) %>%
  ungroup() %>%
  arrange(desc(episode_damage)) %>%
  mutate(episode_damage = formatC(episode_damage, format = "d", big.mark = ","),
         episode_damage = paste0("$", episode_damage)) %>%
  slice(1:10) %>%
  select(event_types, episode_damage, location, episode_narrative) %>%
  rename(`Event types` = event_types,
         `Total damage` = episode_damage,
         Counties = location,
         Narrative = episode_narrative) %>% 
  as.data.frame() %>% 
  pander::pander(split.cells = c(5, 5, 15, 35), justify = "lcll") 
```


Here is a table of the total number of events with non-zero damage (either for property damage or crop damage), the total damage costs across all events of that type in 2015, and the median and maximum damage per event for events of that type with non-zero damage costs. Event types are ordered by total damage costs summed across all events. 

```{r echo = FALSE}
county_damage_2015 %>% 
  group_by(event_type) %>%
  dplyr::mutate(tot_amount = sum(amount)) %>%
  ungroup() %>%
  filter(amount > 0) %>%
  dplyr::mutate(event_type = forcats::fct_reorder(event_type, tot_amount, .desc = TRUE)) %>%
  group_by(event_type) %>%
  dplyr::summarize(n_non_zero = n(),
                   total = sum(amount),
                   median = median(amount),
                   max = max(amount)) %>%
  arrange(desc(total)) %>% 
  mutate(n_non_zero = prettyNum(n_non_zero, big.mark = ","),
         total = paste0("$", formatC(total, format = "d", big.mark = ",")),
         median = paste0("$", formatC(median, format = "d", big.mark = ",")),
         max = paste0("$", formatC(max, format = "d", big.mark = ","))) %>%
  knitr::kable(col.names = c("Event type", "# Events", "Total damage", "Median damage", 
                             "Maximum damage"))
```

The following graph shows the distribution of events between those with listed damage costs greater than $0 and those without.

```{r echo = FALSE, fig.width = 6, fig.height = 7}
county_damage_2015 %>% 
  filter(!is.na(amount)) %>% 
  group_by(event_id) %>% 
  mutate(amount = sum(amount)) %>% 
  ungroup() %>% 
  mutate(no_damage = amount == 0,
         no_damage = factor(no_damage, 
                            levels = c(FALSE, TRUE),
                            labels = c("Damage", "No damage"))) %>% 
  group_by(event_type, no_damage) %>% 
  count() %>% 
  mutate(event_type = as_factor(event_type)) %>% 
  ggplot(aes(x = fct_reorder(event_type, n, .fun = sum), y = n, fill = no_damage)) + 
  geom_bar(stat = "identity") + 
  coord_flip() + 
  theme_classic() + 
  theme(legend.position = "top") + 
  scale_fill_viridis(name = " ", discrete = TRUE, option = "D", begin = 0.3, end = 0.8,
                     direction = -1) + 
  labs(x = "", y = "Number of events")
```

The following figure shows the top 50 episodes in terms of damages (property and crop damage combined) in 2015:

```{r fig.height = 10, fig.width = 7, fig.align = "center", echo = FALSE}
library(maps)
data(state.fips)
ex <- county_events_2015 %>% 
  filter(str_to_lower(state) %in% as.character(state.fips$polyname)) %>% 
  unite(begin_date, begin_yearmonth, begin_day, sep = "-") %>% 
  mutate(begin_date = ymd(begin_date)) %>% 
  unite(end_date, end_yearmonth, end_day, sep = "-") %>% 
  mutate(end_date = ymd(end_date)) %>% 
  mutate(damage_property = parse_damage(damage_property),
         damage_crops = parse_damage(damage_crops),
         state = str_to_title(state)) %>% 
  select(begin_date, end_date, episode_id, state, damage_property, damage_crops, event_type) %>% 
  gather(key = damage_type, value = damage, contains("damage")) %>% 
  group_by(episode_id) %>% 
  summarize(begin_date = min(begin_date),
            end_date = max(end_date),
            state = first(state),
            event_type = paste(unique(event_type), collapse = ", "),
            damage = sum(damage)) %>% 
  arrange(desc(damage)) %>% 
  slice(1:50) %>% 
  arrange(state, begin_date) %>% 
  mutate(n = 1:n(),
         damage_text = ifelse(damage == 0, "", 
                              paste0("$", formatC(damage, format = "d", big.mark = ","))))

ggplot(ex, aes(x = begin_date, y = n, color = state)) + 
  geom_segment(aes(xend = end_date, yend = n)) + 
  geom_point(size = 0.2) + 
  geom_point(aes(x = end_date), size = 0.2) +
  theme_classic() + 
  theme(legend.position = "bottom",
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.line.y = element_blank()) + 
  geom_text(aes(label = event_type, x = begin_date - 2), 
            hjust = "right", size = 2, color = "black") + 
  geom_text(aes(label = damage_text, x = end_date + 2), 
            hjust = "left", size = 2, color = "black") +  
  labs(color = "") + 
  scale_x_date(name = "Episode date", 
               breaks = ymd(c(20150101, 20150301, 20150501, 
                                   20150701, 20150901, 20151101,
                              20160101)), 
               date_labels = "%b",
               limits = c(ymd(20141001), ymd(20160201)))
```

