#####################################################
# Marco FW Gauger
# 2021/02/23
#
# This is the script to obtain fig 4 of the manuscript submitted to PeerJ Gauger et al 2021 Diel influences on bottlenose dolphin acoustic detection in a coastal lagoon in the southwestern Gulf of California
#####################################################

ifelse("mgcv" %in% rownames(installed.packages()) == FALSE, 
       install.packages("mgcv", dependencies = T),
       ("mgcv is already installed"))
ifelse("fields" %in% rownames(installed.packages()) == FALSE, 
       install.packages("fields", dependencies = T),
       ("fields is already installed"))
ifelse("scales" %in% rownames(installed.packages()) == FALSE, 
       install.packages("scales", dependencies = T),
       ("scales is already installed"))
require(mgcv)
require(fields)
require(scales)


Sys.setlocale("LC_ALL","English")
#Sys.setenv(TZ='UTC')

analysed_data_hour <- read.csv("Supplemental Data S1.csv", sep=";")
analysed_data_hour$cluster_hcpc <- as.factor(analysed_data_hour$cluster_hcpc)



b <-
  gam(dp10m_ELAP1~ s(lunar_phase, bs="cc", k=6) +
        te(derivate_tide, tide, bs=c("tp","tp"),
           k=c(10,10)) + 
        s(SST, k=10, bs="tp", m=2) +
        s(h, bs="cc") + 
        effort +
        s(depth, k=4) +
        s(distance_BLAP, k=4),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))
AIC(b)
summary(b)


c <- #dp10m_mod_cluster_S_h <-
  gam(dp10m_ELAP1~ 
        s(lunar_phase, bs="cc", k=6) + 
        te(derivate_tide, tide, bs=c("tp","tp"),
           k=c(10,10)) + 
        s(SST, k=10, bs="tp") + effort + 
        s(h, cluster_hcpc, bs="fs", k=10, xt=list(bs="cc")) + 
        s(depth, k=4) + 
        s(distance_BLAP, k=4), 
      data=analysed_data_hour, method="REML",family="poisson", 
      knots=list(h = c(0,24)))
summary(c)
plot(c, page=1, scale=0, shade=T)

d<-#dp10m_mod_cluster_I_h <-
  gam(dp10m_ELAP1~ 
        s(lunar_phase, bs="cc") + 
        te(derivate_tide, tide, bs=c("tp","tp"),
           k=c(10,10)) + 
        s(SST, k=10, bs="tp") +
        s(h, by= cluster_hcpc, k=10,  bs="cc") + 
        effort +
        s(depth, k=4) +
        s(distance_BLAP, k=4)+
        s(cluster_hcpc, bs="re"),
      data=analysed_data_hour, method="REML",family="poisson",
      knots=list(h = c(0,24)))

summary(d)
plot(d, page=1, scale=0, shade=T)



yb1 <- plot.gam(b, pages=1, plot=F)
qb1 <- yb1[[2]] 
yc1 <- plot.gam(c, pages=1, plot=F)
qc1 <- yc1[[2]] 
yd1 <- plot.gam(d, pages=1, plot=F)
qd1 <- yd1[[2]] 




#x11()
png("Figure 5.png", width = 170, height = 280,     units = "mm", pointsize = 12, res=300)
par(mfrow=c(8,3), oma=c(4,4,1,1), mar=c(1,2,1,1))
###### moon ######
plot(b, select=1, shade=T, ylim=c(-1,1), yaxt="n", xaxt="n"); abline(h=0, lwd=2);axis(side=2, las=2);axis(side=1, at=seq(0,1,.25), labels=c("NM","Wax","FM","Wan","NM"))
mtext(side=2, text="moon", line=3)
mtext(side=3, text="global", line=.5)
plot(c, select=1, shade=T, ylim=c(-1,1), yaxt="n", xaxt="n"); abline(h=0, lwd=2);axis(side=2, las=2);axis(side=1, at=seq(0,1,.25), labels=c("NM","Wax","FM","Wan","NM"))
mtext(side=3, text="model S", line=.5)
plot(d, select=1, shade=T, ylim=c(-1,1), yaxt="n", xaxt="n"); abline(h=0, lwd=2);axis(side=2, las=2);axis(side=1, at=seq(0,1,.25), labels=c("NM","Wax","FM","Wan","NM"))
mtext(side=3, text="model I", line=.5)

###### SST ######
plot(b, select=3, shade=T, ylim=c(-2,3), yaxt="n", ylab=""); abline(h=0, lwd=2, ylab="");axis(side=2, las=2)
mtext(side=2, text="SST [°C]", line=3)
plot(c, select=3, shade=T, ylim=c(-2,3), yaxt="n", ylab="", lty=5); abline(h=0, lwd=2);axis(side=2, las=2)
plot(d, select=3, shade=T, ylim=c(-2,3), yaxt="n", ylab="", lty=5); abline(h=0, lwd=2);axis(side=2, las=2)
###### hour ######
plot(b, select=4, shade=T, ylim=c(-1,1), yaxt="n", xaxt="n", ylab="", xlab="hour"); axis(side=2, las=2, at = seq(-1,1,.5))
axis(side=1, at=seq(0,24,4))
mtext(side=2, text="hour", line=3)
plot(c, select=4, shade=T, ylim=c(-1,1), yaxt="n", xaxt="n", ylab="", xlab="", col=c("white"), lwd=.003); abline(h=0, lwd=2);axis(side=2, las=2, at=seq(-1,1,.5))

lines(qc4$x,qc4$fit[seq(1,length(qc4$fit)/3)+length(qc4$fit)/3*0], lwd=3, lty=1, col="grey35")
lines(qc4$x,qc4$fit[seq(1,length(qc4$fit)/3)+length(qc4$fit)/3*1], lwd=2, lty=1, col="white")
lines(qc4$x,qc4$fit[seq(1,length(qc4$fit)/3)+length(qc4$fit)/3*2], lwd=2, lty=1, col="white")
lines(qc4$x,qc4$fit[seq(1,length(qc4$fit)/3)+length(qc4$fit)/3*1], lwd=3, lty=2, col="magenta")
lines(qc4$x,qc4$fit[seq(1,length(qc4$fit)/3)+length(qc4$fit)/3*2], lwd=3, lty=3, col="#1B9E77")
axis(side=1, at=seq(0,24,4))
abline(h=0, lwd=2)
legend("topleft", legend=c("Aug-Oct","Oct-Mar","May-Jul"), col=c("grey35","magenta","#1B9E77"), lty=c(1,2,3), lwd=3, bty="n", cex=.8)

plot(d, select=5, ylim=c(-1,1), xaxt="n", yaxt="n", ylab="", xlim=c(0,24), shade=T, shade.col=alpha("magenta",.5), lty=2)
axis(side=1, at=seq(0,24,4))
abline(h=0, col="white", lwd=3)
par(new=T)
plot(d, select=4, ylim=c(-1,1), xaxt="n", yaxt="n", ylab="", xlim=c(0,24), shade=T, lty=4, lwd=1)
par(new=T)
plot(d, select=6, ylim=c(-1,1), xaxt="n", yaxt="n", ylab="", xlim=c(0,24), shade=T, shade.col=alpha("#1B9E77",.5),lty=3)
axis(side=2, las=2, at=seq(-1,1,.5))
legend("bottomright", cex=.8, bty="n", legend=c("Aug-Oct", "Oct-Mar","May-Jul"), pch=15, col=alpha(c(1,"magenta","#1B9E77"),.75))
###### depth ######
plot(b, select=5, shade=T, ylim=c(-3,9), yaxt="n", ylab=""); abline(h=0, lwd=2, ylab="");axis(side=2, las=2)
mtext(side=2, text="depth [m]", line=3)
plot(c, select=5, shade=T, ylim=c(-3,9), yaxt="n", ylab="", xlab=""); abline(h=0, lwd=2, ylab="");axis(side=2, las=2)
plot(d, select=7, shade=T, ylim=c(-3,9), yaxt="n", ylab="", xlab=""); abline(h=0, lwd=2, ylab="");axis(side=2, las=2)
###### distance ######
plot(b, select=6, shade=T, ylim=c(-13,3), yaxt="n", ylab="", xlab="distance"); abline(h=0, lwd=2, ylab="");axis(side=2, las=2)
mtext(side=2, text="distance \n[km]", line=3)
plot(c, select=6, shade=T, ylim=c(-13,3), yaxt="n", ylab="", xlab=""); abline(h=0, lwd=2, ylab="");axis(side=2, las=2)
plot(d, select=8, shade=T, ylim=c(-13,3), yaxt="n", ylab="", xlab=""); abline(h=0, lwd=2, ylab="");axis(side=2, las=2)

###### effort ######
plot(b, shade=T, select=7, all.terms=T,  ylim=c(-13,3), yaxt="n", ylab="", xlab="effort"); axis(side=2, las=2)
mtext(side=2, text="effort \n[min]", line=3)
plot(c, shade=T, select=7, all.terms=T,  ylim=c(-13,3), yaxt="n", ylab="", xlab=""); axis(side=2, las=2)
plot(d, shade=T, select=10, all.terms=T,  ylim=c(-13,3), yaxt="n", ylab="", xlab=""); axis(side=2, las=2)
###### season ######
plot(c(0,1), c(0,1), col="white", axes=F, ylab="", xlab="")
mtext(side=2, text="random \neffect", line=3)
plot(c(0,1), c(0,1), col="white", axes=F, ylab="", xlab="")
plot(d, shade=T, select=9, all.terms=T,  ylim=c(-13,3), yaxt="n", ylab="", main="", xlab="",xaxt="n"); axis(side=2, las=2) 
axis(side=1, at=c(-.87,0,.87), labels=c("Aug-Oct","Oct-Mar","May-Jul"))
###### tide ######
bk = seq(-1, 1, by=1/10)
mycols <- c("darkblue",colorRampPalette(colors = c("blue", "cyan", "white", "yellow","red"))(length(bk)-3),"darkred")
for(k in 1:3){
  if(k==1) q1 = qb1; D=b
  if(k==2) q1 = qc1; D=c
  if(k==2) q1 = qd1; D=d
  image(q1$x,q1$y,matrix(q1$fit,nrow=length(q1$x),
                         ncol=length(q1$y)),col=mycols,
        breaks=bk, main="",lwd=1,cex.lab=1,cex.axis=1, 
        font.axis=1, xlab="",ylab="",xaxt="n", yaxt="n", 
        axis.args=list(lwd=2, cex.axis=1, font.axis=1))
  par(new=T);plot(D, select= 2, rug=F, 
                  ylab="", xlab="", 
                  xaxt="n", yaxt="n", main="",
                  cex.axis=.001, col.lab="white" ,
                  col=scales::alpha("white",0), lwd=0.15,lty=1)
  par(new=T);plot(D, select= 2, scheme=0, rug=F, 
                  ylab="", xlab="", 
                  xaxt="n", yaxt="n",
                  col=c("darkgrey","darkgrey","darkgrey",1,"darkgrey",
                        "darkgrey","darkgrey","darkgrey"),
                  cex.axis=.001, col.lab="white" ,main="", too.far=.15, 
                  lwd=c(.5,1.75,.5))
  axis(1)
  axis(2, las=2)
  if(k ==1) mtext(side=2, line=3.5,adj=.5, "tidal\nheight [m]", cex=0.8, outer=F)
  if(k ==2) mtext(side=1, line=2.0,adj=.5, expression(paste("tidal change [m"^2,"/hour]")), cex=0.8, outer=T)
}
###### ######

dev.off()