##Code for manuscript:
#"short term fire frequency manipulation alters community and population responses
#in a Southeastern pine savanna"
#Simha, A and Wright, JP (2023)
#
#Submitted: 30 June 2023
#
#Code written by A. Simha#


#######################
######read in data ####
######################
require(openxlsx)
require(data.table)
require(ggplot2)
require(dplyr)
require(lubridate)
require(tidyr)
require(vegan)
require(lme4)

##########read in cover data from 2011-2017 and add summary variables##############
data = read.xlsx("Cover_2011-2017.xlsx")

#prep data - only want late coversets
data = data[!CoverSet %like% "^E",]

##to make everything a relative cover estimate
cols = colnames(data[,6:221])
data[, pctcover := rowSums(.SD, na.rm = T), .SDcols = cols]
#data[, (cols) := lapply(.SD, function(d) d/pctcover), .SDcols = cols]

#add other summary variables to data
data[, CoverDate := as.Date(CoverDate, origin = "1899-12-30")]
data[, UID := paste0("S", Site, "P", Plot)]
data[, H := diversity(data[,6:220])]
data[, richness := rowSums(data[,6:221] > 0, na.rm = T)]


#############read in burn data from 2008-2018#######################
burns = data.table(read.xlsx(
  "Cover Data/2020-07-13_burn-data-2008-2018.xlsx"))[,1:3]
#burns = data.table(read.xlsx(
#  "Cover Data/2019-11-05_burn-data-2011-2018.xlsx"))[,1:3]
names(burns) = c("Plot", "BurnDate", "BurnYear")
burns[, BurnDate := as.Date(BurnDate, origin = "1899-12-30")]

#############merge burn and cover data, add DSB and FRI variables###

#DSB: take each cover date, identify most recent burn, and subtract to get days
##where UID and Burn year and cover year match...
#cover for year of burn should be taken after burn
##if BurnDate > CoverDate, BurnYear becomes BurnYear + 1
data1 = merge(data, burns, by.x = c("UID", "Year"), 
              by.y = c("Plot", "BurnYear"), all.x = T, all.y = F)
data1 = data1[!duplicated(data1)]
data1[, BurnYear := Year]
data1[, Plot := UID]
for (i in 1:dim(data1)[1]){
  if(!is.na(data1[i, BurnDate]) & data1[i,BurnDate] > data1[i, CoverDate]){
    data1[i, BurnYear := BurnYear + 1]
  }
}
Year_orig = data1[,Year]

#remove burn data with adjusted year column and re-merge on UID and year
burnnames = c("UID", "BurnYear", "BurnDate")
test = data1[, c("Year", "BurnYear", "CoverDate", "BurnDate")]
burns1 = data1[, burnnames, with = FALSE]
burns1[, Year_orig := Year_orig]
burns1 = burns1[!is.na(BurnDate)]
names(burns1) = c("Plot", "BurnYear", "BurnDate", "Year_orig")

##gather by UID, split and add fire return interval by getting prior burn date
data1 = merge(data, burns1, by.x = c("UID", "Year"), 
              by.y = c("Plot", "BurnYear"), all.x = T, all.y = F)

data1 = data1[order(CoverDate, decreasing = FALSE)]
setkey(data1, "UID")
data1[, LastBurnDate := shift(BurnDate, 1L, type = "lag"), by = UID]
#datan = select(data1, UID, Year, BurnDate, LastBurnDate)
#carry last observation forward to replace NAs
setkey(data1, "UID")
library(zoo)
data1[, LastBurnDate := na.locf(LastBurnDate, na.rm = F), by = UID]
#datan = select(data1, UID, Year, BurnDate, LastBurnDate)

#where is.na LastBurnDate, look at 2008-2010 burn data and pull the 
#latest date with that UID
burns = burns[BurnYear < 2011,]
burns = burns[order(BurnYear, decreasing = F)]
#table(table(burns$Plot)) #each plot is only showing up once here
for (i in 1:dim(burns)[1]){
  data1[is.na(LastBurnDate) & UID == burns[i, Plot], LastBurnDate := burns[i, BurnDate]]
}

#for NMDS analysis - number of burns in experimental period
#numburns is the number of times UID is the same and BurnDate is not NA, data since 2011
#for each row, n is a data frame with all the same UID data where year is less than current row coverdate
data1[, numburns := as.integer(NA_integer_)]
for (i in 1:dim(data1)[1]){
  n = data1[data1[,UID] == data1[i, UID],]
  n = n[BurnDate < data1[i, CoverDate],]
  if (dim(n)[1] > 0){
    data1[i, numburns := dim(n)[1]]
  } else{
    data1[i, numburns := 0]
  }
}
#add numburns value even if there wasn't a burn that year by carrying forward most recent number
setkey(data1, "CoverSet")
data1[, numburns := na.locf(numburns, na.rm = F), by = UID]

##for two burns between cover samples, current code adds 2 burns for both
#need to decrease numburns by 1 for earlier burn
data1 = data1[order(Year, decreasing = T)]
data1[duplicated(data1[, c("UID", "Year")]), numburns := numburns-1]
#x = data1[duplicated(data1[, c("UID", "Year")]), UID]
#x = data1[UID %in% x,] #check looks good

##merge on UID and Burn Year and subtract BurnDate from CoverDate
#FRI only calculated for fire years; DSB calculated for all
data1[, FRI := as.numeric(BurnDate - LastBurnDate)]
data1[,DSB := as.numeric(CoverDate-BurnDate)]
data1[is.na(DSB), DSB:=as.numeric(CoverDate-LastBurnDate)]
#test = data1[, c("UID", "Year", "CoverDate", "BurnDate", "LastBurnDate", "DSB", "FRI")]

#define this function for following scripts
isEmpty <- function(x) {
  return(length(x)==0)
}