seasonal plot with trend line

seaPlot.R

# Load Library
library(mblm); # For trend line

# Commented out - added function in script
# devtools::install_github("SwampThingPaul/AnalystHelper")
# library(AnalystHelper); # some plotting helper functions like axis_fun(...), pt_lines(...), etc

## axis function from AnalystHelper package
axis_fun=function (side, at, at2, labels, cex.axis = 1, line = -0.25, 
                   lwd = 1, maj.tcl = -0.6, min.tcl = -0.3, las = 1, axisLine = 0, 
                   ...) 
{
  axis(side, line = line, at = at, labels = labels, las = las, 
       tcl = maj.tcl, lty = 0, cex.axis = cex.axis, ...)
  axis(side, at = at, labels = F, las = las, tcl = maj.tcl, 
       lwd = lwd, line = axisLine)
  axis(side, at = at2, labels = F, tcl = min.tcl, lwd = lwd, 
       line = axisLine)
}

## Some Data
dat=data.frame(CY=c(2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2008, 2008, 
                    2008, 2008, 2008, 2008, 2008, 2008, 2008, 2008, 2008, 2008, 2009, 
                    2009, 2009, 2009, 2009, 2009, 2009, 2009, 2009, 2009, 2009, 2009, 
                    2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 2010, 
                    2010, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 2011, 
                    2011, 2011, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 2012, 
                    2012, 2012, 2012, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 2013, 
                    2013, 2013, 2013, 2013, 2014, 2014, 2014, 2014, 2014, 2014, 2014, 
                    2014, 2014, 2014, 2014, 2014, 2015, 2015, 2015, 2015, 2015, 2015, 
                    2015, 2015, 2015, 2015, 2015, 2015, 2016, 2016, 2016, 2016, 2016, 
                    2016, 2016, 2016, 2016, 2016, 2016, 2016, 2017, 2017, 2017, 2017, 
                    2017, 2017, 2017, 2017, 2017, 2017, 2017, 2017, 2018, 2018, 2018, 
                    2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2018, 2019, 2019, 
                    2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019, 2020, 
                    2020, 2020, 2020, 2020, 2020, 2020, 2020, 2020, 2020, 2020, 2020, 
                    2021, 2021, 2021, 2021),
               month=c(5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 
                       12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 
                       7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 
                       2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 
                       10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 
                       5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 
                       12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 
                       7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 
                       2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4),
               TRF.cm=c(6.3, 17, 23.5, 9.7, 15.8, 9.5, 0.9, 5.8, 2.1, 5.9, 9.3, 6.8, 
                        6.8, 26, 18.2, 35.7, 10.4, 7, 1.5, 2.4, 1.2, 1, 2.1, 1.3, 13.4, 
                        17.1, 13.5, 14.2, 11.5, 1.7, 1.8, 7.8, 5.4, 5, 20.2, 18.5, 16, 
                        10.5, 16, 18, 11.2, 0.3, 3.2, 1.6, 5.3, 0.8, 7.6, 2.2, 3.1, 11.8, 
                        13.5, 18.9, 17.4, 18.8, 6.3, 3.5, 1.4, 4.3, 10.5, 5.2, 9.6, 19.6, 
                        13.1, 28.4, 13.7, 9.2, 0.5, 3.4, 1.7, 4.8, 2.9, 8.2, 14.2, 25, 
                        24.9, 16.5, 16.3, 6.2, 2.1, 0.9, 6.7, 5.4, 5.7, 4.7, 10.1, 22.2, 
                        18.7, 12.8, 17.6, 8.9, 6.2, 1.5, 2.5, 6.8, 2.5, 11.8, 4, 14.2, 
                        14.5, 17.1, 15.9, 3.7, 8.1, 10.8, 23.5, 7.3, 6, 2.4, 17.7, 21.4, 
                        10.9, 19.2, 18.9, 6.2, 0.3, 1.4, 2.3, 1.7, 1.4, 2.9, 9.2, 29.9, 
                        14.8, 13.6, 28.8, 19, 3.6, 3.4, 2.8, 1.9, 1.6, 7.7, 26.3, 16.6, 
                        14.5, 12.6, 10, 2.6, 5.5, 4.3, 9.5, 3.8, 3.6, 4.1, 9.7, 19.6, 
                        16.6, 22.5, 4.5, 8.4, 3.2, 10.1, 3, 5.3, 0.4, 9.4, 15.9, 17.4, 
                        18.7, 21.4, 18.7, 14.9, 9.4, 4.3, 0.5, 6.6, 1.8, 9.5))

## Set up x and y limits and tick marks
ylim.val=c(0,40);by.y=10;ymaj=seq(ylim.val[1],ylim.val[2],by.y);ymin=seq(ylim.val[1],ylim.val[2],by.y/2)
xlim.val=c(2008,2021);by.x=8;xmaj=seq(xlim.val[1],xlim.val[2],by.x);xmin=seq(xlim.val[1],xlim.val[2],by.x/by.x)

# Base plotting setup
# To write to a png
png(filename="S79MonthQ.png",width=8.5,height=2.5,units="in",res=200,type="windows",bg="white")

# sets up spacing and layout
par(family="serif",mar=c(0.5,0.75,0.25,0.1),oma=c(2.5,3.5,1,2));
layout(matrix(1:12,1,12,byrow=T))

# loop through each month to plot lines, points and trends (with 95% CI)
for(i in 1:12){
  tmp=subset(dat,month==i)
  plot(TRF.cm~CY,tmp,ylim=ylim.val,xlim=xlim.val,type="n",axes=F,ann=F)
  abline(h=ymaj,v=xmaj,lty=3,col=adjustcolor("grey",0.5))
  lines(TRF.cm~CY,tmp,lty=2,col="dodgerblue1")
  points(TRF.cm~CY,tmp,pch=21,bg="dodgerblue1",col="grey",lwd=0.1)
  mod=mblm::mblm(TRF.cm~CY,na.omit(tmp[,c("TRF.cm","CY")]),repeated = F)
  x.val=seq(min(tmp$CY),max(tmp$CY),length.out=20)
  mod.pred=predict(mod,data.frame(CY=x.val),interval="confidence")
  lines(x.val,mod.pred[,1],lwd=1,col="red")
  lines(x.val,mod.pred[,2],lwd=0.5,col="red",lty=2)
  lines(x.val,mod.pred[,3],lwd=0.5,col="red",lty=2)
  if(i==1){axis_fun(2,ymaj,ymin,format(ymaj))}else{axis_fun(2,ymaj,ymin,NA)}
  axis_fun(1,xmaj,xmin,xmaj,line=-0.75,cex=0.75)
  box(lwd=1)
  if(i==1){mtext(side=2,line=2.5,"Rainfall (cm month\u207B\u00B9)")}
  mtext(side=3,adj=0,month.abb[i],cex=0.75,col="grey50")
}
mtext(side=1,line=1,outer=T,"Calendar Year")

dev.off()