ecdf-density.Rmd

---
title: "R Notebook"
output:
  html_notebook: default
  pdf_document: default
---
```{r include=FALSE}
knitr::opts_chunk$set(fig.width = 7, fig.height = 4, fig.align = "center", out.width = '80%')
```

# Preamble

First we load the dataset

```{r message=FALSE}
# Source base code -----------------------------------------
source("hadisst_base.R")
source("hurdat2_pdi_base.R")
source("analysis_base.R")

# Create PDI data frame ------------------------------------

# Get hurricanes data frames
hurr.natl.obs <- get_hurr_obs("hurdat2-1851-2016-apr2017.txt")
hurr.epac.obs <- get_hurr_obs("hurdat2-nepac-1949-2016-apr2017.txt")
hurr.all.obs <- rbind(hurr.natl.obs, hurr.epac.obs)

# Create data frame with PDI and year of the storm
hurr.natl.pdi <- get_pdis(hurr.natl.obs) %>%
	filter(storm.id != "AL171988")
attr(hurr.natl.pdi, "title") <- "N. Atl."
hurr.epac.pdi <- get_pdis(hurr.epac.obs) %>%
	filter(storm.id != "EP231989")
attr(hurr.epac.pdi, "title") <- "E. Pac."
hurr.all.pdi <- get_pdis(hurr.all.obs)

# Create SST data frames -----------------------------------

hadsst.raster <- load_hadsst(file = "data/HadISST_sst.nc")

# Windows of activity
years.natl <- 1966:2016
coords.natl <- c("90W", "20E", "5N", "25N")
coords.natl.map <- c("100W", "20E", "0N", "60N")

years.epac <- 1966:2016
coords.epac <- c("120W", "90W", "5N", "20N")
coords.epac.map <- c("160W", "90W", "5N", "35N")

# Construct SST data frames
ssts.natl <- get_mean_ssts(years = years.natl, coords = coords.natl)
attr(ssts.natl, "title") <- "N. Atl."

ssts.epac <- get_mean_ssts(years = years.epac, coords = coords.epac)
attr(ssts.epac, "title") <- "E. Pac."
```

# ECDF & Density functions (PDI / energy)

```{r}
plot(density(hurr.natl.pdi$storm.pdi))
ggplot(hurr.natl.pdi, aes(storm.pdi)) +
	stat_density(geom = "line")
```

```{r}
plot(ecdf(hurr.all.pdi$storm.pdi))
ggplot(hurr.natl.pdi, aes(storm.pdi)) +
	stat_ecdf(geom = "line")
```

# ECDF & Density functions (duration)

```{r}
plot(density(hurr.natl.pdi$storm.duration))
ggplot(hurr.natl.pdi, aes(storm.duration)) +
	stat_density(geom = "line")
```

```{r}
plot(ecdf(hurr.natl.pdi$storm.duration))
ggplot(hurr.natl.pdi, aes(storm.duration)) +
	stat_ecdf(geom = "line")
```