small improvements

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: esd
-Version: 1.10.86
-Date: 2024-07-10
+Version: 1.10.87
+Date: 2024-08-09
 Title: Climate analysis and empirical-statistical downscaling (ESD) package for monthly and daily data
 Author: Rasmus E. Benestad, Abdelkader Mezghani, Kajsa M. Parding, Helene B. Erlandsen, Ketil Tunheim, and Cristian Lussana    
 Maintainer: Rasmus E. Benestad <rasmus.benestad@met.no>

R/diagnose.R

@@ -526,6 +526,7 @@ diagnose.dsensemble.list <- function(x,...,plot=FALSE,is=NULL,ip=NULL,
   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))
@@ -590,8 +591,9 @@ diagnose.dsensemble.list <- function(x,...,plot=FALSE,is=NULL,ip=NULL,
     ## 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,]), sd=sd(deltagcm[i,])))))
+    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)

R/map.R

@@ -413,6 +413,9 @@ map.eof <- function(x,...,it=NULL,is=NULL,new=FALSE,projection="lonlat",what="eo
         z <- map2sphere(X,it=it,lonR=lonR,latR=-90,axiR=axiR,lab=lab,
                         xlim=xlim,ylim=ylim,type=type,gridlines=gridlines,
                         colbar=colbar,new=new,verbose=verbose,...)
+      } else if (length(grep('moll|aea|utm|stere|robin',projection))>0) {
+        z <- map.sf(X,projection=projection,xlim=xlim,ylim=ylim,type=type,
+                    gridlines=gridlines,colbar=colbar,...)
       }
     } else z <- X
   }

R/plot.R

@@ -403,6 +403,7 @@ plot.station <- function(x,...,plot.type="single",new=TRUE,
       }
     }  
   }
+  if (verbose) print(paste('main:',main))
   #if (is.null(main)) main <- attr(x,'longname')[1] 
   if (is.null(col)) {
     if (is.null(dim(x))) {
@@ -452,9 +453,10 @@ plot.station <- function(x,...,plot.type="single",new=TRUE,
   if("seasonalcycle" %in% cls) {
     axis(1,at=seq(1,12),labels=month.abb,cex.axis=cex.axis,las=2)
   }
-
+  if (verbose) print(plot.type)
   if (plot.type=="single") {
-    if (errorbar) {
+    ## REB: added test: 2024-08-06
+    if ( errorbar & (length(x)==length(err(x))) ) {
       segments(index(x),x-err(x),index(x),x+err(x),
                lwd=3,col=rgb(0.5,0.5,0.5,0.25))
     }

demo/illustrate.DS.R

@@ -1,13 +1,30 @@
 ## Illustrate downscaling with common EOFs using sample EOFs from the esd-package
 
+warrow <- function(col='grey',border='darkgrey',x=0.5,y=0.1,n=100) {
+  par(mar=c(1,1,1,1),bty='n',xaxt='n',yaxt='n')
+  plot(c(0,1.5),c(0,1),type='n',xlab='',ylab='')
+  s <- 0.05*sin(seq(0,2*pi*3,length=n))
+  x <- c(0,0.1,seq(0.1,0.9,length=n),0.9,1,
+           0.9,seq(0.9,0.1,length=n),0.1,0)
+
+  y <- c(0.5,1,rep(0.7,n)+s,1,0.5,
+             0,rep(0.3,n)+s,0,0.5)
+  polygon(x,y,col=col,border=border)
+  #rect(0,0,1,1,lty=2)
+}
 
 frame <- function(x,ip,col='black') {
+  data(geoborders)
   d <- dim(x)
   nx <- d[1]+2*d[3]; ny <- d[2]+2*d[3]
   X <- matrix(rep(NA,nx*ny),nx,ny)
-  X[2*ip + 1:d[1],2*ip + 1:d[2]] <- x[,,ip]
+  X[2*ip + (1:d[1]) -1, 2*ip + (1:d[2])-1] <- x[,,ip]
+  #print(c(nx,ny,range(2*ip + (1:d[1]) -1),range(2*ip + (1:d[2])-1)))
   image(X,add=TRUE)
-  rect(2*ip,2*ip,2*ip+d[1],2*ip+d[2],border=col)
+  xcoords <- approx(1:nx,seq(0,1,length=nx),xout=c(2*ip,2*ip+d[1]-1))$y
+  ycoords <- approx(1:ny,seq(0,1,length=ny),xout=c(2*ip+1,2*ip+d[2]-1))$y
+  rect(xcoords[1],ycoords[1],xcoords[2],ycoords[2],border=col)
+  lines(geoborders,col='grey')
 }
 
 ## Function that draws a wavy arrow
@@ -18,35 +35,69 @@ similar <- function(x) {
 ## Get some sample data for illustration purposes  
 library(esd)
 print('NCEP')
-t2m.1 <- annual(subset(t2m.NCEP(),it=c(1991,2020),is=list(lon=c(-180,360),lat=c(-90,90))))
+t2m.1 <- annual(subset(t2m.NCEP(),it=c(1961,1990),is=list(lon=c(-180,360),lat=c(-90,90))))
 print('NorESM.M')
-t2m.2 <- annual(subset(t2m.NorESM.M(),it=c(1991,2020),is=list(lon=c(-180,360),lat=c(-90,90))))
+t2m.2 <- annual(subset(t2m.NorESM.M(),it=c(1961,1990),is=list(lon=c(-180,360),lat=c(-90,90))))
 print('combine')
 t2m <- combine(t2m.1,t2m.2)
 print('EOF')
-eof.t2m.NCEP <- EOF(t2m,n=5)
-x <- attr(eof.t2m.NCEP,'pattern')
+eof.t2m <- EOF(t2m,n=3)
+x <- t(rbind(coredata(t2m),coredata(attr(t2m,'appendix.1'))))
+dim(x) <- c(attr(t2m,'dimensions')[1:2],60)
 ## Dimensions of the spatial patterns
 d <- dim(x)
 print(d)
-nx <- d[1]; ny <- d[2]; np <- d[3]
-col <- c(rep('black',60),rep('grey',60))
+nx <- d[1]; ny <- d[2]; nt <- d[3]
+col <- c(rep('black',30),rep('grey',30))
 par(xaxt='n',yaxt='n',bty='n',xpd=NA)
 X <- matrix(rep(NA,nx*ny),nx,ny)
 
 ## First figure: Illustrate the EOFs.
 print('Fig. A1')
-layout(matrix(c(rep(0,5),c(1,1,0,2,2),c(1,1,3,2,2),c(1,1,0,4,4),c(0,0,5,4,4)),5,5,byrow = TRUE))
-image(X)
-z <- t(as.matrix(t2m))
-dim(z) <- attr(t2m,'dimensions')
-nt <- dim(z)[3]
-for (it in seq(nt,1,by=-1)) frame(z,it,col=col[it])
+layout(matrix(c(c(1,1,0,2,2),c(1,1,0,2,2),c(1,1,3,4,4),c(1,1,0,4,4),c(1,1,0,5,0)),5,5,byrow = TRUE))
+image(X,xlab='',ylab='',main=paste('Maps: ',d[1],'x',d[2],'x',d[3]))
+
+for (it in seq(nt,1,by=-1)) frame(x,it,col=col[it])
+
+image(X,xlab='',ylab='',main=paste('Modes: ',d[1],'x',d[2],'x 5 + (5 x',d[3],'+ 5)'))
+z <- attr(eof.t2m,'pattern')
+np <- dim(z)[3]
+for (ip in seq(np,1,by=-1)) frame(z,ip)
+
+warrow()
+
+par(mar=c(1,1,1,1),xaxt='s',yaxt='s',bty='n',xpd=NA)
+PCs <- zoo(rbind(coredata(eof.t2m),coredata(attr(eof.t2m,"appendix.1"))))
+plot(PCs,plot.type='single',xlab='')
+rect(1,-0.3,30.5,0.3,col=rgb(1,0.8,0.8),border=rgb(1,0.9,0.9))
+rect(30.5,-0.3,60,0.3,col=rgb(0.8,0.8,1),border=rgb(0.9,0.9,1))
+for (ip in 1:np) {
+  lines(PCs[1:30,ip],col='black',lwd=2,lty=ip)
+  lines(PCs[31:60,ip],col='grey',lwd=2,lty=ip)
+}
+
+plot(100*attr(eof.t2m,'eigenvalues')^2/attr(eof.t2m,'tot.var'),
+     type='b',ylab='Variance (%)',xlab='Mode',ylim=c(0,100),lwd=2)
 
 ## Second figure: illustrate the downscaling
 print('Fig. A2')
+par(xaxt='n',yaxt='n',bty='n',xpd=NA)
 layout(matrix(1:4,2,2))
-## Get spatial patterns 
+## Get a time series that represents the predictand (local information)
+
+nino3.4 <- NINO3.4()
+
+par(xaxt='s',yaxt='s',bty='n',xpd=NA)
+plot(PCs,plot.type='single',xlab='')
+rect(1,-0.3,30.5,0.3,col=rgb(1,0.8,0.8),border=rgb(1,0.9,0.9))
+rect(30.5,-0.3,60,0.3,col=rgb(0.8,0.8,1),border=rgb(0.9,0.9,1))
+for (ip in 1:np) {
+  lines(PCs[1:30,ip],col='black',lwd=2,lty=ip)
+  lines(PCs[31:60,ip],col='grey',lwd=2,lty=ip)
+}
+
+plot(nino3.4)
+
+plot(zoo(DS(nino3.4,eof.t2m)))
+
 
-image(X)
-for (ip in seq(np,1,by=-1)) frame(x,ip)