improvement of plot.dsensemble.one - more robust to missing data

metno · Mar 8, 2024 · 1726389 · 1726389
1 parent 63819c8
commit 1726389
Show file tree

Hide file tree

Showing 4 changed files with 51 additions and 31 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: esd
-Version: 1.10.76
-Date: 2024-03-06
+Version: 1.10.77
+Date: 2024-03-08
 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 <[email protected]>

diff --git a/R/DS.R b/R/DS.R
@@ -417,7 +417,8 @@ DS.default <- function(y,X,verbose=FALSE,plot=FALSE,it=NULL,
   attr(r2,'description') <- 'var(fitted.values))/var(y)'
   attr(ds,'quality') <- r2
   attr(ds,'variable') <- attr(y0,'variable')
-  attr(ds,'model') <- model
+  #attr(ds,'model') <- model ## REB 2024-03-08: Save the model summary rather than the model itself
+  attr(ds,'model') <- fsum
   attr(ds,'mean') <- offset
   attr(ds,'method') <- method
   attr(ds,'eof') <- X0
@@ -1025,7 +1026,8 @@ DS.pca <- function(y, X, verbose=FALSE, plot=FALSE, it=NULL, method="lm",
     ##                              # addition to offset
     model$calibration.data <- X
     class(model$fitted.values) <- class(y0)
-    attr(ds,'model') <- model
+    #attr(ds,'model') <- model ## REB 2024-03-08: Save the model summary rather than the model itself
+    attr(ds,'model') <- fsum
     attr(ds,'quality') <- var(coredata(model$fitted.values))/var(y,na.rm=TRUE)
     attr(ds,'eigenvalues') <- pca$d
     attr(ds,'sum.eigenv') <- sum(pca$d)
@@ -1153,7 +1155,8 @@ DS.pca <- function(y, X, verbose=FALSE, plot=FALSE, it=NULL, method="lm",
   attr(ds,'calibration_data') <- attr(z,'calibration_data')
   attr(ds,'fitted_values') <- zoo(fit.val,order.by=index(attr(z,'fitted_values')))
   class(attr(ds,'fitted_values')) <- class(y0)
-  attr(ds,'model') <- model
+  #attr(ds,'model') <- model ## REB 2024-03-08: Save the model summary rather than the model itself
+  attr(ds,'model') <- fsum
   attr(ds,'eof') <- eof
   attr(ds,'original_data') <- y0
   attr(ds,'variable') <- varid(y0)

diff --git a/R/DSensemble.R b/R/DSensemble.R
@@ -2143,25 +2143,26 @@ DSensemble.pca <- function(y,...,plot=TRUE,path="CMIP5.monthly/",rcp="rcp45",bia
       z <- attr(ds,'appendix.1') ## KMP 09.08.2015
 
       ## REB: 2016-11-29
-      if (test) {
-        ## model takes up too much space! can it be stored more efficiently?
-        ## REB 2016-11-29: remove most of the contents and keep only a small part
-        if (verbose) print('Add reduced model information')
-        for (iii in 1:dim(ds)[2]) {
-          print(names(attr(ds,'model')[[iii]]))
-          attr(ds,'model')[[iii]]$residuals <- NULL
-          attr(ds,'model')[[iii]]$effects <- NULL
-          attr(ds,'model')[[iii]]$rank <- NULL
-          attr(ds,'model')[[iii]]$fitted.values <- NULL
-          attr(ds,'model')[[iii]]$assign <- NULL
-          attr(ds,'model')[[iii]]$qr <- NULL
-          attr(ds,'model')[[iii]]$df.residual <- NULL
-          attr(ds,'model')[[iii]]$xlevels <- NULL
-          attr(ds,'model')[[iii]]$model <- NULL
-          attr(ds,'model')[[iii]]$terms <- NULL
-          print(names(attr(ds,'model')[[iii]]))
-        }
-      }
+      ## REB: 2024-03-08: Only save model summary in DS.
+      # if (test) {
+      #   ## model takes up too much space! can it be stored more efficiently?
+      #   ## REB 2016-11-29: remove most of the contents and keep only a small part
+      #   if (verbose) print('Add reduced model information')
+      #   for (iii in 1:dim(ds)[2]) {
+      #     print(names(attr(ds,'model')[[iii]]))
+      #     attr(ds,'model')[[iii]]$residuals <- NULL
+      #     attr(ds,'model')[[iii]]$effects <- NULL
+      #     attr(ds,'model')[[iii]]$rank <- NULL
+      #     attr(ds,'model')[[iii]]$fitted.values <- NULL
+      #     attr(ds,'model')[[iii]]$assign <- NULL
+      #     attr(ds,'model')[[iii]]$qr <- NULL
+      #     attr(ds,'model')[[iii]]$df.residual <- NULL
+      #     attr(ds,'model')[[iii]]$xlevels <- NULL
+      #     attr(ds,'model')[[iii]]$model <- NULL
+      #     attr(ds,'model')[[iii]]$terms <- NULL
+      #     print(names(attr(ds,'model')[[iii]]))
+      #   }
+      # }
 
       attr(z,'predictor.pattern') <- attr(ds,'predictor.pattern')
       ## REB: 2023-12-08 - add information about the regression/calibration coefficients.

diff --git a/R/plot.R b/R/plot.R
@@ -2252,7 +2252,7 @@ plot.dsensemble.one <-  function(x,pts=FALSE,it=0,
   args <- list(...)
   if (verbose) print(names(args))
   ixl <- grep('xlim',names(args))
-  if (length(ixl)==0) xlim <- range(year(z)) else
+  if (length(ixl)==0) xlim <- range(c(year(z),year(attr(x,'station')))) else
     xlim <- args[[ixl]]
   iyl <- grep('ylim',names(args))
   if (length(iyl)==0) ylim <- pscl*range(coredata(z),na.rm=TRUE) else
@@ -2289,20 +2289,36 @@ plot.dsensemble.one <-  function(x,pts=FALSE,it=0,
   for (ii in 1:49) {
     qp1 <- qnorm(1-ii/50,mean=coredata(mu),sd=coredata(si))
     qp2 <- qnorm(ii/50,mean=coredata(mu),sd=coredata(si))
+    ## REB 2024-03-08: Extra instructions dealing with missing data/NA
+    valid.data <- is.finite(qp1) &  is.finite(qp2)
     if(smooth) {
-      qp1 <- smooth.spline(year(z), qp1)$y
-      qp2 <- smooth.spline(year(z), qp2)$y
+      qp1 <- smooth.spline(year(z)[valid.data], qp1[valid.data])$y
+      qp2 <- smooth.spline(year(z)[valid.data], qp2[valid.data])$y
+      if (sum(valid.data)!=length(valid.data)) {
+        if (verbose) print('<smoothing dealing with missing data [qp1-qp2]>')
+        ## REB 2024-03-08: Need vectors with same length
+        qp1 <- approx(year(z)[valid.data],qp1,xout =year(z))$y
+        qp2 <- approx(year(z)[valid.data],qp2,xout =year(z))$y
+      }
     }
     ci <- c(qp1,rev(qp2))
     polygon(t2[!is.na(ci)],ci[!is.na(ci)], col= envcol, border=NA)
   }
   q05 <- qnorm(0.05,mean=mu,sd=si)
   q95 <- qnorm(0.95,mean=mu,sd=si)
-
+  ## REB 2024-03-08: Extra instructions dealing with missing data/NA
+  valid.data <- is.finite(mu) &  is.finite(q05) &  is.finite(q95)
   if(smooth) {
-    mu <- smooth.spline(year(z), mu)$y
-    q05 <- smooth.spline(year(z), q05)$y
-    q95 <- smooth.spline(year(z), q95)$y
+    mu <- smooth.spline(year(z)[valid.data], mu[valid.data])$y
+    q05 <- smooth.spline(year(z)[valid.data], q05[valid.data])$y
+    q95 <- smooth.spline(year(z)[valid.data], q95[valid.data])$y
+    if (sum(valid.data)!=length(valid.data)) {
+      if (verbose) print('<smoothing dealing with missing data [q05-mu-q95]>')
+      ## REB 2024-03-08: Need vectors with same length
+      mu <- approx(year(z)[valid.data],mu,xout =year(z))$y
+      q05 <- approx(year(z)[valid.data],q05,xout =year(z))$y
+      q95 <- approx(year(z)[valid.data],q95,xout =year(z))$y
+    }
   }
 
   lcol <- adjustcolor(envcol,offset=c(0.5,0.5,0.5,0.2))