# North-Aral

library(data.table)
library(ggplot2)
library(ggpubr)
library(lubridate)
library(seefo)
library(geomtextpath)


rm(list=ls())

zs <- 2.77 # depth of surface sensor
zb <- 11.77 # depth of bottom sensor (10.02)

# Modelled data
out <- fread(file = "aral_sensitivity_out.csv")

# observed data
dat <- fread(file = "north-aral-temp-daily.csv")


sdepths <- unique(out[,depth])
sdepths <- sdepths[order(sdepths)]

sexts <- unique(out[,ext])
sexts <- sexts[order(sexts)]

ssalins <- unique(out[,sal])
ssalins <- ssalins[order(ssalins)]


# evaluation --------------------------------------------------------------

out[,month:=as.POSIXlt(date)$mon+1]

july <- out[month==7,.(Ts_Tb=mean(Ts_Tb)),.(depth,ext,sal,month)]


strat <- out[,lapply(.SD, function(x) {
    o2 <- strat_phenology(Ts=Ts, Tb=Tb, dates=date)
    return(mean(as.numeric(o2[,"MaxStratDur"]), na.rm=TRUE))
  }), by=.(flake,depth,ext,sal), 
  .SDcols=c("No")]

# This calculation of stratification duration generates some NAs for 
# shallow and clear combinations. I am pretty sure that in these cases
# the lake never stratifies, so the strat_phenology function returns 
# an NA or logical value. I set these to zero because they correspond to 
# complete polymixis and zero stratification.

strat[,depth:=as.numeric(depth)]
strat[,ext:=as.numeric(ext)]
strat[,sal:=as.numeric(sal)]

strat[is.na(No),No:=0]

setkey(x = july,cols="depth")
setkey(x = strat,cols="depth")

# plot and evaluate sensitivity -------------------------------------------


# 
# library(plot3D)
# library(ggplot2)

## Try to smooth stratification durations for contour plot
d1 <- strat[sal==11 & depth<16 & ext < 1.2]

lomod <- loess(No ~ depth + ext, data = d1, span=0.25)

df <- with(d1, 
           expand.grid(ext = seq(min(ext), max(ext), by=0.01),
                       depth = seq(min(depth), max(depth), by=0.1)))

df$No <- c(predict(lomod, newdata = df, se = FALSE))
setDT(df, key = c("depth","ext"))

p1 <- ggplot(july[sal==11 & depth<16 & ext < 1.2]) +
  geom_raster(aes(depth, ext, fill=Ts_Tb), interpolate=TRUE) +
  geom_textcontour(data = df,
               aes(depth,ext,z=No), colour="black", 
               breaks = c(0,50,100,150, 180),
               straight=TRUE,upright=TRUE,
               size=3) +
  coord_cartesian(xlim=c(6,13), ylim=c(0.3,0.9)) +
  theme_classic() +
  xlab("Depth (m)") + 
  ylab(expression(Extinction~(m^-1))) +
  scale_fill_distiller(palette="Spectral", 
                       name=expression(italic(T[s])-italic(T[b])~(degree*C)))

p1
# ggsave(filename = "../plots/aral-july-Ts-Tb-sensitivity.png", 
#        plot = p1, width = 14, height = 12, units = "cm", dpi = 300)

# create a plot with 3 panels from TWS model: 
# 1) Tb vs time for different depths at obs extinction
# 2) sames as 1) but different extinction
# 3) sensitivity plot
# or aa 2 panel plot with just 1) and 3) 

out[,depth:=as.factor(depth)]
plotdepths <- 6:16 #c(6,8,10,12,14,16,18)

p2 <- ggplot() +
  xlab("") +
  ylab(expression(Temperature~(degree*C))) +
  theme_classic() +
  geom_line(data=out[flake=="flake2TWS" & ext==0.5 & sal==11 &
                       date>=dat[1,date] & date <= as.POSIXct("2019-10-01",tz="UTC") &
                       depth %in% plotdepths],
            aes(x=date, y=Tb, colour=depth), size=1.0) +
  geom_point(data=dat[DEPTH==10.02], # 11.77
             aes(x=date, y=WTEMP), 
             colour="blue", size=1.2,
             alpha=0.7) +
  geom_point(data=dat[DEPTH==zb], # 11.77
             aes(x=date, y=WTEMP), 
             colour="black", size=1.2,
             alpha=0.7) +
  
  scale_colour_brewer(palette = "Spectral", name="Depth (m)")  # "PuOr" "RdYlBu" "Spectral"

p2

p3 <- ggplot(strat[sal==11 & depth<16 & ext < 1.2]) +
  geom_raster(aes(depth, ext, fill=No), interpolate=TRUE) +
  coord_cartesian(xlim=c(6,13), ylim=c(0.3,0.9)) +
   theme_classic() +
  xlab("Depth (m)") + 
  ylab(expression(Extinction~(m^-1))) +
  scale_fill_distiller(palette="Spectral", 
                       name=expression(Strat~(d)))
p3

p4 <- ggplot() +
  xlab("") +
  ylab(expression(Temperature~(degree*C))) +
  theme_classic() +
  geom_line(data=out[flake=="flake2TWS" & ext==0.5 & sal==11 &
                       date>=dat[1,date] & date <= as.POSIXct("2019-10-01",tz="UTC") &
                       depth %in% plotdepths],
            aes(x=date, y=Ts, colour=depth), size=1.0) +
  geom_point(data=dat[DEPTH==zs], # 11.77
             aes(x=date, y=WTEMP), 
             colour="black", size=1.5,
             alpha=0.7) +
  scale_colour_brewer(palette = "Spectral", name="Depth (m)")  # "PuOr" "RdYlBu" "Spectral"

p4


figure <- ggarrange(p2, p1, labels=c("a","b"), ncol=1, nrow=2)
figure

# ggsave(filename = "../plots/aral-sensitivity.png", 
#        plot = figure, width = 12, height = 14, 
#        units = "cm", dpi = 300)
# ggsave(filename = "../plots/aral-sensitivity.eps", 
#        plot = figure, width = 12, height = 14, 
#        units = "cm", dpi = 300)

figure2 <- ggarrange(p2, p3, labels=c("a","b"), ncol=1, nrow=2)
figure2

# ggsave(filename = "../plots/aral-sensitivity-strat.png", 
#        plot = figure2, width = 12, height = 14, units = "cm", dpi = 300)

p4 <- ggplot(data=dat[DEPTH %in% c(zb)]) + # 2.22, 10.02, 11.77
  geom_point(aes(x=date, y=WTEMP, col=DEPTH)) +
  geom_line(aes(x=date, y=WTEMP, col=DEPTH)) +
  scale_colour_distiller(palette = "Spectral", name="Depth (m)")  # "PuOr" "RdYlBu" "Spectral" + 


p4

# ice
out[,ice:=H_ice>0]
modon <- out[diff(c(0,ice))==1 & sal==11 & ext ==0.5 & depth==11,date][1]
modoff <- out[diff(c(0,ice))==-1 & sal==11 & ext ==0.5 & depth==11,date][2]


iceon <- ISOdate(2018,12,11,0, tz="UTC") # from Georgiy  ISOdate(2018,12,9,0, tz="UTC")
iceoff <- ISOdate(2019,3,31,0, tz="UTC")  # from Georgiy

modon - iceon
modoff - iceoff
iceoff - iceon
modoff - modon

depths <- c(2.22,5.82,10.02,11.77)

cols <- c("black","green","red","blue")
CEX <- 0.4

plot_ice <- function() {
  plot(WTEMP~date,data=dat, 
       col=cols[1], type="n", xlab="", ylab="Temperature")
  for(i in 1:length(depths)) {
    points(WTEMP~date,data=dat,subset=DEPTH==depths[i], 
           col=cols[i], type="o", cex=CEX)
  }
  abline(v=c(iceon,iceoff))
  
  lines((H_ice*10)~date,out[depth==10 & sal==11 & ext==0.5], 
        lty=2)
  
  legend("top", c(paste(depths,"m"),"ice"), 
         col=c(cols[1:length(depths)],"black"), 
         pch=c(rep(1,length(depths)),NA), pt.cex=CEX,
         lty=c(rep(1,length(depths)),2), bty="n")
}

plot_ice()

# Goodness of fit ---------------------------------------------------------


Ts <- out[depth==11 & sal==11 & ext==0.5,.(date,Ts)][
  dat[DEPTH==zs,.(date,Ts.obs=WTEMP)], on=.(date=date)
]

Tb <- out[depth==11 & sal==11 & ext==0.5,.(date,Tb)][
  dat[DEPTH==zb,.(date,Tb.obs=WTEMP)], on=.(date=date)
]

Tb7 <- out[depth==7 & sal==11 & ext==0.5,.(date,Tb)][
  dat[DEPTH==zb,.(date,Tb.obs=WTEMP)], on=.(date=date)
]

Tb7a <- out[depth==7 & sal==11 & ext==0.5,.(date,Tb)][
  dat[DEPTH==6.92,.(date,Tb.obs=WTEMP)], on=.(date=date)
]


rmse <- function(mod,obs) {
  sqrt(sum((mod-obs)^2)/length(obs))
}

rmse(Ts$Ts,Ts$Ts.obs) # Ts error
rmse(Tb$Tb,Tb$Tb.obs) # Tb error
rmse(Tb7$Tb,Tb7$Tb.obs) # Tb error with model 7m deep compared to actual lake temp bottom at zb=11.77m
rmse(Tb7a$Tb,Tb7a$Tb.obs) # Tb error with model 7m deep compared to actual lake temp bottom at DEPTH=6.92m


# calculations for manuscript ---------------------------------------------

dat[,month:=month(date)]
dat[,year:=year(date)]
dat[DEPTH==11.77,mean(WTEMP),by=.(year,month)] # mean monthly observed Tb at real bottom
dat[DEPTH==10.02,mean(WTEMP),by=.(year,month)] # mean monthly observed Tb at bottom of homothermal winter layer
out[depth==11&sal==11&ext==0.5,mean(Tb),by=month] # mean monthly modelled Tb 
out[depth==11&sal==11&ext==0.5,min(Tb),by=month] # min monthly modelled Tb 

# stratification duration

xx<-out[depth==11&sal==11&ext==0.5,.(date,Ts,Tb)] # modelled data for strat
strat_phenology(Ts = xx$Ts, Tb = xx$Tb, dates = xx$date) # modelled stratification duration
yy<-dat[DEPTH==2.77, .(date,Ts=WTEMP)] # observed data for strat
yy <- yy[dat[DEPTH==11.77, .(date,Tb=WTEMP)],on=.(date)]
yy <- yy[!is.na(Ts)]
strat_phenology(Ts = yy$Ts, Tb = yy$Tb, dates = yy$date)
yy[,Ts_Tb:=.(Ts - Tb)]
yy[,strat:=.(Ts_Tb>1)]
yy[,month:=month(date)]
yy[,year:=year(date)]
yy[,swtch:=c(0,diff(strat))]
straton <- yy[swtch==1,date]
stratoff <- yy[swtch==-1,date]
stratobs <- data.table(straton,stratoff=c(stratoff,NA))
stratobs[,dur:=.(stratoff - straton)]
stratobs
mean(stratobs[,dur], na.rm=TRUE)
median(stratobs[,dur], na.rm=TRUE)

cbind(july[ext==0.5 & depth%in% 5:14],
 diff=c(NA,diff(july[ext==0.5 & depth%in% 5:14,Ts_Tb])))
cbind(strat[ext==0.5 & depth%in% 5:14], # stratificaiton duration
  diff=c(NA,diff(strat[ext==0.5 & depth%in% 5:14,No]))) # change in strat per 1 m depth change
cbind(df[ext==0.5 & depth%in% 5:14], # smoothed stratification duration
  diff=c(NA,diff(df[ext==0.5 & depth%in% 5:14,No]))) # change in strat per 1 m depth change