fixed conflicts and roxygenized

metno · Sep 3, 2024 · 8be7e28 · 8be7e28
2 parents 789664c + 16daa55
commit 8be7e28
Show file tree

Hide file tree

Showing 7 changed files with 86 additions and 70 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: esd
-Version: 1.10.88
-Date: 2024-08-19
+Version: 1.10.89
+Date: 2024-09-03
 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>

diff --git a/R/as.annual.R b/R/as.annual.R
@@ -149,7 +149,8 @@ annual.default <- function(x,FUN='mean',na.rm=TRUE, nmin=NULL,start=NULL,...,
   ## Check how many valid data points)
   nok <- aggregate(X,year,FUN='nv')
 
-  if (FUN == 'sum') na.rm <- FALSE ## AM
+  ## REB 2024-08-27: The argument 'nmin' takes care of NAs or by setting 'na.rm=FALSE' as argument.
+  #if (FUN == 'sum') na.rm <- FALSE ## AM
   if (verbose) print(paste('aggregate: FUN=',FUN))
 
   if (verbose) str(X)
@@ -163,13 +164,14 @@ annual.default <- function(x,FUN='mean',na.rm=TRUE, nmin=NULL,start=NULL,...,
     y <- aggregate(X,year,FUN=FUN,...,threshold=threshold) ## AM 20-05-2015
   } else if ((sum(is.element(names(formals(FUN)),'na.rm')==1)) |
              (sum(is.element(FUN,c('mean','min','max','sum','quantile')))>0)) {
-    if (verbose) print('Function has na.rm-argument')
+    if (verbose) print(paste(FUN,'has na.rm-argument:',na.rm))
     y <- aggregate(X,year,FUN=FUN,...,na.rm=na.rm)
   } else {
     if (verbose) print('Function has no treshold nor na.rm arguments')
     y <- aggregate(X,year,FUN=FUN,...) # REB
   }
   y[!is.finite(y)] <- NA ## AM
+  if (verbose) print(y)
 
   if (verbose) print('check for incomplete sampling')
   ## Flag the data with incomplete sampling as NA
@@ -187,7 +189,7 @@ annual.default <- function(x,FUN='mean',na.rm=TRUE, nmin=NULL,start=NULL,...,
   ## Check if the series is gappy:
   if (verbose) print('Fill in gaps')
   it <- seq(min(index(y)),max(index(y)),by=1)
-  if (verbose) print(c(range(it),length(it)))
+  if (verbose) print(c(range(it),length(it),sum(is.finite(y))))
   if(is.null(dim(y))) {
     z <- zoo(rep(NA,length(it)),order.by=it)
     z[is.element(it,index(y))] <- y

diff --git a/R/cumugram.R b/R/cumugram.R
@@ -14,14 +14,15 @@
 #' @param FUN a function
 #' @param verbose a boolean; if TRUE print information about progress
 #' @param main main title
+#' @PARAM new TRUE to open a new device/window
 #' @param \dots additional arguments
 #'
 #' @examples
 #' data(bjornholt)
 #' cumugram(bjornholt)
 #'
 #' @export
-cumugram <- function(x,it=NULL,start='-01-01',prog=FALSE,plot=TRUE,verbose=FALSE,FUN='mean',main=NULL,...) {
+cumugram <- function(x,it=NULL,start='-01-01',prog=FALSE,plot=TRUE,verbose=FALSE,FUN='mean',main=NULL,new=FALSE,...) {
   stopifnot(!missing(x),inherits(x,"station"))
 
   if(verbose) print("cumugram")
@@ -42,7 +43,7 @@ cumugram <- function(x,it=NULL,start='-01-01',prog=FALSE,plot=TRUE,verbose=FALSE
     eval(parse(text=paste("main <- paste('",titletext,"',
                           tolower(attr(x,'longname')),sep=' ')")))
 
-   if (plot) {dev.new(); par(bty="n")}
+   if (plot) {if (new) dev.new(); par(bty="n")}
 
   z <- coredata(x)
   ylim <- c(NA,NA)

diff --git a/R/g2dl.R b/R/g2dl.R
@@ -83,8 +83,9 @@ g2dl.field <- function(x,greenwich=TRUE,verbose=FALSE,...) {
   if(verbose) print("g2dl.field")
   attr(x,'longitude') -> lon
   attr(x,'latitude') -> lat
-  d <- attr(x,'dimensions')
-
+  ## REB 2024-08-26 - more robust code using the actual object size
+  #d <- attr(x,'dimensions')
+  d <- c(length(index(x)),length(lon),length(lat))
   if (greenwich) {
     wh <- lon < 0
     lon[wh] <- lon[wh] + 360
@@ -95,10 +96,12 @@ g2dl.field <- function(x,greenwich=TRUE,verbose=FALSE,...) {
 
   xsrt <- order(lon)
   xsrt <- xsrt[!duplicated(lon)]
+  if (verbose) print(paste('g2dl: length(lon)=',length(lon),length(xsrt)))
   X <- t(coredata(x))
-  dim(X) <- d
+  dim(X) <-  c(d[2],d[3],d[1])
   X <- X[xsrt,,]
-  dim(X) <- c(length(lon)*d[2],d[3])
+  if (verbose) {print(dim(X)); print(d)}
+  dim(X) <- c(d[2]*d[3],d[1])
   y <- zoo(t(X),index(x))
   lon <- sort(lon)
 

diff --git a/R/lonlatprojection.R b/R/lonlatprojection.R
@@ -8,7 +8,7 @@ lonlatprojection <- function(x,it=NULL,is=NULL,new=FALSE,projection="lonlat",
                              type=c("fill","contour"),gridlines=FALSE,
                              verbose=FALSE,geography=TRUE,fancy=TRUE,
                              main=NA,cex.sub=0.8,cex.axis=0.8,
-                             fig=NULL,add=FALSE,plot=TRUE,...) {
+                             fig=NULL,add=FALSE,plot=TRUE,useRaster=TRUE,...) {
 
   if (verbose) {print('lonlatprojection'); str(x)}
   ## REB 2024-04-29
@@ -79,13 +79,16 @@ lonlatprojection <- function(x,it=NULL,is=NULL,new=FALSE,projection="lonlat",
     greenwich <- FALSE
   }
   ## Make sure to use the right arrangement: from dateline or Greenwich 
+  if (verbose) {print(dim(x)); print(attr(x,'greenwich')); print(greenwich)}
   if(inherits(x,"matrix") & is.null(attr(x,"dimensions"))) {
     x <- g2dl(x,d=c(length(lon(x)),length(lat(x)),1),
               greenwich=greenwich,verbose=verbose)
   } else {
     x <- g2dl(x,greenwich=greenwich,verbose=verbose)
   }
-  if (verbose) print(paste('dimensions of x:',dim(x),'lon=',length(lon(x)),'lat=',length(lat(x)),collapse=' - '))
+  if (verbose) print(paste('dimensions of x:',paste(dim(x),collapse=' - '),
+                           'tim=',length(index(x)),'lon=',length(lon(x)),
+                           'lat=',length(lat(x))))
   dim(x) <- c(length(lon(x)),length(lat(x)))
   ## Make sure the longitudes are ordered correctly
   srtx <- order(lon(x)); lon <- lon(x)[srtx]
@@ -184,12 +187,15 @@ lonlatprojection <- function(x,it=NULL,is=NULL,new=FALSE,projection="lonlat",
          xlim=xlim,ylim=ylim,main=main,
          xaxt="n",yaxt="n") # AM 17.06.2015
 
-    if (sum(is.element(tolower(type),'fill'))>0)  
-      ## KMP 2024-08-12: Setting useRaster=TRUE only if the grid is regular. I added this because 
-      ##   the example data slp.ERA5 has a slightly uneven grid, for some reason, 
-      ##   which results in an error in one of the examples where it is used (see the function CCI).
-      useRaster <- length(unique(diff(lon)))==1 & length(unique(diff(lat)))==1 
-      image(lon,lat,x,xlab="",ylab="", add=TRUE, useRaster = useRaster, #useRaster = TRUE,
+    if (useRaster) {
+      ## ‘useRaster = TRUE’ can only be used with a regular grid
+      ## Force the coordinates to be evenly spaced
+      if (verbose) print('ensure regular grid')
+      lon <- seq(min(lon),max(lon),length=length(lon))
+      lat <- seq(min(lat),max(lat),length=length(lat))
+    }
+    if (sum(is.element(tolower(type),'fill'))>0)   
+      image(lon,lat,x,xlab="",ylab="", add=TRUE,useRaster = useRaster,
             col=colbar$col,breaks=colbar$breaks)
 
     if (geography) {

diff --git a/R/map.R b/R/map.R
@@ -136,7 +136,7 @@ map.default <- function(x,...,FUN='mean',it=NULL,is=NULL,new=FALSE,
                                      show=TRUE,type="p",cex=2,h=0.6,v=1),
                         type=c("fill","contour"),gridlines=FALSE,cex=2,
                         lonR=NULL,latR=NULL,axiR=NULL,style='plain',
-                        verbose=FALSE,plot=TRUE,add=FALSE) {
+                        verbose=FALSE,plot=TRUE,add=FALSE,useRaster=TRUE) {
 
 
   ## default with no arguments will produce a map showing available station
@@ -147,57 +147,59 @@ map.default <- function(x,...,FUN='mean',it=NULL,is=NULL,new=FALSE,
     z <- map.station(x,FUN=FUN,it=it,is=is,new=new,projection=projection,
                      xlim=xlim,ylim=ylim,zlim=zlim,n=15,
                      colbar= colbar,gridlines=gridlines,verbose=verbose,
-                     plot=plot,...)
+                     plot=plot,useRaster=useRaster,...)
     invisible(z)
-  }
-  par0 <- par(no.readonly = TRUE) # save default, for resetting...
-  if (is.logical(colbar)) colbar <- NULL
-  ## If only a few items are provided in colbar - then set the rest to the default
-  if (!is.null(colbar)) {
-    colbar <- colbar.ini(x,FUN=FUN,colbar=colbar,verbose=FALSE)
-  }
-  x <- subset(x,it=it,is=is)
-  X <- attr(x,'pattern')
-
-  ## if zlim is specified, then mask data outside this range
-  if (!is.null(zlim)) {
-    d <- dim(X)
-    mask <- (X < min(zlim)) | (X > max(zlim))
-    X[mask] <- NA
-    dim(X) <- d
-    if (verbose) {print(zlim); print(dim(X)); print(sum(mask))}
-  }
-  attr(X,'longitude') <- lon(x)
-  attr(X,'latitude') <- lat(x)
-  attr(X,'variable') <- attr(x,'variable')
-  attr(X,'unit') <- attr(x,'unit')[1]
-  if (attr(X,'unit') =='%') attr(X,'unit') <- "'%'"
-  attr(X,'source') <- attr(x,'source')
-  attr(X,'variable') <- varid(x)
-  if (inherits(X,'zoo')) {
-    attr(X,'time') <- range(index(x))
-  } else if (!is.null(attr(x,'time'))) {
-    attr(X,'time') <- attr(x,'time')
-  }
-  if (plot) {
-    if (projection=="lonlat") {
-      z <- lonlatprojection(x=X,xlim=xlim,ylim=ylim,colbar=colbar,verbose=verbose,
-                            lab=lab,type=type,new=new,gridlines=gridlines,...)
-    } else if (projection=="sphere") {
-      z <- map2sphere(x=X,lonR=lonR,latR=latR,axiR=axiR,xlim=xlim,ylim=ylim,
-                      lab=lab,type=type,gridlines=gridlines,colbar=colbar,new=new,...)
-    } else if (projection=="np") {
-      z <- map2sphere(X,lonR=lonR,latR=90,axiR=axiR,xlim=xlim,ylim=ylim,
-                      lab=lab,type=type,gridlines=gridlines,colbar=colbar,new=new,...)
-    } else if (projection=="sp") {
-      z <- map2sphere(X,lonR=lonR,latR=-90,axiR=axiR,new=new,xlim=xlim,ylim=ylim,
-                      lab=lab,type=type,gridlines=gridlines,colbar=colbar,...)
-    } else if (length(grep('+proj=|moll|aea|utm|stere|robin',projection))>0) {
-      z <- map.sf(X,projection=projection,xlim=xlim,ylim=ylim,type=type,
-                  gridlines=gridlines,colbar=colbar,...)
+  } else { 
+    par0 <- par(no.readonly = TRUE) # save default, for resetting...
+    if (is.logical(colbar)) colbar <- NULL
+    ## If only a few items are provided in colbar - then set the rest to the default
+    if (!is.null(colbar)) {
+      colbar <- colbar.ini(x,FUN=FUN,colbar=colbar,verbose=FALSE)
     }
-  } else z <- X
-  invisible(z)
+    x <- subset(x,it=it,is=is)
+    X <- attr(x,'pattern')
+
+    ## if zlim is specified, then mask data outside this range
+    if (!is.null(zlim)) {
+      d <- dim(X)
+      mask <- (X < min(zlim)) | (X > max(zlim))
+      X[mask] <- NA
+      dim(X) <- d
+      if (verbose) {print(zlim); print(dim(X)); print(sum(mask))}
+    }
+    if (verbose) print('map.default: Copy metadata')
+    attr(X,'longitude') <- lon(x)
+    attr(X,'latitude') <- lat(x)
+    attr(X,'variable') <- attr(x,'variable')
+    attr(X,'unit') <- attr(x,'unit')[1]
+    if (attr(X,'unit') =='%') attr(X,'unit') <- "'%'"
+    attr(X,'source') <- attr(x,'source')
+    attr(X,'variable') <- varid(x)
+    if (inherits(X,'zoo')) {
+      attr(X,'time') <- range(index(x))
+    } else if (!is.null(attr(x,'time'))) {
+      attr(X,'time') <- attr(x,'time')
+    }
+    if (plot) {
+      if (projection=="lonlat") {
+        z <- lonlatprojection(x=X,xlim=xlim,ylim=ylim,colbar=colbar,verbose=verbose,
+                              lab=lab,type=type,new=new,gridlines=gridlines,useRaster=useRaster,...)
+      } else if (projection=="sphere") {
+        z <- map2sphere(x=X,lonR=lonR,latR=latR,axiR=axiR,xlim=xlim,ylim=ylim,
+                        lab=lab,type=type,gridlines=gridlines,colbar=colbar,new=new,...)
+      } else if (projection=="np") {
+        z <- map2sphere(X,lonR=lonR,latR=90,axiR=axiR,xlim=xlim,ylim=ylim,
+                        lab=lab,type=type,gridlines=gridlines,colbar=colbar,new=new,...)
+      } else if (projection=="sp") {
+        z <- map2sphere(X,lonR=lonR,latR=-90,axiR=axiR,new=new,xlim=xlim,ylim=ylim,
+                        lab=lab,type=type,gridlines=gridlines,colbar=colbar,...)
+      } else if (length(grep('+proj=|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
+    invisible(z)
+  }
 }
 
 #' @exportS3Method
@@ -780,7 +782,8 @@ map.pca <- function(x,...,it=NULL,is=NULL,ip=1,new=FALSE,projection="lonlat",
       X[mask] <- NA
       dim(X) <- d
       if (verbose) {print(zlim); print(dim(X)); print(sum(mask))}
-    }    
+    }
+    if (verbose) print('map.pca: copy metadata')
     attr(X,'longitude') <- lon(x)
     attr(X,'latitude') <- lat(x)
     attr(X,'mean') <- NULL

diff --git a/man/cumugram.Rd b/man/cumugram.Rd