Zach_Thesis.Rmd

---
title: "Zach_Thesis"
author: "Zach"
date: "10/2/2021"
output: html_document
editor_options: 
  chunk_output_type: console
---

This file will document/organize work flow for the Master's thesis


##Abstract
Wildfires in the Apalachicola National Forest, USA From 1992-2018:
Relationships to Soil Dryness and Lightning

Anthropogenic climate change is shifting the risk of wildfires worldwide. The influence of climate change on fire hazards is diverse and complex and varies by region. Regional studies provide a more local perspective on the changing climate. In this study, the fire threat in the Southeast United States is investigated and the association between soil dryness and lightning-ignited wildfires in the Apalachicola National Forest (ANF) is examined. The relationship between pre-fire season dryness and the number of fires during the fire season is quantified with a statistical model. Soil dryness in the duff layer is estimated using the Keetch-Byram Drought Index (KBDI) which is a measure of moisture deficit based on temperature and net rainfall computed using daily summary-of-the-day values from the nearby first-order weather station. The probability of daily wildfire occurrence in the ANF is found to increase by 27% for every centimeter increase in daily soil moisture deficit and by 14% for every 100 daily lightning strikes. Given the strong statistical relationship between the KBDI during April and wildfires in the ANF from May to July, a seasonal prediction model is constructed to skillfully portend the occurrence of wildfires. Results show nearly a 20% increase in the average rate of wildfires from May through July for every one-centimeter increase in soil moisture deficit during April. Long-term upward trends in soil dryness are also shown, with the trend most pronounced during the driest months. This analysis indicates the potential for a greater risk of fires in the ANF in the future. Finally, the opportunity for future work is highlighted, and another model is proposed to further distinguish between a season with many small fires versus a season with a few big fires.


##Import libraries and data sets
```{r}
library(lubridate)
library(tidyverse)
library(ggpp)
library(ggpubr)
library(scales)
library(patchwork)
library(tidycensus)
library(sf)
library(tmap)
library(rgeos)
library(RColorBrewer)
library(rgdal)
library(USAboundaries)
library(units)
library(equatiomatic)
#library(MASS)
```

Daily summary of the day values from Tallahassee (TLH)
March 1,1940 - April 20, 2020
Includes observed 24hr precip, tmax, and tmin
```{r}
TLH.df <- read_csv(file = 'Data/TLH_Daily1940.csv') %>%
  rename(Date = DATE) %>%
  mutate(Year = year(Date), 
         month = month(Date, label = TRUE, abbr = TRUE),
         doy = yday(Date),
         MaxTemp = TMAX,
         MinTemp = TMIN,
         Rainfall24 = PRCP,
         Rainfall24 = replace_na(Rainfall24, 0),
         Rainfall24mm = Rainfall24 * 25.4)

#Fill in the missing temperature data
TLH.df$MaxTemp[TLH.df$Date == "2005-07-08"] <- 96
```

Calculate KBDI referencing functions that have been built. Append to TLH.df
```{r}
source("netRainfall.R")
TLH.df$NetR <- netRainfall(TLH.df$Rainfall24)

source("droughtIndex.R")
Q <- 269
R <- 59.23 #Average annual rainfall for TLH in inches
TLH.df$Ql <- droughtIndex(Q,R,TLH.df$MaxTemp,TLH.df$NetR)$Ql
TLH.df$Qlm <- TLH.df$Ql * .254  # tenth of an inch to mm
TLH.df$DroughtIndex <- droughtIndex(Q,R,TLH.df$MaxTemp,TLH.df$NetR)$DroughtIndex
```

Include a Ql centimeters column (drough index in cm) in the TLH.df data frame
```{r}
TLH.df <- TLH.df %>%
  mutate(Qlcm = Qlm/10)
```

ANF Boundaries data set:
```{r}
if(!"S_USA.NFSLandUnit" %in% list.files()){
  download.file("https://data.fs.usda.gov/geodata/edw/edw_resources/shp/S_USA.NFSLandUnit.zip",
                "S_USA.NFSLandUnit.zip")
unzip("S_USA.NFSLandUnit.zip")
}

NF_Bounds.sf <- st_read(dsn = "s_USA.NFSLandUnit.shp") %>%
  st_transform(crs = 3086)

anfbounds.sf <- NF_Bounds.sf %>%
  filter(NFSLANDU_2 == "Apalachicola National Forest")
```

Import County Boundaries of the forest. This will make the bounds consistent across both the fire data set and the lightning data set.
```{r}
ANFcountyBounds.sf <- st_read(dsn = "ANFCounties.shp") %>%
  st_transform(crs = 3086)
```

Fires data set:
bounded on counties to match the lightning data boundary. Includes Liberty, Wakulla, Franklin, and Leon. ANF is within these county boundaries.Includes all fire types.
```{r}
#if(!"Fires2018" %in% list.files()){
  #download.file("https://www.fs.usda.gov/rds/archive/products/RDS-2013-0009.5/RDS-2013-0009.5_GPKG.zip",
                #"Data/updatedFireData/Fires2018.zip")
  #unzip("Data/updatedFireData/Fires2018.zip",
        #exdir = "Data/updatedFireData")
#}

county_Fires.sf <- st_read(dsn = "Data/updatedFireData/Data/FPA_FOD_20210617.gpkg", layer = "Fires") %>%
  filter(STATE == "FL") %>%
  st_transform(crs = st_crs(ANFcountyBounds.sf)) %>%
  st_intersection(ANFcountyBounds.sf) %>%
  select(FOD_ID, FIRE_NAME, FIRE_YEAR, DISCOVERY_DATE, NWCG_CAUSE_CLASSIFICATION, NWCG_GENERAL_CAUSE, FIRE_SIZE, FIRE_SIZE_CLASS, LATITUDE, LONGITUDE, NAME, FIPS_CODE, Shape)

#Reformat date column to not include time. Rename to match other merged columns
county_Fires.sf$DISCOVERY_DATE <- as.Date(county_Fires.sf$DISCOVERY_DATE)
county_Fires.sf <- county_Fires.sf %>%
  rename(DISCOVERY_ = DISCOVERY_DATE)

#county lightning fires
county_LF.sf <- county_Fires.sf %>%
  filter(NWCG_GENERAL_CAUSE == "Natural")
```

Bounded to the ANF
```{r}
#all fires
anf_Fires.sf <- county_Fires.sf %>%
  st_transform(crs = st_crs(anfbounds.sf)) %>%
  st_intersection(anfbounds.sf)

#filtered to lighting fires
anf_LF.sf <- anf_Fires.sf %>%
  filter(NWCG_GENERAL_CAUSE == "Natural")
```


lightning data set:
```{r}
#gives lighting data for each individual county based on fips code
lightning.df <- list.files(path = "C:/Users/zlaw9/OneDrive/GITHUB/KDBI code/FIPS_LightningData",
                  pattern = "*.csv", full.names = TRUE) %>%
  lapply(read_csv) %>%
  bind_rows

lightning.df <- lightning.df %>%
  mutate(SEQDAY = as.Date(SEQDAY),
         DAY = day(SEQDAY),
         MONTH = month(SEQDAY),
         YEAR = year(SEQDAY)) %>%
#missing data beyond 5/20/2013 for FIPS #12037. Remove last 5 rows from data set. This addresses the parsing failure and removes data that is NA.
  filter(SEQDAY <=  "2013-05-20")

#make FCOUNT_NLDN column numeric
lightning.df$FCOUNT_NLDN <- as.numeric(lightning.df$FCOUNT_NLDN)

#Combine total lightning strikes in ANF counties by date.
lightning.df <- lightning.df %>%
  group_by(SEQDAY) %>%
  summarise(LightningCount = sum(FCOUNT_NLDN))
```

Viewing lightning fires within the Apalachicola National Forest (437 observations)
```{r}
tmap_mode("view")

tm_shape(anfbounds.sf) +
  tm_borders()

tm_shape(anf_LF.sf) +
  tm_dots(col = "orange")
```

##Exploratory Analysis

Cause of Fires in the Apalachicola National Forest
Natural fires account for nearly 40% of fires in the ANF
```{r}
#table(anf_Fires.sf$NWCG_GENERAL_CAUSE)

df <- anf_Fires.sf %>%
  st_drop_geometry() %>%
  group_by(NWCG_GENERAL_CAUSE) %>%
  summarize(nF = n(),
            perF = nF/nrow(anf_Fires.sf))

ggplot(df,
       mapping = aes(y = reorder(NWCG_GENERAL_CAUSE, perF),
                     x = perF,
                     fill = perF)) +
  geom_col() +
  scale_fill_distiller(palette = "Oranges",
                       direction = 1,
                       guide = "none") +
  scale_x_continuous(labels = percent) +
  ylab("") + xlab("") +
  labs(title = "Lightning is the predominant spark for wildfires\nin the Apalachicola National Forest",
       subtitle = "Based on data from 1992-2018",
       caption = "Data source: Short, Karen (2021)")
  #theme_minimal()
  #theme(plot.title = element_text(size = 15),
        #axis.text = element_text(size = 12))

#dev.off()
```

Locations of lightning fires in the ANF
```{r}
FL_Counties.sf <- us_counties(states = "FL",
                              resolution = "high") %>%
  st_transform(crs = st_crs(anfbounds.sf)) %>%
  st_crop(st_bbox(st_buffer(anfbounds.sf, dist = .25)) )

ggplot() +
  geom_sf(data = FL_Counties.sf, fill = "transparent", col = "gray80") +
  geom_sf(data = anfbounds.sf, fill = "transparent") +
  geom_sf(data = anf_LF.sf,
          mapping = aes(col = FIRE_SIZE_CLASS)) +
  scale_color_brewer(palette = "Oranges",
                     direction = 1,
                     name = "Size Class") +
  theme_bw() 
  #labs(title = "Location of natural-caused wildfires in the Apalachicola National Forest (1992-2018)",
  #     subtitle = "Darker color points indicates the fire resulted in a larger burn area",
  #     caption = "Data source: Short, Karen (2021)")
```

Plot the number of lightning fires per month
Over 80% of lightning-sparked wildfires in the ANF occur during May-July
```{r}
wfplot.df <- anf_LF.sf %>%
  st_drop_geometry() %>%
  #mutate(MonthF = factor(month.name[month(DISCOVERY_)],
  mutate(MonthF = factor(month.abb[month(DISCOVERY_)],
                         #levels = rev(month.abb),
                         #levels = rev(month.name), 
                         levels = month.abb,
                         ordered = TRUE)) %>%
  group_by(MonthF, .drop = FALSE) %>%
  summarize(nF = n(),
            perF = nF/nrow(anf_LF.sf))

ggplot(data = wfplot.df,
       mapping = aes(y = MonthF, 
                     x = perF,
                     fill = perF)) +
  geom_col() +
  scale_fill_distiller(palette = "Oranges",
                       direction = 1,
                       guide = "none") +
  scale_x_continuous(labels = percent) +
  xlab("") + ylab("") +
  labs(title = "Over 80% of lightning-sparked wildfires in the\nApalachicola National Forest occur during May-July",
       subtitle = "Percentage of all lightning-sparked wildfires by month",
       caption = "Period of record: 1992-2018, Data source: Short, Karen (2021)") +
  theme_minimal() +
  theme(plot.title = element_text(size = 15),
        axis.text = element_text(size = 10)) 
  #theme_dark()
```

Plot average monthly rainfall and max temperatures
```{r}
df <- TLH.df %>%
  filter(Year >= 1943 & Year <= 2020) %>%
  group_by(month) %>%
  summarize(totalRain = sum(Rainfall24mm / 10),
            avgMonthlyTotalRain = totalRain/(2020 - 1943 + 1),
            avgDailyHighTemp = mean(MaxTemp))

p1 <- ggplot(data = df,
             mapping = aes(x = month, y = avgMonthlyTotalRain, fill = avgMonthlyTotalRain)) +
        geom_col() +
        scale_fill_distiller(palette = "Greens",
                         direction = 1, 
                         guide = 'none') +
        theme_minimal() +
        labs(x = "", y = "Rainfall (cm)", title = "A")

p2 <- ggplot(data = df,
             mapping = aes(x = month, y = avgDailyHighTemp, color = avgDailyHighTemp)) +
        geom_point(size = 4) +
        scale_color_distiller(palette = "Reds",
                              direction = 1, 
                              guide = 'none') +
        scale_y_continuous(limits = c(60, 100)) +
        theme_minimal() +
        labs(x = "", y = "Temperature (°F)", title = "B")

library(patchwork)
p1 / p2
```


Average number of lightning strikes in the ANF
```{r}
lightningAvg.df <- lightning.df %>%
  filter(year(SEQDAY) <= 2012) %>%
  #mutate(MonthF = factor(month.name[month(SEQDAY)],
  mutate(MonthF = factor(month.abb[month(SEQDAY)],
                         #levels = rev(month.name),
                         levels = month.abb,
                         ordered = TRUE)) %>%
  group_by(MonthF, .drop = FALSE) %>%
  summarize(avgLightningCount = mean(LightningCount, na.rm = TRUE))
```

Plot average lightning
```{r}
lightningAvg.df %>%
  ggplot(mapping = aes(y = MonthF,
                       x = avgLightningCount,
                       fill = avgLightningCount)) +
  geom_col() +
  #scale_fill_gradientn(colors = brewer.pal(4, "OrRd"),
  scale_fill_gradientn(colors = brewer.pal(4, "Blues"),                     
                       guide = "none") +
  labs(x = "", y = "",
       title = "Lightning in the Apalachicola National Forest peaks from June-August",
       subtitle = "Daily average number of cloud-to-ground lightning strikes by month",
       caption = "Period of record 1986 - 2012, Data Source: NCEI") +
  #theme_dark()
  theme_minimal()
```

```{r}
p3 <- ggplot(data = wfplot.df,
       mapping = aes(x = MonthF, 
                     y = perF,
                     fill = perF)) +
       #mapping = aes(x = perF,
                     #y = MonthF,
                     #fill = perF))+
  geom_col() +
  scale_fill_distiller(palette = "Oranges",
                       direction = 1,
                       guide = "none") +
  scale_y_continuous(labels = scales::percent_format(1L)) +
  #scale_x_continuous(labels = scales::percent_format(1L)) +
  labs(x = "", y = "Wildfire frequency",
       title = "A") +
  #labs(x = "wildfire frequency", y = "") +
  #labs(title = "Over 80% of natural wildfires in the Apalachicola National Forest\noccur during May-July",
  #labs(title = "Percentage of all natural\nwildfires by month",
       #subtitle = "Percentage of all natural wildfires by month",
       #caption = "Period of record: 1992-2018, Data source: Short, Karen (2021)") +
  theme_minimal()

p4 <- ggplot(data = lightningAvg.df,
       mapping = aes(x = MonthF, 
                     y = avgLightningCount,
                     fill = avgLightningCount)) +
  geom_col() +
  scale_fill_distiller(palette = "Blues",
                       direction = 1,
                       guide = "none") +
  labs(x = "", y = "Cloud-to-ground strikes",
       title = "B") +
       #title = "Lightning in the Apalachicola National Forest peaks from June-August",
       #title = "Daily average cloud-to-ground\nlightning strikes",
       #subtitle = "Daily average cloud-to-ground lightning strikes by month",
       #caption = "Period of record: 1986-2012, Data source: NCEI") +
  theme_minimal() 

p3/p4
```

Create three panel plot with precip, lightning, and wildfires
```{r}
(p1/p4) | p3
```


Monthly average soil moisture deficit.
```{r}
TLH.df %>%
  filter(Year >= 1941 & Year <= 2019) %>%
  mutate(MonthF = factor(month.name[month(Date)], 
                         levels = rev(month.name), 
                         ordered = TRUE)) %>%
  group_by(MonthF, .drop = FALSE) %>%
  summarize(AvgSoilMoistureDeficit = mean(Qlm)) %>%
ggplot(mapping = aes(y = MonthF, 
                     x = AvgSoilMoistureDeficit,
                     fill = AvgSoilMoistureDeficit)) +
  geom_col() +
  scale_fill_gradientn(colors = terrain.colors(5),
                       guide = 'none') +
  labs(x = "", y = "",
       title = "May through November is the dry season in the Apalachicola National Forest",
       subtitle = "Average dryness (mm)",
       caption = "Period of record: 1941-2019, Data source: NWSFO Tallahassee") +
  theme_dark() 
```

Another interpretation of soil moisture plot
```{r}
TLH.df %>%
  filter(Year >= 1943) %>%
  group_by(month, Year) %>%
  #summarise(Avg = mean(Qlm)) %>%
  summarise(Avg = mean(Qlcm)) %>%
  ggplot(aes(x = Year, y = Avg, color = Avg)) +
  geom_smooth(method = lm, se = FALSE, color = "gray70") +
  geom_point() +  
  scale_color_gradientn(colors = terrain.colors(5), guide = "none") +
  #scale_y_continuous(limits = c(0, 201)) +
  scale_y_continuous(limits = c(0, 20)) +
  scale_x_continuous(limits = c(1940, 2020), breaks = c(1950, 1980, 2010)) +
  ylab("") + xlab("") +
  facet_wrap(~ month, ncol = 12) +
  #theme_dark() +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Risk of wildfires is increasing in the Apalachicola National Forest.",
       subtitle = "Monthly average dryness (cm) by year with trend line (gray)", 
       caption = "Period of record: 1943-2020, Data source: NWSFO Tallahassee") 
```

##General linear regression - Study how KBDI impacts lightning wildfire occurrence.

Filter the KBDI data set to match the time frame of the fire data set (1992 - 2018)
```{r}
TLH_fireYrs.df <- TLH.df %>%
  filter(Year >= 1992 & Year <= 2018)
```

create a data frame for the number of fires that occurred on a given day
```{r}
countfires.df <- anf_LF.sf %>%
  count(DISCOVERY_) %>%
  rename(FireCount = n, Date = DISCOVERY_)
```

Merge the KBDI and fire count data sets
```{r}
KBDI_fire.df <- left_join(x = TLH_fireYrs.df, y = countfires.df, by = c("Date" = "Date")) %>%
  select(Date, FireCount, Ql, Qlm, Qlcm, DroughtIndex)

KBDI_fire.df$fireBool <- !is.na(KBDI_fire.df$FireCount)
```


Binomial Regression Model. Predicting the probability of wildfire occurrence using soil dryness as a predictor.
Ql units = 100th inches
```{r}
glmQl <- glm(formula = fireBool~Ql, family = binomial, data = KBDI_fire.df)

pred_glmQl <- predict(object = glmQl,
                      type = "response",
                      se.fit = TRUE)

low_bound <- pred_glmQl$fit - (1.96*pred_glmQl$se.fit)
up_bound <- pred_glmQl$fit + (1.96*pred_glmQl$se.fit)

ggplot(mapping = aes(x = KBDI_fire.df$Ql, y = pred_glmQl$fit)) +
  geom_ribbon(aes(ymin = low_bound, ymax = up_bound), fill = "grey") +
  geom_line(color = "blue") +
  #stat_regline_equation() +
  #stat_cor(aes(label = ..rr.label..), label.x = 0, label.y = .155) +
  ylab("Predicited Probablility") +
  xlab("Daily KBDI") +
  labs(
    title = "Predicted Probability of Lightning Fire Occurrence Based on Daily KBDI Values",
    subtitle = "Line of Fitted Values With a 95% Confidence Interval",
    caption = "Period of Record: 1992 to 2018")
```

```{r}
extract_eq(glmQl, use_coefs = TRUE, coef_digits = 4, fix_signs = FALSE)
```
$$
\log\left[ \frac { \widehat{P( \operatorname{fireBool} = \operatorname{TRUE} )} }{ 1 - \widehat{P( \operatorname{fireBool} = \operatorname{TRUE} )} } \right] = -5.5013 + 0.0049(\operatorname{Ql})
$$

Ql is statistically significant
```{r}
summary(glmQl)
```

Filtered to fire season months
```{r}
KBDI_fireFS.df <- KBDI_fire.df %>%
  filter(month(Date) == 05 | month(Date) == 06 | month(Date) == 07)

glmQlFS <- glm(formula = fireBool ~ Ql, family = binomial, data = KBDI_fireFS.df)

pred_glmQlFS <- predict(object = glmQlFS,
                      type = "response",
                      se.fit = TRUE)

low_bound <- pred_glmQlFS$fit - (1.96*pred_glmQlFS$se.fit)
up_bound <- pred_glmQlFS$fit + (1.96*pred_glmQlFS$se.fit)

ggplot(mapping = aes(x = KBDI_fireFS.df$Ql, y = pred_glmQlFS$fit)) +
  geom_ribbon(aes(ymin = low_bound, ymax = up_bound), fill = "grey") +
  geom_line(color = "blue") +
  ylab("Predicited Probablility") +
  xlab("Daily Ql") +
  labs(
    title = "Predicted Probability of Fire Occurrence Based on Daily KBDI Values\nFiltered By Fire Season Months",
    subtitle = "With 95% Confidence Interval",
    caption = "Fire Season Months (May - July) 1992 - 2018")
```

Ql is statistically significant
```{r}
summary(glmQlFS)
```


Adjust Binomial Regression model to represent KBDI in cm (SI units)
```{r}
glmQlcm <- glm(formula = fireBool~Qlcm, family = binomial, data = KBDI_fire.df)

pred_glmQlcm <- predict(object = glmQlcm,
                      type = "response",
                      se.fit = TRUE)

low_bound <- pred_glmQlcm$fit - (1.96*pred_glmQlcm$se.fit)
up_bound <- pred_glmQlcm$fit + (1.96*pred_glmQlcm$se.fit)

ggplot(mapping = aes(x = KBDI_fire.df$Qlcm, y = pred_glmQlcm$fit)) +
  geom_ribbon(aes(ymin = low_bound, ymax = up_bound), fill = "grey") +
  geom_line(color = "blue") +
  ylab("Predicited Probablility") +
  xlab("Soil Moisture Deficit (cm)") +
  #labs(
    #title = "Predicted Probability of Lightning Fire Occurrence\nBased on Daily Soil Moisture Deficit Values",
    #subtitle = "Line of Fitted Values With a 95% Confidence Interval",
    #caption = "Period of Record: 1992 to 2018")
  theme_minimal()
```

```{r}
extract_eq(glmQlcm, use_coefs = TRUE, coef_digits = 4, fix_signs = FALSE)
```
$$
\log\left[ \frac { \widehat{P( \operatorname{fireBool} = \operatorname{TRUE} )} }{ 1 - \widehat{P( \operatorname{fireBool} = \operatorname{TRUE} )} } \right] = -5.5013 + 0.1937(\operatorname{Qlcm})
$$

predicting the probability of fire occurrance based on Qlm values
```{r}
temp.df <- data.frame(Qlcm = 19.60179)
predict.glm(glmQlcm, newdata = temp.df, type = "response", se.fit = TRUE)
```


##General linear regression - Study how the number of lightning strikes on a given day impacts lightning wildfire occurrence

In modeling with lightning, we need to use the county_LF.sf data set because the lighting data set is bounded to the county boarders
```{r}
count_cFires.df <- county_LF.sf %>%
  count(DISCOVERY_) %>%
  rename(FireCount = n, Date = DISCOVERY_)
```

Each data set must have the same temporal bounds
```{r}
count_cFiresYrs.df <- count_cFires.df %>%
  filter(year(Date) <= 2012)

lightningYrs.df <- lightning.df %>%
  filter(year(SEQDAY) >= 1992 & year(SEQDAY) <= 2012)
```

merge lightning count data and counted fire data
```{r}
light_fires.df <- left_join(x = lightningYrs.df, y = count_cFiresYrs.df, by = c("SEQDAY" = "Date")) %>%
  rename(Date = SEQDAY) %>%
  select(Date, LightningCount, FireCount)

light_fires.df$FireBool <- !is.na(light_fires.df$FireCount)
```

binomial regression
```{r}
glmLightning <- glm(formula = FireBool ~ LightningCount, family = binomial, data = light_fires.df)

#probabilities based on glm
pred_glmLightning <- predict(object = glmLightning,
                             type = "response",
                             se.fit = TRUE)

#creating confidence interval
low_bound <- pred_glmLightning$fit - (1.96*pred_glmLightning$se.fit)
up_bound <- pred_glmLightning$fit + (1.96*pred_glmLightning$se.fit)

ggplot(mapping = aes(x = light_fires.df$LightningCount, y = pred_glmLightning$fit)) +
  geom_ribbon(aes(ymin = low_bound, ymax = up_bound), fill = "grey") +
  geom_line(color = "blue") +
  ylab("Predicited Probablility") +
  xlab("Number of Lightning Strikes") +
  #labs(
    #title = "Predicted Probability of a Fire Occurrence Based on the\nNumber of Lightning Strikes in a Day",
    #subtitle = "With a 95% Confidence Interval",
    #caption = "Period of Record: 1992 to 2012")
  theme_minimal()
```

lightning count is statistically significant to occurrence of wildfires.
```{r}
summary(glmLightning)
```

Include fitted values in the table
```{r}
light_fires.df$fit <- glmLightning$fitted.values
```

Exploring the confidence intervals
```{r}
temp.df <- data.frame(LightningCount = 2794)
predict.glm(glmLightning, newdata = temp.df, type = "response", se.fit = TRUE)
```


```{r}
extract_eq(glmLightning, use_coefs = TRUE, coef_digits = 4, fix_signs = FALSE)
```
$$
\log\left[ \frac { \widehat{P( \operatorname{FireBool} = \operatorname{TRUE} )} }{ 1 - \widehat{P( \operatorname{FireBool} = \operatorname{TRUE} )} } \right] = -3.2102 + 0.0012(\operatorname{LightningCount})
$$

Filtered to fire season months
```{r}
light_firesFS.df <- light_fires.df %>%
  filter(month(Date) == 05 | month(Date) == 06 | month(Date) == 07)

glmLightningFS <- glm(formula = FireBool ~ LightningCount, family = binomial, data = light_firesFS.df)

#probabilities based on glm
pred_glmLightning <- predict(object = glmLightningFS,
                             type = "response",
                             se.fit = TRUE)

#creating confidence interval
low_bound <- pred_glmLightning$fit - (1.96*pred_glmLightning$se.fit)
up_bound <- pred_glmLightning$fit + (1.96*pred_glmLightning$se.fit)

ggplot(mapping = aes(x = light_firesFS.df$LightningCount, y = pred_glmLightning$fit)) +
  geom_ribbon(aes(ymin = low_bound, ymax = up_bound), fill = "grey") +
  geom_line(color = "blue") +
  ylab("Predicited Probablility") +
  xlab("Number of Lightning Strikes") +
  labs(
    title = "Predicted Probability of a Fire Occurrence Based on the Number of Lightning Strikes in a Day During the Fire Season",
    subtitle = "With a 95% Confidence Interval",
    caption = "Period of Record: 1992 to 2012 Fire Seasons (May - July)")
```



Include daily soil dryness and daily lightning count in one model.

First, create one data frame which includes both data sets.
combine KBDI_fire.df and light_fires.df which were used in building each individual model.
Both data frames start on the same date, left join on light_fires.df because data period stops in 2012.
Important to note in the merge. The bounds of the forest are not quite the same. The lightning data frame bounds the data by county and the dryness data frame bounds the data by forest boundaries. We will assume the lightning count data is an approximation for the forest count here and will reference the fire count and firebool from the soil dryness data frame.
```{r}
logisticModel.df <- left_join(light_fires.df, KBDI_fire.df, by = c("Date" = "Date"))
logisticModel.df <- logisticModel.df %>%
  dplyr::select("Date", "FireCount.y", "fireBool", "Qlcm", "LightningCount") %>%
  dplyr::rename(FireCount = FireCount.y)
```

Create logistic regression with both predictors
```{r}
logistic_DL <- glm(formula = fireBool ~ Qlcm + LightningCount, family = binomial(link = "logit"), data = logisticModel.df)
```

Both daily soil dryness and daily lightning count are statistically significant
```{r}
summary(logistic_DL)
```

add fitted values within a column in the data frame
```{r}
logisticModel.df$fit <- logistic_DL$fitted.values
```


Equation of model
```{r}
extract_eq(logistic_DL, use_coefs = TRUE, coef_digits = 4, fix_signs = FALSE)
```
$$
\log\left[ \frac { \widehat{P( \operatorname{fireBool} = \operatorname{TRUE} )} }{ 1 - \widehat{P( \operatorname{fireBool} = \operatorname{TRUE} )} } \right] = -6.1967 + 0.2368(\operatorname{Qlcm}) + 0.0013(\operatorname{LightningCount})
$$

Further visualize and plot the binomial model with two predictors (Soil Dryness and Lightning count).
Set lightning count to its mean and plot predicted probabilities across the range of soil moisture deficit
Set soil moisture deficit to its mean and plot predicted probabilities across the range of lightning count 

```{r}
logisticModel_Avg.df <- logisticModel.df %>%
  select("Date", "fireBool", "Qlcm", "LightningCount") %>%
  mutate(LightningAvg = mean(LightningCount), QlcmAvg = mean(Qlcm))

logistic_SM <- glm(formula = fireBool ~ Qlcm + LightningAvg, family = binomial(link = "logit"), data = logisticModel_Avg.df)

pred_SM <- predict(object = logistic_SM,
                         type = "response",
                         se.fit = TRUE)

#creating confidence interval
low_boundSM <- pred_SM$fit - (1.96*pred_SM$se.fit)
up_boundSM <- pred_SM$fit + (1.96*pred_SM$se.fit)

SM_plot <- ggplot(mapping = aes(x = logisticModel_Avg.df$Qlcm, y = logistic_SM$fit)) +
  geom_ribbon(aes(ymin = low_boundSM, ymax = up_boundSM), fill = "grey") +
  geom_line(color = "blue") +
  theme_minimal() +
  #ylab("Predicited Probablility") +
  #theme(axis.title.y = element_text(hjust = .8)) +
  ylab("Daily Fire Probability") +
  xlab("Soil Dryness (cm)") +
  theme(aspect.ratio = 6/5)


logistic_L <- glm(formula = fireBool ~ QlcmAvg + LightningCount, family = binomial(link = "logit"), data = logisticModel_Avg.df)

pred_L <- predict(object = logistic_L,
                         type = "response",
                         se.fit = TRUE)

#creating confidence interval
low_boundL <- pred_L$fit - (1.96*pred_L$se.fit)
up_boundL <- pred_L$fit + (1.96*pred_L$se.fit)

L_plot <- ggplot(mapping = aes(x = logisticModel_Avg.df$LightningCount, y = logistic_L$fit)) +
  geom_ribbon(aes(ymin = low_boundL, ymax = up_boundL), fill = "grey") +
  geom_line(color = "blue") +
  #ylab("Predicited Probablility") +
  ylab("Daily Fire Probability") +
  xlab("Daily Lightning Strike Count") +
  theme_minimal() +
  theme(aspect.ratio = 6/5)

patchworkplot <- SM_plot+L_plot

patchworkplot + plot_annotation(
    title = "Predicted probability of daily fire occurrence for each predictor",
    subtitle = "Assuming the second predictor is equal to its mean")
```

Add fitted values to the logisticModel_Avg.df dataframe
```{r}
logisticModel_Avg.df$Qlcm_fit <- logistic_SM$fitted.values
logisticModel_Avg.df$LightningCount_fit <- logistic_L$fitted.values
```

adjust logistic model data frame so it can be used as a table snapshot in a PowerPoint slide
```{r}
Logistic_PPointTable <- logisticModel.df %>%
  select("Date":"LightningCount") %>%
  rename(KBDI_cm = Qlcm)
```

Model vizualization
```{r}
newdat <- expand.grid(Qlcm = seq(0, 20, 0.5),
                      LightningCount = seq(0, 5600, 100))
z <- predict(logistic_DL, newdata = newdat, type = "response")
newdat$br <- z

ggplot(newdat, aes(x = Qlcm, y = LightningCount, fill = br)) + 
        geom_tile() +
        scale_fill_distiller(palette = "Oranges", direction = 1) +
        labs(fill = "Predicted Probability") +
        xlab("Soil Dryness (cm)") + ylab("Lightning Strike Count") +
        theme_minimal() +
        theme(aspect.ratio = 1) +
  labs(title = "Predicted probability of daily fire occurrence")
```



##Study how the burned area over the past year may impact fire occurance.

Create a data frame to include the area that was burned over the past year.
```{r}
anf_Fires.sf <- anf_Fires.sf %>%
  mutate(InSeason = DISCOVERY_ >= as.Date(paste0(as.character(FIRE_YEAR), "-04", "-30")) & DISCOVERY_ <= as.Date(paste0(as.character(FIRE_YEAR+1), "-04", "-29"))) %>%
   mutate(burnedGroupStart = as.Date(paste0(as.character(FIRE_YEAR), "-04", "-30"))) %>%
   mutate(burnedGroupEnd = as.Date(paste0(as.character(FIRE_YEAR+1), "-04", "-29")))

for(i in 1:length(anf_Fires.sf$InSeason)){
  if (anf_Fires.sf$InSeason[i] == FALSE){
    anf_Fires.sf$burnedGroupStart[i] <- as.Date(paste0(as.character(anf_Fires.sf$FIRE_YEAR[i]-1), "-04", "-30"))
    anf_Fires.sf$burnedGroupEnd[i] <- as.Date(paste0(as.character(anf_Fires.sf$FIRE_YEAR[i]), "-04", "-29"))
  }
}

anf_Fires.sf <- anf_Fires.sf %>%
  unite(col = "burnedGroupInterval", c("burnedGroupStart", "burnedGroupEnd"), sep = " -- ") %>%
  #rearrange order of columns
  select(FOD_ID:DISCOVERY_, burnedGroupInterval, NWCG_CAUSE_CLASSIFICATION:Shape)
```

Group data by burned periods and sum the total acres burned over that period
Note this does not include fires in the first 3 months of the data set
```{r}
burnedArea.sf <- anf_Fires.sf %>%
  group_by(burnedGroupInterval) %>%
  summarise(acresBurned = sum(FIRE_SIZE)) %>%
  #create column to define the year the data will be used to forecast for
  mutate(forecastYear = substr(burnedGroupInterval, 15, 18)) %>%
  #reorder columns
  select(forecastYear, burnedGroupInterval:acresBurned) %>%
  #convert forecastYear to an int/ double so it will merge
  transform(forecastYear = as.numeric(forecastYear)) %>%
  filter(forecastYear <= 2018)
```

Create a count of natural wildfires occurring in the anf during the fire season. This will be merged with the burned acres data frame.
```{r}
seasonFires.sf <-  anf_LF.sf %>%
  filter(month(DISCOVERY_) == 05 | month(DISCOVERY_) == 06 | month(DISCOVERY_) == 07) %>%
  count(FIRE_YEAR) %>%
  rename(Year = FIRE_YEAR, nFIRES = n) %>%
  st_set_geometry(NULL)
```

Merge burned data frame with the number of lightning wildfires that occurred in the Apalachicola National Forest for each year.
```{r}
burnedArea.sf <- left_join(x = burnedArea.sf, y = seasonFires.sf, by = c("forecastYear" = "Year"))

burnedArea.sf$nFIRES <- replace_na(burnedArea.sf$nFIRES, 0)
```

The model finds that the number of acres burned over the past year is not statistically significant when modeled by itself.
p-value = 0.0525
```{r}
burnedGLM <- glm(formula = nFIRES ~ acresBurned, family = poisson, data = burnedArea.sf)

pred_burnedGLM <- predict(object = burnedGLM,
                         type = "response",
                         se.fit = TRUE)

summary(burnedGLM)

cor(burnedArea.sf$acresBurned, burnedArea.sf$nFIRES)

#correlation of model residuals
cor(burnedArea.sf$acresBurned, resid(burnedGLM))
```

Explore area as a percentage of the forest that was burned.
```{r}
st_area(anfbounds.sf) #2564052729 m^2
drop_units(st_area(anfbounds.sf))/1000000 #2564.053 kilometers^2
drop_units(st_area(anfbounds.sf))/4047 #approximately 633,568.7 acres
```

Create a column within the burnedArea.sf data frame to represent the percentage of the forest burned prior to the fire season
```{r}
burnedArea.sf <- burnedArea.sf %>%
  mutate(percentBurned = (acresBurned/(drop_units(st_area(anfbounds.sf))/4047))*100) %>%
  #reorder columns
  select(forecastYear:acresBurned, percentBurned, nFIRES:Shape)
```

create model based on the percent of forest that was burned.
We find the percent of forest burned is not statistically significant (p-value = 0.0525, same value)
```{r}
percentBurnedGLM <- glm(formula = nFIRES ~ percentBurned, family = poisson, data = burnedArea.sf)

pred_percentBurnedGLM <- predict(object = percentBurnedGLM,
                         type = "response",
                         se.fit = TRUE)

summary(percentBurnedGLM)

cor(burnedArea.sf$percentBurned, burnedArea.sf$nFIRES)

#correlation of model residuals
cor(burnedArea.sf$percentBurned, resid(percentBurnedGLM))
```


##Model Preseason Soil Dryness as a Predictor

Create a bar chart showing the number of fires per year over the data period
```{r}
seasonFires.sf %>%
  ggplot(mapping = aes(x = Year, 
                       y = nFIRES, 
                       fill = nFIRES)) +
  geom_col() +
  scale_fill_gradientn(colors = brewer.pal(4, "OrRd"),
                       guide = "none") +
  scale_x_continuous(breaks = seq(1995, 2015, 5)) +
  ylab("Number of Fires") +
  xlab("Year") +
  labs(
    title = "Total Number of Fires During the Fire Season",
    caption = "Period of Record: 1992 to 2018") +
  theme_minimal()
```


Add April 30th KBDI to the seasonFires.sf data frame. This will represent preliminary season dryness
```{r}
April30.df <- KBDI_fire.df %>%
  filter(month(Date) == 04 & day(Date) == 30) %>%
  select("Date", "Ql", "Qlm", "Qlcm") %>%
  rename(KBDI_April30 = Ql) %>%
  mutate(Year = year(Date))

#make April30.df left table so that years 1997 and 2005 are filled as zeros
seasonFires.sf <- left_join(x = April30.df, y = seasonFires.sf, by = c("Year" = "Year")) %>%
  select("Year", "nFIRES", "KBDI_April30")
seasonFires.sf$nFIRES <- replace_na(seasonFires.sf$nFIRES, 0)
```

add the area burned over the prior season to seasonFires.sf
```{r}
seasonFires.sf <- left_join(x = seasonFires.sf, y = burnedArea.sf, by = c("Year" = "forecastYear")) %>%
  select("Year":"percentBurned") %>%
  rename(nFIRES = nFIRES.x)
```

Create a temp data frame to start in 1993 because we do not have all acres burned prior to 1992 to be included in the 1992 preseason count.
```{r}
temp <- seasonFires.sf %>%
  filter(Year >= 1993)
```



1) set soil moisture deficit units to cm (SI units)
2) fit a negative binomial regression model

Adjust data frames so that units are in cm
```{r}
seasonFires_cm.sf <- left_join(x = April30.df, y = seasonFires.sf, by = c("Year" = "Year")) %>%
  dplyr::select("Year", "nFIRES", "Qlcm")
seasonFires.sf$nFIRES <- replace_na(seasonFires.sf$nFIRES, 0)
```

Model reflecting only soil moisture deficit
negative binomial regression
```{r}
library(MASS)

nbr1 <- glm.nb(nFIRES ~ Qlcm, data = seasonFires_cm.sf)

#summary(nbr1)

detach(package:MASS, unload=TRUE)
```

```{r}
extract_eq(nbr1, use_coefs = TRUE, coef_digits = 4, fix_signs = FALSE)
```

$$
\log ({ \widehat{E( \operatorname{nFIRES} )} })  = 0.7515 + 0.1773(\operatorname{Qlcm})
$$

visualize the negative binomial regression model
```{r}
library(MASS)

ggplot(seasonFires_cm.sf, aes(x = Qlcm, y = nFIRES)) +
  geom_point() +
  geom_smooth(method = "glm.nb") + #, formula = y ~ x,
              #method.args = list(family = poisson(link = "log"))) +
  ylab("Number of Fires in the ANF (May - July)") +
  xlab("Soil Dryness on April 30th (cm)") +
  theme_minimal() +
  theme(aspect.ratio = 1) +
  labs(
    title = "Predicted number of fires during\nthe fire season")

detach(package:MASS, unload=TRUE)
```

Exploring predicted values and confidence intervals
```{r}
temp.df <- data.frame(Qlcm = 14.768622)
predict.glm(nbr1, newdata = temp.df, type = "response", se.fit = TRUE)
```

Explore the Predicted Values
```{r}
seasonFires_cm.sf$nbr1_fit <- nbr1$fitted.values
```



Setting up data to carry out a multiple regression model which will include both soil moisture deficit and prior area burned
```{r}
seasonFires_cm.sf <- left_join(x = seasonFires_cm.sf, y = burnedArea.sf, by = c("Year" = "forecastYear"))

seasonFires_cm.sf <- seasonFires_cm.sf %>%
  dplyr::select("Year":"percentBurned") %>%
  rename(nFIRES = nFIRES.x)
```

Multiple negative binomial regression model. Soil moisture deficit and prior acres burned.
Use temp data frame because 1992 does not store all acres burned data.
```{r}
library(MASS)

temp_cm <- seasonFires_cm.sf %>%
  filter(Year >= 1993)
nbr2 <- glm.nb(nFIRES ~  Qlcm + acresBurned, data = temp_cm)

detach(package:MASS, unload=TRUE)
```

add nbr2 fitted values to dataframe to compare models
```{r}
temp_cm$nbr2_fit <- nbr2$fitted.values
```

explore the change in fitted values
```{r}
temp_cm$change_fit <- temp_cm$nbr1_fit - temp_cm$nbr2_fit
```


```{r}
extract_eq(nbr2, use_coefs = TRUE, coef_digits = 8, fix_signs = FALSE)
```

8 digits are used to show extent of acres burned beta coefficent
$$
\log ({ \widehat{E( \operatorname{nFIRES} )} })  = 0.75384833 + 0.18424609(\operatorname{Qlcm}) + -1.026e-05(\operatorname{acresBurned})
$$

Acres burned is not statistically significant to the model
```{r}
summary(nbr2)
```

Exploring nrb2 predicted values and confidence intervals
```{r}
temp.df <- data.frame(Qlcm = 0, acresBurned = 518.07)
predict.glm(nbr2, newdata = temp.df, type = "response", se.fit = TRUE)
```



##editing Thesis plots and figures
Create a filtered data frame from 1943 to 2020 to omit missing data
```{r}
TLH_1943.df <- TLH.df %>%
  filter(Year >= 1943)
```

Exploring paper plots in centimeters
```{r}
ggplot(TLH_1943.df, aes(x = Date, y = Qlcm, color = Qlcm)) +
  geom_line() +
#  scale_color_gradientn(colors = terrain.colors(5), guide = 'none') +
  scale_color_distiller(palette = "BrBG",
                        direction = -1,
                        guide = 'none') +
  scale_x_date(date_breaks = "20 years", date_labels = "%Y") +
  ylab("Soil moisture deficit (cm)") + xlab("") +
  geom_smooth(method = lm, se = FALSE, color = "white") +
  facet_wrap(~ month, ncol = 12) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 70, hjust = 1)) 
#  labs(title = "Amount of water needed to bring soil moisture to full capacity in the Apalachicola National Forest",
#       subtitle = "Daily soil moisture deficit (mm)", 
#       caption = "Period of record: 1943-2020, Data source: NWSFO Tallahassee") 
```

Number of days per year with KBDI above 650 (hundredths of an inch) (Qlcm = 16.51)
```{r}
library(MASS)

TLH_1943.df %>%
  group_by(Year) %>%
  summarize(nD = sum(Qlcm > 16)) %>%
  ggplot(aes(Year, nD)) +
  geom_point() +
  labs(x = "", y = "Number of days\ndryness exceeded 16 cm") +
#  labs(title = "Number of days KBDI exceeds 160 mm",
#       x = "", y = "") +
#  stat_smooth(method = "glm.nb",
#              formula = y ~ x, 
#              se = FALSE,
#              col = "red") +
  labs(title = "Prior to 1980, the annual number of days\nexceeding 16 cm KBDI was rarely more than 25",
  subtitle = "Since then, nearly half of the years have had at least 25 days above 16cm KBDI", 
  caption = "Period of record: 1943-2020, Data source: NWSFO Tallahassee") +
  theme_minimal()
detach(package:MASS, unload=TRUE)
```


```{r}
library(ggrepel)

ggplot(data = seasonFires_cm.sf,
       mapping = aes(x = Qlcm, y = nFIRES, label = Year)) +
  geom_point() +
  geom_label_repel() +
  labs(x = "Soil moisture deficit (cm) on April 30th", y = "Number of fires in the ANF (May-Jul)") +
  theme_minimal()
```


##exploring potential clustering
Are the locations of natural-caused wildfires during May, June, and July within the ANF clustered?

Create a {spatstat} object containing the locations of lightning-kindled fires in the ANF during May-July. Get a planar projection suggested from `suggest_crs()`. Here projected 2779 NAD83(HARN) / Florida North
```{r}
library(spatstat)
library(maptools)
#remotes::install_github("walkerke/crsuggest")
library(crsuggest)

suggest_crs(anfbounds.sf)

W <- anfbounds.sf %>% 
  st_transform(crs = 2779) %>%
  as_Spatial() %>% 
  as.owin()

NF.ppp <- anf_LF.sf %>% 
  #filter(month(DISCOVERY_) %in% c(5, 6, 7)) %>%
  #filter(month(DISCOVERY_) == 05 | month(DISCOVERY_) == 06 | month(DISCOVERY_) == 07) %>%
  st_geometry() %>%
  st_transform(crs = 2779) %>%
  as_Spatial() %>% 
  as.ppp()

NF.ppp <- NF.ppp[W] %>%
  rescale(s = 1000, 
          unitname = "km")
summary(NF.ppp)
```

Area is 2,564 square km. There are 437 natural wildfires for an average spatial intensity of 17.05 fires per 10 square kilometer during May-July over the 27-year period.
```{r}
NF.ppp %>%
  density() %>%
  plot()

NF.ppp %>%
  dirichlet() %>%
  plot()

NF.ppp %>%
  quadratcount(nx = 20, ny = 15) %>%
  plot()

NF.ppp %>%
  quadrat.test(nx = 20, ny = 15,
               method = "M",
               nsim = 9999,
               alternative = "clustered")
```

There is some weak evidence that the locations are clustered. But the evidence is not strong.

Compute Ripley K using the `Kstat()` function
```{r}
df <- NF.ppp %>%
  envelope(fun = Kest,
           nsim = 100,
           global = TRUE,
           nrank = 1) %>%
  as.data.frame()

df <- df %>%
  mutate(theo = theo * intensity(NF.ppp),
         obs = obs * intensity(NF.ppp),
         lo = lo * intensity(NF.ppp),
         hi = hi * intensity(NF.ppp))

ggplot(data = df,
       mapping = aes(x = r, y = obs)) +
  geom_ribbon(aes(ymin = lo, 
                  ymax = hi), fill = "gray80") +
  geom_line(col = "orange") +
  geom_line(aes(y = theo), color = "white") +
  xlab("Distance (km)") + ylab("K(r)") +
  labs(title = "Locations of natural wildfires in the ANF are inconsistent with a random process",
       subtitle = "Average observed (orange) versus modeled (white) number of additional wildfires within a distance r of a random wildfire.\nThe gray ribbon is the 99% uncertainy band based on 100 simulations from the spatial random process model.") +
  theme_minimal()
```

At longer lags there are features that make natural fires somewhat more or less likely (e.g., Bradwell Bay Wilderness--swamp) 

Correlation by day of month during April.
```{r}
library(xtable)

TLH.df %>%
  dplyr::filter(Year >= 1992 & Year <= 2018) %>%
  dplyr::filter(month == "Apr") %>%
  group_by(doy, Year) %>%
  mutate(Qlm = Qlm) %>%
  dplyr::select(Year, doy, Qlm) %>%
  filter(doy != 91 & doy != 121) %>%
  group_by(doy) %>%
  summarize(r = cor.test(Qlm, seasonFires_cm.sf$nFIRES)$estimate,
            rLo = cor.test(Qlm, seasonFires_cm.sf$nFIRES, conf.level = .9)$conf.int[1],
            rHi = cor.test(Qlm, seasonFires_cm.sf$nFIRES, conf.level = .9)$conf.int[2],
            rs = cor(Qlm, seasonFires_cm.sf$nFIRES, method = "s")) %>%
  xtable()
```

auto correlation - number of fires in each fire season
```{r}
library(tseries)

acf(seasonFires_cm.sf$nFIRES, pl = FALSE)
```


##Thesis Apendix C
Diminishing effect of the start value
```{r}
Q_value.df <- TLH.df %>%
  select("Date", "Rainfall24", "MaxTemp", "NetR")

Q0 <- 0
Q200 <- 200
Q400 <- 400
Q600 <- 600
Q800 <- 800

R <- 59.23 #Average annual rainfall for TLH in inches

source("droughtIndex.R")
Q_value.df$Q0_Ql <- droughtIndex(Q0, R, Q_value.df$MaxTemp, Q_value.df$NetR)$Ql
Q_value.df$Q0_Qlcm <- Q_value.df$Q0_Ql * .0254

Q_value.df$Q200_Ql <- droughtIndex(Q200, R, Q_value.df$MaxTemp, Q_value.df$NetR)$Ql
Q_value.df$Q200_Qlcm <- Q_value.df$Q200_Ql * .0254

Q_value.df$Q400_Ql <- droughtIndex(Q400, R, Q_value.df$MaxTemp, Q_value.df$NetR)$Ql
Q_value.df$Q400_Qlcm <- Q_value.df$Q400_Ql * .0254

Q_value.df$Q600_Ql <- droughtIndex(Q600, R, Q_value.df$MaxTemp, Q_value.df$NetR)$Ql
Q_value.df$Q600_Qlcm <- Q_value.df$Q600_Ql * .0254

Q_value.df$Q800_Ql <- droughtIndex(Q800, R, Q_value.df$MaxTemp, Q_value.df$NetR)$Ql
Q_value.df$Q800_Qlcm <- Q_value.df$Q800_Ql * .0254
```

plot first year of data
```{r}
Q_plot1 <- Q_value.df %>%
  filter(year(Date) < 1941) %>%
  filter(month(Date) < 10) %>%
  ggplot(mapping = aes(x = Date)) +
    geom_line(aes(y = Q0_Qlcm, color = "0"), size = 1.5) +
    geom_line(aes(y = Q400_Qlcm, color = "400"), size = 1.5) +
    geom_line(aes(y = Q800_Qlcm, color = "800"), size = 1.5) +
    scale_color_manual("Initial Q-values",
                       breaks = c("0", "400", "800"),
                       values = c("wheat2", "tan", "goldenrod4")) +
  ylab("Soil Moisture\nDeficit (cm)") +
  scale_x_date(breaks = as.Date(c("1940-03-01", "1940-03-15", "1940-04-01", "1940-04-15", "1940-05-01", "1940-05-15", "1940-06-01", "1940-06-15", "1940-07-01", "1940-07-15", "1940-08-01", "1940-08-15", "1940-09-01", "1940-09-15")),
               date_labels = "%b %e") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Q_plot2 <- Q_value.df %>%
  filter(year(Date) < 1941) %>%
  filter(month(Date) < 10) %>%
  ggplot(mapping = aes(x = Date)) +
    geom_line(mapping = aes(y = Rainfall24), color = "blue3", size = 1.2) +
    ylab("24Hr Net Rainfall (in)") +
    scale_x_date(breaks = as.Date(c("1940-03-01", "1940-03-15", "1940-04-01", "1940-04-15", "1940-05-01", "1940-05-15", "1940-06-01", "1940-06-15", "1940-07-01", "1940-07-15", "1940-08-01", "1940-08-15", "1940-09-01", "1940-09-15")),
               date_labels = "%b %e") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Q_plot1/Q_plot2

```