########Chapter 2
library(tidyverse)
library(splitstackshape)
library(data.table)
library(lme4)
library(MuMIn)
library(gstat)  
library(sp)
library(spatstat)
library(data.table)
library(viridis)
library(hrbrthemes)
library(plotly)
library(jtools)
library(lmerTest)
rm(list=ls())


#################    Read in Data        ####################

data = fread("H:/My Drive/Pipelines/Chapter3/dataWYgroups.csv")
data2 = fread("H:/My Drive/Pipelines/Chapter3/dataWYgroups1_19.csv")
data3 = fread("H:/My Drive/Pipelines/Chapter3/dataWYgroups3_30.csv")
data2$YSD = as.numeric(data2$YSD)
data4 = data2
data4$YSD = ifelse(data4$pipeline == "NW",(data2$YSD -6), data2$YSD)

#--------------------Cleaning---------------------------
data10 = subset(data4, compnpp > 100 & pipenpp > 100 &YSD>0) # subset removing pixels with <10g*C*m/yr
data10$pipewoody = data10$pipeshr+data10$pipetree
data10$compwoody = data10$compshr+data10$comptree
data10$YSD2 = 1/sqrt(data10$YSD)
data10$D = 1
data10comp = data10
data10comp$D = 0
data10comp$pipenpp = data10comp$compnpp
data10comp$pipeherb = data10comp$compherb
data10comp$pipewoody = data10comp$compwoody
data10comp$YSD2 = 0
datalong = subset(rbind(data10,data10comp))

#----------------Models-----------------------------

initial2 = lm(log(pipenpp/10) ~ meanWY + Pdev + Pdev:YSD2 +meanWY:YSD2 + meanWY:Pdev + meanWY:Pdev:YSD2, data = datalong)
initial2pipe = lmer(log(pipenpp/10) ~ meanWY + Pdev + Pdev:YSD2 +meanWY:YSD2 + meanWY:Pdev + meanWY:Pdev:YSD2 + (-1 + YSD2|pipeline), data = datalong)



initial22 = lm(log(pipenpp/10) ~ meanWY + Pdev + Pdev:YSD2 +meanWY:YSD2 + meanWY:Pdev + meanWY:Pdev:YSD2, data = datalong)
initial22$coefficients[1:7] = fixef(initial2pipe)


  
r.squaredGLMM(initial2)
summary(initial2pipe)
plot(initial2pipe)

sensyy3 = lm(log(pipenpp/10) ~     meanWY + Pdev +  Pdev:pipeherb +Pdev:pipewoody + meanWY:pipeherb + meanWY:pipewoody , data =  datalong)
sensyy33 = lmer(log(pipenpp/10) ~  meanWY + Pdev +  Pdev:pipeherb +Pdev:pipewoody + meanWY:pipeherb + meanWY:pipewoody + (1|pipeline), data =  datalong)
sensyy4 = lm(log(pipenpp/10) ~     meanWY + Pdev +  Pdev:pipeherb +Pdev:pipewoody + meanWY:pipeherb + meanWY:pipewoody + YSD2:Pdev:pipeherb + YSD2:Pdev:pipewoody + YSD2:meanWY:pipeherb + YSD2:meanWY:pipewoody, data =  datalong)
sensyy44 = lmer(log(pipenpp/10) ~  meanWY + Pdev +  Pdev:pipeherb +Pdev:pipewoody + meanWY:pipeherb + meanWY:pipewoody + YSD2:Pdev:pipeherb + YSD2:Pdev:pipewoody + YSD2:meanWY:pipeherb + YSD2:meanWY:pipewoody + (1|pipeline), data =  datalong)

summary(initial2pipe)
model.sel(initial2pipe,sensyy33,sensyy44)
summary(sensyy4)
summary(initial2pipe)
r.squaredGLMM(initial2)
r.squaredGLMM(initial2pipe)
r.squaredGLMM(sensyy3)
r.squaredGLMM(sensyy33)
r.squaredGLMM(sensyy4)
r.squaredGLMM(sensyy44)



#---------------Raw Data Graphs-----------------------
ud = subset(data10comp, YSD>0 & YSD < 10)
ud$YSD2 =0
ud2 = ud
ud3 = ud
ud2$Pdev = -100
ud3$Pdev = 100

d = subset(data10, YSD>0 & YSD < 10)
d$YSD2 = 1/sqrt(5)
d2 = d
d3 = d
d2$Pdev = -100
d3$Pdev = 100

dd = subset(data10, YSD>0 & YSD < 10)
dd$YSD2 = 0
dd2 = dd
dd3 = dd
dd2$Pdev = -100
dd3$Pdev = 100
ud$pipewoody

#Woody Cover Change
ggplot(ud, aes(y = pipewoody, x = meanWY, color = "Undisturbed")) + 
  geom_point( alpha = .1, size = .5) + 
  geom_point(data =d, aes(y = pipewoody, x = meanWY,color = "Disturbed"), alpha =.1, size = .5) + 
  geom_smooth(method = "lm", size = 1.2) + 
  geom_smooth(data =d, aes(y = pipewoody, x = meanWY, color = "Disturbed"), method = "lm", size = 1.2) +
  ylab("Woody Cover (%)") + 
  xlab("Mean Annual Precipitation (mm)")+
  theme_bw() + 
  theme(
    plot.title = element_text(face = "bold", size = 12),
    legend.background = element_rect(fill = "white", size = 1, colour = "grey70"),
    legend.justification = c(0, 1),
    legend.position = c(0.02, .98),
    #legend.box = element_rect(color = "black"),
    legend.title = element_blank(),
    axis.ticks = element_line(colour = "grey80", size = 0.2),
    panel.grid.major = element_line(colour = "grey80", size = 0.2),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.x = element_text(size = 15),
    axis.title.y = element_text(size = 15),
    legend.text = element_text(size = 13)
  ) +
  scale_color_manual(breaks=c('Undisturbed', 'Disturbed'),
                     values=c('Undisturbed'='black', 'Disturbed'='red'))+
  guides(color=guide_legend(override.aes=list(fill=NA)))+
  ylim(0,65)

#ggsave("H:/My Drive/Pipelines/Chapter3/More Figures/Woodydist2.png", dpi = 700)


#Herbaceous Cover Change
ggplot(ud, aes(y = pipeherb, x = meanWY, color = "Undisturbed")) + 
  geom_point( alpha = .1, size = .5) + 
  geom_point(data =d, aes(y = pipeherb, x = meanWY,color = "Disturbed"), alpha =.1, size = .5) + 
  geom_smooth(method = "lm", size = 1.2) + 
  geom_smooth(data =d, aes(y = pipeherb, x = meanWY, color = "Disturbed"), method = "lm", size = 1.2) +
  ylab("Herbaceous Cover (%)") + 
  xlab("Mean Annual Precipitation (mm)")+
  theme_bw() + 
  theme(
    plot.title = element_text(face = "bold", size = 12),
    legend.background = element_rect(fill = "white", size = 1, colour = "grey70"),
    legend.justification = c(0, 1),
    legend.position = "none",#c(0.02, .98),
    #legend.box = element_rect(color = "black"),
    legend.title = element_blank(),
    axis.ticks = element_line(colour = "grey80", size = 0.2),
    panel.grid.major = element_line(colour = "grey80", size = 0.2),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.x = element_text(size = 15),
    axis.title.y = element_text(size = 15),
    legend.text = element_text(size = 13)
  ) +
  scale_color_manual(breaks=c('Undisturbed', 'Disturbed'),
                     values=c('Undisturbed'='black', 'Disturbed'='red'))+
  guides(color=guide_legend(override.aes=list(fill=NA)))+
  ylim(0, 70)

#ggsave("H:/My Drive/Pipelines/Chapter3/More Figures/Herbdist3.png", dpi = 700)

ggplot(ud, aes(y = compnpp/10, x = meanWY, color = "Undisturbed")) + 
  geom_point(size = .9, alpha = .2) + 
  geom_point(data =d, aes(y = pipenpp/10, x = meanWY),color = "#f24068",   size = .9, alpha = .2) +  
  geom_smooth(data =ud, aes(y = compnpp/10, x = meanWY, color = "Undisturbed"), method = "loess", span =1.3, linewidth = 1.3) +
  geom_smooth(data =d, aes(y = pipenpp/10, x = meanWY, color = "Disturbed"), method = "loess", span = .9, linewidth = 1.3) +
  ylab("NPP (g*C*yr/m^2)") + 
  xlab("Mean Annual Precipitation (mm)")+
  theme_bw() + 
  theme(
    plot.title = element_text(face = "bold", size = 12),
    legend.background = element_rect(fill = "white", size = 1, colour = "grey70"),
    legend.justification = c(0, 1),
    legend.position = c(0.02, .98),
    #legend.box = element_rect(color = "black"),
    legend.title = element_text(size = 12),
    axis.ticks = element_line(colour = "grey80", size = 0.2),
    panel.grid.major = element_line(colour = "grey80", size = 0.2),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.x = element_text(size = 15),
    axis.title.y = element_text(size = 15),
    legend.text = element_text(size = 11)
  ) +
  scale_color_manual(name = 'Actual Data',
                     breaks=c('Undisturbed', 'Disturbed'),
                     values=c('Undisturbed'='black', 'Disturbed'='red'))+
  guides(color=guide_legend(override.aes=list(fill=NA)))+
  ylim(0,800)

#ggsave("H:/My Drive/Pipelines/Chapter3/More Figures/NPPdist2.png", dpi = 700)

ud$predictedud = predict(sensyy33,ud)
dd$predictedd = predict(sensyy33, dd)
d$predictedYSD = predict(sensyy44,d)
d$predictedYSDonly = exp(predict(initial2,d))

rubyud = subset(ud) %>%
  group_by(meanWY) %>%
  summarise(pipesd = sd(pipenpp, na.rm = T),
            compsd = sd(compnpp, na.rm = T),
            predsd = sd(predictedud, na.rm = T),
            difsd = (sd(pipenpp, na.rm = T)-sd(compnpp, na.rm = T))/sd(compnpp, na.rm = T),
            upperpipe = quantile(pipenpp, .8),
            meanpipe = mean(pipenpp, na.rm = T),
            meanpredpipe = mean(exp(predictedud), na.rm = T),
            lowerpipe = quantile(pipenpp,.2),
            uppercomp = quantile(compnpp, .8),
            lowercomp = quantile(compnpp,.2),
            meancomp = mean(compnpp, na.rm =T),
            Pdevsd = sd(Pdev, na.rm =T))
rubyd = subset(dd) %>%
  group_by(meanWY) %>%
  summarise(pipesd = sd(pipenpp, na.rm = T),
            predsd = sd(predictedd, na.rm = T),
            compsd = sd(compnpp, na.rm = T),
            difsd = (sd(pipenpp, na.rm = T)-sd(compnpp, na.rm = T))/sd(compnpp, na.rm = T),
            upperpipe = quantile(pipenpp, .8),
            meanpipe = mean(pipenpp, na.rm = T),
            meanpredpipe = mean(exp(predictedd), na.rm = T),
            lowerpipe = quantile(pipenpp,.2),
            uppercomp = quantile(compnpp, .8),
            lowercomp = quantile(compnpp,.2),
            meancomp = mean(compnpp, na.rm =T),
            Pdevsd = sd(Pdev, na.rm =T))
rubyysd = subset(d) %>%
  group_by(meanWY) %>%
  summarise(pipesd = sd(pipenpp, na.rm = T),
            predsd = sd(predictedYSD, na.rm = T),
            compsd = sd(compnpp, na.rm = T),
            difsd = (sd(pipenpp, na.rm = T)-sd(compnpp, na.rm = T))/sd(compnpp, na.rm = T),
            upperpipe = quantile(pipenpp, .8),
            meanpipe = mean(pipenpp, na.rm = T),
            meanpredpipe = mean(exp(predictedYSD), na.rm = T),
            meanpredpipe2 = mean(exp(predictedYSDonly), na.rm = T),
            lowerpipe = quantile(pipenpp,.2),
            uppercomp = quantile(compnpp, .8),
            lowercomp = quantile(compnpp,.2),
            meancomp = mean(compnpp, na.rm =T),
            Pdevsd = sd(Pdev, na.rm =T))

ggplot(rubyud,aes(y = (meanpredpipe), x = meanWY,color = "Undisturbed"))+ geom_point(size = .9, alpha = .2)+ 
  geom_point(data = rubyd,aes(y = (meanpredpipe), color = "Disturbed (composition effects alone)"),   size = .9, alpha = .2) + 
  #geom_point(data = rubyysd,aes(y = (meanpredpipe), color = "Disturbed (Dist. can change PUE)"),size = .9, alpha = .2) + 
  geom_smooth( method = "loess", span = 1, size = 1.3) +  
  geom_smooth(data = rubyd,aes(y = (meanpredpipe), color = "Disturbed (composition effects alone)"), method = "loess", span = 1, size = 1.3) +
  geom_smooth(data = rubyysd,aes(y = (meanpredpipe), color = "Disturbed (comp + other dist effects)"), linetype = "dashed", method = "loess", span = 1, size = 1.3) +
  #geom_smooth(data = rubyysd,aes(y = (meanpredpipe2), color = "Post-Dist Composition & Dist Effects"), linetype = "solid", method = "loess", span = 1, size = 1.3) +
  ylab("Average NPP")+
  xlab("Mean Annual Precipitation (mm)")+
  theme_bw() + 
  theme(
    plot.title = element_text(face = "bold", size = 12),
    legend.background = element_rect(fill = "white", size = 1, colour = "grey70"),
    legend.justification = c(0, 1),
    legend.position = c(0.01, .99),
    #legend.box = element_rect(color = "black"),
    legend.title = element_text(size = 12),
    axis.ticks = element_line(colour = "grey80", size = 0.2),
    panel.grid.major = element_line(colour = "grey80", size = 0.2),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.x = element_text(size = 15),
    axis.title.y = element_text(size = 15),
    legend.text = element_text(size = 11.7)  ) +
  scale_color_manual(name = 'Model Predictions',
                     breaks=c('Undisturbed', 'Disturbed (composition effects alone)', 'Disturbed (comp + other dist effects)'),
                     values=c('Undisturbed'='black', 'Disturbed (composition effects alone)'='red', 'Disturbed (comp + other dist effects)'='darkred'))+
  guides(color=guide_legend(override.aes=list(fill=NA)))#+
  

#ggsave("H:/My Drive/Pipelines/Chapter3/More Figures/CompositionModelsAvgNPP4.png", dpi = 700)

  
ggplot(rubyud,aes(y = (predsd), x = meanWY,color = "Undisturbed"))+ 
  geom_point(data = rubyud,aes(y = (predsd), color = "Undisturbed"),   size = .9, alpha = .2)+ 
  geom_point(data = rubyd,aes(y = (predsd), color = "Disturbed (composition effects alone)"),   size = .9, alpha = .2) + 
  #geom_point(data = rubyysd,aes(y = (meanpredpipe), color = "Disturbed (Dist. can change Sensitivity)"),size = .9, alpha = .2) + 
  geom_smooth(aes(y = predsd), method = "loess", span = 1, size = 1.3) +  
  geom_smooth(data = rubyd,aes(y = (predsd), color = "Disturbed (composition effects alone)"), method = "loess", span = 1, size = 1.3) +
  geom_smooth(data = rubyysd,aes(y = (predsd), color = "Disturbed (comp + other dist effects)"), linetype = "dashed", method = "loess", span = 1, size = 1.3) +
  ylab("Temporal Production Variability") +
  xlab("Mean Annual Precipitation (mm)") +
  theme_bw() + 
  theme(
    plot.title = element_blank(),#text(face = "bold", size = 12),
    legend.background = element_rect(fill = "white", size = 1, colour = "grey80"),
    legend.justification = c(0, 1),
    legend.position = c(0.01, .99),
    #legend.box = element_rect(color = "black"),
    legend.title = element_blank(),#_text("Years Since Disturbance", size =12),
    axis.ticks = element_line(colour = "grey80", size = 0.2),
    panel.grid.major = element_line(colour = "grey90", size = 0.2),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.x = element_text(size = 15),
    axis.title.y = element_text(size = 14.7),
    legend.text = element_text(size = 11.7),
    legend.margin=margin(c(.5,5,5,5))
  ) +
  scale_color_manual(name = 'Model Predictions',
                     breaks=c('Undisturbed', 'Disturbed (composition effects alone)', 'Disturbed (comp + other dist effects)'),
                     values=c('Undisturbed'='black', 'Disturbed (composition effects alone)'='red', 'Disturbed (comp + other dist effects)'='darkred'))+
  guides(color=guide_legend(override.aes=list(fill=NA)))+
  ylim(0,.8) 

ggsave("H:/My Drive/Pipelines/Chapter3/More Figures/CompositionModelsTempNPP2.png", dpi = 700)##


#---------------PUE Graphs--------------------------


meanWY <- seq(from =100, to =1000, by =50)

Pdev <- c(100,0)
YSD2 <- c(0,1/sqrt(5), 1/sqrt(60))

# Calling expand.grid() Function
superset = expand.grid(meanWY,Pdev,YSD2)
colnames(superset) <- c("meanWY","Pdev", "YSD2")

#Get predictions and standard error of predictions
superset$predicted = predict(initial22, superset)
superset$upper95 = predict(initial22, superset, interval = "confidence")[,3]
superset$lower95 = predict(initial22, superset, interval = "confidence")[,2]


#Format for graphing
graph <- superset %>%
  group_by(meanWY,YSD2) %>%
  summarise(PUE = mean(predicted),
            up = max(predicted),
            down = min(predicted),
            range= ((max(predicted)-min(predicted))/min(predicted))*mean(exp(predicted/100)),
            upper = ((max(upper95)-min(predicted))/min(predicted))*mean(exp(upper95/100)),
            lower = ((max(lower95)-min(predicted))/min(predicted))*mean(exp(lower95/100)))
graph$YSD2= as.factor(graph$YSD2)
graph$Years_Since_Disturbance = rep(c("Undisturbed","65","5"),19)
graph$Years_Since_Disturbance <- factor(graph$Years_Since_Disturbance, c("5","65", "Undisturbed"))




#ggplot(datalong, aes(x = WYprcp, y = compnpp, color = meanWY2)) + geom_smooth(method = "lm", se = F)

datalong$meanWY = as.numeric(datalong$meanWY)
summary(datalong$yeargroup)

#Sensitivity Plot
ggplot(subset(graph, Years_Since_Disturbance == "Undisturbed"), aes(x = meanWY, y = range*10, ymin = upper*10,ymax=lower*10)) + 
  geom_smooth(aes(color = "Undisturbed"), size = 1.5, se = F)+
  geom_smooth(data =subset(graph, Years_Since_Disturbance == "5"), aes(y = range*10, ymin = upper*10,ymax=lower*10,color = "0-5 YSD"), size = 1.5, method = "loess", span = 10) +
  geom_smooth(data =subset(graph, Years_Since_Disturbance == "65"), aes(y = range*10, ymin = upper*10,ymax=lower*10,color = "60-70 YSD"), size = 1.5, method = "loess", span = 10) +
  geom_ribbon(data =subset(graph, Years_Since_Disturbance == "65"), aes(y = range*10, ymin = upper*10,ymax=lower*10),fill = "#00a6fb", alpha = .1) +
  geom_ribbon(data =subset(graph, Years_Since_Disturbance == "5"), aes(y = range*10, ymin = upper*10,ymax=lower*10),fill = "#FF0054", alpha = .1) +
  geom_ribbon(data =subset(graph, Years_Since_Disturbance == "Undisturbed"), aes(y = range*10, ymin = upper*10,ymax=lower*10),fill = "black", alpha = .1) +
  xlim(100,850)+
  ylab("Sensitivity g*C/mm") +
  xlab("Mean Annual Precipitation ") +
  #theme_classic()+
  #xlim(100,1000)+
  theme_bw() + 
  theme(
    plot.title = element_text(face = "bold", size = 12),
    legend.background = element_rect(fill = "white", size = 1, colour = "grey70"),
    legend.justification = c(0, 1),
    legend.position = c(0.01, .99),
    #legend.box = element_rect(color = "black"),
    legend.title = element_blank(),#text(size = 14.3),
    axis.ticks = element_line(colour = "grey80", size = 0.2),
    panel.grid.major = element_line(colour = "grey80", size = 0.2),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 15),
    axis.text.y = element_text(size = 15),
    axis.title.x = element_text(size = 17),
    axis.title.y = element_text(size = 17),
    legend.text = element_text(size = 13),
    plot.background = element_rect(fill = "white")
  ) + scale_color_manual(breaks=c( '0-5 YSD', '60-70 YSD','Undisturbed'),
                     values=c('Undisturbed'='black', '0-5 YSD'='#FF0054', '60-70 YSD' = "#00a6fb"))+
  guides(color=guide_legend(override.aes=list(fill=NA))) 

#ggsave("H:/My Drive/Pipelines/Chapter3/More Figures/SensitivitySimplified.png", dpi = 700)

meanWY2 <- seq(from =100, to =1000, by =10)
Pdev2 <- 0
YSD12 <- c(0,1/sqrt(2), 1/sqrt(60))

# Calling expand.grid() Function
superset2 = expand.grid(meanWY2,Pdev2,YSD12)
colnames(superset2) <- c("meanWY","Pdev", "YSD2")

#Get predictions and standard error of predictions
superset2$predicted = predict(initial22, superset2)
superset2$predselower = predict(initial22, superset2, interval = "confidence")[,2]
superset2$predseupper = predict(initial22, superset2, interval = "confidence")[,3]

#Format for graphing
superset2$YSD2 = as.factor(superset2$YSD2)
superset2$Years_Since_Disturbance = c(rep("Undisturbed",91),rep("5",91),rep("65",91))
superset2$Years_Since_Disturbance <- factor(superset2$Years_Since_Disturbance, c("5","65", "Undisturbed"))
superset2$PUE = exp(superset2$predicted)
superset2$PUEseupper = exp(superset2$predseupper)*.98
superset2$PUEselower = exp(superset2$predselower)*.97

#, ymax = PUEseupper, ymin = PUEselower
d$PUE = (d$pipenpp/10)
ud$PUE = (ud$pipenpp/10)

#PUE Plot
ggplot(subset(superset2, Years_Since_Disturbance == "Undisturbed"), aes(x = meanWY, y = PUE,  color = "Undisturbed")) + 
  geom_point(data = subset(data10, YSD <6), aes(x = meanWY, y = pipenpp/10, color = "0-5 YSD"), size = .7, alpha =.3)+
  geom_point(data = subset(data10comp, YSD <6), aes(x = meanWY, y = pipenpp/10, color = "Undisturbed"), size = .7, alpha = .3)+
  #geom_point(data = subset(data10, YSD >60 & YSD <70), aes(x = meanWY, y = pipenpp/10, color = "60-70"), size = .8)+
  geom_smooth(aes(color = "Undisturbed"), method = "gam", size = 1.3) +
  geom_smooth(data = subset(superset2, Years_Since_Disturbance == "5"), aes(x = meanWY, y = PUE,  color = "0-5 YSD"), method = "gam",size = 1.3)+
  #geom_ribbon(alpha=.4) +
  ylab("NPP") +
  xlab("Mean Annual Precipitation") +
  ylab("Average NPP") +
  #theme_classic()+
  xlim(100,850)+
  #ylim(0,950)+
  theme_bw() + 
  theme(
    plot.title = element_text(face = "bold", size = 12),
    legend.background = element_rect(fill = "white", size = 1, colour = "grey70"),
    legend.justification = c(0, 1),
    legend.position = c(0.02, .99),
    #legend.box = element_rect(color = "black"),
    legend.title = element_blank(),#text(size = 14.3),
    axis.ticks = element_line(colour = "grey80", size = 0.2),
    panel.grid.major = element_line(colour = "grey80", size = 0.2),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 15),
    axis.text.y = element_text(size = 15),
    axis.title.x = element_text(size = 17),
    axis.title.y = element_text(size = 17),
    legend.text = element_text(size = 13),
    plot.background = element_rect(fill = "white")
  )+
  scale_color_manual(breaks=c('Undisturbed', '0-5 YSD'),
                     values=c('Undisturbed'='black', '0-5 YSD'='#FF0054'))+
  guides(color=guide_legend(override.aes=list(fill=NA))) +
  ylim(0,980)

#ggsave("H:/My Drive/Pipelines/Chapter3/More Figures/ModelLines_RawPoints0_5.png", dpi = 700)


ggplot(subset(superset2, Years_Since_Disturbance == "Undisturbed"), aes(x = meanWY, y = PUE,  color = "Undisturbed")) + 
  geom_point(data = subset(data10, YSD >55 & YSD <65), aes(x = meanWY, y = pipenpp/10, color = "55-65 YSD"), size = .7, alpha = .5)+
   geom_point(data = subset(data10comp, YSD >55 & YSD <65), aes(x = meanWY, y = pipenpp/10, color = "Undisturbed"), size = .7, alpha = .4)+
  geom_smooth(data = subset(superset2, Years_Since_Disturbance == "65"), aes(x = meanWY, y = PUE,  color = "55-65 YSD"), method = "gam", size = 1.3)+
  geom_smooth(aes(color = "Undisturbed"), method = "gam", size = 1.3) +
   #geom_ribbon(alpha=.4) +
  ylab("NPP") +
  xlab("Mean Annual Precipitation") +
  ylab("Average NPP") +
  #theme_classic()+
  xlim(100,850)+
  #ylim(0,950)+
  theme_bw() + 
  theme(
    plot.title = element_text(face = "bold", size = 12),
    legend.background = element_rect(fill = "white", size = 1, colour = "grey70"),
    legend.justification = c(0, 1),
    legend.position = c(0.02, .99),
    #legend.box = element_rect(color = "black"),
    legend.title = element_blank(),#(size = 14.3),
    axis.ticks = element_line(colour = "grey80", size = 0.2),
    panel.grid.major = element_line(colour = "grey80", size = 0.2),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 15),
    axis.text.y = element_text(size = 15),
    axis.title.x = element_text(size = 17),
    axis.title.y = element_text(size = 17),
    legend.text = element_text(size = 13),
    plot.background = element_rect(fill = "white")
  )+
  scale_color_manual(breaks=c('Undisturbed', '55-65 YSD'),
                     values=c('Undisturbed'='black', '55-65 YSD'='#00a6fb'))+
  guides(color=guide_legend(override.aes=list(fill=NA))) + 
  ylim(0,980)

#ggsave("H:/My Drive/Pipelines/Chapter3/More Figures/ModelLines_RawPoints60_70.png", dpi = 700)


#calculate differences with time
ysd5 = subset(datalong, YSD < 10 & YSD >-1)
ysd5$MAP = round(ysd5$meanWY/10, 0)*10

comparison5 <- ysd5 %>%
  group_by(WYCoords) %>%
  summarize(dist_total = sum(pipenpp/10),
            comp_total = sum(compnpp/10),
            dist_sd = sd(pipenpp/10/10),
            comp_sd = sd(compnpp),
            meanWY = mean(meanWY))

ysd45 = subset(datalong, YSD < 50 & YSD >40)
ysd45$MAP = round(ysd45$meanWY/10, 0)*10

comparison45 <- ysd45 %>%
  group_by(WYCoords) %>%
  summarize(dist_total = sum(pipenpp/10),
            comp_total = sum(compnpp/10),
            dist_sd = sd(pipenpp/10/10),
            comp_sd = sd(compnpp),
            meanWY = mean(meanWY))


comparison5$percentdrop = comparison5$dist_total/comparison5$comp_total
comparison5$cummulativetotaldif = comparison5$comp_total-comparison5$dist_total
comparison45$percentdrop = comparison45$dist_total/comparison45$comp_total
comparison45$cummulativetotaldif = comparison45$comp_total-comparison45$dist_total

#YSD 0-10
below350 = subset(comparison5, meanWY <350)
above600 = subset(comparison5, meanWY >600)

mean(below350$percentdrop)
mean(above600$percentdrop)

mean(below350$cummulativetotaldif)
mean(above600$cummulativetotaldif)

#convert from g/m to kg/hectare
mean((below350$cummulativetotaldif)/1000)*10000
mean((above600$cummulativetotaldif)/1000)*10000

#--------------------------------------------
#YSD 40-50
below350_2 = subset(comparison45, meanWY <350)
above600_2 = subset(comparison45, meanWY >600)

mean(below350_2$percentdrop)
mean(above600_2$percentdrop)

mean(below350_2$cummulativetotaldif)
mean(above600_2$cummulativetotaldif)

#convert from g/m to kg/hectare
mean((below350_2$cummulativetotaldif)/1000)*10000
mean((above600_2$cummulativetotaldif)/1000)*10000


#---------------Response to drought-----------------
ysddrought5 = subset(datalong, YSD < 10 & YSD >-1 & WYprcp < meanWY*.66)
ysddrought5$percentdroughtdif = ysddrought5$pipenpp/ysddrought5$compnpp
ggplot(ysddrought5, aes(y= percentdroughtdif, x = meanWY)) + geom_smooth(method = 'loess')

ysddrought45 = subset(datalong, YSD < 50 & YSD >40 & WYprcp < meanWY*.66)
ysddrought45$percentdroughtdif = ysddrought45$pipenpp/ysddrought45$compnpp
ggplot(ysddrought45, aes(y= percentdroughtdif, x = meanWY)) + geom_smooth(method = 'loess')

ysdabove5 = subset(datalong, YSD < 10 & YSD >-1 & WYprcp < meanWY*1.33)
ysdabove5$percentabovedif = ysdabove5$pipenpp/ysdabove5$compnpp
ggplot(ysdabove5, aes(y= percentabovedif, x = meanWY)) + geom_smooth(method = 'loess')

ysdabove45 = subset(datalong, YSD < 50 & YSD >40 & WYprcp < meanWY*1.33)
ysdabove45$percentabovedif = ysdabove45$pipenpp/ysdabove45$compnpp
ggplot(ysdabove45, aes(y= percentabovedif, x = meanWY)) + geom_smooth(method = 'loess')




h = superset2 %>%
  group_by(YSD2, meanWY) %>%
  summarize(NPP = mean(PUE))
hh = pivot_wider(h, names_from = YSD2, values_from = NPP)

jj=data.frame(ggplot_build(a)[[1]])
jj$MAP = round(jj$x)
jj$y2 = exp(jj$y)-1
dif2 = pivot_wider(jj[,c(10,15,16)],names_from = group, values_from = y2)



data$lat = as.numeric(substr(data$WYCoords,12,19))/1000000
data$long = as.numeric(substr(data$WYCoords,1,10))/1000000

ruby = subset(data10,  YSD>-10 & YSD<10)
ruby$D = ifelse(ruby$YSD >-1,1,0)

rnum <- ruby %>%
  group_by(lat,long) %>%
  summarize(map = mean(meanWY),
            deserttype = mean(bio8))
rnum$number = 1:998

ruby2 = merge(ruby,rnum, by = c("lat","long"))

for (k in 1:998) {
  bb = subset(ruby2, number == k)
  b = lm((bb$pipenpp)~bb$Pdev*bb$D)
  rnum$disteffect[k] = b$coefficients[4]
  rnum$predistsensitivity[k] = b$coefficients[2]
  rnum$predistPUE[k] = b$coefficients[1]
  rnum$PdevUD[k] = b$coefficients[2]
  rnum$disteffectPUE[k] = b$coefficients[3]
  rnum$rsqr[k] = r.squaredGLMM(b)[1]
  #rnum$pvalue[k] = summary(b)$coefficients[4,4]
  #rnum$pvaluePUE[k] = summary(b)$coefficients[3,4]
  rnum$woodycover = mean(subset(bb, D == 0)$pipetree)
  c = subset(bb, D == 0)
  cc = subset(bb,D ==1)
  rnum$changeinwoodycover[k] = (mean(c$pipetree +c$pipeshr)-mean(cc$pipetree +cc$pipeshr))#/mean(c$pipetree +c$pipeshr)
  rnum$changeinshrubcover[k] = (mean(cc$pipeshr)-mean(c$pipeshr))#/mean(c$pipetree +c$pipeshr)
  rnum$UDshrubcover[k] = (mean(c$pipeshr))#/mean(c$pipetree +c$pipeshr)
  rnum$Dshrubcover[k] = (mean(cc$pipeshr))#/mean(c$pipetree +c$pipeshr)
  rnum$UDafgccover[k] = (mean(c$pipeafgc))#/mean(c$pipetree +c$pipeshr)
  rnum$Dafgccover[k] = (mean(cc$pipeafgc))#/mean(c$pipetree +c$pipeshr)
  rnum$UDtreecover[k] = (mean(c$pipetree))#/mean(c$pipetree +c$pipeshr)
  rnum$Dtreecover[k] = (mean(cc$pipetree))#/mean(c$pipetree +c$pipeshr)
  rnum$UDpfgccover[k] = (mean(c$pipepfgc))#/mean(c$pipetree +c$pipeshr)
  rnum$Dpfgccover[k] = (mean(cc$pipepfgc))#/mean(c$pipetree +c$pipeshr)
  rnum$proportionwoody[k] = mean(c$pipetree +c$pipeshr)/mean(c$totalplantcover)
  rnum$proportionwoodyD[k] = mean(cc$pipetree +cc$pipeshr)/mean(cc$totalplantcover)
  rnum$proportionherb[k] = mean(c$pipeafgc +c$pipepfgc)/mean(c$totalplantcover)
  rnum$proportionherbD[k] = mean(cc$pipeafgc +cc$pipepfgc)/mean(cc$totalplantcover)
  rnum$proportionshrub[k] = mean(c$pipeshr)/mean(c$totalplantcover)
  rnum$proportionshrubD[k] = mean(cc$pipeshr)/mean(cc$totalplantcover)
  rnum$proportionafgc[k] = mean(c$pipeafgc )/mean(c$totalplantcover)
  rnum$proportionafgcD[k] = mean(cc$pipeafgc )/mean(cc$totalplantcover)
  rnum$proportionpfgc[k] = mean(c$pipepfgc)/mean(c$totalplantcover)
  rnum$proportionpfgcD[k] = mean(cc$pipepfgc)/mean(cc$totalplantcover)
  rnum$changeinherbcover[k] = mean(cc$pipeafgc +cc$pipepfgc)-(mean(c$pipeafgc +c$pipepfgc))#/mean(c$pipetree +c$pipeshr)
  rnum$changeintotalcover[k] = mean(cc$totalplantcover)/(mean(c$totalplantcover))#/mean(c$pipetree +c$pipeshr)
  rnum$changeinafgccover[k] = (mean(cc$pipeafgc)-mean(c$pipeafgc))#/mean(c$pipetree +c$pipeshr)
  rnum$changeinpfgccover[k] = mean(cc$pipepfgc)-mean(c$pipepfgc)#/mean(c$pipetree +c$pipeshr)
  rnum$changeintreecover[k] = (mean(cc$pipetree )-mean(c$pipetree ))#/mean(c$pipetree +c$pipeshr)
  rnum$totalcover[k] = mean(c$totalplantcover)
  rnum$woodycover2[k] = mean(c$pipewoody)
  rnum$predistNPP[k] = mean(c$pipenpp)
}

ggplot(rnum, aes((y=predistPUE+disteffectPUE)/10,x=map, color = "Disturbed")) + 
  geom_point(alpha =.3)+ 
  geom_point(aes(y = predistPUE/10, color = "Undisturbed"),alpha =.3, , span = 40)+
  geom_smooth(aes(y = predistPUE/10, color = "Undisturbed"),method = 'loess')+
  geom_smooth(method = 'gam', span = 1)+
  xlab("Mean Annual Precipitation")+ 
  ylab("Average NPP") + 
  theme_bw()+
  theme(
    plot.title = element_text(face = "bold", size = 12),
    legend.background = element_rect(fill = "white", size = 1, colour = "grey70"),
    legend.justification = c(0, 1),
    legend.position = c(0.6, .2),
    #legend.box = element_rect(color = "black"),
    legend.title = element_blank(),
    axis.ticks = element_line(colour = "grey80", size = 0.2),
    panel.grid.major = element_line(colour = "grey80", size = 0.2),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.x = element_text(size = 15),
    axis.title.y = element_text(size = 15),
    legend.text = element_text(size = 13)
  ) +
  scale_color_manual(breaks=c('Undisturbed', 'Disturbed'),
                     values=c('Undisturbed'='black', 'Disturbed'='red'))+
  guides(color=guide_legend(override.aes=list(fill=NA))) 

ggsave("H:/My Drive/Pipelines/Chapter3/More Figures/PrePostProduction.png", dpi = 700)

ggplot(rnum, aes((y=predistsensitivity+disteffect)/10,x=map, color = "Disturbed")) + 
  ylim(-.3,.4) +geom_point(alpha =.3)+ 
  geom_smooth(method = 'lm', span = 1)+
  geom_smooth(aes(y = predistsensitivity/10, color = "Undisturbed"),method = 'lm')+
  geom_point(aes(y = predistsensitivity/10, color = "Undisturbed"),alpha =.3, , span = 40)+
  xlab("Mean Annual Precipitation")+ 
  ylab("Sensitivity") + 
  theme_bw()+
  theme(
    plot.title = element_text(face = "bold", size = 12),
    legend.background = element_rect(fill = "white", size = 1, colour = "grey70"),
    legend.justification = c(0, 1),
    legend.position = c(0.6, .2),
    #legend.box = element_rect(color = "black"),
    legend.title = element_blank(),
    axis.ticks = element_line(colour = "grey80", size = 0.2),
    panel.grid.major = element_line(colour = "grey80", size = 0.2),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.x = element_text(size = 15),
    axis.title.y = element_text(size = 15),
    legend.text = element_text(size = 13)
  ) +
  scale_color_manual(breaks=c('Undisturbed', 'Disturbed'),
                     values=c('Undisturbed'='black', 'Disturbed'='red'))+
  guides(color=guide_legend(override.aes=list(fill=NA)))


ggplot(rnum, aes((y=(disteffect-predistsensitivity)/predistsensitivity),x=map)) +geom_point(alpha =.3)+ geom_smooth(method = 'lm')+scale_color_gradient(low="blue", high="red") + xlab("Pre-Disturbance NPP")+ ylab("Effect of Disturbance on PUE") + theme_classic() 



## ------------------------------------Raw Data Plots--------------------------------------------------------

ysd10 = subset(data10, YSD>-1 &YSD<11)
ysd60 = subset(data10, YSD>60 &YSD<70)
ysd30 = subset(data10, YSD>15 &YSD<30)

ggplot(ysd10, aes(y = pipenpp/10, x = meanWY)) + 
  geom_jitter(aes(color = "Disturbed"), alpha = .3, size = .4)+ 
  geom_point(aes(y = compnpp/10, color = "Undisturbed"), alpha = .3, size = .4)+
  geom_smooth(aes(color = "Disturbed"), method = "loess", size = 1.3)+
  geom_smooth(aes(y = compnpp/10),color = "grey20", method = "loess",size = 1.3)+
  xlab("Mean Annual Precipitation")+ 
  ylab("NPP") + 
  theme_bw()+
  ggtitle("0-10 Years After Disturbance")+
  theme(
    plot.title = element_text(face = "bold", size = 15, hjust = 0.5),
    legend.background = element_rect(fill = "white", size = 1, colour = "grey70"),
    legend.justification = c(0, 1),
    legend.position = c(0.01, .99),
    #legend.box = element_rect(color = "black"),
    legend.title = element_blank(),
    axis.ticks = element_line(colour = "grey80", size = 0.2),
    panel.grid.major = element_line(colour = "grey80", size = 0.2),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.x = element_text(size = 15),
    axis.title.y = element_text(size = 15),
    legend.text = element_text(size = 13)
  ) +
  scale_color_manual(breaks=c('Undisturbed', 'Disturbed'),
                     values=c('Undisturbed'='black', 'Disturbed'='red'))+
  guides(color=guide_legend(override.aes=list(fill=NA)))+
  ylim(0,1000)+
  xlim(100,850)

#ggsave("H:/My Drive/Pipelines/Chapter3/More Figures/0-10NPPraw.png", dpi = 700)


ggplot(ysd30, aes(y = pipenpp/10, x = meanWY)) + 
  geom_jitter(color = "red", alpha = .3, size = .4)+ 
  geom_point(aes(y = compnpp/10), color = "black", alpha = .3, size = .4)+
  geom_smooth(color = "red", method = "loess", size = 1.3)+
  geom_smooth(aes(y = compnpp/10), color = "grey20", method = "loess", size = 1.3)+
  xlab("Mean Annual Precipitation")+ 
  ylab("NPP") + 
  theme_bw()+
  ggtitle("20-30 Years After Disturbance")+
  theme(
    plot.title = element_text(face = "bold", size = 15, hjust = .5),
    legend.background = element_rect(fill = "white", size = 1, colour = "grey70"),
    legend.justification = c(0, 1),
    legend.position = c(0.6, .2),
    #legend.box = element_rect(color = "black"),
    legend.title = element_text(size = 18),
    axis.ticks = element_line(colour = "grey80", size = 0.2),
    panel.grid.major = element_line(colour = "grey80", size = 0.2),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.x = element_text(size = 15),
    axis.title.y = element_text(size = 15),
    legend.text = element_text(size = 13)
  ) +
  scale_color_manual(breaks=c('Undisturbed', 'Disturbed'),
                     values=c('Undisturbed'='black', 'Disturbed'='red'))+
  guides(color=guide_legend(override.aes=list(fill=NA)))+
  ylim(0,500)
#ggsave("H:/My Drive/Pipelines/Chapter3/More Figures/20-30NPPraw.png", dpi = 700)


ggplot(ysd60, aes(y = pipenpp/10, x = meanWY)) + 
  geom_jitter(color = "red", alpha = .3, size = .4)+ 
  geom_point(aes(y = compnpp/10), color = "black", alpha = .3, size = .4)+
  geom_smooth(color = "red", method = "loess", span = 3, size = 1.3)+
  geom_smooth(aes(y = compnpp/10), color = "grey20", method = "loess", span = 1.1, size = 1.3)+
  xlab("Mean Annual Precipitation")+ 
  ylab("NPP") + 
  theme_bw()+
  ggtitle("60-70 Years After Disturbance")+
  theme(
    plot.title = element_text(face = "bold", size = 15, hjust = 0.5),
    legend.background = element_rect(fill = "white", size = 1, colour = "grey70"),
    legend.justification = c(0, 1),
    legend.position = c(0.6, .2),
    #legend.box = element_rect(color = "black"),
    legend.title = element_text(size = 18),
    axis.ticks = element_line(colour = "grey80", size = 0.2),
    panel.grid.major = element_line(colour = "grey80", size = 0.2),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 13),
    axis.text.y = element_text(size = 13),
    axis.title.x = element_text(size = 15),
    axis.title.y = element_text(size = 15),
    legend.text = element_text(size = 13)
  ) +
  scale_color_manual(breaks=c('Undisturbed', 'Disturbed'),
                     values=c('Undisturbed'='black', 'Disturbed'='red'))+
  guides(color=guide_legend(override.aes=list(fill=NA)))+
  ylim(0,1000)+
  xlim(100,850)
ggsave("H:/My Drive/Pipelines/Chapter3/More Figures/60-70NPPraw.png", dpi = 700)


ggplot(ysd30, aes(y = pipenpp/10, x = meanWY)) + geom_point(color = "red", alpha = .2)+ geom_point(aes(y = compnpp/10), color = "black", alpha = .2)+ 
  geom_smooth(color = "red", method = "loess")+ geom_smooth(aes(y = compnpp/10), color = "black", method = "loess")

library(ggridges)
datalong$D= as.factor(datalong$D)
lowlong = subset(datalong, WYCoords == "-119866700_41821068")

ggplot(subset(lowlong), aes(x = pipenpp/10, y = D)) +
  geom_density_ridges() +
  theme_ridges() + 
  theme(legend.position = "none")

midlong = subset(datalong, WYCoords == "-120235009_42162428")
ggplot(subset(midlong), aes(x = pipenpp/10, y = D)) +
  geom_density_ridges() +
  theme_ridges() + 
  theme(legend.position = "none")


b = hist(low$pipenpp/10,breaks = 40, plot = FALSE)
a = hist(low$compnpp/10,breaks = 40, plot = FALSE)
plot(a)
plot(b, col = "red", add = TRUE)

b = hist(high$pipenpp/10,breaks = 20, plot = FALSE)
a = hist(high$compnpp/10,breaks = 40, plot = FALSE)
plot(a)
plot(b, col = "red", add = TRUE)


corly = datalong %>%
  group_by(D,meanWY,WYCoords)%>%
  summarize(PdevCor = cor(Pdev,pipenpp))

ggplot(corly, aes(y = PdevCor, x = meanWY, color = D)) + geom_point()+ geom_smooth()

#=========================Test Autocorrelation===============================================

#-----------------Variogram Using geoR-------------------
library(geoR)
#subset data
datalong$distresiduals = (predict(initial2pipe, datalong)-log(datalong$pipenpp))
datalong$compresiduals = sensyy3$residuals
datalong$compdistresiduals = sensyy4$residuals

datalong$lat = as.numeric(substr(datalong$WYCoords,12,19))/1008
datalong$long = as.numeric(substr(datalong$WYCoords,1,10))/1008

sub = na.omit(subset(datalong, year == 2013 ))
sub2 = na.omit(subset(datalong, year == 2013 & YSD2 == 0))

EC971 <- as.geodata(sub, coords.col = 46:47, data.col = 53)
EC972 <- as.geodata(sub, coords.col = 46:47, data.col = 54)
EC973 <- as.geodata(sub, coords.col = 46:47, data.col = 55)

dup.coords(EC971)
Var1 <- variog(EC971, trend = "1st", max.dist = 55)
Var2 <- variog(EC972, trend = "1st", max.dist = 55)
Var3 <- variog(EC973, trend = "1st", max.dist = 55)

plot.new()
par(mfrow = c(3, 1))
plot(Var1, pch = 19, col = "blue", main = "Disturbance Only Model", xlab = "distance (km)")
plot(Var2, pch = 19, col = "red", main = "Composition Only Model", xlab = "distance (km)")
plot(Var3, pch = 19, col = "black", main = "Composition + Disturbance Model", xlab = "distance (km)")

ini.vals <- expand.grid(seq(0, .5, by = .01), seq(0, 1, by = .01))
ols1 <- variofit(Var1, ini = ini.vals, fix.nug = TRUE, wei = "equal", max.dist = .1)
wls1 <- variofit(Var1, ini = ini.vals, fix.nug = FALSE, fix.kappa = F, wei = "npairs", max.dist = .1)
summary(ols1)
summary(wls1)

plot.new()
par(mfrow = c(1, 1))
plot(Var1, pch = 19, col = "blue", main = "Variogram GSL")
par(new = TRUE)
plot(Var2, pch = 19, col = "red", xaxt = "n", yaxt = "n")
lines(wls1)
lines(ols1, lty = 2, col = "blue")



subb = subset(sub, select = c("lat", "long"))
subb$lat = subb$lat*111
subb$long = subb$long*111
dd = dist(subb)

quantile(dd, c(.01,.02,.03,.04,.05,.06))
hist(dd*111, breaks = 50, main = "Histogram of Distances")
abline(v =.6)

dd2 = as.numeric(toString(dd))
hist(dd2)
dist.mat <- dist(sub2, long, lat)

#Load datasets from Earth Engine
allNPP = fread("ALLpipeNPP.csv", header = T)
allcompNPP = fread("ALLcompNPP.csv", header = T)

#Create spatial object to calculate closest neighbor
allNPP$lat = allNPP$lat/1008
allNPP$long = allNPP$long/1008
allcompNPP$lat = allcompNPP$lat/1008
allcompNPP$long = allcompNPP$long/1008
X = ppp(allNPP$lat, allNPP$long, c(min(allNPP$lat)-10000,max(allNPP$lat)+10000), c(min(allNPP$long)-10000,max(allNPP$long)+10000))

Y = ppp(allcompNPP$lat, allcompNPP$long, c(min(allcompNPP$lat)-10000,max(allcompNPP$lat)+10000), c(min(allcompNPP$long)-10000,max(allcompNPP$long)+10000))

#Create vector contaitning which row of Y contains the closest undisturbed pixel
N <- nncross(X,Y, what = "which") #closest
N2 <- nncross(X,Y, k=2, what = "which") #Second closest
N3 <- nncross(X,Y, k=3, what = "which") #Third closest

NN <- nncross(X,Y, what = "dist") # distance to closest



source("G:/My Drive/Pipelines/Functions.R")


list.coefficients.comp<-list()
list.coefficients.pipe<-list()
list.r.squared.comp<-list()
list.r.squared.pipe<-list()

for(i in 1:100)
{
  
  # Stratify by region and the MAP of the region #
  stratified_start<-stratified(datalong, "meanWY", .01)
  stratified_final<-datalong[datalong$WYCoords %in% stratified_start$WYCoords,]
  
  # Climate interaction with annual precipitation deviation, varies by ecoregion #
  boot_pipe = lmer(log(pipenpp/10) ~ meanWY + Pdev + Pdev:YSD2 +meanWY:YSD2 + meanWY:Pdev + meanWY:Pdev:YSD2+ (-1 + YSD2|pipeline) , data = stratified_final)
  
  #Now get and store information from the models in the lists 
  
  # Coefficients #
  #list.coefficients.map_npp[[i]] <- df.coef.create2(stratified_final_npp)
  list.coefficients.pipe[[i]] <- df.coef.create(boot_pipe)
  
  # r-squared: total #
  #list.r.squared.map_npp[[1]] <- df.r.squared.create(stratified_final_npp)
  list.r.squared.pipe[[i]] <- r.squaredGLMM(boot_pipe)[1]
}

# Get Coefficients and 95% confidence intervals    
Get_Coef(list.coefficients.comp)
Get_Coef(list.coefficients.pipe)

r2_totalcomp<- bind.r.squared2(list.r.squared.comp)
r2_totalpipe<- bind.r.squared2(list.r.squared.pipe)

mean(r2_totalcomp)
mean(r2_totalpipe)



summary(boot_pipe)
& YSD2 == 0

datalong2 = subset(stratified_final, year == 2015)
library(ade4)
station.dists <- dist(cbind(datalong2$long, datalong2$lat))
ozone.dists <- dist(datalong2$distresiduals)
ozone.dists <- dist(datalong2$compnpp)


as.matrix(station.dists)
as.matrix(ozone.dists)

plot(mantel.rtest(station.dists, ozone.dists, nrepet = 99))





breaks = seq(0, 90, by =3)
dists <- dist(ozone[,3:4]) 