Merge branch 'master' of https://github.com/metno/esd

Author: abdelkaderm

NAMESPACE

@@ -539,10 +539,11 @@ export(subset.trajectory)
 export(corfield.trajectory)
 export(DS.trajectory)
 
-#export(dX)
-#export(dY)
-#export(dT)
-#export(regfit)
+export(CCI)
+export(dX)
+export(dY)
+export(dT)
+export(regfit)
 
 export(histwet)
 export(visprob)

R/DS.R

@@ -264,12 +264,12 @@ DS.station <- function(y,X,biascorrect=FALSE,mon=NULL,
 
     if ( (!inherits(X,'eof')) & (inherits(X,'field')) ) {
                                         #print("HERE")
-        ds <- DS.field(y=y,X=X,biascorrect=biascorrect,mon=mon,
-                       method=method,swsm=swsm,m=m,
-                       rmtrend=rmtrend,eofs=eofs,
-                       area.mean.expl=area.mean.expl,verbose=verbose,
-                       weighted=TRUE,pca=FALSE,npca=20,...) 
-        return(ds)
+      ds <- DS.field(y=y,X=X,biascorrect=biascorrect,mon=mon,
+                     method=method,swsm=swsm,m=m,
+                     rmtrend=rmtrend,eofs=eofs,
+                     area.mean.expl=area.mean.expl,verbose=verbose,
+                     weighted=TRUE,pca=FALSE,npca=20,...) 
+      return(ds)
     } else if (is.list(X)) {
                                         # REB 2014-10-08
         if (verbose) print("The predictor is a list")
@@ -426,7 +426,7 @@ DS.comb <- function(y,X,biascorrect=FALSE,mon=NULL,
         Y <- zoo(z,order.by=index(Z))
         Y <- attrcp(Z,Y)
 
-                                        # Store in a similar structure as in the common EOF:
+        # Store in a similar structure as in the common EOF:
         eval(parse(text=paste("attr(ds,'appendix.",i,"') <- Y",sep="")))
     }
     rm("X0"); gc(reset=TRUE)
@@ -665,6 +665,25 @@ DS.pca <- function(y,X,biascorrect=FALSE,mon=NULL,
                    verbose=verbose,weighted=weighted,pca=pca,npca=npca,...)
       return(z)
     }
+
+    if (!inherits(X,"eof") & inherits(X,"zoo")) {
+      ## If the predictor is an index, then use the same code to
+      ## estimate teleconnection patterns
+      if (verbose) print('Predictor is a zoo object')
+      attr(X,'history') <- history.stamp()
+      attr(X,'pattern') <- attr(y,'pattern')
+      attr(X,'eigenvalues') <- attr(y,'eigenvalues')
+      attr(X,'longitude') <- lon(y)
+      attr(X,'latitude') <- lat(y)
+      attr(X,'variable') <- varid(y)
+      attr(X,'unit') <- unit(y)
+  
+      class(X) <- c('eof',class(X))
+      z <- DS.pca(y,X,method=method,swsm=swsm,m=m,
+                  eofs=eofs,rmtrend=rmtrend,verbose=verbose,
+                  weighted=weighted,...)
+      return(z)
+    }
     stopifnot(!missing(y),!missing(X),
               inherits(X,"eof"),inherits(y,"station"))
 
@@ -675,8 +694,8 @@ DS.pca <- function(y,X,biascorrect=FALSE,mon=NULL,
                                         # synchronise the two zoo objects through 'merge' (zoo)
     y <- matchdate(y,it=X,verbose=verbose) # REB: 2014-12-16
     X <- matchdate(X,it=y,verbose=verbose) # REB: 2014-12-16
-    dy <- dim(y)
-    dx <- dim(X)
+    dy <- dim(y); if (is.null(dy)) dy <- c(length(y),1)
+    dx <- dim(X); if (is.null(dx)) dx <- c(length(X),1)
 
                                         # Use method for downscaling
                                         #str(y); str(X)
@@ -729,9 +748,12 @@ DS.pca <- function(y,X,biascorrect=FALSE,mon=NULL,
           if (verbose) print(attr(X0,'n.apps'))
           Xp <- attr(X0,'appendix.1')
         } else Xp <- X
+
+        
         y.out <- matrix(rep(NA,dx[1]*dy[2]),dx[1],dy[2])
         fit.val <- y.out
-        yp.out <- matrix(rep(NA,dim(Xp)[1]*dy[2]),dim(Xp)[1],dy[2])
+        dxp <- dim(Xp); if (is.null(dxp)) dxp <- c(length(Xp),1)
+        yp.out <- matrix(rep(NA,dxp[1]*dy[2]),dxp[1],dy[2])
 
         if (verbose) print(paste('PC',eofs,collapse=' '))
                                         # Loop over the PCs...
@@ -740,16 +762,16 @@ DS.pca <- function(y,X,biascorrect=FALSE,mon=NULL,
         ## pattern for each PC. Combine into one matrix. The predictor pattern
         ## for each station can be recovered by multiplying with the PCA pattern
         if (verbose) print('Predictor pattern')
-        #browser()
-        x0p <- attr(X0,'pattern') 
+        x0p <- attr(X0,'pattern')
         dp <- dim(x0p)
         if (is.null(dp)) dp <- c(length(x0p),1,1)  # list combining EOFs
+        if (length(dp)==2) dp <- c(dp,1)[c(1,3,2)] # if PCA rather than EOF
         if (verbose) {str(x0p); print(dp); print(dy)}
         predpatt <- rep(NA,dp[1]*dp[2]*dy[2])
         dim(predpatt) <- c(dp[1]*dp[2],dy[2])
         dim(x0p) <- c(dp[1]*dp[2],dp[3])
         ## This works for single predictors, but is more complicated for
-        #3 multiple predictors.
+        ## multiple predictors.
         if (dp[3] == length(attr(X0,'eigenvalues'))) x0p <- x0p %*% diag(attr(X0,'eigenvalues'))
         model <- list(); eof <- list()
         for (i in 1:dy[2]) {

R/DSensemble.R

@@ -374,7 +374,7 @@ DSensemble.t2m <- function(y,plot=TRUE,path="CMIP5.monthly/",
                       attr(X,'area.mean.expl') <- FALSE
   class(X) <- c("dsensemble","zoo")
   if (!is.null(path.ds)) file.ds <- file.path(path.ds,file.ds)
-  save(file=file.ds,X)#"DSensemble.rda",X)
+  #save(file=file.ds,X)#"DSensemble.rda",X)
   print("---")
   invisible(X)
 } 
@@ -1028,7 +1028,7 @@ DSensemble.pca <- function(y,plot=TRUE,path="CMIP5.monthly/",
     # Recursive: do each season seperately if there are more than one season
     if (length(table(season(y)))>1) {
       if (verbose) print('--- Apply DS to seasons seperately ---')
-      Z <- list(info='DSensembe.pca for different seasons')
+      Z <- list(info='DSensemble.pca for different seasons')
       for (season in names(table(season(y)))) {
         if (verbose) print(paste('Select',season))
         z <- DSensemble.pca(subset(y,it=season),plot=plot,path=path,
@@ -1080,8 +1080,8 @@ DSensemble.pca <- function(y,plot=TRUE,path="CMIP5.monthly/",
   flog <- file("DSensemble.pca-log.txt","at")
 
   ## Set up a list environment to keep all the results
-  Q <- list(info='DSensemble.pca',pca=y)
-  
+  Q <- list(info='DSensemble.pca',pca=y) ## KMP 06-08-2015
+ 
   if (verbose) print("loop...") 
   for (i in 1:N) {
     if (verbose) print(ncfiles[select[i]])
@@ -1092,7 +1092,7 @@ DSensemble.pca <- function(y,plot=TRUE,path="CMIP5.monthly/",
     gcmnm[i] <- paste(attr(gcm,'model_id'),attr(gcm,'realization'),sep="-r")
     if (inherits(y,'season')) {
       GCM <- as.4seasons(gcm,FUN=FUNX)
-      GCM <- subset(gcm,it=season(T2M)[1])
+      GCM <- subset(GCM,it=season(T2M)[1])
     } else if (inherits(y,'annual')) {
       GCM <- annual(gcm,FUN=FUNX)
     } else if (inherits(y,'month')) {
@@ -1102,6 +1102,7 @@ DSensemble.pca <- function(y,plot=TRUE,path="CMIP5.monthly/",
     }
     rm("gcm"); gc(reset=TRUE)
     T2MGCM <- combine(T2M,GCM)
+    
     if (verbose) print("- - - > EOFs")
     Z <- try(EOF(T2MGCM))
 
@@ -1118,36 +1119,26 @@ DSensemble.pca <- function(y,plot=TRUE,path="CMIP5.monthly/",
     if (biascorrect) Z <- biasfix(Z)
     ds <- try(DS(y,Z,eofs=eofs,verbose=verbose))
     if (verbose) print("post-processing")
-
+    
     ## Keep the results for the projections:
     if (verbose) print('Extract the downscaled projection')
-    z <- attr(ds,'appendix.1')
+    z <- attr(ds,'appendix.1') ## KMP 09.08.2015
+    ##attr(z,'model') <- attr(ds,'model') ## KMP 09-08-2015
+    ## model takes up too much space! can it be stored more efficiently?
+    attr(z,'predictor.pattern') <- attr(ds,'predictor.pattern')
+    attr(z,'evaluation') <- attr(ds,'evaluation')
+    ##class(z) <- class(ds)
+    
     cl <- paste('Q$i',i,'_',gsub('-','.',gcmnm[i]),' <- z',sep='')
     eval(parse(text=cl))
-
+    ##browser()
     if (verbose) {
-        print('Test to see if as.station has all information needed')
-        test.stations <- as.station(z)
-    }
-    
-    # The test lines are included to assess for non-stationarity
-    if (non.stationarity.check) {
-        if (verbose) print('non.stationarity.check')
-        testds <- DS(testy,testZ,biascorrect=biascorrect,eofs=eofs)
-                                                                # REB 29.04.2014
-        testz <- attr(testds,'appendix.1')                      # REB 29.04.2014
-        difference.z <- testy - testz                           # REB 29.04.2014
+       print('Test to see if as.station has all information needed')
+       test.stations.ds <- as.station(ds)
+       a <- attrcp(y,z);  class(a) <- c("ds",class(y))
+       test.stations.z <- as.station(a)
     }
-        
-    i1 <- is.element(paste(years,months,sep='-'),
-                     paste(year(z),month(z),sep='-'))
-    i2 <- is.element(paste(year(z),month(z),sep='-'),
-                     paste(years,months,sep='-'))
-    
-    if (verbose) print(paste('i=',i,gcmnm[i],'data points',
-                             sum(i1),'=',sum(i2)))
-    X[,i,i1] <- z[i2,]
-    
+       
     # Diagnose the residual: ACF, pdf, trend. These will together with the
     # cross-validation and the common EOF diagnostics provide a set of
     # quality indicators.
@@ -1220,51 +1211,33 @@ DSensemble.pca <- function(y,plot=TRUE,path="CMIP5.monthly/",
                     round(quality)))
   }
 
-  if (verbose) print('Downscaling finished - need to convert PCAs to stations')
+  if (verbose) print('Downscaling finished')
 
   ## Unpacking the information tangled up in GCMs, PCs and stations:
   ## Save GCM-wise in the form of PCAs
   gcmnm <- gsub('-','.',gcmnm)
-  Y <- as.station(y)
-
-
-  ## Save the information station-wise 
-#  if (verbose) print('Save the results station-wise')
-#  x <- matrix(rep(NA,dim(X)[3]*dim(X)[2]),dim(X)[3],dim(X)[2])
-#  ns <- length(loc(y))
-#  if (verbose) print(loc(y))
-#  for (i in 1:(ns)) {
-#    x[,] <- NA
-#    ## Loop over the GCMs and pick the selected downscaled station series:
-#    for (j in 1:N){
-#      x[,j] <- as.station(Z[[j+1]])[,i]
-#    }
-#    yloc <- tolower(loc(y)[i])
-#    yloc <- gsub('-','.',yloc)
-#    yloc <- gsub(' ','.',yloc)
-#    yloc <- gsub('/','.',yloc)
-#    if (verbose) {print(yloc); str(x)}
-#    ds.x <- zoo(x,order.by=index(Z[[2]]))
-#    attr(ds.x,'station') <- subset(Y,is=i)
-#    class(ds.x) <- c('dsensemble','zoo')
-#    eval(parse(text=paste('Z$',yloc,' <- ds.x',sep='')))
-#  }
+  ##Y <- as.station(y)
 
   #Z <- attrcp(y,Z)
   if (verbose) print('Set attributes')
-  attr(Z,'predictor') <- attr(T2M,'source')
-  attr(Z,'domain') <- list(lon=lon,lat=lat)
-  attr(Z,'scorestats') <- scorestats
-  attr(Z,'path') <- path
-  attr(Z,'scenario') <- rcp
-  attr(Z,'history') <- history.stamp(y)
+  attr(Q,'predictor') <- attr(T2M,'source')
+  attr(Q,'domain') <- list(lon=lon,lat=lat)
+  attr(Q,'scorestats') <- scorestats
+  attr(Q,'path') <- path
+  attr(Q,'scenario') <- rcp
+  attr(Q,'variable') <- attr(y,"variable")[1]
+   attr(Q,'unit') <- attr(y,"unit")[1]
+  attr(Q,'history') <- history.stamp(y)
   if (non.stationarity.check)
-    attr(Z,'on.stationarity.check') <- difference.z else
-    attr(Z,'on.stationarity.check') <- NULL
-  attr(Z,'area.mean.expl') <- FALSE
+    attr(Q,'on.stationarity.check') <- difference.z else
+    attr(Q,'on.stationarity.check') <- NULL
+  attr(Q,'area.mean.expl') <- FALSE
+  class(Q) <- c("dsensemble","pca","list")
 
+  ##browser()
+  dse.pca <- Q
   if(!is.null(path.ds)) file.ds <- file.path(path.ds,file.ds) ## KMP 2015-08-05
-  save(file=file.ds,Z)
+  save(file=file.ds,dse.pca)
   if (verbose) print("---")
-  invisible(Z)
+  invisible(dse.pca)
 }

R/diagnose.R

@@ -405,8 +405,8 @@ diagnose.dsensemble <- function(x,plot=TRUE,type='target',...) {
   N <- sum(i1,na.rm=TRUE)
   ## browser()
   if (plot) {
-    x <- -round(200*(0.5-pbinom(outside,size=N,prob=0.1)),2)
-    y <- -round(200*(0.5-pnorm(deltaobs,mean=mean(deltagcm),sd=sd(deltagcm))),2)
+    x <- -round(200*(0.5-pnorm(deltaobs,mean=mean(deltagcm),sd=sd(deltagcm))),2)
+    y <- -round(200*(0.5-pbinom(outside,size=N,prob=0.1)),2)
     #print(c(x,y))
 
     par(bty="n",xaxt="n",yaxt="n")

R/infographics.R

@@ -744,7 +744,7 @@ visprob.station.precip <- function(x,y=NULL,is=1,threshold=1,dy=0.005,
   if (verbose) print('visprob.station.precip')
   ## If y is provided, synchronise the two time series:
   if (!is.null(y)) {
-    y <- subset(y,,it=c(start(x),end(x)))
+    y <- subset(y,it=c(start(x),end(x)))
     x <- subset(x,it=c(start(y),end(y)))
     y <- annual(y)
   } else y <- year(annual(x))

R/plot.R

@@ -15,11 +15,9 @@ plot.station <- function(x,plot.type="single",new=TRUE,
   ##     else              
   ##         ylim <- apply(x,2,pretty,n=5)
 
-  ##<<<<<<< HEAD
-  ##=======
   if (is.null(ylim))
       ylim <- pretty(as.numeric(x))
-  ##>>>>>>> d6d9c84656c9b9b73e711c2e7ee2c8d0fb230980
+  
   unit <- attr(x,'unit')[1]
   for (i in 1:length(unit)) {
     if ( (is.na(unit[i]) | is.null(unit[i])) ) unit[i] <- " "

R/wheel.R

@@ -4,13 +4,23 @@
 
 wheel <- function(x,...) UseMethod("wheel")
 
-wheel.station <- function(x,new=TRUE,lwd=2,col=NULL,bg="grey90",...) {
+wheel.station <- function(x,y=NULL,new=TRUE,lwd=2,col=NULL,
+                          bg="grey90",verbose=FALSE,...) {
+
+  if (verbose) print('wheel.station')
+  ## Copied from visprob.station.precip:
+  ## If y is provided, synchronise the two time series:
+  if (!is.null(y)) {
+    y <- subset(y,it=c(start(x),end(x)))
+    x <- subset(x,it=c(start(y),end(y)))
+    years <- annual(y)
+  } else years <- year(annual(x))
   mx <- max(abs(coredata(x)),na.rm=TRUE)
 #  r <- mean(coredata(x),na.rm=TRUE)
   r <- mx
   #jday <- julian(zoo(x))
   MD <- month(x)*100 + day(x)
-  years <- as.integer(rownames(table(year(x))))
+  #REB2015-08-20 years <- as.integer(rownames(table(year(x))))
   #print(years)
   ny <- length(years)
   md <- as.integer(rownames(table(MD)))
@@ -33,6 +43,12 @@ wheel.station <- function(x,new=TRUE,lwd=2,col=NULL,bg="grey90",...) {
   w <- seq(0,2*pi,length=m)
   j <- 1:ny
   col <- rgb(j/ny,abs(sin(pi*j/ny)),(1-j/ny))
+  ## REB 2015-08-20: if y is given, use it for setting the colours.
+  ## e.g. according to an index like NINO3.4
+  srtc <- order(years)
+  col <- col[srtc]
+  years <- year(annual(x))
+  if (verbose) {print(srtc); print(col)}
 
   for (i in 1:m) {
     wj <- -w[i]
@@ -62,13 +78,22 @@ wheel.station <- function(x,new=TRUE,lwd=2,col=NULL,bg="grey90",...) {
   image(1:2,years,colbar,col=col)
 }
 
-wheel.spell <- function(x,new=TRUE,lwd=2,col=NULL...) {
+wheel.spell <- function(x,y=NULL,new=TRUE,lwd=2,col=NULL,verbose=FALSE,...) {
+
+    if (verbose) print('wheel.spell')
+    ## Copied from visprob.station.precip:
+  ## If y is provided, synchronise the two time series:
+  if (!is.null(y)) {
+    y <- subset(y,it=c(start(x),end(x)))
+    x <- subset(x,it=c(start(y),end(y)))
+    years <- annual(y)
+  } else years <- year(annual(x))
   mx <- max(abs(coredata(x)),na.rm=TRUE)
 #  r <- mean(coredata(x),na.rm=TRUE)
   r <- mx
   #jday <- julian(zoo(x))
   MD <- month(x)*100 + day(x)
-  years <- as.integer(rownames(table(year(x))))
+  #REB2015-08-20 years years <- as.integer(rownames(table(year(x))))
   #print(years)
   ny <- length(years)
   md <- as.integer(rownames(table(MD)))
@@ -91,6 +116,12 @@ wheel.spell <- function(x,new=TRUE,lwd=2,col=NULL...) {
   w <- seq(0,2*pi,length=m)
   j <- 1:ny
   col <- rgb(j/ny,abs(sin(pi*j/ny)),(1-j/ny))
+  ## REB 2015-08-20: if y is given, use it for setting the colours.
+  ## e.g. according to an index like NINO3.4
+  srtc <- order(years)
+  col <- col[srtc]
+  years <- year(annual(x))
+  if (verbose) {print(srtc); print(col)}
 
   for (i in 1:m) {
     wj <- -w[i]