Created
September 12, 2024 15:08
-
-
Save SwampThingPaul/539220ae2d25354eefe820d3d8a6fc4e to your computer and use it in GitHub Desktop.
Ice phenology plots
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # 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() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment