R/diagnose.R

Defines functions diagnose.matrix diagnose.dsensemble.list diagnose.dsensemble diagnose.distr diagnose.ds.pca diagnose.ds diagnose.cca diagnose.mvr diagnose.station diagnose.cca diagnose.comb.eof diagnose.eof diagnose.comb diagnose.default diagnose

Documented in diagnose diagnose.cca diagnose.comb diagnose.comb.eof diagnose.distr diagnose.ds diagnose.dsensemble diagnose.ds.pca diagnose.eof diagnose.matrix diagnose.mvr diagnose.station

#' Diagnose
#' 
#' Diagnose and examine combined fields, MVR, and CCA results. applies some
#' tests to check for consistency.
#' 
#' The method \code{diagnose.comb.eof} which estimates the difference in the
#' mean for the PCs of the calibration data and GCMs over a common period in
#' addition to the ratio of standard deviations and lag-one autocorrelation.
#' The x-axis shows the difference in the mean of the segments in the PCs
#' representing the different data souces, the y-axis shows difference in
#' standard deviation and the size of the symbols the difference in the
#' autocorrelation (open symbols if the autocorrelation have different signs).
#' 
#' \code{climvar} estimates the climatological variance, e.g. how the
#' inter-annual variance varies with seasons.
#'
#' @importFrom stats spectrum qqline na.omit
#' @importFrom graphics matplot arrows
#' 
#' @aliases diagnose diagnose.comb diagnose.eof
#' diagnose.comb.eof diagnose.mvr diagnose.cca diagnose.ds diagnose.station
#' diagnose.distr diagnose.matrix diagnose.dsensemble diagnose.ds.pca
#' @param x data object
#' @param it teporal selection - see \code{\link{subset}}
#' @param plot if TRUE, plot results
#' @param plot.type type of plot
#' @param verbose Logical value defaulting to FALSE. If FALSE, do not display
#' comments (silent mode). If TRUE, displays extra information on progress.
#' @param new if TRUE plot in new window
#' @param xlim range of x-axis
#' @param xlim range of y-axis
#' @param alpha factor modifying the opacity alpha; typically in [0,1]
#' @param map.show if TRUE show map
#' @param xrange longitude range to display in map
#' @param yrange latitude range to display in map
#' @param main main label in plot
#' @param sub smaller label (subtitle) in plot
#' @param xlab label of x-axis
#' @param ylab label of y-axis
#' @param probs quantile to display in plot, e.g., probs=0.95 gives a diagnosis
#' of the 95th percentile of the data.
#' @param \dots additional arguments
#'
#' @return A 'diag' object containing test results
#' @author R.E. Benestad
#' @keywords utilities
#' @examples
#' 
#' t2m <- t2m.DNMI(lon=c(-40,40),lat=c(30,70))
#' T2m <- t2m.NorESM.M(lon=c(-40,40),lat=c(30,70))
#' # Combine in time to compute common EOFs:
#' X <- combine(t2m,T2m)
#' diagnose(X)
#' 
#' ceof <- EOF(X,it="jan")
#' plot(diagnose(ceof))
#' 
#' slp <- slp.NCEP(lat=c(-40,40),anomaly=TRUE)
#' sst <- sst.NCEP(lat=c(-40,40),anomaly=TRUE)
#' eof.1 <- EOF(slp,it="jan")
#' eof.2 <- EOF(sst,it="jan")
#' cca <- CCA(eof.1,eof.2)
#' diagnose(cca)
#' 
#' 
#' @export diagnose
diagnose <-function(x,...) UseMethod("diagnose")

#' @exportS3Method
#' @export
diagnose.default <- function(x,...) {
}

#' @exportS3Method
#' @export diagnose.comb
diagnose.comb <- function(x,...) {
  n.app <- attr(x,'n.apps')
  cols <- c("black","red","blue","darkgreen","darkred","darblue",
            "grey","green","mangenta","cyan")
  par(bty="n")
  plot(colMeans(x),cex=0.5,pch=19,
       main="Grid box mean value for combined fields")
  for (i in 1:n.app) {
    col <- cols[i%%10+1]
    y <- eval(parse(text=paste("attr(x,'appendix.",i,"')",sep="")))
    points(colMeans(y,na.rm=TRUE),col=col,cex=0.3)
  }
}

#' @exportS3Method
#' @export diagnose.eof
diagnose.eof <- function(x,...) {
  if (inherits(x,'comb')) y <- diagnose.comb.eof(x,...) else
    y <- x
  return(y)
}

#' @exportS3Method
#' @export diagnose.comb.eof
diagnose.comb.eof <- function(x,...,verbose=FALSE) {
  if (verbose) print("diagnose.comb.eof")
  ACF <- function(x) acf(x,plot=FALSE,na.action=na.omit)$acf[2]
  sign <- function(x,y) {z<-x*y; z[z<0] <- -1; z[z>0] <- 1; z}
  
  stopifnot(!missing(x), inherits(x,"eof"),inherits(x,"comb"))
  #print("diagnose.comb.eof")
  
  # The original field, e.g. reanalyses
  Y <- zoo(coredata(x),order.by=index(x))
  n <- attr(x,'n.apps')
  m <- length(attr(x,'eigenvalues'))
  dm <- rep(NA,n*m); dim(dm) <- c(n,m)
  sr <- dm; ar <- sr
  # The appended fields, e.g. GCM results
  rowname <- rep("GCM",n)
  z <- NULL
  for ( i in 1:n ) {
    eval(parse(text=paste("z <- attr(x,'appendix.",i,"')",sep="")))
    y <- zoo(coredata(z),order.by=index(z))
    
    # Extract a common period:
    X <- merge(Y,y,all=FALSE)
    Ym <- apply(coredata(X),2,mean,na.rm=TRUE)
    #print(Ym)
    #plot(Ym)
    Ys <- apply(coredata(X),2,sd,na.rm=TRUE)
    AR <- apply(coredata(X),2,ACF)
    ## KMP 2016-11-02 normalise the difference with the standard deviation
    #dm[i,] <- (Ym[1:m] - Ym[(m+1):(2*m)])
    dm[i,] <- (Ym[1:m] - Ym[(m+1):(2*m)])/Ys[1:m]
    # ratio: GCM/original
    # The problem is when the denominator is close to zero...
    #sr.test <- abs((Ys[(m+1):(2*m)] - Ys[1:m])/Ys[(m+1):(2*m)] )
    ## KMP 2016-11-02 removed abs and made the reference data the denominator
    sr.test <- (Ys[1:m] - Ys[(m+1):(2*m)] )/Ys[1:m]
    sr.test[!is.finite(sr.test)] <- 0
    sr[i,] <- sr.test
    ar[i,] <- 0.5*( 2- abs(AR[(m+1):(2*m)] - AR[1:m]) )*sign(AR[(m+1):(2*m)],AR[1:m])
    if (!is.null(attr(z,'source'))) rowname[i] <- attr(z,'source') else
      if (!is.null(attr(z,'model_id'))) rowname[i] <- attr(z,'model_id') 
  }
  rownames(dm) <- rowname
  rownames(sr) <- rowname
  rownames(ar) <- rowname
  ## KMP 19-11-2015: added abs(dm) because the sign is arbitrary
  diag <- list(mean.diff=dm,sd.ratio=sr,autocorr.ratio=ar,
               common.period=range(index(Y)),sd0=Ys,
               calibrationdata=attr(x,'source'))
  attr(diag,'variable') <- attr(x,'variable')
  attr(diag,'evaluation_period') <- paste(start(Y),end(Y),sep='-')
  #print(summary(diag))
  attr(diag,'history') <- history.stamp(x)
  class(diag) <- c("diagnose","comb","eof","list")
  invisible(diag)
}

#' @exportS3Method
#' @export diagnose.cca
diagnose.cca <- function(x,...) {
  par(bty="n")
  plot(x$r,pch=19,cex=1.5,main="Canonical correlations",
       ylim=c(-1,1),ylab="correlation",xlab="pattern number")
  lines(c(0,length(x$r)),rep(0,2),col="grey")
  grid()
}

#' @exportS3Method
#' @export diagnose.station
diagnose.station <- function(x,...,main='Data availability',
                             xlab='',ylab='station',nletters=6,
                             sub=src(x),verbose=FALSE) {
  if(verbose) print("diagnose.station")
  d <- dim(x)
  if (is.null(d)) {
    z <- diagnose.distr(x,...)
    return(z)
  }
  par(mar=c(5, 4, 4, 5),las=1,xpd=TRUE,cex.lab=0.5,cex.axis=0.5,xaxt='n')

  image(1:d[1],1:d[2],coredata(x),
        main=main,xlab=xlab,ylab=ylab,
        sub=sub,col=rev(rainbow(21)),useRaster = TRUE,...)
  par(xaxt='s')
  axis(1,at=1:d[1],labels=year(x),cex.lab=0.5,col='grey')
  axis(4,at=1:d[2],labels=substr(loc(x),1,nletters),
       cex.lab=0.5,col='grey')
  par(xpd=FALSE)
  nyrs <- length(rownames(table(year(x))))
  grid(nx=nyrs,ny=d[2])
}

diagnose.mvr <- function(x,...) {
  print("Not finished")
}

#' @exportS3Method
#' @export diagnose.cca
diagnose.cca <- function(x,...) {
  par(bty="n")
  plot(x$r,pch=19,cex=1.5,main="Canonical correlations",
       ylim=c(-1,1),ylab="correlation",xlab="pattern number")
  lines(c(0,length(x$r)),rep(0,2),col="grey")
  grid()
}

# Display cross-validation and statistics on the residual
#' @exportS3Method
#' @export diagnose.ds
diagnose.ds <- function(x,...,ip=1,plot=FALSE,verbose=FALSE,new=TRUE) {
  
  ## the attribute 'evaluation' contains cross-validation
  if (verbose) print("diagnose.ds")
  stopifnot(inherits(x,"ds"))
  if (inherits(x,"pca")) {
    diagnostics <- diagnose.ds.pca(x,plot=plot,verbose=verbose)
    return(diagnostics)
  }
  
  if (!is.null(attr(x,'evaluation'))) xval <- attr(x,'evaluation') else
    xval <- crossval(x)
  ## Check the residuals
  if (verbose) print('residuals')
  y <- zoo(x) - zoo(attr(x,'eof'))
  z <- as.original.data(x)
  anova <- summary(attr(x,'model'))
  eof <- attr(x,'eof')
  if (inherits(eof,'comb')) bias.diag <- diagnose(eof) else
    bias.diag <- NULL
  
  if (verbose) print('spectrum')
  spectrum(coredata(y)[,ip],plot=FALSE) -> s
  sp <- data.frame(y=log(s$spec),x=log(s$freq))
  ## beta is used to estimate the Hurst coefficient 
  if (length(dim(y))==0) {
    beta <- -summary(lm(y ~ x, data=sp))$coefficient[2]
    beta.error <- summary(lm(y ~ x, data=sp))$coefficient[4]
    ar1 <- acf(y,plot=FALSE)$acf[2]
  } else {beta <- NA; beta.error <- NA; ar1 <- NA}
  
  ## Auto-correlation of the residual
  ar <- acf(coredata(y)[,ip],plot=FALSE)
  
  if (plot) {
    ## Timer series of the residual
    if (new) dev.new()
    par(bty="n",mfcol=c(3,2))
    plot(xval,plot.type='single',col=c("blue","red"),
         main='cross-validation',
         sub=paste('correlation=',round(cor(xval)[2,1],2)))
    
    
    for (i in 1:dim(y)[2]) {
      plot(y[,i],main='contains a trend?',lwd=2,col='grey')
      lines(trend(y[,i]),lty=2)
    }
    
    plot(ar$lag,ar$acf,type='b',main='Residual ACF?')
    
    ## Rsidual correlated with original data?
    plot(coredata(z)[,ip],coredata(y)[,ip],
         main='Residual correlated with original data?')
    
    #sp <- spectrum(y,plot=FALSE)
    plot(s$freq,s$spec,type='l',main='Residual power-spectrum',log='xy')
    
    ## Residual normally distributed?
    qqnorm(coredata(y)[,ip],main='Residual normally distributed?')
    qqline(y)
    
    if  (!is.null(attr(x,'diagnose'))) plot(attr(x,'diagnose'))
  }
  
  ii1 <- is.element(names(xval),paste('X.PCA',ip,sep='.'))
  ii2 <- is.element(names(xval),paste('Z.PCA',ip,sep='.'))
  x.corr <- cor(xval)[ii1,ii2]
  R2 <- summary(attr(x,'model')[[ip]])$r.squared
  
  diagnostics <- list(residual=y,anova=anova,xval=xval,bias.diag=bias.diag,
                      ar1=ar1,beta=beta, H=(beta+1)/2, beta.error=beta.error,
                      x.corr = x.corr,R2=R2)
  return(diagnostics)
}

# Display cross-validation and statistics on the residual
#' @exportS3Method
#' @export diagnose.ds.pca
diagnose.ds.pca <- function(x,...,plot=FALSE,verbose=FALSE,new=TRUE) {
  
  ## the attribute 'evaluation' contains cross-validation
  if (verbose) print("diagnose.ds.pca")
  if (!is.null(attr(x,'evaluation'))) xval <- attr(x,'evaluation') else
    xval <- crossval(x)
  ## Check the residuals
  y <- as.residual(x)
  y <- as.station(y)
  z <- matchdate(as.original.data(x),y)
  anova <- summary(attr(x,'model'))
  eof <- attr(x,'eof')
  if (inherits(eof,'comb')) bias.diag <- diagnose(eof) else
    bias.diag <- NULL
  
  w <- coredata(y)
  ok <- apply(w,1,function(x) !any(is.na(x)))
  w <- w[ok,]
  spectrum(w,plot=FALSE) -> s
  sp <- data.frame(y=log(s$spec),x=log(s$freq))
  if (length(dim(y))==0) {
    beta <- -summary(lm(y ~ x, data=sp))$coefficient[2]
    beta.error <- summary(lm(y ~ x, data=sp))$coefficient[4]
    ar1 <- acf(y,plot=FALSE)$acf[2]
  } else {beta <- NA; beta.error <- NA; ar1 <- NA}
  
  if (plot) {
    ## Timer series of the residual
    if (new) dev.new()
    #par(bty="n",mfcol=c(3,2))
    plot(xval,plot.type='single',col=c("blue","red"),
         main='cross-validation',
         sub=paste('correlation=',round(cor(xval)[2,1],2)))
    
    if (new) dev.new()
    matplot(y,type="p",pch=1,main='contains a trend?')
    matplot(trend(y),type="l",add=TRUE)
    
    ## Auto-correlation of the residual
    if (new) dev.new()
    ar <- acf(w,plot=FALSE)
    matplot(ar$lag,ar$acf,type='b',main='Residual ACF?')
    
    ## Residual correlated with original data?
    if (new) dev.new()
    matplot(coredata(z),coredata(y),main='Residual correlated with original data?')
    
    #sp <- spectrum(y,plot=FALSE)
    if (new) dev.new()
    matplot(s$freq,s$spec,type='l',main='Residual power-spectrum',log='xy')
    
    ## Residual normally distributed?
    if (new) dev.new()
    qqnorm(w,main='Residual normally distributed?')
    qqline(w)
    
    if  (!is.null(attr(x,'diagnose'))) 
      plot(attr(x,'diagnose'))
  }
  
  diagnostics <- list(residual=y,anova=anova,xval=xval,bias.diag=bias.diag,
                      ar1=ar1,beta=beta, H=(beta+1)/2, beta.error=beta.error)
  return(diagnostics)
}

#' @exportS3Method
#' @export diagnose.distr
diagnose.distr <- function(x,...,main=NULL,
                           xlab='mean',ylab=expression(q[p]),
                           sub=src(x),probs=0.95,plot=TRUE,verbose=FALSE) {
  if(verbose) print("diagnose.distr")
  x0 <- x
  if (is.T(x)) {
    y <- anomaly(x)
    djf <- subset(y,it='djf')
    mam <- subset(y,it='mam')
    jja <- subset(y,it='jja')
    son <- subset(y,it='son')
    m.djf <- aggregate(djf,year,FUN='mean',na.rm=TRUE)
    q.djf <- aggregate(djf,year,FUN='quantile',probs=probs,na.rm=TRUE)
    m.mam <- aggregate(mam,year,FUN='mean',na.rm=TRUE)
    q.mam <- aggregate(mam,year,FUN='quantile',probs=probs,na.rm=TRUE)
    m.jja <- aggregate(jja,year,FUN='mean',na.rm=TRUE)
    q.jja <- aggregate(jja,year,FUN='quantile',probs=probs,na.rm=TRUE)
    m.son <- aggregate(son,year,FUN='mean',na.rm=TRUE)
    q.son <- aggregate(son,year,FUN='quantile',probs=probs,na.rm=TRUE)
    
    x <- c(coredata(m.djf),coredata(m.mam),coredata(m.jja),coredata(m.son))
    y <- c(coredata(q.djf),coredata(q.mam),coredata(q.jja),coredata(q.son))
    col <- c(rep(rgb(0.2,0.2,0.6,0.3,length(m.djf))),
             rep(rgb(0.1,0.6,0.1,0.3,length(m.mam))),
             rep(rgb(0.7,0.7,0.1,0.3,length(m.jja))),
             rep(rgb(0.7,0.5,0.5,0.3,length(m.son))))
    par(bty='n')
    r <- cor.test(x,y)
    if (is.null(main)) main <- paste(loc(x0),'T(2m): mean v.s. quantile')
    model.djf <- lm(coredata(q.djf) ~ coredata(m.djf))
    model.mam <- lm(coredata(q.mam) ~ coredata(m.mam))
    model.jja <- lm(coredata(q.jja) ~ coredata(m.jja))
    model.son <- lm(coredata(q.son) ~ coredata(m.son))
    signf <- c(rep(summary(model.djf)$coefficients[8] < 0.05,length(m.djf)),
               rep(summary(model.mam)$coefficients[8] < 0.05,length(m.mam)),
               rep(summary(model.jja)$coefficients[8] < 0.05,length(m.jja)),
               rep(summary(model.son)$coefficients[8] < 0.05,length(m.son)))
    
    if (plot) {
      plot(x,y,main=main,xlab=xlab,ylab=ylab,col=col,pch=19,
           sub=paste('Anomaly: p=',probs,'; r=',round(r$estimate,3),' [',
                     round(r$conf.int[1],3),', ',round(r$conf.int[1],3),']',sep=''))
      
      abline(model.djf,col=rgb(0.2,0.2,0.6))
      abline(model.mam,col=rgb(0.1,0.6,0.1))
      abline(model.jja,col=rgb(0.7,0.7,0.1))
      abline(model.son,col=rgb(0.7,0.5,0.5))
      grid()
      points(x[signf],y[signf])
      legend(min(x),max(y),c(paste('DJF',round(model.djf$coefficients[2],3)),
                             paste('MAM',round(model.mam$coefficients[2],3)),
                             paste('JJA',round(model.jja$coefficients[2],3)),
                             paste('SON',round(model.son$coefficients[2],3))),pch=19,
             col=c(rgb(0.2,0.2,0.6),rgb(0.1,0.6,0.1),
                   rgb(0.7,0.7,0.1),rgb(0.7,0.5,0.5)),bty='n')
      text(max(x),min(y),expression(q[p]==mu + sigma*sqrt(2)*erf^-1 *(2*p - 1)),
           pos=2,col='grey')
    }
    invisible(list(DJF=model.djf,MAM=model.mam,JJA=model.jja,SON=model.son))
  }
}

#' @exportS3Method
#' @export diagnose.dsensemble
diagnose.dsensemble <- function(x,...,plot=TRUE,type='target',xrange=NULL,
                                yrange=NULL,main=NULL,map.show=TRUE,
                                map.type="points",verbose=FALSE) {
  if (verbose) print('diagnose.dsensemble')
  # Trend-evaluation: rank
  # Counts outside 90% confidence: binomial distrib. & prob.
  stopifnot(!missing(x),inherits(x,"dsensemble"))
  
  if (inherits(x,"pca") | inherits(x,"list")) {
    diag <- diagnose.dsensemble.list(x,plot=plot,map.show=map.show,
                                     main=main,xrange=xrange,yrange=yrange,
                                     verbose=verbose,...)
    invisible(diag)
  } else {
    if (is.null(attr(x,'station'))) {
      if (verbose) print('Found no station data - premature exit')
      return()
    }
    
    z <- x
    d <- dim(z)
    t <- index(z)
    y <- attr(x,'station')
    
    ## REB 2017-11-07 Set the index to the year if the data is annual or only one per year
    if ( inherits(y,'annual') | (length(rownames(table(month(x))))==1) ){
      if (verbose) print('set the time index to years')
      index(z) <- year(z)
      index(y) <- year(y)
    }
    
    if (sum(is.finite(y))==0) {
      print('diagnose.dsensemble: problem detected - no valid station data'); print(match.call())
      return(NULL)
    }
    
    ## Use the same dates
    yz <- merge( zoo(y), zoo(z),all=FALSE )
    #plot(yz)
    
    # statistics: past trends
    ## KMP 2020-08-05: Making code more efficient
    delta <- apply(coredata(yz),2,FUN='trend.coef')
    deltaobs <- delta[1]
    deltagcm <- delta[2:length(delta)]
    robs <- round(100*sum(deltaobs < deltagcm)/d[2])
    if(verbose) {print(deltaobs); print(deltagcm); print(order(c(deltaobs,deltagcm))[1])}
    
    #  i1 <- is.element(year(y)*100 + month(y),year(z)*100 + month(z))
    #  i2 <- is.element(year(z)*100 + month(z),year(y)*100 + month(y))
    #  obs <- data.frame(y=y[i1],t=year(y)[i1])
    #  obs <- data.frame(y=c(yz[,1]),t=year(yz))
    #  if (verbose) print(summary(obs))
    #  if (sum(is.finite(obs$y))==0) {
    #    print('diagnose.dsensemble: problem detected - no valid station data'); print(match.call())
    #    return(NULL)
    #  }
    #deltaobs <- round(lm(y ~ t,data=obs)$coefficients[2]*10,4)  # deg C/decade
    #  deltagcm <- rep(NA,d[2])
    #  if (verbose) print(dim(deltagcm))
    #  for (j in 1:d[2]) {
    #    gcm <- data.frame(y=c(yz[,j+1]),t=year(yz))
    #    deltagcm[j] <- round(lm(y ~ t,data=gcm)$coefficients[2]*10,2)  # deg C/decade
    #  }
    ## REB 2016-11-07: faster and more efficient code than for-loop.
    #  deltagcm <- c(apply(coredata(yz)[,-1],2,FUN='trend.coef'))
    #  robs <- round(100*sum(deltaobs < deltagcm)/d[2])
    #  if(verbose) {print(deltaobs); print(deltagcm); print(order(c(deltaobs,deltagcm))[1])}
    # apply to extract mean and sd from the selected objects:
    mu <- apply(coredata(yz[,-1]),1,mean,na.rm=TRUE)
    si <- apply(coredata(yz[,-1]),1,sd,na.rm=TRUE)
    q05 <- qnorm(0.05,mean=mu,sd=si)
    q95 <- qnorm(0.95,mean=mu,sd=si)
    # number of points outside conf. int. (binom)
    above <- yz[,1] > q95
    below <- yz[,1] < q05
    outside <- sum(above,na.rm=TRUE) + sum(below,na.rm=TRUE)
    N <- sum(is.finite(yz[,1]),na.rm=TRUE)
    
    # add plot title
    if(is.null(main)) main <- attr(x,"variable")
    
    if (verbose) print('make list object')
    diag <- list(z=z,robs=robs,deltaobs=deltaobs,deltagcm=deltagcm,
                 outside=outside,above=above,below=below,
                 y=yz[,1],N=N,xrange=xrange,yrange=yrange,
                 mu=zoo(mu,order.by=index(x)),
                 si=zoo(si,order.by=index(x)),
                 q05=zoo(q05,order.by=index(x)),
                 q95=zoo(q95,order.by=index(x)))
    attr(diag,'history') <- history.stamp(x)
    class(diag) <- c('diagnose','dsensembles','list')
    if (plot) plot(diag,map.show=map.show,map.type=map.type,
                   main=main,verbose=verbose)
    if (verbose) print(paste('exit diagnose.dsensemble: N=',N,'obs.change=',deltaobs,
                             'simulated change in [',min(deltagcm),',',max(deltagcm),']'))
    invisible(diag)
  }
}

#' @exportS3Method
#' @export
diagnose.dsensemble.list <- function(x,...,plot=FALSE,is=NULL,ip=NULL,
                                     map.show=TRUE,alpha=0.6,xrange=NULL,yrange=NULL,
                                     main=NULL,verbose=FALSE,new=TRUE) {
  X <- x
  if (verbose) print('diagnose.dsensemble.list')
  stopifnot(inherits(X,"dsensemble") & inherits(X,"list"))
  locations <- loc(X)
  if (inherits(X,"pca")) X <- as.station(X,is=is,ip=ip,verbose=verbose)
  gcms <- attr(X[[1]],"model_id")
  if (length(gcms)==0) gcms <- 1:dim(X[[1]])[2]
  if (verbose) print("Compare variance and trends")
  outside <- matrix(NA,length(X))
  deltagcm <- matrix(NA,length(X),length(gcms))
  deltaobs <- matrix(NA,length(X))
  N <- matrix(NA,length(X))
  for (i in 1:length(X)) {
    if (verbose) print(loc(X[[i]]))
    di <- diagnose(X[[i]],plot=FALSE,verbose=verbose)
    outside[i] <- di$outside
    deltaobs[i] <- di$deltaobs
    deltagcm[i,] <- di$deltagcm
    N[i] <- di$N
  }
  d <- list(outside=outside,deltaobs=deltaobs,deltagcm=deltagcm,
            N=di$N,location=locations)
  
  if(is.null(main)) main <- attr(X,"variable")[1]
  if(plot) {
    if(verbose) print("target plot") 
    par0 <- par()
    if (new) dev.new()
    par(bty="n",fig=c(0.05,0.95,0,0.95),mgp=c(2,1,.5),xpd=TRUE)
    plot(c(-100,100),c(-100,100),type="n",
         axes=FALSE,ylab="",xlab="",main=main)
    mtext("trend",side=1,line=1.5,cex=par("cex"))
    mtext("interannual variability",side=2,line=2,cex=par("cex"))
    u <- par("usr")
    dx <- (u[2]-u[1])/20
    dy <- (u[4]-u[3])/20
    arrows(u[1]+dx,u[3],u[2]-dx,u[3],
           lwd=0.75,length=0.1,angle=20,code=2,xpd=NA)
    arrows(u[2]-dx,u[3],u[1]+dx,u[3],
           lwd=0.75,length=0.1,angle=20,code=2,xpd=NA)
    arrows(u[1],u[4]-dy,u[1],u[3]+dy,
           lwd=0.75,length=0.1,angle=20,code=2,xpd=NA)
    arrows(u[1],u[3]+dy,u[1],u[4]-dy,
           lwd=0.75,length=0.1,angle=20,code=2,xpd=NA)
    mtext("ensemble > obs",side=1,line=0,adj=0,cex=par("cex")*0.75)
    mtext("ensemble < obs",side=1,line=0,adj=1,cex=par("cex")*0.75)
    mtext("ensemble > obs",side=2,line=0.5,adj=0,cex=par("cex")*0.75)
    mtext("ensemble < obs",side=2,line=0.5,adj=1,cex=par("cex")*0.75)  
    bcol=c("grey95","grey40")
    for (i in 1:10) {
      ri <- (11-i)*10
      polygon(ri*cos(pi*seq(0,2,length=360)),
              ri*sin(pi*seq(0,2,length=360)),
              col=bcol[i %% 2 + 1],border="grey15")
    }
    xlon <- sapply(X,lon)
    xlat <- sapply(X,lat)
    nx <- (xlon-min(xlon))/diff(range(xlon))
    ny <- (xlat-min(xlat))/diff(range(xlat))
    if(all(is.finite(nx)) & all(is.finite(ny))) {
      col <- rgb(1-ny,nx,ny,alpha)
    } else {
      col <- adjustcolor(col="blue",alpha.f=alpha)
    }
    if(is.null(xrange)) xrange <- c(min(xlon)-1,max(xlon)+1)
    if(is.null(yrange)) yrange <- c(min(xlat)-1,max(xlat)+1)
    ## KMP 2020-08-04: Changed x (probability of observed trends) because the 
    ## calculations did not seem right. Trend distributions for different stations 
    ## were mixed up when applying mean(deltagcm) and sd(deltagcm). 
    #x <- -round(200*(0.5-pnorm(deltaobs,mean=mean(deltagcm),
    #                           sd=sd(deltagcm))),2)
    x <- -round(200*(0.5-sapply(seq_along(deltaobs), 
                                function(i) pnorm(deltaobs[i], 
                                        mean=mean(deltagcm[i,],na.rm=TRUE), sd=sd(deltagcm[i,],na.rm=TRUE)))))
    y <- -round(200*(0.5-pbinom(outside,size=N,prob=0.1)),2)
    d$x <- x; d$y <- y
    points(x,y,pch=21,cex=2*par("cex"),col='black',bg=col)
    if(map.show) {
      geoborders <- NULL
      data(geoborders, envir = environment())
      lon <- geoborders$x
      lat <- geoborders$y
      ok <- lon>min(xrange) & lon<max(xrange) &
        lat>min(yrange) & lat<max(yrange)
      lon2 <- attr(geoborders,"borders")$x
      lat2 <- attr(geoborders,"borders")$y
      ok2 <- lon2>min(xrange) & lon2<max(xrange) &
        lat2>min(yrange) & lat2<max(yrange)
      par(fig=c(0.8,0.97,0.8,0.97),new=TRUE, mar=c(0,0,0,0),
          cex.axis=0.75*par("cex"),xpd=NA,col.main="grey",bty="n")
      plot(lon[ok],lat[ok],lwd=1,col="black",type='l',
           xlab=NA,ylab=NA,axes=FALSE,xlim=xrange,ylim=yrange)
      axis(1,mgp=c(3,.5,0))
      axis(2,mgp=c(2,.5,0))
      lines(lon2[ok2],lat2[ok2],col = "pink",lwd=1)
      points(xlon,xlat,pch=21,cex=1,col='Grey',bg=col,lwd=0.5)
    }
    par(bty=par0$bty,fig=par0$fig,mgp=par0$mgp,xpd=par0$xpd)
  }
  invisible(d)
}

#' @exportS3Method
#' @export diagnose.matrix
diagnose.matrix <- function(x,...,xlim=NULL,ylim=NULL,verbose=FALSE) {
  if (verbose) print('diagnose.matrix')
  plot.diagnose.matrix(x,xlim=xlim,ylim=ylim,verbose=verbose)
}
metno/esd documentation built on Sept. 9, 2024, 5:07 a.m.