ANF_Fires.Rmd

---
title: "Natural Wildfires in the Apalachicola National Forest (ANF)"
output: html_document
editor_options: 
  chunk_output_type: console
---

Check to see if the wind speeds have increased during the fire season. Maybe not the average wind speed but the strongest wind speeds.

R version 4.1.2 (2021-11-01) -- "Bird Hippie"

## Get required packages
```{r}
library(sf)
library(tidyverse)
library(tidycensus)
library(lubridate)
library(ggplot2)
library(ggrepel)
library(tmap)
library(scales)
library(USAboundaries)
```

## Get a boundary file for the ANF

https://data.fs.usda.gov/geodata/edw/edw_resources/shp/S_USA.NFSLandUnit.zip
```{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",
      files = "S_USA.NFSLandUnit")
}

ANF_Boundary.sf <- st_read(dsn = "S_USA.NFSLandUnit") |>
  filter(NFSLANDU_2 == "Apalachicola National Forest")

st_area(ANF_Boundary.sf) # 2,567 sq. km (~634,000 acre)
```

Make a map showing the ANF boundary and the airport location. Export > Save as Image > tiff increasing width & height then open in Preview and convert to pdf.
```{r}
Airport.sf <- data.frame(Name = c("Regional Airport"), 
                         Latitude = c(30.39306),
                         Longitude = c(-84.35333)) |>
  st_as_sf(coords = c("Longitude", "Latitude"),
           crs = 4326) |>
  st_transform(crs = st_crs(ANF_Boundary.sf))

tmap_mode("view")
tm_shape(ANF_Boundary.sf) +
  tm_borders() +
tm_shape(Airport.sf) +
  tm_dots(size = .03, 
          col = "darkgreen")
```

## Get the wildfire data

Severity of fires: Area burned 
Frequency of fires: Number of fires

```{r}
#unzip("RDS-2013-0009.5_GPKG.zip")
#st_layers(dsn = "Data/FPA_FOD_20210617.gpkg")

t0 <- Sys.time()
Fires.sf <- st_read(dsn = "Data/FPA_FOD_20210617.gpkg",
                    layer = "Fires") |>
  filter(STATE == "FL") |>
  st_transform(crs = st_crs(ANF_Boundary.sf)) |>
  st_intersection(ANF_Boundary.sf)
Sys.time() - t0  # ~ 4 min on desktop computer, ~ 2.1 min on macbook
```

NOT USED: See also NASA burned-area data https://firms.modaps.eosdis.nasa.gov/ and https://www.mtbs.gov/
```{r, eval=FALSE}
Fires2.sf <- st_read(dsn = "NIFC_ANF_firesshapefile",
                     layer = "NIFCboundedfires")

table(year(Fires2.sf$FireDiscov))

Fires3.sf <- st_read(dsn = "mtbs_fod_pts_data",
                     layer = "mtbs_FODpoints_DD") |>
  st_intersection(ANF_Boundary.sf) |>
  mutate(Date = as.Date(Ig_Date))

tmap_mode("view")
tm_shape(ANF_Boundary.sf) +
  tm_borders() +
tm_shape(Fires3.sf) +
  tm_dots(size = .03, 
          col = "orange")
```

Make a bar chart of frequency by cause.
```{r}
library(xtable)
table(Fires.sf$NWCG_GENERAL_CAUSE) |>
  xtable()

df <- Fires.sf |>
  st_drop_geometry() |>
  group_by(NWCG_GENERAL_CAUSE) |>
  summarize(nF = n(),
            perF = nF/nrow(Fires.sf)) |>
  xtable()

ggplot(data = 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 in the Apalachicola National Forest",
       subtitle = "Based on data from 1992-2018",
       caption = "Data source: Short, Karen (2021)") +
  theme_minimal()
```

Create an interactive map showing the locations of all fires over the period of record.
```{r}
tm_shape(ANF_Boundary.sf) +
  tm_borders() +
tm_shape(Fires.sf) +
  tm_bubbles(size = "FIRE_SIZE", col = "orange") +
tm_shape(Airport.sf) +
  tm_dots(size = .1, col = "black") 
```

Create a new data frame containing only the natural-caused wildfires and one containing all fires.
```{r}
NaturalFires.sf <- Fires.sf |>
  filter(NWCG_GENERAL_CAUSE == "Natural") |>
  mutate(Year = year(DISCOVERY_DATE),
         Month = month(DISCOVERY_DATE),
         Day = day(DISCOVERY_DATE),
         YDay = yday(DISCOVERY_DATE),
         YDayF = factor(YDay, levels = as.character(1:366)),
         SizeF = factor(FIRE_SIZE_CLASS, ordered = TRUE)) |>
  dplyr::select(Year, Month, Day, YDay, YDayF, Size = FIRE_SIZE, SizeF)

AllFires.sf <- Fires.sf |>
  mutate(Year = year(DISCOVERY_DATE),
         Month = month(DISCOVERY_DATE),
         Day = day(DISCOVERY_DATE),
         YDay = yday(DISCOVERY_DATE),
         YDayF = factor(YDay, levels = as.character(1:366)),
         SizeF = factor(FIRE_SIZE_CLASS, ordered = TRUE)) |>
  dplyr::select(Year, Month, Day, YDay, YDayF, Size = FIRE_SIZE, SizeF)
```

Draw a static map showing the location of the natural-caused wildfires by size of area burned.
```{r}
FL_Counties.sf <- us_counties(states = "FL",
                              resolution = "high") %>%
  st_transform(crs = st_crs(ANF_Boundary.sf)) %>%
  st_crop(st_bbox(st_buffer(ANF_Boundary.sf, dist = .25)) )

ggplot() +
  geom_sf(data = FL_Counties.sf, fill = "transparent", col = "gray80") +
  geom_sf(data = ANF_Boundary.sf, fill = "transparent") +
  geom_sf(data = NaturalFires.sf,
          mapping = aes(col = SizeF)) +
  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)")
```

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(ANF_Boundary.sf)

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

NF.ppp <- NaturalFires.sf %>% 
#  filter(Month %in% c(5, 6, 7)) %>%
  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 356 natural wildfires for an average spatial intensity of 13.9 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) 

Monthly occurrence (relative frequency) of natural wildfires in the ANF.
```{r}
df <- NaturalFires.sf |>
  st_drop_geometry() |>
  mutate(MonthF = factor(month.abb[Month], 
                         levels = month.abb, 
                         ordered = TRUE)) |>
  group_by(MonthF, .drop = FALSE) |>
  summarize(nF = n(),
            perF = nF/nrow(NaturalFires.sf))

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

Graph showing the seasonal number of natural fires and seasonal total acres burned.
```{r}
df <- NaturalFires.sf %>%
  st_drop_geometry() %>%
  filter(Month %in% c(5, 6, 7)) %>%
  mutate(YearF = factor(Year, levels = as.character(1992:2018), ordered = TRUE)) %>%
  group_by(YearF, .drop = FALSE) %>%
  summarize(nFires = n(),
            totalAcresBurned = sum(Size)) %>%
  mutate(Year = 1992:2018)

ggplot(data = df,
       mapping = aes(x = Year, y = nFires, color = totalAcresBurned)) +
  geom_line(col = "gray80") +
  geom_point(mapping = aes(size = totalAcresBurned)) +
  scale_size_continuous(guide = "none") +
  scale_color_distiller(palette = "Oranges",
                        direction = 1,
                        guide = "none") +
  scale_x_continuous(breaks = seq(1992, 2018, by = 4)) +
  theme_minimal() +
  labs(x = "", y = "Number of wildfires")
#       title = "Large inter-seasonal variability in the number of natural wildfires in the ANF",
#       subtitle = "Larger and darker points indicate larger burn area")

ggplot(data = df,
       mapping = aes(x = Year, y = totalAcresBurned, color = nFires)) +
  geom_line(col = "gray80") +
  geom_point(mapping = aes(size = nFires)) +
  scale_size_continuous(guide = "none") +
  scale_color_distiller(palette = "Oranges",
                        direction = 1,
                        guide = "none") +
  scale_x_continuous(breaks = seq(1992, 2018, by = 3)) +
  theme_minimal() +
  labs(x = "", y = "",
       title = "Large inter-seasonal variability in acres burned from natural wildfires in the ANF",
       subtitle = "Larger and darker color points indicate more fires")
```

Annual statistics
```{r}
df$Year[df$nFires == 0]
max(df$nFires)
df$Year[df$nFires == max(df$nFires)]
mean(df$nFires); var(df$nFires)
```

Seasonal temporal autocorrelation in frequency and total area burned.
```{r}
plot(acf(df$nFires))
plot(acf(df$totalAcresBurned))
cor.test(df$nFires, df$totalAcresBurned)
cor.test(df$nFires[-1], df$totalAcresBurned[-24])
```

No significant inter-seasonal autocorrelation. Significant correlation between number of fires and burned area as expected.

Future work: Make fire size an ordered factor and create a model that predicts the probability by size category similar to what we did for predicting tornado EF rating. See `PredictTorFreqByEF.Rmd` in the `EF-dist` distribution.
```{r}
NaturalFires.sf |>
  group_by(SizeF) |>
  summarize(nFires = n(),
            Range1 = range(Size)[1],
            Range2 = range(Size)[2])
```

## Get lightning data

Daily county-level counts 1986-2013. Data location: https://www1.ncdc.noaa.gov/pub/data/swdi/reports/county/byFips/

First get the Florida FIPS codes for Liberty, Wakulla, and Franklin counties, then get the data.
```{r}
FLfips <- fips_codes %>%
  filter(state == "FL") %>%
  filter(county %in% c("Liberty County", "Wakulla County")) %>%
  pull(county_code)

fn <- paste0("https://www1.ncdc.noaa.gov/pub/data/swdi/reports/county/byFips/swdireport-12",
             FLfips, "-BETA.csv")

df <- data.frame()
for(i in 1:length(fn)){
  X <- read.csv(fn[i], na.strings = "NULL", header = TRUE, stringsAsFactors = FALSE)
  df <- rbind(df, X)
}

df <- df %>%
  mutate(SEQDAY = as.Date(SEQDAY),
         DAY = day(SEQDAY),
         MONTH = month(SEQDAY),
         YEAR = year(SEQDAY))
```

Seasonal cycle.
```{r}
sc <- df %>%
  mutate(MonthF = factor(month.abb[MONTH], 
                         levels = month.abb, 
                         ordered = TRUE)) %>%
  group_by(MonthF, .drop = FALSE) %>%
  summarize(AvgDailyCount = mean(FCOUNT_NLDN, na.rm = TRUE))

p2 <- ggplot(data = sc,
       mapping = aes(x = MonthF, 
                     y = AvgDailyCount,
                     fill = AvgDailyCount)) +
  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",
#       subtitle = "Daily average cloud-to-ground lightning strikes by month",
#       caption = "Period of record: 1986-2013, Data source: NCEI") +
  theme_minimal() 

p1 / p2
```

## Get TLH airport daily weather data and compute the Keetch & Byram drought index

Original paper outlining the rationale and how to create it: https://www.srs.fs.usda.gov/pubs/rp/rp_se038.pdf

NWS Tallahassee daily weather data:
https://www.ncdc.noaa.gov/cdo-web/datasets/GHCND/stations/GHCND:USW00093805/detail
https://www1.ncdc.noaa.gov/pub/data/ghcn/daily/readme.txt
Column explanations: https://docs.google.com/document/d/1q2WEpXndpMx9lUq-0ON63GojkaOzixrGyYEhI_xkPtw/edit?usp=sharing
AWND = Average daily wind speed (tenths of meters per second) 

Import the summary of the day data and add columns to the data frame.
```{r}
TLH.df <- read.csv(file = 'Data/TLH_Daily1940-2020.csv',
                   stringsAsFactors = FALSE,
                   header = TRUE) |>
  mutate(Date = as.Date(DATE)) |>
  mutate(Year = year(Date), 
         month = month(Date, label = TRUE, abbr = TRUE),
         doy = yday(Date),
         MaxTemp = TMAX,
         MinTemp = TMIN,
         Rainfall24 = PRCP,
         Rainfall24mm = Rainfall24 * 25.4) |>
  filter(Year <= 2020)

sum(is.na(TLH.df$PRCP))
sum(is.na(TLH.df$PRCP)) / nrow(TLH.df) * 100  # only 5 missing values [~ .017% of all days]
sum(is.na(TLH.df$TMAX))
sum(is.na(TLH.df$TMAX)) / nrow(TLH.df) * 100  # only 1 missing value [~ .0034% of all days]

# Missing temperature and rainfall filled in from NCDC daily summaries
# https://gis.ncdc.noaa.gov/maps/ncei/summaries/daily

TLH.df$MaxTemp[TLH.df$Date == "2005-07-08"] <- 95
TLH.df$Rainfall24[is.na(TLH.df$PRCP)] <- c(.1, .3, 0, .1, 0)
TLH.df$Rainfall24mm[is.na(TLH.df$PRCP)] <- c(3, 8.9, 0, 2.5, 0)

AnnualRainfall <- TLH.df |>
  filter(Year >= 1943) |>  # complete data started in March 1943
  group_by(Year) |>
  summarize(totalRainfall = sum(Rainfall24),
            totalRainfallcm = totalRainfall * 2.54)

AnnualRainfall |> 
  summarize(avgRainfall = mean(totalRainfallcm),
            sdRainfall = sd(totalRainfallcm),
            varRainfall = var(totalRainfallcm))
AnnualRainfall |>
  arrange(totalRainfallcm)
```

Plot average daily mean wind speeds during May-July by year.
```{r}
TLH.df |>
  filter(Year >= 1984) |>
  filter(month %in% c("May", "Jun", "Jul")) |>
  group_by(Year) |>
  summarize(AvgWindSpeed = mean(AWND) / 10) |>
ggplot(mapping = aes(x = Year, y = AvgWindSpeed)) +
  geom_point() +
  geom_smooth(method = lm)
```

Plot average daily minimum during May-July by year.
```{r}
TLH.df |>
  filter(month %in% c("May", "Jun", "Jul")) |>
  group_by(Year) |>
  summarize(AvgMinTemp = mean(MinTemp, na.rm = TRUE)) |>
ggplot(mapping = aes(x = Year, y = AvgMinTemp)) +
  geom_point() +
  geom_smooth(method = lm) +
  labs(y = "Average minimum temperature (F)") +
  theme_minimal()
```

Plot monthly values of average daily high temperature and average month rainfall.
```{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
```

Compute daily values of the KBDI. Original paper outlining the rationale and how to create it: https://www.srs.fs.usda.gov/pubs/rp/rp_se038.pdf with a minor correction identified in Alexander1992.pdf (paper found in Dropbox/Literature).

Step one: Compute net rainfall on each day. Units are inches.
```{r}
source("netRainfall.R")
TLH.df$NetR <- netRainfall(Rainfall24 = TLH.df$Rainfall24)
```

Step two: Compute daily drought index. Units of temperature are degrees F. Q is set initially to 400. The influence of this start value diminishes to zero after a few months R is the annual average rainfall for TLH in inches
```{r}
source("droughtIndex.R")
DI <- droughtIndex(Q = 400, R = 60, MaxTemp = TLH.df$MaxTemp, NetR = TLH.df$NetR)
TLH.df$Ql <- DI$Ql
TLH.df$Qlm <- TLH.df$Ql * .254  # tenth of an inch to mm
TLH.df$DroughtIndex <- DI$DroughtIndex
```

Note: There is a package to compute the standardized precip-evapotranspiration index {SPEI}. Only good for monthly data. Same for the Palmer Drought Severity Index (PDSI).

Daily
```{r}
TLH.df <- TLH.df |>
  filter(Year >= 1943) |>
  mutate(Date2 = as.Date(doy - 1, origin = "2020-01-01"))

p1 <- ggplot(data = TLH.df, 
       mapping = aes(x = Date2, y = Year, fill = Qlm)) +
  geom_tile() +
  scale_y_reverse(breaks = seq(1945, 2020, by = 10)) +
  scale_x_date(date_breaks = "month", date_labels = "%b", position = "top") +
  scale_fill_distiller(palette = "BrBG",
                       direction = -1,
                       name = "mm") +
  labs(title = "A", x = "", y = "") +
  theme_minimal() 
```

Monthly average soil moisture deficit.
```{r}
p2 <- TLH.df |>
  mutate(MonthF = factor(month.abb[month], 
                         levels = month.abb, 
                         ordered = TRUE)) |>
  group_by(MonthF, .drop = FALSE) |>
  summarize(AvgSoilMoistureDeficit = mean(Qlm)) |>
ggplot(mapping = aes(x = MonthF, 
                     y = AvgSoilMoistureDeficit,
                     fill = AvgSoilMoistureDeficit)) +
  geom_col() +
  scale_fill_distiller(palette = "BrBG",
                       direction = -1,
                       guide = 'none') +
  labs(title = "B", y = "Soil Moisture Deficit (mm)", x = "") +
  theme_minimal() 

p1 / p2
```

Make some trend plots. This plot is statistically more realistic because it removes the temporal autocorrelations.
```{r}
ggplot(TLH.df, aes(x = Date, y = Qlm, color = Qlm)) +
  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 (mm)") + xlab("") +
  geom_smooth(method = lm, se = FALSE, color = "white") +
  facet_wrap(~ month, ncol = 12) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, 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") 
```

Repeat using data from Apalachicola (see `compare.rmd`). https://www.ncdc.noaa.gov/  > I want to search for data at a particular location > DATA TOOLS > Final a Station

Linear trend model. Statistical significance of the trend by month.
```{r}
summary(lm(Qlm ~ Date, data = TLH.df)) # this model assumes daily Qlm are independent (they are certainly not)

library(broom)

( Trends.df <- TLH.df %>% 
  group_by(Year, month) %>%
  summarize(AvgQlm = mean(Qlm)) %>%
  group_by(month) %>%
  do(tidy(lm(AvgQlm ~ Year, data = .))) %>%
  dplyr::filter(term == "Year") )

Avg <- TLH.df %>%
  group_by(month) %>%
  summarize(AvgQlm = mean(Qlm)) %>%
  pull(AvgQlm)

cor(Trends.df$estimate, Avg) # positive, indicating upward trends are occurring in months with higher drought index values

Avg <- TLH.df %>%
  group_by(month) %>%
  summarize(nYear = n_distinct(Year), 
            TotalPRCP = sum(PRCP, na.rm = TRUE),
            AvgMonthlyPRCP = TotalPRCP/nYear) %>%
  pull(AvgMonthlyPRCP)

cor(Trends.df$estimate, Avg)  # negative as expected since upward trends are occurring in months with less rainfall
```

Soil moisture deficits are increasing at .3 mm/year on average during April. Which amounts to 3 mm/decade.

What is the correlation between PRCP, TMAX and Qlm grouped by month?
```{r}
TLH.df %>%
  group_by(month) %>%
  summarize(rPRCP = cor(Qlm, PRCP, use = "complete"),
            rTMAX = cor(Qlm, TMAX, use = "complete"),
            rTMIN = cor(Qlm, TMIN, use = "complete"),
            rTAVG = cor(Qlm, TAVG, use = "complete"))
```

Seasonal (fire season) temperature vs KDBI scatter plot. Create a character variable called `FireSeason` (May, June, July). Then group by that variable computing the average temp and total precipitation (or avg KDBI).
```{r}
df <- TLH.df %>%
  mutate(FireSeason = month %in% c("May", "Jun", "Jul")) %>%
  filter(FireSeason) %>%
  group_by(Year) %>%
  summarize(AvgHiTemp = mean(MaxTemp),
            AvgLoTemp = mean(MinTemp),
            TotalPrcp = sum(Rainfall24mm),
            TotalNetR = sum(NetR),
            AvgQlm = mean(Qlm))

ggplot(df, 
       mapping = aes(y = TotalNetR, x = AvgHiTemp, label = Year, color = factor(Year))) +
  geom_point() +
  geom_label_repel(size = 2.5) +
  scale_y_continuous(limits = c(0, NA)) +
  scale_color_grey(start = .8, end = .2, guide = "none") +
#  geom_smooth(method = lm, se = FALSE, col = "black") +
  ylab("Total Precipitation (mm)") + xlab("Average Daily High Temperature (°F)") +
  ggtitle(label = "Fire seasons in the Apalachicola National Forest (May-Jul) that are drier tend to be hotter",
          subtitle = "Dry days are sunnier and thus they tend to get hotter. Later years have darker labels.") +
  theme_minimal() 
```

Get SOI data and join with TLH.df. LINK NO LONGER WORKS
```{r}
SOI.df <- read.csv(file = "https://www.ncdc.noaa.gov/teleconnections/enso/indicators/soi/data.csv",
                   skip = 1, header = TRUE) %>%
  mutate(Date = parse_date_time(as.character(Date), "ym"),
         Year = year(Date),
         Month = month(Date),
         month = month(Date, label = TRUE, abbr = TRUE),
         SOI = Value) %>%
  dplyr::select(Year, Month, month, SOI)

Joined.df <- TLH.df %>% 
  dplyr::filter(Year >= 1951) %>%
  group_by(Year, month) %>%
  summarize(AvgQlm = mean(Qlm)) %>%
  dplyr::select(Year, month, AvgQlm) %>%
  left_join(SOI.df, by = c("Year", "month"))

Joined.df %>%
  group_by(month) %>%
  do(tidy(lm(AvgQlm ~ Year + SOI, data = .))) %>%
  dplyr::filter(term == "Year")

Joined.df %>%
  group_by(month) %>%
  do(tidy(lm(AvgQlm ~ Year + SOI, data = .))) %>%
  dplyr::filter(term == "SOI")
```

SOI is positively correlated with moisture deficit. La Nina conditions (positive SOI) during winter into spring leads to larger moisture deficits (higher Qlm). The upward trend in Qlm remains significant after accounting for SOI.

In terms of frequency by Drought Index category. Proportion of the year in each drought category.
```{r}
TLH.df %>%
  group_by(Year, DroughtIndex) %>%
  summarize(N = n()) %>%
  ggplot(aes(x = Year, y = N)) +
  geom_col(aes(fill = as.factor(DroughtIndex))) +
  scale_fill_manual(values = terrain.colors(8), name = "KDBI") +
  scale_x_continuous(name = "", breaks = seq(1941, 2020, 4)) +
  scale_y_reverse() +
  theme_minimal() +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank()) +
  guides(fill = guide_legend(reverse = TRUE)) +
  labs(title = "The proportion of the year in drought (KDBI > 5) in the ANF has increased since 1960",
       subtitle = "KDBI = 7-6: Drought, 4-5: Dry, 1-3: Moist, 0: Saturated",
       caption = "Period of record: 1941-2020, Data source: NWSFO Tallahassee") 
```

Daily
```{r}
ggplot(data = TLH.df, 
       mapping = aes(x = Date2, y = Year, fill = factor(DroughtIndex))) +
  geom_tile() +
  scale_y_reverse(breaks = seq(1943, 2020, 4)) +
  scale_x_date(date_breaks = "month", date_labels = "%b", position = "top") +
#  scale_fill_ordinal(name = "", direction = 1, alpha = .6) +
   scale_fill_brewer(palette = "BrBG",
                       direction = -1,
                       name = "") +
  ylab("") + xlab("") +
  theme_minimal() +
  guides(fill = guide_legend(reverse = TRUE)) +
  labs(title = "Daily KBDI category in the Apalachicola National Forest",
       subtitle = "KDBI = 7-6: Drought, 4-5: Dry, 1-3: Moist, 0: Saturated",
       caption = "Period of record: 1949-2019, Data source: NWSFO Tallahassee")
```

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

TLH.df |>
  group_by(Year) |>
  summarize(nD = sum(Qlm > 160)) |>
  ggplot(aes(Year, nD)) +
  geom_point() +
  labs(x = "", y = "Number of days\nSoil moisture deficit exceeded 160 mm") +
#  labs(title = "Number of days KBDI exceeds 160 mm",
#       x = "", y = "") +
#  stat_smooth(method = "glm.nb",
#              formula = y ~ x, 
#              se = FALSE,
#              col = "red") +
  theme_minimal()
```

Combine seasonal fire data with dryness data. Include a column with previous year's total acres burned.
```{r}
SeasonalFires <- NaturalFires.sf %>%
  st_drop_geometry() %>%
  filter(Month %in% c(5, 6, 7)) %>%
  mutate(YearF = factor(Year, levels = as.character(1992:2018), ordered = TRUE)) %>%
  group_by(YearF, .drop = FALSE) %>%
  summarize(nFires = n(),
            totalAcresBurned = sum(Size)) %>%
  mutate(Year = 1992:2018)

SeasonalFires2 <- AllFires.sf %>%
  st_drop_geometry() %>%
  filter(Month %in% c(5, 6, 7)) %>%
  mutate(YearF = factor(Year, levels = as.character(1992:2018), ordered = TRUE)) %>%
  group_by(YearF, .drop = FALSE) %>%
  summarize(nFires = n(),
            totalAcresBurned = sum(Size)) %>%
  mutate(Year = 1992:2018)

Dryness <- TLH.df %>%
  dplyr::filter(Year >= 1992 & Year <= 2018) %>%
  dplyr::filter(month == "Apr") %>%
  group_by(Year) %>%
  summarize(QlmLastDay = last(Qlm),
            AvgQlm = mean(Qlm))

#AprilSOI <- SOI.df %>%
#  dplyr::filter(Year >= 1992 & Year <= 2018) %>%
#  dplyr::filter(month == "Apr") %>%
#  pull(SOI)

PreviousYearBurnArea.df <- Fires.sf |>
  st_drop_geometry() |>
  mutate(Date = DISCOVERY_DATE,
         Year = year(DISCOVERY_DATE),
         Month = month(DISCOVERY_DATE),
         BurnYear = Year + (Month >= 5))|>
  group_by(BurnYear) |>
  mutate(BurnArea = cumsum(FIRE_SIZE)) |>
  summarize(BurnArea = last(BurnArea)) |>
  filter(BurnYear < 2019)

df <- data.frame(Year = Dryness$Year, 
                 QlmLastDay = Dryness$QlmLastDay,
                 AvgQlm = Dryness$AvgQlm,
                 nFires = SeasonalFires$nFires, 
                 AcresBurned = SeasonalFires$totalAcresBurned,
                 nFiresAll = SeasonalFires2$nFires,
                 AcresBurnedAll = SeasonalFires2$totalAcresBurned,
                 preYearAcresBurned = PreviousYearBurnArea.df$BurnArea)
#                 AprilSOI = AprilSOI)

write_csv(df, file = "SeasonalForecastModelData.csv")
```

Correlations and scatter plots.
```{r}
cor(df$AvgQlm, df$AcresBurned) # .54
cor(df$QlmLastDay, df$AcresBurned) #.48
cor(df$QlmLastDay, df$nFires) # .61
cor(df$AvgQlm, df$nFires) # .56

ggplot(data = df,
       mapping = aes(x = QlmLastDay, y = nFires, label = Year)) +
  geom_point() +
  geom_label_repel() +
  labs(x = "Soil moisture deficit (mm) on April 30th", y = "Number of fires in the ANF (May-Jul)") +
  theme_minimal()
```

Correlation by day of month during April.
```{r}
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, SeasonalFires$nFires)$estimate,
            rLo = cor.test(Qlm, SeasonalFires$nFires, conf.level = .9)$conf.int[1],
            rHi = cor.test(Qlm, SeasonalFires$nFires, conf.level = .9)$conf.int[2],
            rs = cor(Qlm, SeasonalFires$nFires, method = "s")) |>
  xtable()
```

Make a plot
```{r}
  ggplot(data = df, 
         mapping = aes(x = Date, y = r)) +
  geom_ribbon(aes(ymin = rLo, ymax = rHi), fill = "gray90") +
  geom_point() +
  geom_line() +
  scale_y_continuous(limits = c(0, 1)) +
#  scale_x_continuous(breaks = c(1,	32,	60,	91, 121, 152, 182, 213, 244, 274, 305, 335),
#                     labels = month.abb, position = "top") +
  theme_minimal()
```

A slight upward trend in the correlation from early to late April but no significant difference.

Continue with `PredictionModel.Rmd` and `GriddedClimateData.Rmd`.