Zach_Fires.Rmd

---
title: "Zach_Fires"
author: "Zach Law"
date: "3/15/2021"
output: html_document
editor_options: 
  chunk_output_type: console
---

This file will be used to further analyze fire data bounded by the Apalachicola National Forest.


##Importing/creating boundaries of Apalachicola National Forest

Link to National Forest Boundary File: https://data.fs.usda.gov/geodata/edw/datasets.php?xmlKeyword=National+Forest+System+Land+Units

Link to zipfile containing a shapefile for all forest boundaries: https://data.fs.usda.gov/geodata/edw/edw_resources/shp/S_USA.NFSLandUnit.zip 

*Note: This link provides boundaries for every national forest boundary in the U.S.
Because we only want the boundary for the Apalachicola National Forest, the file was edited in esri arcpro and a shapefile containing only the Apalachicola National Forest boundaries was exported for use here.

```{r}
library(tidyverse)
library(sf)

anfbounds.sf <- read_sf(dsn = "forest boundaries",
                     layer = "ANFboundaries") %>%
  st_transform(crs = 3086)

st_crs(anfbounds.sf)
  
```
create a zoomable map
```{r}
library(tmap)

tmap_mode("view")

tm_shape(anfbounds.sf) +
  tm_borders()

```


Create window for apalacicola national forest boundary

```{r}
library(spatstat)
library(maptools)

anf.win <- anfbounds.sf %>%
  as_Spatial() %>%
  as.owin()

plot(anf.win)
```

Alternative method for importing forest shapefile using only R.

```{r}
library(tidyverse)
library(sf)

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)
```

Filter shapefile for just Apalachicola National Forest. Table won't load and we don't know how to isolate the apalachicola forest boundary.


Window for all forest boundaries
```{r}
library(spatstat)
library(maptools)

nf.win <- NF_Bounds.sf %>%
  as_Spatial() %>%
  as.owin()

plot(nf.win)
```

##Importing Fire Data

```{r}
if(!"FL_Fires" %in% list.files()){
download.file("http://myweb.fsu.edu/jelsner/temp/data/FL_Fires.zip",
"FL_Fires.zip")
unzip("FL_Fires.zip")
}
FL_Fires.sf <- st_read(dsn = "FL_Fires") %>%
st_transform(crs = 3086)
```

Filter FL_Fires simple feature dataframe for only fires bounded by anf
```{r}
anf_fires.sf <- st_join(FL_Fires.sf, anfbounds.sf, join = st_within) %>%
  filter(NFSLANDU_2 == 'Apalachicola National Forest')
```


Create fire data into ppp object
```{r}
FL_fires.ppp <- FL_Fires.sf %>%
  as_Spatial() %>%
  as.ppp() %>%
  unmark()
```


Filter for fires only in the Apalachicola National Forest Boundary
```{r}
anf_fires.ppp <- FL_fires.ppp[anf.win]

anf_fires.ppp

anf_firesS.ppp <- rescale(anf_fires.ppp, s = 1000, unitname = "kilometer")

summary(anf_firesS.ppp)
plot(anf_firesS.ppp)
```
It is noted that a large amount of fires are clustered in the northeast boundary of the forest.

Note: We want to rescale the ppp object find the units.

```{r}
intensity(anf_firesS.ppp)

anf_firesS.ppp %>%
  density() %>%
  plot()
```

```{r}
library(tmap)

tmap_mode("view")

tm_shape(anfbounds.sf) +
  tm_borders()

tm_shape(anf_fires.sf) +
  tm_dots(col = "orange")

```


Look at fire data for the apalachicola national forest by month and year
```{r}
library(lubridate)
library(ggplot2)

anf_fires.sf %>% 
  ggplot(mapping = aes(x = month(DISCOVERY_))) +
    geom_bar() +
    facet_wrap(~ FIRE_Y)
``` 

View all apalachicola fires by month
Running this code we see that fires are most active during May and June and least active during December
```{r}
anf_fires.sf %>% 
  ggplot(mapping = aes(x = month(DISCOVERY_))) +
    geom_bar()
```

plot all may fires in the apalachicola boundary
```{r}
june_fires <- anf_fires.sf %>%
  filter(month(DISCOVERY_) == '6')
  
june_fires.ppp <- june_fires %>%
  as_Spatial() %>%
  as.ppp() %>%
  unmark()

june_fires.ppp <- june_fires.ppp[anf.win]

plot(june_fires.ppp)

```

plot apalachicola fires for the active fire period of June 2011

```{r}
june2011fire <- anf_fires.sf %>%
  filter(FIRE_Y == "2011") %>%
  filter(month(DISCOVERY_) == "6")

june2011fire.ppp <- june2011fire %>%
  as_Spatial() %>%
  as.ppp() %>%
  unmark()

june2011fire.ppp <- june2011fire.ppp[anf.win]

plot(june2011fire.ppp)

```



Plots based on fire size to note any trends
```{r}
anf_fires.sf %>% 
  ggplot(mapping = aes(x = FIRE_SIZE_)) +
    geom_bar() +
    facet_wrap(~ FIRE_Y)
```



All Florida fires visualized by year

Was 1996 a la nina year? Was 2010 an el nino year? 
```{r}
FL_Fires.sf %>% 
  ggplot(mapping = aes(x = month(DISCOVERY_))) +
    geom_bar() +
    facet_wrap(~ FIRE_Y)
```

##Plot Florida fires and drought for June 2011


Create window for Florida

```{r}
library(USAboundaries)

FL.sf <- us_states(states = "Florida") %>%
  st_transform(crs = st_crs(FL_Fires.sf))

FL.win <- FL.sf %>%
  as_Spatial() %>%
  as.owin()
```



All Florida Fires June 2011

https://droughtmonitor.unl.edu/Data/GISData.aspx

```{r}
june2011.sf <- FL_Fires.sf %>%
  filter(FIRE_Y == '2011') %>%
  filter(month(DISCOVERY_) == '6')

june2011.ppp <- june2011.sf %>%
  as_Spatial() %>%
  as.ppp() %>%
  unmark()

june2011.ppp <- june2011.ppp[FL.win]

plot(june2011.ppp)
```

Download drought shapefile for june
US Drought Monitor link to a list of shape files https://droughtmonitor.unl.edu/Data/GISData.aspx

Here the shape file for May 31, 2011 was used to represent drought for the period of June 2011

```{r}
if(!"USDM_20110531" %in% list.files()){
download.file("https://droughtmonitor.unl.edu/data/shapefiles_m/USDM_20110531_M.zip",
"USDM_20110531_M.zip")
unzip("USDM_20110531_M.zip")
}
drought062011.sf <- st_read(dsn = "USDM_20110531.shp") %>%
st_transform(crs = st_crs(FL_Fires.sf))
```

Subset on the FL polygons
```{r}
FL_Fires2.sf <- june2011.sf %>%
  st_intersection(FL.sf)

drought062011.sf <- drought062011.sf %>%
    st_intersection(FL.sf)
```


```{r}
library(tmap)

tm_shape(drought062011.sf) +
  tm_fill(col = "DM") +
  tm_shape(june2011.sf) +
  tm_dots()
```


plot fires within the drought shapefile
why aren't colors showing up?

```{r}
drought062011.ppp <- drought062011.sf %>%
  as_Spatial() %>%
  as.owin()

drought062011.ppp <- drought062011.ppp[FL.win]

plot(june2011.ppp[drought062011.ppp])

```
```{r}
plot(drought062011.sf)
```



###Copied Code from ANF_Fires.Rmd


## Get the required packages
```{r}
library(sf)
library(tidyverse)
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.shp")
}

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

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

Make a map showing the ANF boundary and the airport location.
```{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(anfbounds.sf))

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



Make a bar chart of frequency by cause.

Cause of Fires in the Apalachicola National Forest
```{r}
#table(Fires.sf$STAT_CAU_1)
table(anf_fires.sf$STAT_CAU_1)

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

ggplot(df,
       mapping = aes(y = reorder(STAT_CAU_1, perF),
                     x = perF,
                     fill = perF)) +
  geom_col() +
  scale_fill_distiller(palette = "Oranges",
                       direction = 1,
                       guide = FALSE) +
  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-2015",
       caption = "Data source: Short, Karen (2017)") +
  theme_minimal()
```


Lightning-sparked fires.
```{r}
LightningFires.sf <- anf_fires.sf %>%
  filter(STAT_CAU_1 == "Lightning") %>%
  mutate(Year = year(DISCOVERY_),
         Month = month(DISCOVERY_),
         Day = day(DISCOVERY_),
         YDay = yday(DISCOVERY_),
         YDayF = factor(YDay, levels = as.character(1:366)),
         SizeF = factor(FIRE_SIZE_, ordered = TRUE)) %>%
  dplyr::select(Year, Month, Day, YDay, YDayF, FIRE_SIZE, SizeF)
```

Static map of lightning-sparked fires by size of the area burned.
```{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 = LightningFires.sf,
          mapping = aes(col = SizeF), show.legend = FALSE) +
  scale_color_brewer(palette = "Oranges",
                     direction = 1,
                     guide = FALSE) +
  theme_bw() +
  labs(title = "Location of lightning-sparked wildfires in the Apalachicola National Forest (1992-2015)",
       subtitle = "Darker color points indicates the fire resulted in a larger burn area",
       caption = "Data source: Short, Karen (2017)")
```

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

ggplot(data = df,
       mapping = aes(y = MonthF, 
                     x = perF,
                     fill = perF)) +
  geom_col() +
  scale_fill_distiller(palette = "Oranges",
                       direction = 1,
                       guide = FALSE) +
  scale_x_continuous(labels = percent) +
  xlab("") + ylab("") +
  labs(title = "Over 80% of lightning-sparked wildfires in the Apalachicola National Forest occur during May-July",
       subtitle = "Percentage of all lightning-sparked wildfires by month",
       caption = "Period of record: 1992-2015, Data source: Short, Karen (2017)") +
  theme_minimal() 
```

## 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

The Keetch-Byram Drought Index assesses the risk of fire by representing the net effect of evapotranspiration and precipitation in producing cumulative moisture deficiency in deep duff and upper soil layers. The index ranges from zero, the point of no moisture deficiency, to 800, the maximum drought that is possible. Units are 100's of inches. KBDI relates current and recent weather conditions to potential or expected fire behavior. It was originally developed for the Southeast United States and is based primarily on recent rainfall patterns. It is one of the only drought indexes specifically developed to equate the effects of drought with potential fire behavior. 

For different soil types, the depth of soil required to hold 8 inches of moisture varies (loam = 30", clay = 25" and sand = 80"). A prolonged drought (high KBDI) influences fire intensity largely because more fuel is available for combustion (i.e. fuels have a lower moisture content). In addition, the drying of organic material in the soil can make it harder to suppress fires.

High values of the KBDI indicate conditions favorable for the occurrence and spread of wildfires, but drought is not by itself a prerequisite for wildfires. Other weather factors, such as wind, temperature, relative humidity and atmospheric stability, play a major role in determining the actual fire danger.

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. Assign no rainfall on days with missing values.
```{r}
TLH.df <- read.csv(file = 'Data/TLH_Daily1940.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,
         Rainfall24 = replace_na(Rainfall24, 0),
         Rainfall24mm = Rainfall24 * 25.4)

TLH.df$MaxTemp[TLH.df$Date == "2005-07-08"] <- 96 # missing max temperature. Filled from weather underground data

sum(is.na(TLH.df$PRCP))
nrow(TLH.df)
sum(is.na(TLH.df$PRCP)) / nrow(TLH.df) * 100  # only 5 missing values [< .02% of all days]
```

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}
Rainfall24 <- TLH.df$Rainfall24
PR <- dplyr::lag(Rainfall24)
PR[1] <- 0

CumR <- 0
NetR <- numeric()

for(i in 1:length(Rainfall24)) {
  R24 <- Rainfall24[i]
  if (R24 == 0) {
    NetR[i] <- 0
    CumR <- 0
  } 
  else if(R24 > 0 & R24 <= .2) {
      CumR <- CumR + R24
      if (PR[i] > .2 | CumR > .2) NetR[i] <- R24
      else if (CumR > .2) NetR[i] <- CumR - .2
      else NetR[i] <- 0
    }
  else if (R24 > .2) {
      if (CumR <= .2) {
      NetR[i] <- CumR + R24 - .2
      CumR <- CumR + R24
      }
      else {
      NetR[i] <- R24
      CumR <- CumR + R24
      }
  }
}

TLH.df$NetR <- NetR
```

Step two: Compute daily drought index. Units of temperature are degrees F.
```{r}
Q <- 269  
R <- 59.23 # average annual rainfall for TLH in inches

MaxTemp <- TLH.df$MaxTemp

Ql <- numeric()
DeltaQl <- numeric()
for(i in 1:length(Rainfall24)){
  DeltaQ <- (800 - Q) * (.968 * exp(.0486 * MaxTemp[i]) - 8.3) /(1 + 10.88 * exp(-.0441 * R)) * .001 
  Q <- ifelse(NetR[i] == 0,  Q + DeltaQ,  (Q + DeltaQ) - NetR[i] * 100)
  Q <- ifelse(Q < 0, 0, Q) 
  Ql <- c(Ql, Q)
  DeltaQl <- c(DeltaQl, DeltaQ)
}

TLH.df$Ql <- Ql
TLH.df$Qlm <- Ql * .254  # tenth of an inch to mm
TLH.df$DeltaQl <- DeltaQl
TLH.df$DroughtIndex <- floor(Ql/100)

TLH.df <- TLH.df %>%
 dplyr::filter(Year >= 1946 & Year <= 2019) # Only full years
range(TLH.df$Ql)
```

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).

Monthly average soil moisture deficit.
```{r}
TLH.df %>%
  mutate(MonthF = factor(month.name[month], 
                         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_distiller(palette = "BrBG",
#                       direction = -1,
#                       guide = FALSE) +
  scale_fill_gradientn(colors = terrain.colors(5),
                       guide = FALSE) +
#  scale_fill_gradient2(low = "#FC8B93", mid = "#965784", high = "#D9DE6E", guide = FALSE) +
  labs(x = "", y = "",
       title = "May through November is the dry season in the Apalachicola National Forest",
       subtitle = "Average soil moisture deficit (mm)",
       caption = "Period of record: 1949-2019, Data source: NWSFO Tallahassee") +
  theme_minimal() 
```



Interactive map showing the locations of all fires over the period of record.
```{r}
tm_shape(anfbounds.sf) +
  tm_borders() +
tm_shape(LightningFires.sf) +
  tm_bubbles(size = .1, col = "Month")
#tm_shape(Airport.sf) +
  #tm_dots(size = .1, col = "black")
```



lightning wildfires by year

```{r}
df.year <- LightningFires.sf %>%
  st_drop_geometry() %>%
  mutate(YearF = factor(Year, 
                        ordered = TRUE)) %>%
  group_by(YearF, .drop = FALSE) %>%
  summarize(nF = n())
```

```{r}
ggplot(data = df.year, 
       mapping = aes(x = YearF,
                     y = nF,
                     fill = nF)) +
  geom_col() +
  scale_fill_distiller(palette = "Oranges",
                       direction = 1,
                       guide = FALSE) +
  xlab("Year") + ylab("Number of Fires") +
  labs(title = "Number of Lighting-Sparked Wildfires in the Apalachicola National Forest",
       subtitle = "Could years of with high fire counts correspond to La Nina events?",
       caption = "Period of record: 1992-2015, Data source: Short, Karen (2017)") +
  theme_minimal()
  

```


KBDI for year 2000
```{r}
TLH.df %>%
  filter(year(Date) == '2000') %>%
ggplot(mapping = aes(x = Date, y = Ql)) +
  geom_line() +
  xlab("Date") + ylab("KBDI") +
  labs(title = "In 2000, wildfire occurance peaked with high values of KBDI",
       caption = "Period of record: 1949-2019, Data source: NWSFO Tallahassee")
  
```


Fires in year 2000


```{r}
df2000 <- LightningFires.sf %>%
  filter(Year == '2000') %>%
  st_drop_geometry() %>%
  mutate(MonthF = factor(month.name[Month], 
                         levels = month.name, 
                         ordered = TRUE)) %>%
  group_by(MonthF, .drop = FALSE) %>%
  summarize(nF = n())
```

```{r}

ggplot(data = df2000, 
       mapping = aes(x = MonthF,
                     y = nF,
                     fill = nF)) +
  geom_col() +
  scale_fill_distiller(palette = "Oranges",
                       direction = 1,
                       guide = FALSE) +
  xlab("Month") + ylab("Number of Fires") +
  labs(title = "The 2000 fire season peaked in July",
       subtitle = "How does this correspond to KBDI values?",
       caption = "Period of record: 2000, Data source: Short, Karen (2017)") +
  theme_minimal()
```