#####
#Load packages

library(dplyr)
library(ggplot2)
library(readxl)
library(stringr)
library(lubridate)
library(tidyr)
library(ggalluvial)
library(broom)
#####
#Plot themes
theme_scatter <- function() {
  theme_classic(base_size = 12, base_family = "serif")%+replace%
    theme(axis.line = element_line(size = 1.1),
          axis.ticks = element_line(size = 1.1),
          axis.text = element_text(color = "black"),
          legend.position = c(0.80,0.80))
}

theme_afs_point <- function(base_size = 14, base_family = "Times"){
  theme_classic(base_size = base_size, base_family = base_family) %+replace%
    theme(
      axis.title.y = element_text(size = 14, angle = 90, margin = margin(t=0, r=10, b=0.1, l=0), colour = "black"),
      axis.title.x = element_text(size = 14, margin = margin(t=10, r=0, b=0, l=0),colour = "black"),
      plot.title = element_text(face = "bold", family = "Arial"),
      legend.position = c(0.3,0.85),
      legend.text = element_text(size = 12),
      plot.margin = unit(c(0.5,0.5,0.5,0.5), "cm"),
      axis.ticks.y = element_line(size = 0.5),
      axis.ticks.x = element_line(size = 0.5),
      axis.ticks.length = unit(0.2, "cm"),
      axis.text.y = element_text(colour = "black", size = 14, angle = 0, vjust = 0.5, hjust = 1, margin = margin(t=0, r=5,b=0,l=0)),
      axis.text.x = element_text(colour = "black", size = 14, angle = 0, vjust = 0, hjust = 0.5, margin = margin(t=5,r=0,b=0,l=0)),
      axis.line = element_line(colour = "black", size = 0.5, lineend = "square")
    )
}


#####
#Habitat data at 30m resolution in each pool per year for each habitat type (i.e., main channel, Off main channel, dyke field)
hab_data <- read.csv(file = "Desktop/Arkansas_River_Habitat/Data/JRC_yearly_ArkChan_Pool_Hab2.csv", header = TRUE, sep = ",") #Must change directory to read in the data
names(hab_data)
colnames(hab_data)<-c("Year", "Pool", "Hab_Pool", "Habitat", "No_Observ", "Not_Water", "Permanent", "Seasonal")
hab_data<-hab_data %>% filter(!Habitat == "Distant")
hab_data$Pool<-as.factor(hab_data$Pool)

#convert pixel counts to area - each pixel is 30 by 30 m
hab_data$Permanent <- ((hab_data$Permanent)*(900))*(1/10000)
hab_data$Seasonal <- ((hab_data$Seasonal)*(900))*(1/10000)


#Plots of raw data for each habitat type
permanent_raw<-ggplot(data = hab_data, aes(y = Permanent, x = Year, alluvium = Habitat)) +
  geom_alluvium(aes(fill = Habitat, colour = Habitat))+
  facet_wrap(~Pool, scales = "free_y", ncol = 3)+
  coord_cartesian(expand = FALSE)+
  labs(x = "Year", y = "Permanent Water (ha)")+
  theme_bw()+
  theme(panel.grid = element_blank(),
        axis.text.x = element_text(angle = 90),
        axis.text.y = element_text(angle = 90),
        legend.position = c(0.80,0.07))

seasonal_raw<-ggplot(data = hab_data, aes(y = Seasonal, x = Year, alluvium = Habitat)) +
  geom_alluvium(aes(fill = Habitat, colour = Habitat))+
  facet_wrap(~Pool, scales = "free_y", ncol = 3)+
  coord_cartesian(expand = FALSE)+
  labs(x = "Year", y = "Episodic Water (ha)")+
  theme_bw()+
  theme(panel.grid = element_blank(),
        axis.text.x = element_text(angle = 90),
        axis.text.y = element_text(angle = 90),
        legend.position = c(0.80,0.07))

ggsave(permanent_raw, file = "permanent_raw.png")
ggsave(seasonal_raw, file = "seasonal_raw.png")

#Stage data for each pool
Stage_Data <- read.csv(file = "Desktop/Arkansas_River_Habitat/Data/Hydro_Data_AR.csv", header = TRUE, sep = ",")#Must change directory to read in data
Stage_Data$Date <- as.character(Stage_Data$Date)
Stage_Data$Date <- dmy(Stage_Data$Date)
Stage_Data$Year <- year(Stage_Data$Date)
str(Stage_Data)
colnames(Stage_Data[,3:13])<-c("2","3","4","5","6","7","8","9","10","12","13")

Stage_Sum <- Stage_Data %>% 
  gather(key = "Pool", value = "Stage", -Obs, -Date, -Year) %>% 
  separate(Pool,into = c(NA, "Pool"), sep = "(?<=[A-Za-z])(?=[0-9])") %>% 
  filter(Year > 1983 & Year < 2016)

Stage_Year_Sum <- Stage_Data %>% 
  gather(key = "Pool", value = "Stage", -Obs, -Date, -Year) %>% 
  separate(Pool,into = c(NA, "Pool"), sep = "(?<=[A-Za-z])(?=[0-9])") %>% 
  group_by(Year,Pool) %>% 
  mutate(Stage_2 = Stage/3.281) %>%
  summarize(Mean_Stage = mean(Stage_2,na.rm = TRUE),
            SD_Stage = sd(Stage_2, na.rm = TRUE),
            Test_stage = mean(Stage,na.rm = TRUE),
            test_sd_stage = sd(Stage, na.rm = TRUE),
            Min_Stage = min(Stage_2,na.rm = TRUE),
            Max_Stage = max(Stage_2,na.rm = TRUE),
            CV_Stage = (sd(Stage_2, na.rm = TRUE)/mean(Stage_2, na.rm = TRUE))*100) %>% 
  filter(Year > 1983 & Year < 2016)


#Stage data are in ft; convert to m
Stage_Year_Sum$Mean_Stage_2<-Stage_Year_Sum$Mean_Stage/3.281

stage_plot <- ggplot()+
geom_line(data = Stage_Year_Sum, aes(x=Year,y=Mean_Stage_2))+
  facet_wrap(~factor(Stage_Year_Sum$Pool, levels = c("2","3","4","5","6","7","8","9","10","12","13")), scales = "free_y", ncol = 3)+
  labs(y="Mean Annual River Stage (m)")+
  theme_bw()+
  theme(panel.grid = element_blank(),
        axis.text.y = element_text(colour = "black", size = 10, angle = 45, vjust = 0.5, hjust = 1, margin = margin(t=0, r=5,b=0,l=0)),
        axis.text.x = element_text(colour = "black", size = 10, angle = 0, vjust = 0, hjust = 0.5, margin = margin(t=5,r=0,b=0,l=0)))

ggsave(stage_plot, file = "stage_plot.png")

stage_cv_plot<-ggplot()+
  geom_line(data = Stage_Year_Sum, aes(x=Year,y=CV_Stage))+
  facet_wrap(~factor(Stage_Year_Sum$Pool, levels = c("2","3","4","5","6","7","8","9","10","12","13")), scales = "free_y")+
  labs(y = "River Stage Coefficient of Variation")+
  theme_bw()+
  theme(panel.grid = element_blank(),
        axis.text.y = element_text(colour = "black", size = 10, angle = 45, vjust = 0.5, hjust = 1, margin = margin(t=0, r=5,b=0,l=0)),
        axis.text.x = element_text(colour = "black", size = 10, angle = 0, vjust = 0, hjust = 0.5, margin = margin(t=5,r=0,b=0,l=0)))

ggsave(stage_cv_plot, file = "stage_cv_plot.png")

hab_data_2<-inner_join(hab_data, Stage_Year_Sum, by = c("Year","Pool"))

#Permanent habitat models
hab_perm_30m_mod_1 <- glm(Permanent ~ Year, data = hab_data_2, family = "Gamma"(link=log))
hab_perm_30m_mod_2 <- glm(Permanent ~ Pool, data = hab_data_2, family = "Gamma"(link=log))
hab_perm_30m_mod_3 <- glm(Permanent ~ Habitat, data = hab_data_2, family = "Gamma"(link=log))
hab_perm_30m_mod_4 <- glm(Permanent ~ Year + Pool + Habitat, data = hab_data_2, family = "Gamma"(link=log))
hab_perm_30m_mod_5 <- glm(Permanent ~ Year + Pool*Habitat, data = hab_data_2, family = "Gamma"(link=log))
hab_perm_30m_mod_6 <- glm(Permanent ~ Habitat + Year*Pool, data = hab_data_2, family = "Gamma"(link=log))
hab_perm_30m_mod_7 <- glm(Permanent ~ Pool + Year*Habitat, data = hab_data_2, family = "Gamma"(link=log))
hab_perm_30m_mod_8 <- glm(Permanent ~ Pool*Year*Habitat, data = hab_data_2, family = "Gamma"(link=log))
hab_perm_30m_mod_9 <- glm(Permanent ~ Pool + Mean_Stage*Habitat, data = hab_data_2, family = "Gamma"(link=log))
hab_perm_30m_mod_10 <- glm(Permanent ~ Pool*Mean_Stage*Habitat, data = hab_data_2, family = "Gamma"(link=log))
hab_perm_30m_mod_11 <- glm(Permanent ~ Mean_Stage, data = hab_data_2, family = "Gamma"(link=log))
hab_perm_30m_mod_12 <- glm(Permanent ~ Pool*Mean_Stage, data = hab_data_2, family = "Gamma"(link=log))
hab_perm_30m_mod_13 <- glm(Permanent ~ Pool*CV_Stage, data = hab_data_2, family = "Gamma"(link=log))
hab_perm_30m_mod_14 <- glm(Permanent ~ 1, data = hab_data_2, family = "Gamma"(link=log))
hab_perm_30m_mod_15 <- glm(Permanent ~ Pool*CV_Stage*Habitat, data = hab_data_2, family = "Gamma"(link=log))


#Seasonal habitat models
hab_seas_30m_mod_1 <- glm(Seasonal ~ Year, data = hab_data_2, family = "Gamma"(link=log))
hab_seas_30m_mod_2 <- glm(Seasonal ~ Pool, data = hab_data_2, family = "Gamma"(link=log))
hab_seas_30m_mod_3 <- glm(Seasonal ~ Habitat, data = hab_data_2, family = "Gamma"(link=log))
hab_seas_30m_mod_4 <- glm(Seasonal ~ Year + Pool + Habitat, data = hab_data_2, family = "Gamma"(link=log))
hab_seas_30m_mod_5 <- glm(Seasonal ~ Year + Pool*Habitat, data = hab_data_2, family = "Gamma"(link=log))
hab_seas_30m_mod_6 <- glm(Seasonal ~ Habitat + Year*Pool, data = hab_data_2, family = "Gamma"(link=log))
hab_seas_30m_mod_7 <- glm(Seasonal ~ Pool + Year*Habitat, data = hab_data_2, family = "Gamma"(link=log))
hab_seas_30m_mod_8 <- glm(Seasonal ~ Pool*Year*Habitat, data = hab_data_2, family = "Gamma"(link=log))
hab_seas_30m_mod_9 <- glm(Seasonal ~ Pool + Mean_Stage*Habitat, data = hab_data_2, family = "Gamma"(link=log))
hab_seas_30m_mod_10 <- glm(Seasonal ~ Pool*Mean_Stage*Habitat, data = hab_data_2, family = "Gamma"(link=log))
hab_seas_30m_mod_11 <- glm(Seasonal ~ Mean_Stage, data = hab_data_2, family = "Gamma"(link=log))
hab_seas_30m_mod_12 <- glm(Seasonal ~ Pool*Mean_Stage, data = hab_data_2, family = "Gamma"(link=log))
hab_seas_30m_mod_13 <- glm(Seasonal ~ Pool*CV_Stage, data = hab_data_2, family = "Gamma"(link=log))
hab_seas_30m_mod_14 <- glm(Seasonal ~ 1, data = hab_data_2, family = "Gamma"(link=log))
hab_seas_30m_mod_15 <- glm(Seasonal ~ Pool*CV_Stage*Habitat, data = hab_data_2, family = "Gamma"(link=log))

write.csv(tidy(hab_perm_30m_mod_8), "coefs.perm.csv")
write.csv(tidy(hab_seas_30m_mod_8), "coefs.seas.csv")

perm_models<-list(hab_perm_30m_mod_1,hab_perm_30m_mod_2,hab_perm_30m_mod_3,hab_perm_30m_mod_4,hab_perm_30m_mod_5,hab_perm_30m_mod_6,hab_perm_30m_mod_7,hab_perm_30m_mod_8,hab_perm_30m_mod_9,hab_perm_30m_mod_10,hab_perm_30m_mod_11,hab_perm_30m_mod_12,hab_perm_30m_mod_13,hab_perm_30m_mod_14,hab_perm_30m_mod_15)
seas_models<-list(hab_seas_30m_mod_1,hab_seas_30m_mod_2,hab_seas_30m_mod_3,hab_seas_30m_mod_4,hab_seas_30m_mod_5,hab_seas_30m_mod_6,hab_seas_30m_mod_7,hab_seas_30m_mod_8,hab_seas_30m_mod_9,hab_seas_30m_mod_10,hab_seas_30m_mod_11,hab_seas_30m_mod_12,hab_seas_30m_mod_13,hab_seas_30m_mod_14,hab_seas_30m_mod_15)

#create a aic function using the following steps
aic = sapply(perm_models, AIC)#for seasonal models; switch to perm_models for permanent models
which.min(aic)
models[which.min(aic)]
#Obtaining the Del AIC values which are the distance from the best model to each model in the model set#
deltas = aic - min(aic)
weights = exp(-deltas/2)
weights = weights/sum(weights)
#Obtaining the number of parameters and combining the model selection results into one table
npars = sapply(perm_models,function(mm)length(coef(mm)))+1#accounts for the addition of a variance parameter
cbind(AIC=aic,k=npars,deltas=deltas,weights=weights)
aic.table = (cbind(perm_models=1:15, AICc = aic, k = npars, deltas = round(deltas,2), weights = round(weights, 3)))
AICctable <- aic.table[order(aic),]
write.csv(tidy(AICctable), "perm_aic.csv")
write.csv(tidy(AICctable), "seas_aic.csv")


#####
#model 8  permanent predictions
nd_perm_8_hab <- expand.grid(Year=1984:2015,
                           Pool = c("2","3","4","5","6","7","8","9","10","12","13"),
                           Habitat = c("Main Channel", "Dyke Field", "Off Main Channel")) 
nd.predict.8.perm<-predict(hab_perm_30m_mod_8, newdata=nd_perm_8_hab, type = "response", se.fit=TRUE)

pred.vals.perm.8<-cbind(nd_perm_8_hab,nd.predict.8.perm)

perm_hab_sum <-pred.vals.perm.8 %>% 
  select(Year, Pool, Habitat, fit) %>% 
  filter(Year == 1984 | Year == 2015) %>% 
  group_by(Pool, Habitat)  %>% 
  spread(Year, fit) 
perm_hab_sum<-data.frame(perm_hab_sum)
perm_hab_sum$Difference <- ((perm_hab_sum[,4] - perm_hab_sum[,3])/perm_hab_sum[,4])*100


#model 8  seasonal predictions
nd_seas_8_hab <- expand.grid(Year=1984:2015,
                             Pool = c("2","3","4","5","6","7","8","9","10","12","13"),
                             Habitat = c("Main Channel", "Dyke Field", "Off Main Channel")) 
nd.predict.8.seas<-predict(hab_seas_30m_mod_8, newdata=nd_seas_8_hab, type="response", se.fit=TRUE)
pred.vals.seas.8<-cbind(nd_seas_8_hab,nd.predict.8.seas)

seas_hab_sum <-pred.vals.seas.8 %>% 
  select(Year, Pool, Habitat, fit) %>% 
  filter(Year == 1984 | Year == 2015) %>% 
  group_by(Pool, Habitat)  %>% 
  spread(Year, fit) 
seas_hab_sum<-data.frame(seas_hab_sum)
seas_hab_sum$Difference <- ((seas_hab_sum[,4] - seas_hab_sum[,3])/seas_hab_sum[,4])*100