baytrends033: Long Term Water Quality Trend Analysis

Documented in gamPlotDisp

# ####
#' Plot censored gam fits vs. time
#'
#' @param gamResult output from procedure gamTest
#' @param analySpec analytical specifications
#' @param fullModel GAM # for displaying full GAM (e.g., 0, 1, 2)
#' @param seasAvgModel GAM # for displaying seasonally average GAM
#' @param seasonalModel GAM # for displaying seasonal GAM
#' @param diffType plot predicted baseline mean ('regular') or adjusted baseline mean ('adjusted')
#' @param obserPlot logical field indicating whether to plot observations
#' @param interventionPlot logical field indicating whether to plot interventions (e.g., method changes)
#' @param seasAvgPlot logical field indicating whether to plot seasonal average GAM
#' @param seasAvgConfIntPlot logical field indicating whether to plot confidence interval for seasonal average GAM
#' @param seasAvgSigPlot logical field indicating whether to plot significant increasing and decreasing trends for seasonal average GAM
#' @param fullModelPlot logical field indicating whether to plot full GAM
#' @param seasModelPlot logical field indicating whether to plot seasonal GAM
#' @param BaseCurrentMeanPlot logical field indicating whether to plot baseline and current mean
#' @param adjustedPlot logical field indicating whether to plot adjusted model
#' @export
#' @importFrom stats quantile
#'
# ####
gamPlotDisp <- function(gamResult=gamResult, analySpec=analySpec,
                         fullModel=2, seasAvgModel=2, seasonalModel=2, diffType='regular',
                         obserPlot=TRUE, interventionPlot=TRUE,
                         seasAvgPlot=TRUE, seasAvgConfIntPlot=TRUE, seasAvgSigPlot=TRUE,
                         fullModelPlot=TRUE, seasModelPlot=TRUE, BaseCurrentMeanPlot=TRUE,
                         adjustedPlot=FALSE) {

# ----- Change history --------------------------------------------
# 12Jul2018: JBH: add magenta line when additional independent variables are in model
# 05Aug2017: JBH: added some place holder code for flw_sal bar & fitted model
# 14Mar2017: JBH: updated off scale y settings
# 05Jan2017: JBH: modified ConfInt shading to allow for missing values in pdat
#                 modified Derivative line to allow for "0's" in pdat$sa.pred1.sig
# 04Nov2016: JBH: modified to print gamDiff regular and gamDiff adjusted baseline and
#                 current means.
# 25Oct2016: JBH: added statement to reset par settings at end of function
# 20Oct2016: JBH: made adjustment to lower offscale setting calculation, yOffScale
# 17Oct2016: JBH: separated gamPlot to gamPlotCalc (calculations) and gamPlotDisp
#                 (display). This enables gamPlotDisp have a cleaner strategy for
#                 customizing plots.
# 16Jun2016: JBH: re-activated derivative code; updated pdat code to always compute
#                 predicted values for seasons and return pdat to calling function;
#                 added code to tally range of dates for significant increases/decreases
# 04Jun2016: JBH: added argument dayStep to thin plot density and reduce processing time
# 03Jun2016: JBH: Added unpack of stat and layer from iSpec
# 27Apr2016: JBH: Explicit use of "::" for non-base functions added.
# 04Feb2016: JBH: added horizontal grid lines to plots

# Unpack & initialization #####
  tsdat     <- gamResult$data
  tsdat.all <- gamResult$data.all
  iSpec     <- gamResult$iSpec
  date.range <- c(iSpec$dateBegin, iSpec$dateEnd)
  centerYear <- iSpec$centerYear
  formula    <- iSpec$gamForm
  transform  <- iSpec$transform
  stat       <- iSpec$stat
  layer      <- iSpec$layer
  dep        <- iSpec$dep
  stat.layer <- paste(stat,layer,sep="-")
  analySpec$gamLegend$On <-FALSE
  showGamNumOnPlot <- analySpec$showGamNumOnPlot

  gridCol <- 'gray80'; gridlty <- 6
  censorlty <- 3; censorlwd <- 1
  seasAvgSiglwd <- 4

# Review data points and set up x- and y-axis ranges ####
  conc <- if (iSpec$isSurv) {
    data.frame(tsdat.all$date, tsdat.all$recensor, unSurv(tsdat.all[,dep])[,1:2]) 
  } else {
    data.frame(tsdat.all$date, tsdat.all$recensor, tsdat.all[,dep], tsdat.all[,dep]) 
  }
                              
  names(conc)   <- c("date", "recensor", "lower", "upper" )
  conc$point    <- "?"

  yRange <- range(c(conc[is.finite(conc$lower),"lower"],
                    conc[is.finite(conc$upper),"upper"]))
  xRange <- range(conc$date)

  #20Oct2016: switch lower vale from "yRange[1] - abs(yRange[1])"
  #           to "yRange[1] - abs(yRange[2])"
  yOffScale <- c(yRange[1] - 1e3*abs(yRange[2]) , yRange[2] + 1e3*abs(yRange[2]))

  conc[conc$recensor,"point"]  <- "Recen"
  conc[!conc$recensor & conc$lower==conc$upper,"point"]    <- "Uncen"
  conc[!conc$recensor & is.infinite(conc$lower) & is.finite(conc$upper),"point"]     <- "LT"
  conc[!conc$recensor & is.finite(conc$lower)   & is.infinite(conc$upper),"point"]   <- "GT"
  conc[!conc$recensor & is.finite(conc$lower)   & is.finite(conc$upper) & conc$lower!=conc$upper,"point"]    <- "Int"

  conc$lower2   <- conc$lower
  conc$upper2   <- conc$upper
  conc[is.infinite(conc$lower2),"lower2"] <- yOffScale[1]
  conc[is.infinite(conc$upper2),"upper2"] <- yOffScale[2]

# Set up plot, ####
  #opar <- par()
  #par(oma = c(1.5, 0,0,0))
  #par(fig = c(0, 1, 0, 1), oma = c(1.5, 0,0,0), mar = c(5.1, 4.1, 4.1, 2.1), new = FALSE)
  par(fig = c(0, 1, 0, 1), oma = c(1.5, 0,0,0), mar = c(3.5, 4.1, 2.7, 2.1), new = FALSE)

  if(transform) {
    plot(y=exp(yRange[1]), x=xRange[1], log='y',  type='n', ylim=exp(yRange), xlim=xRange,
         cex=.6 , cex.axis=.8 ,ann=FALSE  , par(mar=c(3,3.8,3,1)))
  } else {
    plot(y=(yRange[1]), x=xRange[1], log='',  type='n', ylim=(yRange), xlim=xRange,
         cex=.6 , cex.axis=.8 ,ann=FALSE  , par(mar=c(3,3.8,3,1)))
  }

# add flow or salinity indicators #05Aug2017 ####
  if('flw_sal' %in% names(tsdat)) {
    myPCH <- 124
    myCEX <- 1.0
    HiFlw_LoSal <- "springgreen"
    LoFlw_HiSal <- "darkblue"
    y_hydro <- ifelse(transform, exp(yRange[2]), yRange[2])
    prob_hydro <- 0.80
    if(iSpec$hydroTermSel == 'flow') {
      tmp <-  data.frame(x=tsdat[tsdat$flw_sal > quantile(tsdat$flw_sal,probs=prob_hydro), "date"],
                        y=y_hydro)
      points(y=tmp$y , x=tmp$x, col=HiFlw_LoSal, pch=myPCH, cex=myCEX)
      tmp <-  data.frame(x=tsdat[tsdat$flw_sal < quantile(tsdat$flw_sal,probs=(1-prob_hydro)), "date"],
                        y=y_hydro)
      points(y=tmp$y , x=tmp$x, col=LoFlw_HiSal, pch=myPCH, cex=myCEX)

    } else if(iSpec$hydroTermSel == 'salinity') {
      tmp <-  data.frame(x=tsdat[tsdat$flw_sal < quantile(tsdat$flw_sal,probs=(1-prob_hydro)), "date"],
                        y=y_hydro)
      points(y=tmp$y , x=tmp$x, col=HiFlw_LoSal, pch=myPCH, cex=myCEX)
      tmp <-  data.frame(x=tsdat[tsdat$flw_sal > quantile(tsdat$flw_sal,probs=prob_hydro), "date"],
                        y=y_hydro)
      points(y=tmp$y , x=tmp$x, col=LoFlw_HiSal, pch=myPCH, cex=myCEX)
    }
  }

# add axes, grid lines, and labels  ####
  grid (NA, NULL, lty = gridlty, col = gridCol)
  abline(v=axis.POSIXct(1, x=pretty(tsdat$date), labels = FALSE),lty = gridlty, col = gridCol)
  mtext(side=2,text=paste0(iSpec$parmName," [",iSpec$parmUnits,"]"),line=2.2,cex=.9)
  title(paste0(iSpec$parmName,"-",iSpec$layerName," at ", iSpec$stat),cex.main=.9)
  if (showGamNumOnPlot) mtext(side=3,text=fullModel,line=0.5,cex=.8,adj=1, col='black')   #2019Jun24

# fit.GAM #05Aug2017 ; 01Oct2017 ; 24Nov2017####
  # gamOutput4 <-gamResult$gamOutput4$gamRslt$fitted.values
  # fittedValues <- TRUE
  fit.GAM <- ifelse (length(grep('flw_sal',
                                 gamResult[[paste0("gamOutput",fullModel)]]$gamRslt$formula)) == 0, FALSE, TRUE)
  
  if (!fit.GAM) {
    # identify additional independent variables #12Jul2018
    # excluding dependent variable and "gamK*" 
    indVar <- setdiff (all.vars(gamResult[[paste0("gamOutput",fullModel)]]$gamRslt$formula)
                       , c(iSpec$dep, "cyear", "doy", "intervention", "flw_sal"))
    indVar <- indVar[-grep("^gamK", indVar)]
    fit.GAM <- length(indVar) > 0
  }
  
  if(fit.GAM) {
    pdat <- data.frame(date = gamResult$data$date,
                       pred1 = gamResult[[paste0("gamOutput",fullModel)]]$gamRslt$fitted.values)
    colSelA  <- 'magenta'
    lwdSelA  <- 1
    ltySelA  <- 3

    if(transform) {
      lines(exp(pdat$pred1)~pdat$date,col=colSelA, lwd=lwdSelA,lty=ltySelA)
    } else {
      lines((pdat$pred1)~pdat$date,col=colSelA, lwd=lwdSelA,lty=ltySelA)
    }
  }

# shading for CI ####
  if(seasAvgConfIntPlot) {
    pdat <- gamResult[[paste0("gamOutput",seasAvgModel)]]$predictions

    # 05Jan2017: remove records in prediction data set with NAs
    pdat <- pdat[!is.na(pdat$cilb) & !is.na(pdat$ciub),]

    colSelA <- analySpec$gamLegend$colSel[which(analySpec$gamLegend$descrip =="seasConfInt")]

    analySpec$gamLegend$On[which(analySpec$gamLegend$descrip =="seasConfInt")] <- TRUE

    ciub <- as.matrix(pdat[, c('ndate','ciub')])
    cilb <- as.matrix(pdat[, c('ndate','cilb')])
    cilb <- cilb[order(cilb[,1],decreasing=TRUE),]
    ci.polygon <- rbind(ciub,cilb)

    # polys <- unique(pdat$intervention)
    # table(pdat$intervention)

    if(transform) {
      polygon(ci.polygon[,1], exp(ci.polygon[,2]),col=colSelA,border=NA)
    } else {
      polygon(ci.polygon[,1], (ci.polygon[,2]),col=colSelA,border=NA)
    }
  }

# add intervention dates ####
  if(interventionPlot & nrow(iSpec$intervenList)>1) {

    colSelA <- analySpec$gamLegend$colSel[which(analySpec$gamLegend$descrip =="interventionchange")]
    ltySelA <- analySpec$gamLegend$ltyLegend[which(analySpec$gamLegend$descrip =="interventionchange")]
    lwdSelA <- analySpec$gamLegend$lwdLegend[which(analySpec$gamLegend$descrip =="interventionchange")]

    analySpec$gamLegend$On[which(analySpec$gamLegend$descrip =="interventionchange")] <- TRUE

    abline(v=iSpec$intervenList[2:nrow(iSpec$intervenList),"beginDate"], col=colSelA, lty=ltySelA, lwd=lwdSelA)
  }

# plot points ####
  if(obserPlot) {

    #Recensored points
    colSelA <- analySpec$gamLegend$colSel[which(analySpec$gamLegend$descrip =="recensored")]
    pchSelA <- analySpec$gamLegend$pchLegend[which(analySpec$gamLegend$descrip =="recensored")]

    if(nrow(conc[conc$point=="Recen",])>0) {
      analySpec$gamLegend$On[which(analySpec$gamLegend$descrip =="recensored")] <- TRUE
    }

    if(transform) {
      points(conc[conc$point=="Recen","date"], exp(conc[conc$point=="Recen","upper"]), pch=pchSelA, col=colSelA, cex=0.6)
    } else {
      points(conc[conc$point=="Recen","date"],    (conc[conc$point=="Recen","upper"]), pch=pchSelA, col=colSelA, cex=0.6)
    }

    #Select and graph censored data
    colSelA <- analySpec$gamLegend$colSel[which(analySpec$gamLegend$descrip =="censored")]
    pchSelA <- analySpec$gamLegend$pchLegend[which(analySpec$gamLegend$descrip =="censored")]

    # Subset censored data to tmp data frame
    tmp <- conc [conc$point=="LT" | conc$point=="GT" | conc$point=="Int",]

    if(nrow(tmp)>0) {
      analySpec$gamLegend$On[which(analySpec$gamLegend$descrip =="censored")] <- TRUE
    }

    # do plot function if at least one record in tmp
    if (nrow(tmp)>0) {
      # loop thru each point and plot short line
      for (iRow in 1:nrow(tmp)) {
        x  <- tmp[iRow,"date"]
        y1 <- tmp[iRow,"lower2"]
        y2 <- tmp[iRow,"upper2"]
        if(transform) {
          lines(c(x,x), c(exp(y1),exp(y2)), type="l", lwd=censorlwd, col=colSelA, lty=censorlty)
        } else {
          lines(c(x,x), c(   (y1),   (y2)), type="l", lwd=censorlwd, col=colSelA, lty=censorlty)
        }
      }
      # plot symbols at lower and upper limit
      if(transform) {
        lines(tmp[,"date"],     exp(tmp[,"lower2"])         , type="p", pch=pchSelA,  col=colSelA, cex=0.6)
        lines(tmp[,"date"],     exp(tmp[,"upper2"])         , type="p", pch=pchSelA,  col=colSelA, cex=0.6)
      } else {
        lines(tmp[,"date"],        (tmp[,"lower2"])         , type="p", pch=pchSelA,  col=colSelA, cex=0.6)
        lines(tmp[,"date"],        (tmp[,"upper2"])         , type="p", pch=pchSelA,  col=colSelA, cex=0.6)
      }
    }

    # Uncensored
    colSelA <- analySpec$gamLegend$colSel[which(analySpec$gamLegend$descrip =="uncensored")]
    pchSelA <- analySpec$gamLegend$pchLegend[which(analySpec$gamLegend$descrip =="uncensored")]

    if(nrow(conc[conc$point=="Uncen",])>0) {
      analySpec$gamLegend$On[which(analySpec$gamLegend$descrip =="uncensored")] <- TRUE
    }

    if(transform) {
      points(conc[conc$point=="Uncen","date"], exp(conc[conc$point=="Uncen","upper"]), pch=pchSelA, col=colSelA, cex=0.6)
    } else {
      points(conc[conc$point=="Uncen","date"],    (conc[conc$point=="Uncen","upper"]), pch=pchSelA, col=colSelA, cex=0.6)
    }
  }

# Derivative processing #####
  if(seasAvgSigPlot) {
    pdat <- gamResult[[paste0("gamOutput",seasAvgModel)]]$predictions
    colSelA <- analySpec$gamLegend$colSel[which(analySpec$gamLegend$descrip =="seasSignif")]

    analySpec$gamLegend$On[which(analySpec$gamLegend$descrip =="seasSignif")] <- TRUE

    # 05Jan2017:
    pdat[!is.na(pdat$sa.pred1.sig) & pdat$sa.pred1.sig==0, "sa.pred1.sig"] <- NA

    # derivative plotting
    if(transform) {
      lines(pdat[pdat$sa.pred1.sig==-1,"date"], exp(pdat[pdat$sa.pred1.sig==-1,"sa.pred1"]),  lwd = seasAvgSiglwd, col = colSelA)
      lines(pdat[pdat$sa.pred1.sig== 1,"date"], exp(pdat[pdat$sa.pred1.sig== 1,"sa.pred1"]),  lwd = seasAvgSiglwd, col = colSelA)
    } else {
      lines(pdat[pdat$sa.pred1.sig==-1,"date"],    (pdat[pdat$sa.pred1.sig==-1,"sa.pred1"]),  lwd = seasAvgSiglwd, col = colSelA)
      lines(pdat[pdat$sa.pred1.sig== 1,"date"],    (pdat[pdat$sa.pred1.sig== 1,"sa.pred1"]),  lwd = seasAvgSiglwd, col = colSelA)
    }
  }

# fitted GAM for seasonally-adjusted GAM  ####
  if(seasAvgPlot) {
    pdat <- gamResult[[paste0("gamOutput",seasAvgModel)]]$predictions
    colSelA  <- analySpec$gamLegend$colSel[which(analySpec$gamLegend$descrip =="seasAverage")]
    lwdSelA  <- analySpec$gamLegend$lwdLegend[which(analySpec$gamLegend$descrip =="seasAverage")]
    ltySelA  <- analySpec$gamLegend$ltyLegend[which(analySpec$gamLegend$descrip =="seasAverage")]

    analySpec$gamLegend$On[which(analySpec$gamLegend$descrip =="seasAverage")] <- TRUE

    if(transform) {
      lines(exp(pdat$sa.pred1)~pdat$date,col=colSelA, lwd=lwdSelA, lty=ltySelA)
      if(adjustedPlot & "sa.pred1.adjusted" %in% names(pdat)) {
        lines(exp(pdat$sa.pred1.adjusted)~pdat$date,col=colSelA, lwd=lwdSelA, lty=ltySelA)
      }
    } else {
      lines((pdat$sa.pred1)~pdat$date,col=colSelA, lwd=lwdSelA, lty=ltySelA)
      if(adjustedPlot & "sa.pred1.adjusted" %in% names(pdat)) {
        lines((pdat$sa.pred1.adjusted)~pdat$date,col=colSelA, lwd=lwdSelA, lty=ltySelA)
      }
    }
  }

# fitted GAM for full model GAM  ####
  if(fullModelPlot) {
    pdat <- gamResult[[paste0("gamOutput",fullModel)]]$predictions
    colSelA  <- analySpec$gamLegend$colSel[which(analySpec$gamLegend$descrip =="fullGAM")]
    lwdSelA  <- analySpec$gamLegend$lwdLegend[which(analySpec$gamLegend$descrip =="fullGAM")]
    ltySelA  <- analySpec$gamLegend$ltyLegend[which(analySpec$gamLegend$descrip =="fullGAM")]

    analySpec$gamLegend$On[which(analySpec$gamLegend$descrip =="fullGAM")] <- TRUE

    if(transform) {
      lines(exp(pdat$pred1)~pdat$date,col=colSelA, lwd=lwdSelA,lty=ltySelA)
      if(adjustedPlot & "pred1.adjusted" %in% names(pdat)) {
        lines(exp(pdat$pred1.adjusted)~pdat$date,col=colSelA, lwd=lwdSelA,lty=ltySelA)
      }
    } else {
      lines((pdat$pred1)~pdat$date,col=colSelA, lwd=lwdSelA,lty=ltySelA)
      if(adjustedPlot & "pred1.adjusted" %in% names(pdat)) {
        lines((pdat$pred1.adjusted)~pdat$date,col=colSelA, lwd=lwdSelA,lty=ltySelA)
      }
    }
  }

# fitted GAM for seasonal models  ####
  if(seasModelPlot) {
    pdat <- gamResult[[paste0("gamOutput",seasonalModel)]]$predictions

    for (i in 1:length(analySpec$gamLegend[analySpec$gamLegend$season, "legend"])) {
      descrip  <- paste0("season",i)
      colSelA  <- analySpec$gamLegend$colSel[which(analySpec$gamLegend$descrip ==descrip)]
      lwdSelA  <- analySpec$gamLegend$lwdLegend[which(analySpec$gamLegend$descrip ==descrip)]
      ltySelA  <- analySpec$gamLegend$ltyLegend[which(analySpec$gamLegend$descrip ==descrip)]

      analySpec$gamLegend$On[which(analySpec$gamLegend$descrip ==descrip)] <- TRUE

      if(transform) {
        lines(exp(pdat[, paste0("seas.",i)])~pdat$date,col=colSelA,lty=ltySelA,lwd=lwdSelA)
        if(adjustedPlot & paste0("seas.",i,".adjusted") %in% names(pdat)) {
          lines(exp(pdat[, paste0("seas.",i,".adjusted")])~pdat$date,col=colSelA,lty=ltySelA,lwd=lwdSelA)
        }
      } else {
        lines(   (pdat[, paste0("seas.",i)])~pdat$date,col=colSelA,lty=ltySelA,lwd=lwdSelA)
        if(adjustedPlot & paste0("seas.",i,".adjusted") %in% names(pdat)) {
          lines(   (pdat[, paste0("seas.",i,".adjusted")])~pdat$date,col=colSelA,lty=ltySelA,lwd=lwdSelA)
        }
      }
    }
  }

# Baseline and current mean plot ####
  if(BaseCurrentMeanPlot) {

    if(diffType=='regular' | diffType=='both') {
      porDiff <- gamResult[[paste0("gamOutput",fullModel)]]$porDiff.regular
      colSelA <- analySpec$gamLegend$colSel[which(analySpec$gamLegend$descrip =="BaseCurrentMean")]
      pchSelA <- analySpec$gamLegend$pchLegend[which(analySpec$gamLegend$descrip =="BaseCurrentMean")]

      analySpec$gamLegend$On[which(analySpec$gamLegend$descrip =="BaseCurrentMean")] <- TRUE

      por <-  date.range + c(365*24*60*60,-365*24*60*60)
      if(transform) {
        points(exp(c(porDiff$per.mn))~por,pch=pchSelA,col=colSelA,cex=1.5)
      } else {
        points(   (c(porDiff$per.mn))~por,pch=pchSelA,col=colSelA,cex=1.5)
      }
    }

    if(diffType=='adjusted' | diffType=='both') {
      porDiff <- gamResult[[paste0("gamOutput",fullModel)]]$porDiff.adjusted
      colSelA <- analySpec$gamLegend$colSel[which(analySpec$gamLegend$descrip =="BaseCurrentMeanAdj")]
      pchSelA <- analySpec$gamLegend$pchLegend[which(analySpec$gamLegend$descrip =="BaseCurrentMeanAdj")]

      analySpec$gamLegend$On[which(analySpec$gamLegend$descrip =="BaseCurrentMeanAdj")] <- TRUE

      por <-  date.range + c(365*24*60*60,-365*24*60*60)
      if(transform) {
        points(exp(c(porDiff$per.mn))~por,pch=pchSelA,col=colSelA,cex=1.0)
      } else {
        points(   (c(porDiff$per.mn))~por,pch=pchSelA,col=colSelA,cex=1.0)
      }
    }

  }

# Legend ####

  myLegend <- analySpec$gamLegend[analySpec$gamLegend$On, ]

  # update legend for flow or salinity #05Aug2017
  if('flw_sal' %in% names(tsdat)) {
    if(iSpec$hydroTermSel == 'flow') {
      myLegend <- rbind(myLegend,
                        data.frame(legend     = 'Hi Flw',      colSel       = HiFlw_LoSal,    colLegend = HiFlw_LoSal,
                                   lwdLegend  = NA_real_,      ltyLegend    = NA_real_,      pchLegend = 15,
                                   fillLegend = NA_character_, borderLegend = NA_character_, season    = FALSE,
                                   descrip    = "hydro1",      On           = TRUE),
                        data.frame(legend     = 'Lo Flw',      colSel       = LoFlw_HiSal,    colLegend = LoFlw_HiSal,
                                   lwdLegend  = NA_real_,      ltyLegend    = NA_real_,      pchLegend = 15,
                                   fillLegend = NA_character_, borderLegend = NA_character_, season    = FALSE,
                                   descrip    = "hydro2",      On           = TRUE))


    } else if(iSpec$hydroTermSel == 'salinity') {
      myLegend <- rbind(myLegend,
                        data.frame(legend     = 'Lo Sal',      colSel       = HiFlw_LoSal,    colLegend = HiFlw_LoSal,
                                   lwdLegend  = NA_real_,      ltyLegend    = NA_real_,      pchLegend = 15,
                                   fillLegend = NA_character_, borderLegend = NA_character_, season    = FALSE,
                                   descrip    = "hydro1",      On           = TRUE),
                        data.frame(legend     = 'Hi Sal',      colSel       = LoFlw_HiSal,    colLegend = LoFlw_HiSal,
                                   lwdLegend  = NA_real_,      ltyLegend    = NA_real_,      pchLegend = 15,
                                   fillLegend = NA_character_, borderLegend = NA_character_, season    = FALSE,
                                   descrip    = "hydro2",      On           = TRUE))
    }
  }

  #update legend for fitted model
  if(fit.GAM) {
    myLegend <- rbind(myLegend,
                      data.frame(legend     = 'Fit.GAM',     colSel       = 'magenta',    colLegend = 'magenta',
                                 lwdLegend  = 1 ,            ltyLegend    = 3 ,      pchLegend = NA_real_  ,
                                 fillLegend = NA_character_, borderLegend = NA_character_, season    = FALSE,
                                 descrip    = "fit.gam",      On           = TRUE))
  } else {
    myLegend[myLegend$legend=="Adj.GAM", "legend"] <- "Fit.GAM"
  }

  ncolSel <- 1
  ncolSel <- ifelse(nrow(myLegend)>7, ceiling(nrow(myLegend)/2) , nrow(myLegend))

  par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(.5, 3, 0, 2), new = TRUE)
  plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n")
  legend(x="bottom",
         legend  = myLegend$legend,
         col     = myLegend$colLegend ,
         lwd     = myLegend$lwdLegend ,
         lty     = myLegend$ltyLegend ,
         pch     = myLegend$pchLegend ,
         fill    = myLegend$fillLegend ,
         border  = myLegend$borderLegend , ncol=ncolSel,
         cex=.65, xpd = TRUE, horiz = FALSE,inset = c(0,0), x.intersp=0.4)

# reset graphical parameter settings ####

  # 25Oct2016
  par(fig = c(0, 1, 0, 1), oma= c(0, 0, 0, 0), mar=c(5.1, 4.1, 4.1, 2.1))


}