Ice phenology plots

caze_ice.r

# Libraries
# devtools::install_github("SwampThingPaul/AnalystHelper")
library(AnalystHelper);
library(plyr)
library(reshape2)

# -------------------------------------------------------------------------
# DATA_CSV = data location on PC
caze.ice = read.csv(DATA_CSV,sep="\t",na.strings = "-999")

caze.ice$FirstFreezeDate = with(caze.ice,date.fun(as.character(FirstFreezeDate),form = "%Y%m%d"))
caze.ice$IceOutDate = with(caze.ice,date.fun(as.character(IceOutDate),form = "%Y%m%d"))
caze.ice$year = as.numeric(format(caze.ice$FirstFreezeDate,'%Y'))
caze.ice$cover_dur = with(caze.ice, IceOutDate-FirstFreezeDate); # difference in days

caze.ice$Ice_On.WY=WY(caze.ice$FirstFreezeDate,"Fed"); # "Water Year" - starting in October
caze.ice$Ice_Off.WY=WY(caze.ice$IceOutDate,"Fed"); # "Water Year" - starting in October

caze.ice$Ice_On.DOWY=hydro.day(caze.ice$FirstFreezeDate,"Fed"); # day of water year (based on federal water year)
caze.ice$Ice_Off.DOWY=hydro.day(caze.ice$IceOutDate,"Fed"); # day of water year (based on federal water year)
caze.ice$DOWY.duration=with(caze.ice,Ice_Off.DOWY-Ice_On.DOWY); # another difference calc (difference in days)

with(caze.ice,cor.test(year,DOWY.duration,method="kendall")); # quick trend test 

library(mblm)
mblm(DOWY.duration~year,subset(caze.ice,is.na(DOWY.duration)==F)); # Thiel-Sen slope estimator

# plot.path = full path on PC
# png(filename=paste0(plot.path,"caze_ice.png"),width=5,height=6,units="in",res=200,type="windows",bg="white")
par(family="serif",mar=c(2,2,0.25,1),oma=c(1,2.5,0.75,0.25));
layout(matrix(1:2,2,1))

xlim.val=c(1830,2030);by.x=50;xmaj=seq(xlim.val[1],xlim.val[2],by.x);xmin=seq(xlim.val[1],xlim.val[2],by.x/5)
ylim.val.date=date.fun(c("2023-10-01","2024-07-01"));ymaj.date=seq(ylim.val.date[1],ylim.val.date[2],"3 months");ymin.date=seq(ylim.val.date[1],ylim.val.date[2],"1 months")
ylim.val=hydro.day(ylim.val.date,"Fed");ymaj=hydro.day(ymaj.date,"Fed");ymin=hydro.day(ymin.date,"Fed")

plot(Ice_On.DOWY~year,caze.ice,ylim=ylim.val,xlim=xlim.val,type="n",axes=F,ann=F)
abline(h=ymaj,v=xmaj,lty=3,col="grey",lwd=0.5)
with(caze.ice,segments(year,Ice_Off.DOWY,year,Ice_On.DOWY),lty=2,lwd=1.5)
points(Ice_On.DOWY~year,caze.ice,pch=21,bg="white",lwd=0.01,cex=0.75)
points(Ice_Off.DOWY~year,caze.ice,pch=21,bg="dodgerblue1",lwd=0.01,cex=0.75)
axis_fun(2,ymaj,ymin,format(ymaj.date,"%b"))
axis_fun(1,xmaj,xmin,xmaj,line=-0.5);box(lwd=1)
mtext(side=2,line=2.5,'Month')
legend("topleft",legend=c("Ice-Free","Ice-On"),
       pch=c(21),pt.bg=c("dodgerblue1","white"),
       lty=c(0),lwd=c(0.1),col=c("black"),
       pt.cex=1.25,cex=0.75,ncol=1,bty="n",y.intersp=1.25,x.intersp=0.75,xpd=NA,xjust=0.5,yjust=0.5)

mtext(side=3,adj=0,"Lake Cazenovia, New York")

ylim.val=c(0,180);by.y=60;ymaj=seq(ylim.val[1],ylim.val[2],by.y);ymin=seq(ylim.val[1],ylim.val[2],by.y/2)
plot(DOWY.duration~year,caze.ice,ylim=ylim.val,xlim=xlim.val,type="n",axes=F,ann=F)
abline(h=ymaj,v=xmaj,lty=3,col="grey",lwd=0.5)
with(caze.ice,pt_line(year,DOWY.duration,2,"grey50",1,21,"lightblue",pt.lwd=0.01,cex=1))
axis_fun(2,ymaj,ymin,ymaj)
axis_fun(1,xmaj,xmin,xmaj,line=-0.5);box(lwd=1)
mtext(side=2,line=2.5,'Ice Cover Duration (Days)')
mtext(side=1,line=1.75,'Water Year (Oct - Sept)')

dev.off()