R/vis.trends.R

In metno/esd: Climate analysis and empirical-statistical downscaling (ESD) package for monthly and daily data

#' Visualise trends for multiple overlapping periods
#' 
#' Produce a plot showing trends for multiple periods within a time series. The
#' strength of the trend is represented by the color scale and significant
#' trends are marked with black borders.
#'
#' @importFrom stats anova lm
#' @importFrom grDevices colorRampPalette
#' 
#' @param x the 'x' argument provides the time series for which the trend
#' analysis is performed. Only zoo objects are accepted.
#' @param minlen minimum time interval to calculate trends for in units of
#' years.
#' @param unitlabel unit of x.
#' @param varlabel name of x.
#' @param vmax upper limit of trend scale.
#' @param show.significance TRUE to mark statistically significant trends.
#' @param pmax maximum p-value of trends marked as significant.
#' @param is spatial index for subsetting station data
#' @param lwd width of lines
#' @param verbose TRUE or FALSE.
#' @param new if TRUE plot in new window
#' 
#' @author Kajsa Parding
#' @keywords trend
#' @examples
#' 
#' t <- seq(as.Date("1955-01-01"),as.Date("2004-12-31"),by=1)
#' x <- zoo(sample(seq(-30,30,1e-1),length(t),rep=TRUE),order.by=t)
#' vis.trends(x, show.significance=FALSE, new=FALSE)
#' 
#' data(Oslo)
#' vis.trends(Oslo, unitlabel="oC", varlabel = "Temperature",
#'   pmax = 1e-2, minlen = 40, new=FALSE)
#' vis.trends(subset(Oslo,it='jja'), unitlabel="oC",
#'   varlabel = "Temperature JJA",
#'   pmax = 1e-3, vmax=0.5, minlen = 40, new=FALSE)
#' vis.trends(subset(Oslo,it='mam'), unitlabel="oC",
#'   varlabel = "Temperature MAM",
#'   pmax = 1e-3, vmax=0.5, minlen = 40, new=FALSE)
#' 
#' @export vis.trends
vis.trends <- function(x,...,unitlabel=NULL,varlabel=NULL,is=1,
                       pmax=0.01,minlen=15,lwd=NA,vmax=NA,new=TRUE,
                       show.significance=TRUE,verbose=FALSE) {
  if(verbose) print("vis.trends")
  if(is.null(unitlabel)) unitlabel <- attr(x,"unit")
  if(is.null(varlabel)) varlabel <- attr(x,"variable")
  Tr <- calculate.trends(x,minlen=minlen,is=is,verbose=verbose)
  trends <- Tr$trends*10
  p <- Tr$p
  cols <- as.numeric(colnames(trends))
  rows <- as.numeric(rownames(trends))
  significant <- ifelse(p<pmax,trends,NA)
  
  ticks <- seq(1,length(cols),signif(length(cols)/10,1))
  if (is.na(lwd)) lwd <- max(3-0.05*length(cols),0.2)
  
  if (is.na(vmax) | vmax=="q995") vmax <- q995(abs(trends))
  if (vmax=="max") vmax <- max(abs(trends),na.rm=T)
  if (vmax<1) vmax <- signif(vmax,1)
  if (vmax>1) vmax <- signif(vmax,2)
  dv <- signif(vmax/8,1)
  v0 <- 0#signif(dv/2,1)
  vstep <- seq(v0,vmax,dv)
  vstep <- unique(c(-1*vstep,vstep))
  vstep <- vstep[order(vstep)]
  cticks <- vstep[2:length(vstep)]-dv/2
  
  #cstep <- colscal(n=length(vstep)-1,pal="t2m")
  cmin <- rgb(239,138,98,maxColorValue=255) # blue
  cmid <- rgb(247,247,247,maxColorValue=255) # white
  cmax <- rgb(103,169,207,maxColorValue=255) # red
  rgb.palette <- colorRampPalette(c(cmax,cmid,cmin),space="rgb")
  cstep <- rgb.palette(n=length(vstep)-1)
  # Plot trend as color
  if (new) dev.new()
  image(cols,rows,t(trends),breaks=vstep,col=cstep,
        xlab='start year',ylab='length of period (years)',
        main=paste(c(varlabel," trend (",unitlabel,"/decade)"),collapse=""))
  
  trends.plus <- t(trends)
  trends.plus[trends.plus<max(vstep)] <- NA
  image(cols,rows,trends.plus,col=cstep[length(cstep)],add=TRUE)
  trends.minus <- t(trends)
  trends.minus[trends.minus>min(vstep)] <- NA
  image(cols,rows,trends.minus,col=cstep[1],add=TRUE)
  
  # Mark significant trends with dark borders
  if(show.significance) {
    if(verbose) print(paste("mark significant trends (p<",pmax,")",sep=""))
    i <- which((is.finite(t(p)) & t(p)<pmax))
    x <- array(sapply(cols,function(x) rep(x,nrow(p))),length(p))[i]
    y <- rep(rows,nrow(p))[i]
    matlines(rbind(x-1/2,x+1/2),rbind(y-1/2,y-1/2),col='black',lwd=lwd,lty=1)
    matlines(rbind(x-1/2,x+1/2),rbind(y+1/2,y+1/2),col='black',lwd=lwd,lty=1)
    matlines(rbind(x-1/2,x-1/2),rbind(y-1/2,y+1/2),col='black',lwd=lwd,lty=1)
    matlines(rbind(x+1/2,x+1/2),rbind(y-1/2,y+1/2),col='black',lwd=lwd,lty=1)
  }
  colbar(cticks,cstep,fig=c(0.8,0.85,0.65,0.85))
}

calculate.trends <- function(x,minlen=15,is=1,verbose=FALSE){
  # Calculate trends of time series x
  if(verbose) print("calculate.trends - calculate trends for all subperiods")
  stopifnot(inherits(x,'zoo'))
  if(!is.null(dim(x))) {
    if(ncol(x)>1) {
      x <- subset(x,is=is)
      if(!is.null(dim(x))) {
        x <- apply(x,2,mean,na.rm=TRUE)
      }
    }
    dim(x) <- length(x)
  }

  if(!inherits(x,c("annual","season"))) {
    xm <- aggregate(x,by=as.yearmon(index(x)),FUN="mean",na.rm=TRUE)
    xy <- aggregate(xm,by=year(xm),FUN="mean",na.rm=TRUE)
    ny <- aggregate(xm,by=year(xm),FUN="nv")
    xy <- xy[(ny==max(ny))] # exclude years with missing months
  } else xy <- x
  year <- as.numeric(index(xy))
  firstyear <- min(year):(max(year)-minlen+1)
  trendlen <- minlen:(diff(range(year))+1)
  #lastyear <- firstyear+minlen-1
  n <- length(firstyear)
  trends <- matrix(NA,n,n)
  rownames(trends) <- trendlen#firstyear
  colnames(trends) <- firstyear#lastyear
  p <- trends
  # speed up with apply?
  for (i in firstyear) {
    jvec <- i+trendlen[trendlen<=(max(year)-i+1)]-1#(i+minlen-1):(max(year)+1)
    for (j in jvec) {
      if(all(c(i,j) %in% year)) if(!is.na(xy[year==i]) & !is.na(xy[year==j])) {
        #if(verbose) print(paste(i,j))
        ij <- which(year %in% i:j & !is.na(xy))
        ij.model <- lm(xy[ij]~year[ij])
        #ij.kendall <- Kendall(x[ij],year[ij])
        iout <- firstyear==i#as.numeric(colnames(trends))==j
        jout <- trendlen==(j-i+1)#as.numeric(rownames(trends))==i
        trends[jout,iout] <- ij.model$coefficients[2]
        p[jout,iout] <- anova(ij.model)$Pr[1]#ij.kendall$sl[1]
      }
    }
  }  
  return(list("trends"=trends,"p"=p))
}