diff --git a/DESCRIPTION b/DESCRIPTION index 551e61d1..7db48753 100644 --- 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 diff --git a/R/as.annual.R b/R/as.annual.R index 53a5859e..0b8074c9 100644 --- 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 index 17b6fc91..f9ebc413 100644 --- a/R/cumugram.R +++ b/R/cumugram.R @@ -14,6 +14,7 @@ #' @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 @@ -21,7 +22,7 @@ #' 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 index 413639d6..4a619c32 100644 --- 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 index eb932f5a..59eb6a88 100644 --- 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 index 7745c07f..1c43bbbc 100644 --- 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 index a2911092..47e2a1ba 100644 --- a/man/cumugram.Rd +++ b/man/cumugram.Rd @@ -13,6 +13,7 @@ cumugram( verbose = FALSE, FUN = "mean", main = NULL, + new = FALSE, ... ) }