cape-shear.Rmd

---
title: "CAPE/SHEAR vs Cluster Characteristics"
author: "James Elsner"
output: html_document
editor_options: 
  chunk_output_type: console
---

Load packages.
```{r}
library(dplyr)
library(sf)
library(ggplot2)
library(lme4)
library(lubridate)
library(xtable)
library(USAboundaries)
library(tmap)
library(ggpubr)
```

Load the cluster-level tornado data.
```{r}
load("BigDays.RData")

BigDays.sfdfT <- BigDays.sfdfT %>%
  mutate(A = as.numeric(HullArea)/10^10, #convert square kilometers
         CAPE = maxCAPE/1000,
         CIN = minCIN/100,
         DLBS = maxBS_deep/10,
         SLBS = maxBS_shallow/10,
         P = totalPOP/100000, 
         M = maxMR/10, 
         HLCY = maxHLCY/10,
         LCL =maxLCL/1000,
         CAPE_avg = avgCAPE/1000,
         CIN_avg = avgCIN/100,
         DLBS_avg = avgBS_deep/10,
         SLBS_avg = avgBS_shallow/10,
         HLCY_avg = avgHLCY/10)
dim(BigDays.sfdfT)
```

## Generate 4 figures of example clusters. 

*Create a figure of 4 clusters in the data set.* 
`Generate the state and county borders`
```{r}
sts <- state.name[!state.name %in% c("Alaska", "Hawaii")]
stateBorders <- us_states(states = sts)
counties.sf <- us_counties()
```

`Generate a color ramp that you like`
```{r}
cr <- RColorBrewer::brewer.pal(9, "Greys")
cr <- cr[-c(1:3)]
```

`Create the unique ID for the All_Tornadoes file`
```{r}
All_Tornadoes<- All_Tornadoes %>%
   mutate(ID = paste0(gsub("-", "", cDate), groupNumber))
```

Look at the distribution of CAPE, CIN, DLBS, SLBS, SRH
```{r}
( A <- ggplot(BigDays.sfdfT, (aes(maxCAPE))) +
  #geom_smooth(aes(x = week, y = probs/100), span = .5, se = FALSE, color = "gray70", size = 1) +
  geom_freqpoly((aes(maxCAPE, color = "Maximum CAPE")), lwd = 1, show.legend = TRUE) +
  geom_freqpoly((aes(avgCAPE, color = "Average CAPE")), lwd = 1, show.legend = TRUE) +
  xlab("CAPE Distribution") +
  ylab("") +
  scale_x_continuous(expand = c(0, 0), limits = c(-500, 7000),breaks = c(-500, 0, 1000, 2000, 3000, 4000, 5000, 6000, 7000),labels = c("","0", "1000", "2000", "3000", "4000", "5000", "6000", "7000")) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 125), breaks = c(0, 25, 50, 75, 100, 125),labels = c("0", "25", "50", "75", "100", "125")) +
  theme_bw() )

A <- update_labels(A, list(colour=""))

A <- A + theme(legend.position = c(0.8, 0.8),
          legend.box = "vertical")
```

```{r}
( B <- ggplot(BigDays.sfdfT, (aes(maxBS_shallow))) +
  #geom_smooth(aes(x = week, y = probs/100), span = .5, se = FALSE, color = "gray70", size = 1) +
  geom_freqpoly((aes(maxBS_shallow, color = "Maximum SLBS")), lwd = 1, show.legend = TRUE) +
  geom_freqpoly((aes(avgBS_shallow, color = "Average SLBS")), lwd = 1, show.legend = TRUE) +
  xlab("Shallow-Layer Bulk Shear Distribution") +
  ylab("") +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 5, 10, 15, 20, 25, 30, 35),labels = c("0", "5", "10", "15", "20", "25", "30", "35")) +
  scale_y_continuous(expand = c(0, 0),  limits = c(0, 100), breaks = c(0, 25, 50, 75, 100),labels = c("0", "25", "50", "75", "100")) +
  theme_bw() )

B <- update_labels(B, list(colour=""))

B <- B + theme(legend.position = c(0.8, 0.8),
          legend.box = "vertical")
```

```{r}
( C <- ggplot(BigDays.sfdfT, (aes(maxBS_deep))) +
  #geom_smooth(aes(x = week, y = probs/100), span = .5, se = FALSE, color = "gray70", size = 1) +
  geom_freqpoly((aes(maxBS_deep, color = "Maximum DLBS")), lwd = 1, show.legend = TRUE) +
  geom_freqpoly((aes(avgBS_deep, color = "Average DLBS")), lwd = 1, show.legend = TRUE) +
  xlab("Deep-Layer Bulk Shear Distribution") +
  ylab("") +
  scale_x_continuous(expand = c(0, 0), breaks = c(0, 10, 20, 30, 40, 50),labels = c("0", "10", "20", "30", "40", "50")) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 100), breaks = c(0, 25, 50, 75, 100),labels = c("0", "25", "50", "75", "100")) +
  theme_bw() )

C <- update_labels(C, list(colour=""))

C <- C + theme(legend.position = c(0.8, 0.8),
          legend.box = "vertical")
```

```{r}
( D <- ggplot(BigDays.sfdfT, (aes(maxHLCY))) +
  #geom_smooth(aes(x = week, y = probs/100), span = .5, se = FALSE, color = "gray70", size = 1) +
  geom_freqpoly((aes(maxHLCY, color = "Maximum HLCY")), lwd = 1, show.legend = TRUE) +
  geom_freqpoly((aes(avgBS_deep, color = "Average HLCY")), lwd = 1, show.legend = TRUE) +
  xlab("Storm Relative Helicity Distribution") +
  ylab("") +
  scale_x_continuous(expand = c(0, 0), limits = c(-100,900), breaks = c(0, 100, 200, 300, 400, 500, 600, 700, 800, 900),labels = c("0", "100", "200", "300", "400", "500", "600", "700", "800", "900")) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 600), breaks = c(0, 100, 200, 300, 400, 500, 600),labels = c("0", "100", "200", "300", "400", "500", "600")) +
  theme_bw() )

D <- update_labels(D, list(colour=""))

D <- D + theme(legend.position = c(0.8, 0.8),
          legend.box = "vertical")
```

```{r}
( E <- ggplot(BigDays.sfdfT, (aes(minCIN))) +
  #geom_smooth(aes(x = week, y = probs/100), span = .5, se = FALSE, color = "gray70", size = 1) +
  geom_freqpoly((aes(minCIN, color = "Minimum CIN")), lwd = 1, show.legend = TRUE) +
  geom_freqpoly((aes(avgCIN, color = "Average CIN")), lwd = 1, show.legend = TRUE) +
  xlab("Convective Inhibition Distribution") +
  ylab("") +
  scale_x_continuous(expand = c(0, 0), limits = c(-700,50), breaks = c(-700, -600, -500, -400, -300, -200, -100, 0), labels = c("-700", "-600", "-500", "-400", "-300", "-200", "-100", "0")) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 350), breaks = c(0, 50, 100, 150, 200, 250, 300, 350), labels = c("0", "50", "100", "150", "200", "250", "300", "350")) +
  theme_bw() )

E <- update_labels(E, list(colour=""))

E <- E + theme(legend.position = c(0.2, 0.8),
          legend.box = "vertical")
```

```{r}
ggarrange(A, B, C, D, E)
ggarrange(A, B, C, ncol = 1)
```

```{r}
merc = "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"

All_Tornadoes <- All_Tornadoes %>%
 mutate(ID = paste0(gsub("-", "", cDate), groupNumber))

BigDayTornadoes <- BigDayTornadoes %>%
  mutate(ID = paste0(gsub("-", "", cDate), groupNumber))
```

```{r}
#Extract the big day using the unique ID created.
cluster1 <- BigDays.sfdfT %>%
  filter(ID == 19940719217)

#Generate a convex hull around the cluster.
cluster1 <- st_convex_hull(cluster1)

#Extract all tornadoes that are in the Day 1

cluster1torns <- All_Tornadoes %>%
  filter(ID == 19940719217)

cluster1torns$mag <- as.numeric(cluster1torns$mag)
cluster1torns$mag2 <- cut(cluster1torns$mag, breaks=c(-1, 0, 1, 2, 3, 4, 5))
```

```{r}
A <- tm_shape(stateBorders) + 
  tm_borders(col = "gray70") +
tm_shape(cluster1) + 
  tm_borders(col = "gray15", 
             alpha = 1, 
             lwd = 2) +
  tm_scale_bar(color.dark = "gray70", 
               width = .3, 
               size = 1, 
               lwd = 2, 
               position = c("left","bottom")) +
    tm_compass(color.dark = "gray70", 
             size = 5, 
             lwd = 2, 
             position = c("left","top")) + 
    tm_format("World", 
              attr.position = c("left", "top"),
              legend.frame = FALSE,
              inner.margins = c(.25, .2, .2, .2)) +
    tm_layout(legend.bg.color = "white", 
            legend.text.size = .75) + 
  tm_shape(cluster1torns, 
         projection = merc, 
         is.master = TRUE) +
    tm_symbols(size = 3, 
               col = "mag2", 
               n = 6, 
               palette = cr, 
               alpha = 0.8, 
               border.alpha = 0, 
               labels = c("0", "1", "2", "3", "4", "5"), 
               title.col = "EF Rating") +
    tm_layout(title = "July 19, 1994 \n 10 tornadoes", 
              title.position = c("center", "top"), 
              legend.title.size = 1.4,
              legend.position = c("right", "bottom"), 
              legend.stack = "horizontal",
              legend.frame = FALSE, 
              legend.text.size = 1.2, 
              legend.width = -0.2, 
              title.size = 1.5)
A  
```

```{r}
cluster2 <- BigDays.sfdfT %>%
  filter(ID == 199906061644)

cluster2 <- st_convex_hull(cluster2)

cluster2torns <- BigDayTornadoes %>%
  filter(ID == 199906061644) 

cluster2 %>%
  summarize(Area = HullArea/(10**6),
            Duration = Duration/3600)

cluster2torns$mag <- as.numeric(cluster2torns$mag)
cluster2torns$mag2 <- cut(cluster2torns$mag, breaks=c(-1, 0, 1, 2, 3, 4, 5))
```

```{r}
B <- tm_shape(stateBorders) + 
  tm_borders(col = "gray70") +
tm_shape(cluster2) + 
  tm_borders(col = "gray15", 
             alpha = 1, 
             lwd = 2) +
  tm_scale_bar(color.dark = "gray70", 
               width = .2, 
               size = 1, 
               lwd = 2, 
               position = c("left","bottom")) +
    tm_compass(color.dark = "gray70", 
             size = 5, 
             lwd = 2, 
             position = c("left","top")) + 
    tm_format("World", 
              attr.position = c("left", "top"),
              legend.frame = FALSE,
              inner.margins = c(.1, .1, .1, .1)) +
    tm_layout(legend.bg.color = "white", 
            legend.text.size = .75) + 
  tm_shape(cluster2torns, 
         projection = merc, 
         is.master = TRUE) +
    tm_symbols(size = 3, 
               col = "mag2", 
               n = 6, 
               palette = cr, 
               alpha = 0.8, 
               border.alpha = 0, 
               labels = c("0", "1", "2", "3", "4", "5"), 
               title.col = "EF Rating") +
    tm_layout(title = "June 6, 1999 \n 36 tornadoes", 
              title.position = c("center", "top"), 
              legend.title.size = 1.4,
              legend.position = c("right", "bottom"), 
              legend.stack = "horizontal",
              legend.frame = FALSE, 
              legend.text.size = 1.2, 
              legend.width = -0.15, 
              title.size = 1.5)
  B
```

```{r}
cluster3 <- BigDays.sfdfT %>%
  filter(ID == 200802053876)

cluster3 <- st_convex_hull(cluster3)

cluster3torns <- BigDayTornadoes %>%
  filter(ID == 200802053876) 

cluster3 %>%
  summarize(Area = HullArea/(10**6),
            Duration = Duration/3600)

cluster3torns$mag <- as.numeric(cluster3torns$mag)
cluster3torns$mag2 <- cut(cluster3torns$mag, breaks=c(-1, 0, 1, 2, 3, 4, 5))
```

```{r}
C <- tm_shape(stateBorders) + 
  tm_borders(col = "gray70") +
tm_shape(cluster3) + 
  tm_borders(col = "gray15", 
             alpha = 1, 
             lwd = 2) +
  tm_scale_bar(color.dark = "gray70", 
               width = .3, 
               size = 1, 
               lwd = 2, 
               position = c("left","bottom")) +
    tm_compass(color.dark = "gray70", 
             size = 5, 
             lwd = 2, 
             position = c("left","top")) + 
    tm_format("World", 
              attr.position = c("left", "top"),
              legend.frame = FALSE,
              inner.margins = c(.1, .1, .2, .1)) +
    tm_layout(legend.bg.color = "white", 
            legend.text.size = .75) + 
  tm_shape(cluster3torns, 
         projection = merc, 
         is.master = TRUE) +
    tm_symbols(size = 3, 
               col = "mag2", 
               n = 6, 
               palette = cr, 
               alpha = 0.8, 
               border.alpha = 0, 
               labels = c("0", "1", "2", "3", "4", "5"), 
               title.col = "EF Rating") +
    tm_layout(title = "February 5, 2008 \n 85 tornadoes", 
              title.position = c("center", "top"), 
              legend.title.size = 1.4,
              legend.position = c("right", "bottom"), 
              legend.stack = "horizontal",
              legend.frame = FALSE, 
              legend.text.size = 1.2, 
              legend.width = -0.2, 
              title.size = 1.5)
C  
```

```{r}
#Extract the big day using the unique ID created. 
cluster4 <- BigDays.sfdfT %>%
  filter(ID == 201104274630)

#Generate a convex hull around the big day. 
cluster4 <- st_convex_hull(cluster4)

#Extract all tornadoes that are in the biggestday 
cluster4torns <- BigDayTornadoes %>%
  filter(ID == 201104274630) 

cluster4torns$mag <- as.numeric(cluster4torns$mag)
cluster4torns$mag2 <- cut(cluster4torns$mag, breaks=c(-1, 0, 1, 2, 3, 4, 5))
```

```{r}
D <- tm_shape(stateBorders) + 
  tm_borders(col = "gray70") +
tm_shape(cluster4,projection = merc) + 
  tm_borders(col = "gray15", 
             alpha = 1, 
             lwd = 2) +
  tm_scale_bar(color.dark = "gray70", 
               width = .3, 
               size = 1, 
               lwd = 2, 
               position = c("left","bottom")) +
    tm_compass(color.dark = "gray70", 
             size = 5, 
             lwd = 2, 
             position = c("left","top")) + 
    tm_format("World", 
              attr.position = c("left", "top"),
              legend.frame = FALSE,
              inner.margins = c(.15, .1, .2, .1)) +
    tm_layout(legend.bg.color = "white", 
            legend.text.size = .75) + 
  tm_shape(cluster4torns, 
         projection = merc, 
         is.master = TRUE) +
    tm_symbols(size = 3, 
               col = "mag2", 
               n = 6, 
               palette = cr, 
               alpha = 0.8, 
               border.alpha = 0, 
               labels = c("0", "1", "2", "3", "4", "5"), 
               title.col = "EF Rating") +
    tm_layout(title = "April 27, 2011 \n 173 tornadoes", 
              title.position = c("center", "top"), 
              legend.title.size = 1.4,
              legend.position = c("right", "bottom"), 
              legend.stack = "horizontal",
              legend.frame = FALSE, 
              legend.text.size = 1.2, 
              legend.width = -0.2, 
              title.size = 1.5)
  
D
```

```{r}
tmap_arrange(A, B, C, D)
```
`Figure 1: ` \label{fig:Clusters}

## Cluster statistics by time of day
```{r}
BigDays.sfdfT %>%
  group_by(NARRZtime) %>%
  summarize(clusnum = n(), 
            numtorn = sum(nT),
            tornperclus = numtorn/clusnum,
            avgdur = mean(Duration)/3600)
```

## Distribution of casualties
```{r}
BigDays.sfdfT %>%
  st_drop_geometry() %>%
  group_by(GroupDayCas) %>%
  summarize(nC = n(),
            pC = nC/nrow(BigDays.sfdfT)) %>%
  arrange(desc(nC))
```

## How are casualties changing annually?
```{r}
x <- BigDays.sfdfT %>%
  group_by(Year) %>%
  summarize(totclus = n(), 
            avgcas = (mean(GroupDayCas)/totclus),
            totcas = (sum(GroupDayCas)/totclus))

a <- ggplot(x, (aes(x = Year, y = avgcas))) +
  geom_line(color = "black", lwd = 1) +
  geom_smooth(method = lm, se = FALSE) +
  scale_x_continuous(expand = c(0, 0), breaks = c(seq(1994,2018,1))) +
  theme_bw() +
  xlab("Year") +
  ylab("Average Casualties per Cluster\n ") 

a + theme(axis.title.y = element_text(colour = "black"), axis.text.y = element_text(color = "black"),  axis.text.x = element_text(angle = 45, hjust = 1), axis.text = element_text(size = 12), axis.title=element_text(size=14))
```

Create a time series of average tornadoes per cluster by year. 
```{r}
tornyear <- BigDays.sfdfT %>%
  group_by(Year) %>%
  summarize(AvgTornperYear = mean(nT))
```

```{r}
TornYear <- ggplot(tornyear, (aes(Year, y = AvgTornperYear))) +
  geom_smooth(method = lm, se = FALSE, size = 1, col = "blue") +
  geom_line(color = "black", lwd = 1) +
  scale_x_continuous(expand = c(0, 0), breaks = c(seq(1994,2018,1)), limits = c(1994, 2018)) +
  scale_y_continuous(limits = c(10, 40), sec.axis = sec_axis(~.)) + 
#  geom_smooth(aes(x = week, y = probs/100), span = .5, se = FALSE, color = "gray70", size = 1) +
  theme_bw() +
  xlab("") +
  ylab("Average Tornadoes per Cluster\n ")

TornYear <- TornYear + theme(axis.title.y = element_text(colour = "black"), axis.text.y = element_text(color = "black"),  axis.text.x = element_text(angle = 45, hjust = 1), axis.text = element_text(size = 12), axis.title=element_text(size=14)) + ggtitle("Time Series")
TornYear
```

*Separate the clusters by week of the year.*

```{r}
BigDays.sfdfT <- BigDays.sfdfT %>%
  mutate(TorPerHour = nT/as.numeric(Duration) * 3600,
         TorPerKm = nT/as.numeric(HullArea) * 10^6)
dim(BigDays.sfdfT)
```

```{r}
TornsbyWeek <- BigDays.sfdfT %>% 
  group_by(week = week(cDate)) %>%
  summarize(totalclusters = n(),
            numtorn = sum(nT),
            avgsize = numtorn / totalclusters,
            tornrate = mean(TorPerHour), 
            torndens = mean(TorPerKm), 
            avgarea = mean(HullArea * (10^-6)), 
            avgcas = mean(GroupDayCas),
            avgcasperkm = mean((GroupDayCas/((totalPOP+1)/1000000))))
```

**For each week sum the number of clusters and divide by the number of years to get the rate lambda. Then take 1-dpois(lambda) * 100 to get the probability**

## Empirical Probability of getting a cluster by week`
```{r}
Years <- 2018-1994 + 1

lambdas <- TornsbyWeek$totalclusters / Years
probs <- (1 - dpois(0, lambdas)) * 100

TornsbyWeek <- cbind(TornsbyWeek, probs)
```

```{r}
x.Date <- as.Date(paste(rep(1994:2018, each = 12), rep(1:12, 2), 1, sep = "-"))
library(zoo)
x <- zoo(rnorm(24), x.Date)
times <- time(x)
ticks <- as.data.frame(x = seq(times[1], times[length(times)], by = "weeks"))

week <- as.data.frame(ticks[1:53,])
  
months <- as.data.frame(format(week, "%b"))
Mo <- as.data.frame(format(week, "%m"))
day <- as.data.frame(format(week, "%d"))

dat <- as.data.frame(cbind(week, Mo, months, day))
colnames(dat) <- c("Week", "Mo", "Month", "Day") 

dat <- dat %>%
  group_by(Month, Day, Mo) %>%
  summarize(count = n(),
            Week = paste0(Month, " ", Day))


dat <- dat[order(dat$Mo),]

labels = dat$Week
```

```{r}
A <- ggplot(TornsbyWeek, (aes(week))) +
  #geom_smooth(aes(x = week, y = probs/100), span = .5, se = FALSE, color = "gray70", size = 1) +
  geom_line(aes(y = probs/100), color = "black", lwd = 1) +
  scale_x_continuous(expand = c(0, 0), breaks = c(seq(1,53,3)), limits = c(1, 53), labels = labels[seq(1, length(labels), 3)]) +
  scale_y_continuous(limits = c(0, 1)) + 
#  geom_smooth(aes(x = week, y = probs/100), span = .5, se = FALSE, color = "gray70", size = 1) +
  theme_bw() +
  xlab("") +
  ylab("Probability of a cluster\n ")

A <- A + theme(panel.grid = element_blank(), axis.title.y = element_text(colour = "black"), axis.text.y = element_text(color = "black"),  axis.text.x = element_text(angle = 45, hjust = 1), axis.text = element_text(size = 12), axis.title=element_text(size=14)) + ggtitle("A")
```

```{r}
B <- ggplot(TornsbyWeek, (aes(week))) +
  #geom_smooth(aes(x = week, y = avgsize), span = .5, se = FALSE, color = "gray70", size = 1) +
  geom_line(aes(y = avgsize), color = "black", lwd = 1) +
  scale_x_continuous(expand = c(0, 0), breaks = c(seq(1,53,3)), limits = c(1, 53), labels = labels[seq(1, length(labels), 3)]) +
  scale_y_continuous(breaks = c(seq(0,50,10)), limits = c(0, 50)) + 
  theme_bw() +
  xlab("") +
  ylab("Number of tornadoes\n ")

B <- B + theme(panel.grid = element_blank(), axis.title.y = element_text(colour = "black"), axis.text.y = element_text(color = "black"),  axis.text.x = element_text(angle = 45, hjust = 1), axis.text = element_text(size = 12), axis.title=element_text(size=14)) + ggtitle("B")
```

```{r}
C <- ggplot(TornsbyWeek, (aes(week))) +
  #geom_smooth(aes(x = week, y = as.numeric(avgcasperkm)), span = .5, se = FALSE, color = "gray70", size = 1) +
  geom_line(aes(y = as.numeric(avgcasperkm)), color = "black", lwd = 1) +
  scale_x_continuous(expand = c(0, 0), breaks = c(seq(1,53,3)), limits = c(1, 53), labels = labels[seq(1, length(labels), 3)]) +
  #scale_y_continuous(breaks = c(seq(0,150,25)), limits = c(0, 150)) + 
  theme_bw() +
  xlab("") +
  ylab("Number of casualties\n")

C <- C + theme(panel.grid = element_blank(), axis.title.y = element_text(colour = "black"), axis.text.y = element_text(color = "black"),  axis.text.x = element_text(angle = 45, hjust = 1), axis.text = element_text(size = 12), axis.title=element_text(size=14)) + ggtitle("C")
C
```

## Combine into one figure
```{r}
ggarrange(A, B, C, ncol = 1)
```

## Get the maximum, minimum, and average of the environmental variables. Create a table displaying this information.
```{r}
maximum <- as.matrix(c(max(BigDays.sfdfT$maxCAPE), max(BigDays.sfdfT$maxVSTM), max(BigDays.sfdfT$maxBS_deep), max(BigDays.sfdfT$maxBS_shallow), max(BigDays.sfdfT$maxHLCY), (max(BigDays.sfdfT$HullArea)*1e-6), max(BigDays.sfdfT$minCIN), max(BigDays.sfdfT$maxUSTM)), ncol = 1)

minimum <- as.matrix(c(min(BigDays.sfdfT$maxCAPE), min(BigDays.sfdfT$maxVSTM), min(BigDays.sfdfT$maxBS_deep), min(BigDays.sfdfT$maxBS_shallow), min(BigDays.sfdfT$maxHLCY), (min(BigDays.sfdfT$HullArea)*1e-6), min(BigDays.sfdfT$minCIN), min(BigDays.sfdfT$maxUSTM)), ncol = 1)

average <- as.matrix(c(mean(BigDays.sfdfT$maxCAPE), mean(BigDays.sfdfT$maxVSTM), mean(BigDays.sfdfT$maxBS_deep), mean(BigDays.sfdfT$maxBS_shallow), mean(BigDays.sfdfT$maxHLCY), (mean(BigDays.sfdfT$HullArea)*1e-6), mean(BigDays.sfdfT$minCIN), mean(BigDays.sfdfT$maxUSTM)), ncol = 1)

Variable <- as.matrix(c("Convective Available Potential Energy", "Northward Storm Motion", "Deep-Layer Bulk Shear", "Shallow-Layer Bulk Shear", "Helicity", "Outbreak Area", "CIN", "Eastward Storm Motion"), ncol = 1)

abbr <- as.matrix(c("CAPE", "VSTM", "DLBS","SLBS", "HLCY", "AREA", "CIN", "USTM"), ncol = 1)

env_variation <- cbind(Variable, abbr, maximum, minimum, average)
colnames(env_variation) <- c("Variable Name", "Abbreviation", "Maximum", "Minimum", "Average")

xtable(as.data.frame((env_variation))) #round to sig fig
```
`Table 2: The range and average highest (lowest for CIN) value of the environmental variables across the 767 tornado clusters used in the study.` \label{tab:VarValues}

## Create a correlation matrix of the environmental variables

```{r}
env.dat <- cbind(BigDays.sfdfT$maxCAPE, BigDays.sfdfT$minCIN, BigDays.sfdfT$maxBS_deep, BigDays.sfdfT$maxBS_shallow, BigDays.sfdfT$maxHLCY,BigDays.sfdfT$maxDEW, BigDays.sfdfT$maxMR, BigDays.sfdfT$A, BigDays.sfdfT$Lat, BigDays.sfdfT$Lon, BigDays.sfdfT$P)

dat = na.omit(env.dat)

xtable(cor(dat))

```


## Tornado count model

Consider only CAPE, CIN, BS (deep & shallow) for the scaled maximum values of the variables.
```{r}
library(MASS)


modelInitial <- glm.nb(nT ~ A + Lat + Lon + Year + CAPE + CIN + DLBS + SLBS + HLCY + LCL, data = BigDays.sfdfT)

modelInitial <- glm.nb(nT ~ A + Lat + Lon + Year + CAPE + CIN + DLBS + SLBS + HLCY, data = BigDays.sfdfT)

summary(modelInitial)

modelInitialPoisson <- glm(nT ~ A + Lat + Lon + CAPE + CIN + DLBS + SLBS + HLCY, family = "poisson", data = BigDays.sfdfT)
summary(modelInitialPoisson)

modelFinal <- glm.nb(nT ~ A + CAPE  + DLBS + SLBS , data = BigDays.sfdfT)
summary(modelFinal)

(exp(coef(modelFinal)) - 1) * 100

hist(resid(modelFinal))
range(exp(predict(modelFinal)))
plot(log(BigDays.sfdfT$nT), predict(modelFinal))
cor.test(BigDays.sfdfT$nT, predict(modelFinal, type = "response"))
mean(abs(BigDays.sfdfT$nT - predict(modelFinal, type = "response")))
mean(abs(BigDays.sfdfT$nT - predict(modelFinal, type = "response")))/(max(BigDays.sfdfT$nT) - min(BigDays.sfdfT$nT))
mean(abs(BigDays.sfdfT$nT - predict(modelFinal, type = "response")))/(max(predict(modelFinal, type = "response")) - min(predict(modelFinal, type = "response")))

modelglmer.nb <- glmer.nb(nT ~ scale(A) + scale(CAPE) + scale(DLBS) + scale(SLBS) + (1|Month), data = BigDays.sfdfT)
summary(modelglmer.nb) # mixed effects model does not improve things
```

## Observed versus predicted plot
```{r}
df <- data.frame(Observed = BigDays.sfdfT$nT,
                 Predicted = predict(modelFinal, type = "response"))

ggplot(df, aes(x = Observed, y = Predicted)) +
  geom_abline() +
  geom_point(col = "gray70") +
  geom_smooth(method = lm, size = 1, col = "black") +
  scale_x_log10() +
  scale_y_log10() +
  theme_minimal() 
```

Break them into terciles for a random sample of 99 cases. 
```{r}
obs <- df[sample(nrow(df), 99), ]
pred <- obs$Predicted
obs <- obs$Observed
table(ntile(obs, 3), ntile(pred, 3))
```


Consider only CAPE, CIN, BS (deep & shallow) for the scaled average values of the variables.
```{r, eval = FALSE}
library(MASS)

modelInitial <- glm.nb(nT ~ A + Lat + Lon + Year + CAPE_avg + CIN_avg + DLBS_avg + SLBS_avg + HLCY_avg, data = BigDays.sfdfT)
summary(modelInitial)

modelInitialPoisson <- glm(nT ~ A + Lat + Lon + CAPE_avg + CIN_avg + DLBS_avg + SLBS_avg, family = "poisson", data = BigDays.sfdfT)
summary(modelInitialPoisson)

modelFinal <- glm.nb(nT ~ A + CAPE_avg + DLBS_avg + SLBS_avg, data = BigDays.sfdfT)
summary(modelFinal)

(exp(coef(modelFinal)) - 1) * 100

hist(resid(modelFinal))
range(exp(predict(modelFinal)))
plot(log(BigDays.sfdfT$nT), predict(modelFinal))
cor.test(BigDays.sfdfT$nT, predict(modelFinal, type = "response"))
mean(abs(BigDays.sfdfT$nT - predict(modelFinal, type = "response")))
mean(abs(BigDays.sfdfT$nT - predict(modelFinal, type = "response")))/(max(BigDays.sfdfT$nT) - min(BigDays.sfdfT$nT))
mean(abs(BigDays.sfdfT$nT - predict(modelFinal, type = "response")))/(max(predict(modelFinal, type = "response")) - min(predict(modelFinal, type = "response")))

modelglmer.nb <- glmer.nb(nT ~ scale(A) + scale(CAPE_avg) + scale(DLBS_avg) + scale(SLBS_avg) + (1|Month), data = BigDays.sfdfT)
summary(modelglmer.nb) # mixed effects model does not improve things
```

## Observed versus predicted plot
```{r}
df <- data.frame(Observed = BigDays.sfdfT$nT,
                 Predicted = predict(modelFinal, type = "response"))

ggplot(df, aes(x = Observed, y = Predicted)) +
  geom_abline() +
  geom_point(col = "gray70") +
  geom_smooth(method = lm, size = 1, col = "black") +
  scale_x_log10() +
  scale_y_log10() +
  theme_minimal() 
```

## Casualty model
```{r}
modelInitialC <- glm.nb(GroupDayCas ~ P + Lat + Lon + Year + CAPE + CIN + DLBS + SLBS, data = BigDays.sfdfT)
summary(modelInitialC)

modelFinalC <- glm.nb(GroupDayCas ~ P + Lat + Lon + Year + CAPE + DLBS + SLBS, data = BigDays.sfdfT)
summary(modelFinalC)

hist(resid(modelFinalC))
range(exp(predict(modelFinalC)))
plot(log(BigDays.sfdfT$GroupDayCas), predict(modelFinalC))
cor.test(BigDays.sfdfT$GroupDayCas, exp(predict(modelFinalC)))

mean(abs(BigDays.sfdfT$GroupDayCas - predict(modelFinalC, type = "response")))
mean(abs(BigDays.sfdfT$GroupDayCas - predict(modelFinalC, type = "response")))/(max(BigDays.sfdfT$GroupDayCas) - min(BigDays.sfdfT$GroupDayCas))
mean(abs(BigDays.sfdfT$GroupDayCas - predict(modelFinalC, type = "response")))/(max(predict(modelFinalC, type = "response")) - min(predict(modelFinalC, type = "response")))
```

## Observed versus predicted plot
```{r}
df <- data.frame(Observed = BigDays.sfdfT$GroupDayCas,
                 Predicted = predict(modelFinalC, type = "response"))

ggplot(df, aes(x = Observed +.54, y = Predicted)) +
  geom_abline() +
  geom_point(col = "gray70") +
  geom_smooth(method = lm, size = 1, col = "black") +
  xlab("Observed") +
  scale_x_log10() +
  scale_y_log10() +
  theme_minimal() 
```

## Cross validation
```{r}
formulaT <- nT ~ A + CAPE + DLBS + SLBS
formulaC <- GroupDayCas ~ P + Lat + Lon + Year + CAPE + DLBS + SLBS #+ nT
modelFinalT <- glm.nb(formulaT, data = BigDays.sfdfT)
modelFinalC <- glm.nb(formulaC, data = BigDays.sfdfT)

TCount <- numeric()
CCount <- numeric()

for(i in 1:nrow(BigDays.sfdfT)){
TCount[i] <- predict(glm.nb(formulaT, data = BigDays.sfdfT[-i, ]), 
                       newdata = BigDays.sfdfT[i,],
                       type = "response")
CCount[i] <- predict(glm.nb(formulaC, data = BigDays.sfdfT[-i, ]), 
                       newdata = BigDays.sfdfT[i,],
                       type = "response")
}
cor(BigDays.sfdfT$nT, TCount)
cor(BigDays.sfdfT$GroupDayCas, CCount)

mean(abs(BigDays.sfdfT$nT - TCount))
mean(abs(BigDays.sfdfT$nT - TCount))/(max(BigDays.sfdfT$nT) - min(BigDays.sfdfT$nT))
mean(abs(BigDays.sfdfT$nT - predict(modelFinalT, type = "response")))/(max(predict(modelFinalT, type = "response")) - min(predict(modelFinalT, type = "response")))

mean(abs(BigDays.sfdfT$GroupDayCas - CCount))
mean(abs(BigDays.sfdfT$GroupDayCas - CCount))/(max(BigDays.sfdfT$GroupDayCas) - min(BigDays.sfdfT$GroupDayCas))
mean(abs(BigDays.sfdfT$GroupDayCas - predict(modelFinalC, type = "response")))/(max(predict(modelFinalC, type = "response")) - min(predict(modelFinalC, type = "response")))
```

## Under/over predictions
```{r}
df <- data.frame(Observed = BigDays.sfdfT$GroupDayCas,
                 Predicted = CCount) %>%
  mutate(Diff = Observed - Predicted)
```

## Tornado predictions 1-D
```{r}
theta <- 6.246
nTthresh <- 20
predict(modelFinalT, newdata = data.frame(A = c(mean(BigDays.sfdfT$A), mean(BigDays.sfdfT$A)),
                                         CAPE = c(0, 5),
                                         DLBS = c(mean(BigDays.sfdfT$DLBS), mean(BigDays.sfdfT$DLBS)),
                                         SLBS = c(mean(BigDays.sfdfT$SLBS), mean(BigDays.sfdfT$SLBS))),
        type = "response")

1 - pnbinom(nTthresh, size = theta, mu = 18.725) # chance that the outbreak will have more than 50 tornadoes when CAPE = 0
1 - pnbinom(nTthresh, size = theta, mu = 23.55)  # chance that the outbreak will have more than 50 tornadoes when CAPE = 5000

predictions  <- predict(modelFinalT, 
                        newdata = data.frame(A = rep(mean(BigDays.sfdfT$A), times = 6),
                                             CAPE = 0:5,
                                             DLBS = rep(mean(BigDays.sfdfT$DLBS), times = 6),
                                             SLBS = rep(mean(BigDays.sfdfT$SLBS), times = 6)),
                        type = "response")
```

```{r}
( CAPElevels <- round(quantile(BigDays.sfdfT$maxCAPE, probs = c(.05, .5, .95))) )
CAPElevels <- c(100, 2000, 5000)
theta <- 6.246
nTthresh <- 30
Area <- 40 # Area of 10\% tor probability April 12, 2020 outbreak 12UTC

predictions  <- predict(modelFinalT, 
                        newdata = data.frame(A = Area,
                                             CAPE = CAPElevels[1]/1000,
                                             DLBS = seq(1, 4, by = .25),
                                             SLBS = mean(BigDays.sfdfT$SLBS)),
                        type = "response",
                        se.fit = TRUE)
plot.predictions <- data.frame(CAPE = CAPElevels[1],
                               probability = 1 - pnbinom(nTthresh, size = theta, mu = predictions$fit),
                               probabilityU = 1 - pnbinom(nTthresh, size = theta, mu = predictions$fit + 1 * predictions$se.fit),
                               probabilityL = 1 - pnbinom(nTthresh, size = theta, mu = predictions$fit - 1 * predictions$se.fit),
                               DLBS = seq(1, 4, by = .25) * 10)

predictions  <- predict(modelFinalT, 
                        newdata = data.frame(A = Area,
                                             CAPE = CAPElevels[2]/1000,
                                             DLBS = seq(1, 4, by = .25),
                                             SLBS = mean(BigDays.sfdfT$SLBS)),
                        type = "response",
                        se.fit = TRUE)

plot.predictions2 <- data.frame(CAPE = CAPElevels[2],
                               probability = 1 - pnbinom(nTthresh, size = theta, mu = predictions$fit),
                               probabilityU = 1 - pnbinom(nTthresh, size = theta, mu = predictions$fit + 1 * predictions$se.fit),
                               probabilityL = 1 - pnbinom(nTthresh, size = theta, mu = predictions$fit - 1 * predictions$se.fit),
                               DLBS = seq(1, 4, by = .25) * 10)

predictions  <- predict(modelFinalT, 
                        newdata = data.frame(A = Area,
                                             CAPE = CAPElevels[3]/1000,
                                             DLBS = seq(1, 4, by = .25),
                                             SLBS = mean(BigDays.sfdfT$SLBS)),
                        type = "response",
                        se.fit = TRUE)

plot.predictions3 <- data.frame(CAPE = CAPElevels[3],
                               probability = 1 - pnbinom(nTthresh, size = theta, mu = predictions$fit),
                               probabilityU = 1 - pnbinom(nTthresh, size = theta, mu = predictions$fit + 1 * predictions$se.fit),
                               probabilityL = 1 - pnbinom(nTthresh, size = theta, mu = predictions$fit - 1 * predictions$se.fit),
                               DLBS = seq(1, 4, by = .25) * 10)


plot.predictions <- rbind(plot.predictions, plot.predictions2, plot.predictions3)


library(ggplot2)
ggplot(plot.predictions, aes(x = DLBS,
                             y = probability, 
                             color = factor(CAPE, levels = rev(CAPElevels)),
                             group = factor(CAPE, levels = rev(CAPElevels)))) +
  geom_line(size = 2) +
  geom_ribbon(aes(ymin = probabilityL, ymax = probabilityU), color = NA, fill = "gray80", alpha = .2) +
  scale_y_continuous(limits = c(0, 1), labels = scales::percent) +
  scale_colour_ordinal(name = "CAPE (J/kg)", direction = -1, alpha = .6) +
  ylab(paste("Probability of at least", nTthresh, "tornadoes")) +
  xlab("Deep-layer bulk shear (m/s)") +
  theme_minimal()
```

## Tornado predictions 2-D
```{r}
library(directlabels)

pgrid <- expand.grid(CAPE = seq(1, 5, .05),
                     DLBS = seq(1, 4, .05),
                     A = Area,
                     SLBS = mean(BigDays.sfdfT$SLBS))

predictions <- predict.glm(modelFinalT, newdata = pgrid, type = "response")
pgrid$probability <-1 - pnbinom(nTthresh, size = theta, mu = predictions)

( p <- ggplot(pgrid, aes(x = CAPE * 1000, y = DLBS * 10, fill = probability)) +
          geom_raster() +
          scale_fill_viridis_c(alpha = .6) +
          geom_contour(aes(color = ..level.., z = probability), breaks = seq(.3, .6, by = .1), color = "black") +
          ylab("Deep-layer bulk shear (m/s)") +
          xlab("CAPE (J/kg)") + 
          theme_minimal() )
#    ggtitle("Probability of at least 30 tornadoes in an outbreak")
direct.label(p, list("angled.boxes"))
```

##Model Residual 
Rank the residuals: 
```{r}
#Model Underpredicts
head(sort(Residuals_My))
#row 556 and 381
(resid(modelFinalT)[556])
resid(modelFinalT)[381]

#Model Overpredicts
tail(sort(resid(modelFinalT)))
#row 149 and 279
resid(modelFinalT)[149]
resid(modelFinalT)[279]
```


```{r}
df <- data.frame(ID = BigDays.sfdfT$ID,
                 Observed = BigDays.sfdfT$nT,
                 Predicted = predict(modelFinalT, type = "response"), 
                 Area_sqkm = BigDays.sfdfT$A,
                 CAPE = BigDays.sfdfT$CAPE,
                 DLBS = BigDays.sfdfT$DLBS,
                 SLBS = BigDays.sfdfT$SLBS, 
                 nT = BigDays.sfdfT$nT)

df <- df%>%
  mutate(residuals = Observed - Predicted,
         clusdens = Observed/Area_sqkm) 

Residuals_My <- df$Observed - df$Predicted
mean(Residuals_My)
#-0.06287

Model_lessdens <- df %>%
  arrange(clusdens) %>%
  slice(6:(n()))

formulaT <- nT ~ Area_sqkm + CAPE + DLBS + SLBS
modelFinalT2 <- glm.nb(formulaT, data = Model_lessdens)
summary(modelFinalT2)

df <- Model_lessdens %>%
  mutate(Predicted = predict(modelFinalT2, type = "response"))

#Mean Absolute Error: 
mean(abs(df$nT - df$Predicted))

library(Metrics)
mae(df$nT, df$Predicted)

```


Create a graph of the 2 most overpredicted clusters by the model. 
```{r}
#Extract the big day using the unique ID created. 
clusover1 <- BigDays.sfdfT %>%
  filter(ID == 201105234655)

#Generate a convex hull around the big day. 
clusover1 <- st_convex_hull(clusover1)

clusover1$nT

#Extract all tornadoes that are in the biggestday 
clusover1torns <- BigDayTornadoes %>%
  filter(ID == 201105234655) 

predover1.df <- df %>% filter(ID == 201105234655)
predover1.df

clusover1torns$mag <- as.numeric(clusover1torns$mag)
clusover1torns$mag2 <- cut(clusover1torns$mag, breaks=c(-1, 0, 1, 2, 3, 4, 5))
```

```{r}
Over1 <- tm_shape(stateBorders) + 
  tm_borders(col = "gray70") +
tm_shape(clusover1,projection = merc) + 
  tm_borders(col = "gray15", 
             alpha = 1, 
             lwd = 2) +
  tm_scale_bar(color.dark = "gray70", 
               width = .3, 
               size = 1, 
               lwd = 2, 
               position = c("left","bottom")) +
    tm_compass(color.dark = "gray70", 
             size = 5, 
             lwd = 2, 
             position = c("left","top")) + 
    tm_format("World", 
              attr.position = c("left", "top"),
              legend.frame = FALSE,
              inner.margins = c(.15, .1, .2, .1)) +
    tm_layout(legend.bg.color = "white", 
            legend.text.size = .75) + 
  tm_shape(clusover1torns, 
         projection = merc, 
         is.master = TRUE) +
    tm_symbols(size = 3, 
               col = "mag2", 
               n = 6, 
               palette = cr, 
               alpha = 0.8, 
               border.alpha = 0, 
               labels = c("0", "1", "2", "3", "4", "5"), 
               title.col = "EF Rating") +
    tm_layout(title = "May 23, 2011  (Overpredict) \n 17 tornadoes", 
              title.position = c("center", "top"), 
              legend.title.size = 1.4,
              legend.position = c("right", "bottom"), 
              legend.stack = "horizontal",
              legend.frame = FALSE, 
              legend.text.size = 1.2, 
              legend.width = -0.2, 
              title.size = 1.5)
  
Over1
```

```{r}
#Extract the big day using the unique ID created. 
clusover2 <- BigDays.sfdfT %>%
  filter(ID == 200805253953)

#Generate a convex hull around the big day. 
clusover2 <- st_convex_hull(clusover2)

clusover2$nT
clusover2$ID
#Extract all tornadoes that are in the biggestday 
clusover2torns <- BigDayTornadoes %>%
  filter(ID == 200805253953) 

predover2.df <- df %>% filter(ID == 200805253953)
predover2.df

clusover2torns$mag <- as.numeric(clusover2torns$mag)
clusover2torns$mag2 <- cut(clusover2torns$mag, breaks=c(-1, 0, 1, 2, 3, 4, 5))
```

```{r}
Over2 <- tm_shape(stateBorders) + 
  tm_borders(col = "gray70") +
tm_shape(clusover2,projection = merc) + 
  tm_borders(col = "gray15", 
             alpha = 1, 
             lwd = 2) +
  tm_scale_bar(color.dark = "gray70", 
               width = .3, 
               size = 1, 
               lwd = 2, 
               position = c("left","bottom")) +
    tm_compass(color.dark = "gray70", 
             size = 5, 
             lwd = 2, 
             position = c("left","top")) + 
    tm_format("World", 
              attr.position = c("left", "top"),
              legend.frame = FALSE,
              inner.margins = c(.15, .1, .2, .1)) +
    tm_layout(legend.bg.color = "white", 
            legend.text.size = .75) + 
  tm_shape(clusover2torns, 
         projection = merc, 
         is.master = TRUE) +
    tm_symbols(size = 3, 
               col = "mag2", 
               n = 6, 
               palette = cr, 
               alpha = 0.8, 
               border.alpha = 0, 
               labels = c("0", "1", "2", "3", "4", "5"), 
               title.col = "EF Rating") +
    tm_layout(title = "May 25, 2008  (Overpredict) \n 32 tornadoes", 
              title.position = c("center", "top"), 
              legend.title.size = 1.4,
              legend.position = c("right", "bottom"), 
              legend.stack = "horizontal",
              legend.frame = FALSE, 
              legend.text.size = 1.2, 
              legend.width = -0.2, 
              title.size = 1.5)
  
Over2
```

Create a graph of the 2 most underpredicted clusters by the model. 
```{r}
#Extract the big day using the unique ID created. 
clusunder1 <- BigDays.sfdfT %>%
  filter(ID == 199901211532)

#Generate a convex hull around the big day. 
clusunder1 <- st_convex_hull(clusunder1)

clusunder1$nT
clusunder1$ID

#Extract all tornadoes that are in the biggestday 
clusunder1torns <- BigDayTornadoes %>%
  filter(ID == 199901211532) 

predunder1.df <- df %>% filter(ID == 199901211532)
predunder1.df

clusunder1torns$mag <- as.numeric(clusunder1torns$mag)
clusunder1torns$mag2 <- cut(clusunder1torns$mag, breaks=c(-1, 0, 1, 2, 3, 4, 5))
```

```{r}
Under1 <- tm_shape(stateBorders) + 
  tm_borders(col = "gray70") +
tm_shape(clusunder1,projection = merc) + 
  tm_borders(col = "gray15", 
             alpha = 1, 
             lwd = 2) +
  tm_scale_bar(color.dark = "gray70", 
               width = .3, 
               size = 1, 
               lwd = 2, 
               position = c("left","bottom")) +
    tm_compass(color.dark = "gray70", 
             size = 5, 
             lwd = 2, 
             position = c("left","top")) + 
    tm_format("World", 
              attr.position = c("left", "top"),
              legend.frame = FALSE,
              inner.margins = c(.15, .1, .2, .1)) +
    tm_layout(legend.bg.color = "white", 
            legend.text.size = .75) + 
  tm_shape(clusunder1torns, 
         projection = merc, 
         is.master = TRUE) +
    tm_symbols(size = 3, 
               col = "mag2", 
               n = 6, 
               palette = cr, 
               alpha = 0.8, 
               border.alpha = 0, 
               labels = c("0", "1", "2", "3", "4", "5"), 
               title.col = "EF Rating") +
    tm_layout(title = "January 21, 1999  (Underpredict) \n 99 tornadoes", 
              title.position = c("center", "top"), 
              legend.title.size = 1.4,
              legend.position = c("right", "bottom"), 
              legend.stack = "horizontal",
              legend.frame = FALSE, 
              legend.text.size = 1.2, 
              legend.width = -0.2, 
              title.size = 1.5)
  
Under1
```

```{r}
#Extract the big day using the unique ID created. 
clusunder2 <- BigDays.sfdfT %>%
  filter(ID == 201104274630)

#Generate a convex hull around the big day. 
clusunder2 <- st_convex_hull(clusunder2)

clusunder2$nT
clusunder2$ID
#Extract all tornadoes that are in the biggestday 
clusunder2torns <- BigDayTornadoes %>%
  filter(ID == 201104274630) 

predunder2.df <- df %>% filter(ID == 201104274630)
predunder2.df

clusunder2torns$mag <- as.numeric(clusunder2torns$mag)
clusunder2torns$mag2 <- cut(clusunder2torns$mag, breaks=c(-1, 0, 1, 2, 3, 4, 5))
```

```{r}
Under2 <- tm_shape(stateBorders) + 
  tm_borders(col = "gray70") +
tm_shape(clusunder2,projection = merc) + 
  tm_borders(col = "gray15", 
             alpha = 1, 
             lwd = 2) +
  tm_scale_bar(color.dark = "gray70", 
               width = .3, 
               size = 1, 
               lwd = 2, 
               position = c("left","bottom")) +
    tm_compass(color.dark = "gray70", 
             size = 5, 
             lwd = 2, 
             position = c("left","top")) + 
    tm_format("World", 
              attr.position = c("left", "top"),
              legend.frame = FALSE,
              inner.margins = c(.15, .1, .2, .1)) +
    tm_layout(legend.bg.color = "white", 
            legend.text.size = .75) + 
  tm_shape(clusunder2torns, 
         projection = merc, 
         is.master = TRUE) +
    tm_symbols(size = 3, 
               col = "mag2", 
               n = 6, 
               palette = cr, 
               alpha = 0.8, 
               border.alpha = 0, 
               labels = c("0", "1", "2", "3", "4", "5"), 
               title.col = "EF Rating") +
    tm_layout(title = "April 27, 2011 (Underpredict) \n 173 tornadoes", 
              title.position = c("center", "top"), 
              legend.title.size = 1.4,
              legend.position = c("right", "bottom"), 
              legend.stack = "horizontal",
              legend.frame = FALSE, 
              legend.text.size = 1.2, 
              legend.width = -0.2, 
              title.size = 1.5)
  
Under2
```

Create a 4-panel map of the two  most over predicted and two most under predicted variables. 
```{r}
tmap_arrange(Over1, Over2, Under1, Under2)
```

Compare the correlation between Residuals and the Cluster Density: 
```{r}
cor(df$residuals, df$clusdens)
#0.02132019
```

Compare the residuals for April 27, 2011 and June 6, 1999 to the average residuals from the model. 
```{r}
mean(resid(modelFinalT))

#The average residual is -0.1317874

#April 27, 2011: Row 550

predApr27.df <- df %>% filter(ID == 201104274630)
predApr27.df

predApr27.df$Observed - predApr27.df$Predicted

#June 6, 1999: Row 168
predJun6.df <- df %>% filter(ID == 199906061644)
predJun6.df
predJun6.df$Observed - predJun6.df$Predicted
```

## Casualty predictions 1-D

Start with shear (DLBS)
```{r}
( Plevels <- c(10, 40, 80) ) # 100,000s
theta <- .261
nCthresh <- 50 # arbitrary threshold number of casualties

predictions  <- predict(modelFinalC, 
                        newdata = data.frame(P = Plevels[1],
                                             Lat = mean(BigDays.sfdfT$Lat),
                                             Lon = mean(BigDays.sfdfT$Lon),
                                             Year = 2018,
                                             CAPE = mean(BigDays.sfdfT$CAPE),
                                             DLBS = seq(1, 4, by = .25),
                                             SLBS = mean(BigDays.sfdfT$SLBS)),
                        type = "response",
                        se.fit = TRUE)
plot.predictions <- data.frame(P = Plevels[1],
                               probability = 1 - pnbinom(nCthresh, size = theta, mu = predictions$fit),
                               probabilityU = 1 - pnbinom(nCthresh, size = theta, mu = predictions$fit + 1 * predictions$se.fit),
                               probabilityL = 1 - pnbinom(nCthresh, size = theta, mu = predictions$fit - 1 * predictions$se.fit),
                               DLBS = seq(1, 4, by = .25) * 10)

predictions  <- predict(modelFinalC, 
                        newdata = data.frame(P = Plevels[2],
                                             Lat = mean(BigDays.sfdfT$Lat),
                                             Lon = mean(BigDays.sfdfT$Lon),
                                             Year = 2018,
                                             CAPE = mean(BigDays.sfdfT$CAPE),
                                             DLBS = seq(1, 4, by = .25),
                                             SLBS = mean(BigDays.sfdfT$SLBS)),
                        type = "response",
                        se.fit = TRUE)

plot.predictions2 <- data.frame(P = Plevels[2],
                               probability = 1 - pnbinom(nCthresh, size = theta, mu = predictions$fit),
                               probabilityU = 1 - pnbinom(nCthresh, size = theta, mu = predictions$fit + 1 * predictions$se.fit),
                               probabilityL = 1 - pnbinom(nCthresh, size = theta, mu = predictions$fit - 1 * predictions$se.fit),
                               DLBS = seq(1, 4, by = .25) * 10)

predictions  <- predict(modelFinalC, 
                        newdata = data.frame(P = Plevels[3],
                                             Lat = mean(BigDays.sfdfT$Lat),
                                             Lon = mean(BigDays.sfdfT$Lon),
                                             Year = 2018,
                                             CAPE = mean(BigDays.sfdfT$CAPE),
                                             DLBS = seq(1, 4, by = .25),
                                             SLBS = mean(BigDays.sfdfT$SLBS)),
                        type = "response",
                        se.fit = TRUE)

plot.predictions3 <- data.frame(P = Plevels[3],
                               probability = 1 - pnbinom(nCthresh, size = theta, mu = predictions$fit),
                               probabilityU = 1 - pnbinom(nCthresh, size = theta, mu = predictions$fit + 1 * predictions$se.fit),
                               probabilityL = 1 - pnbinom(nCthresh, size = theta, mu = predictions$fit - 1 * predictions$se.fit),
                               DLBS = seq(1, 4, by = .25) * 10)


plot.predictions <- rbind(plot.predictions, plot.predictions2, plot.predictions3)

cbPalette <- c("#E69F00", "#56B4E9", "#009E73")
(p1 <- ggplot(plot.predictions, aes(x = DLBS, 
                                    y = probability, 
                                    color = factor(P * 100000, levels = rev(Plevels * 100000)),
                                    group = factor(P * 100000, levels = rev(Plevels * 100000)))) +
  geom_line(size = 2) +
  geom_ribbon(aes(ymin = probabilityL, ymax = probabilityU), color = NA, fill = "gray80", alpha = .2) +
  scale_y_continuous(limits = c(0, .3), labels = scales::percent) +
  scale_colour_manual(name = "Population", labels = c("8M", "4M", "1M"), values = cbPalette) + 
  ylab(paste("Probability of at least", nCthresh, "casualties")) +
  xlab("Deep-layer bulk shear (m/s)") +
  theme_minimal() )
```

Repeat for with CAPE
```{r}
( Plevels <- c(10, 40, 80) ) # 100,000s
theta <- .261
nCthresh <- 50 # arbitrary threshold number of casualties
predictions  <- predict(modelFinalC, 
                        newdata = data.frame(P = Plevels[1],
                                             Lat = mean(BigDays.sfdfT$Lat),
                                             Lon = mean(BigDays.sfdfT$Lon),
                                             Year = 2018,
                                             CAPE = seq(0, 5, by = .25),
                                             DLBS = mean(BigDays.sfdfT$DLBS),
                                             SLBS = mean(BigDays.sfdfT$SLBS)),
                        type = "response",
                        se.fit = TRUE)
plot.predictions <- data.frame(P = Plevels[1],
                               probability = 1 - pnbinom(nCthresh, size = theta, mu = predictions$fit),
                               probabilityU = 1 - pnbinom(nCthresh, size = theta, mu = predictions$fit + 1 * predictions$se.fit),
                               probabilityL = 1 - pnbinom(nCthresh, size = theta, mu = predictions$fit - 1 * predictions$se.fit),
                               CAPE = seq(0, 5, by = .25) * 1000)

predictions  <- predict(modelFinalC, 
                        newdata = data.frame(P = Plevels[2],
                                             Lat = mean(BigDays.sfdfT$Lat),
                                             Lon = mean(BigDays.sfdfT$Lon),
                                             Year = 2018,
                                             CAPE = seq(0, 5, by = .25),
                                             DLBS = mean(BigDays.sfdfT$DLBS),
                                             SLBS = mean(BigDays.sfdfT$SLBS)),
                        type = "response",
                        se.fit = TRUE)

plot.predictions2 <- data.frame(P = Plevels[2],
                               probability = 1 - pnbinom(nCthresh, size = theta, mu = predictions$fit),
                               probabilityU = 1 - pnbinom(nCthresh, size = theta, mu = predictions$fit + 1 * predictions$se.fit),
                               probabilityL = 1 - pnbinom(nCthresh, size = theta, mu = predictions$fit - 1 * predictions$se.fit),
                               CAPE = seq(0, 5, by = .25) * 1000)

predictions  <- predict(modelFinalC, 
                        newdata = data.frame(P = Plevels[3],
                                             Lat = mean(BigDays.sfdfT$Lat),
                                             Lon = mean(BigDays.sfdfT$Lon),
                                             Year = 2018,
                                             CAPE = seq(0, 5, by = .25),
                                             DLBS = mean(BigDays.sfdfT$DLBS),
                                             SLBS = mean(BigDays.sfdfT$SLBS)),
                        type = "response",
                        se.fit = TRUE)

plot.predictions3 <- data.frame(P = Plevels[3],
                               probability = 1 - pnbinom(nCthresh, size = theta, mu = predictions$fit),
                               probabilityU = 1 - pnbinom(nCthresh, size = theta, mu = predictions$fit + 1 * predictions$se.fit),
                               probabilityL = 1 - pnbinom(nCthresh, size = theta, mu = predictions$fit - 1 * predictions$se.fit),
                               CAPE = seq(0, 5, by = .25) * 1000)


plot.predictions <- rbind(plot.predictions, plot.predictions2, plot.predictions3)

( p2 <- ggplot(plot.predictions, aes(x = CAPE, 
                                     y = probability, 
                                     color = factor(P * 100000, levels = rev(Plevels * 100000)),
                                     group = factor(P * 100000, levels = rev(Plevels * 100000)))) +
  geom_line(size = 2) +
  geom_ribbon(aes(ymin = probabilityL, ymax = probabilityU), color = NA, fill = "gray80", alpha = .2) +
  scale_y_continuous(limits = c(0, .3), labels = scales::percent) +
   scale_colour_manual(name = "Population", labels = c("8M", "4M", "1M"), values = cbPalette) + 
  ylab(paste("Probability of at least", nCthresh, "casualties")) +
  xlab("CAPE (J/kg)") +
  theme_minimal() )

library(patchwork)

p1 + p2
```

## Casualty predictions 2-D
```{r}
theta <- modelFinalC$theta
nCthresh <- 25

pgrid1 <- expand.grid(CAPE = seq(1, 5, .05),
                     DLBS = seq(1, 4, .05),
                     SLBS = mean(BigDays.sfdfT$SLBS),
                     P = 40, # 4 million
                     Year = 2018,
                     Lat = 42.5, # Sioux City IA and Birmingham AL
                     Lon = -96.4,
                     Location = "Sioux City, IA") 
pgrid2 <- expand.grid(CAPE = seq(1, 5, .05),
                     DLBS = seq(1, 4, .05),
                     SLBS = mean(BigDays.sfdfT$SLBS),
                     P = 40, # 4 million
                     Year = 2018,
                     Lat = 33.5,
                     Lon = -86.8,
                     Location = "Birmingham, AL") 
pgrid <- rbind(pgrid1, pgrid2)


predictions <- predict.glm(modelFinalC, newdata = pgrid, type = "response")
pgrid$probability <- 1 - pnbinom(nCthresh, size = theta, mu = predictions)

ggplot(pgrid, aes(x = CAPE * 1000, y = DLBS * 10, fill = probability)) +
          geom_raster(alpha = .8) +
          scale_fill_gradientn(colours = rev(cbPalette), name = "Probability", labels = scales::percent) +
          facet_wrap(~ Location) + 
          ylab("Deep-layer bulk shear (m/s)") +
          xlab("CAPE (J/kg)") + 
          theme_minimal()
```

Percentage point differences
```{r}
pgrid1 <- expand.grid(CAPE = seq(1, 5, .05),
                     DLBS = seq(1, 4, .05),
                     SLBS = mean(BigDays.sfdfT$SLBS),
                     P = 40, # 4 million
                     Year = 2018,
                     Lat = 42.5, # Sioux Falls IA 
                     Lon = -96.4) 
predictions <- predict.glm(modelFinalC, newdata = pgrid1, type = "response")
probs1 <- 1 - pnbinom(nCthresh, size = theta, mu = predictions)
pgrid1$probability <- probs1

pgrid2 <- expand.grid(CAPE = seq(1, 5, .05),
                     DLBS = seq(1, 4, .05),
                     SLBS = mean(BigDays.sfdfT$SLBS),
                     P = 40, # 4 million
                     Year = 2018,
                     Lat = 33.5, #Birmingham AL
                     Lon = -86.8) 
predictions <- predict.glm(modelFinalC, newdata = pgrid2, type = "response")
probs2 <- 1 - pnbinom(nCthresh, size = theta, mu = predictions)

pgrid1$ppDiff <- (probs2 - probs1) * 100

( p <- ggplot(pgrid1, aes(x = CAPE * 1000, y = DLBS * 10, fill = probability)) +
          geom_raster() +
          scale_fill_gradientn(colours = rev(cbPalette)) +
 #         geom_contour(aes(color = ..level.., z = ppDiff), breaks = seq(2, 12, by = 2), color = "black") +
          ylab("Deep-layer bulk shear (m/s)") +
          xlab("CAPE (J/kg)") + 
          theme_minimal() )
#direct.label(p, list("angled.boxes"))
```

## April 12, 2020 case https://www.spc.noaa.gov/cgi-bin-spc/getacrange.pl?date0=20200412&date1=20200412
```{r}
library(sf)

day1out <- read_sf(dsn = "day1otlk_20200412_1200-shp",
                         layer = "day1otlk_20200412_1200_torn") %>%
  st_transform(crs = 3857) 

#%>%
#  st_area()

quantile(BigDays.sfdfT$HullArea)
```

## Trends
```{r}
BigDays.sfdfT %>%
#  filter(Mo %in% c(3, 4, 5, 6)) %>%
  group_by(Year) %>%
  summarize(Avg = mean(SLBS)) %>%
ggplot(aes(x = Year, y = Avg)) +
  geom_point() +
  geom_smooth(method = lm)
```

Pllot the SPC Outlook
```{r}
tm_shape(day1out) + 
  tm_borders()
```

```{r}
cr <- RColorBrewer::brewer.pal(9, "Greys")
cr <- cr[-c(1:3)]
```

```{r}
Apr12 <- read.csv("April122020_Outbreak.csv", header = TRUE)
```

```{r}
Apr12.sfdfT = st_as_sf(Apr12, coords = c("LON", "LAT"), 
                 crs = 4326)
```


```{r}
day1out <- day1out[1:4,]

mycols <- colors()[c(51, 142, 585, 555)]

tm_shape(stateBorders) + 
  tm_borders(col = "gray70", alpha = 1) +
  tm_compass(size = 3, fontsize = 1, lwd = 2, color.dark = "gray70") + 
  tm_scale_bar(width = .3, size = 0.8, lwd = 1.75, color.dark = "gray70") +
  tm_layout(legend.bg.color = "white", 
            legend.text.size = .75, 
            attr.position = c("left", "bottom"), 
            inner.margins = c(.15, .15, .15, .15)) +
tm_shape(day1out, is.master = TRUE, projection = merc) +
  tm_fill("LABEL",
            title = "Probability",
            palette = mycols, alpha = .6) +
tm_shape(Apr12.sfdfT) +
  tm_bubbles(size = 0.05, col = "black", border.lwd = 0, alpha = 0.6 ) +
  tm_layout(#title = "April 12, 2020 1200 UTC \n Day 1 Convective Outlook", 
              title.position = c("center", "top"), 
              legend.title.size = 1.1,
              legend.position = c("right", "bottom"), 
              legend.stack = "horizontal",
              legend.frame = FALSE, 
              legend.text.size = 0.9, 
              legend.width = -0.1, 
              title.size = 1.2)
  
```