From 5220e2e96dc35aa7096c3087022a97a813ca9a98 Mon Sep 17 00:00:00 2001 From: brasmus Date: Thu, 22 Aug 2024 13:29:31 +0200 Subject: [PATCH 1/4] minor improvements --- DESCRIPTION | 4 ++-- R/cumugram.R | 5 +++-- R/g2dl.R | 8 +++++--- R/lonlatprojection.R | 16 +++++++++++++--- R/map.R | 6 +++--- 5 files changed, 26 insertions(+), 13 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 2b3c8c01..df0906da 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: esd -Version: 1.10.87 -Date: 2024-08-09 +Version: 1.10.88 +Date: 2024-08-22 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/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..ffeb4d4e 100644 --- a/R/g2dl.R +++ b/R/g2dl.R @@ -95,10 +95,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) <- d y <- zoo(t(X),index(x)) lon <- sort(lon) @@ -106,7 +108,7 @@ g2dl.field <- function(x,greenwich=TRUE,verbose=FALSE,...) { #nattr <- softattr(x,ignore=c('greenwich','longitude')) #for (i in 1:length(nattr)) # attr(y,nattr[i]) <- attr(x,nattr[i]) - attr(y,'dimensions') <- attr(x,'dimensions') + attr(y,'dimensions') <- d attr(y,'longitude') <- lon attr(y,'greenwich') <- as.logical(greenwich) attr(y,'history') <- history.stamp(x) diff --git a/R/lonlatprojection.R b/R/lonlatprojection.R index 6629dcb1..43df2113 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 @@ -80,13 +80,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] @@ -185,8 +188,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 (useRaster) { + ## ‘useRaster = TRUE’ can only be used with a regular grid + ## Force the coordinates to be evenly spaced + if (verbose) print('ensure reggular 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 = TRUE, + 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 423f4080..3f65051a 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,7 +147,7 @@ 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... @@ -182,7 +182,7 @@ map.default <- function(x,...,FUN='mean',it=NULL,is=NULL,new=FALSE, 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,...) + 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,...) From 1f6705d07b8addb214416e044e0829092672430d Mon Sep 17 00:00:00 2001 From: brasmus Date: Fri, 23 Aug 2024 14:42:41 +0200 Subject: [PATCH 2/4] finor fix in map.station --- R/map.R | 101 +++++++++++++++++++++++++++++--------------------------- 1 file changed, 52 insertions(+), 49 deletions(-) diff --git a/R/map.R b/R/map.R index 3f65051a..02e496fa 100644 --- a/R/map.R +++ b/R/map.R @@ -149,55 +149,57 @@ map.default <- function(x,...,FUN='mean',it=NULL,is=NULL,new=FALSE, colbar= colbar,gridlines=gridlines,verbose=verbose, 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,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('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('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 @@ -781,7 +783,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 From 810b481d986f4606c33dce95c08219da38d0f283 Mon Sep 17 00:00:00 2001 From: brasmus Date: Mon, 26 Aug 2024 12:34:40 +0200 Subject: [PATCH 3/4] bugfix in g2dl.field --- R/g2dl.R | 9 +++++---- R/regrid.R | 2 +- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/R/g2dl.R b/R/g2dl.R index ffeb4d4e..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 @@ -100,7 +101,7 @@ g2dl.field <- function(x,greenwich=TRUE,verbose=FALSE,...) { dim(X) <- c(d[2],d[3],d[1]) X <- X[xsrt,,] if (verbose) {print(dim(X)); print(d)} - dim(X) <- d + dim(X) <- c(d[2]*d[3],d[1]) y <- zoo(t(X),index(x)) lon <- sort(lon) @@ -108,7 +109,7 @@ g2dl.field <- function(x,greenwich=TRUE,verbose=FALSE,...) { #nattr <- softattr(x,ignore=c('greenwich','longitude')) #for (i in 1:length(nattr)) # attr(y,nattr[i]) <- attr(x,nattr[i]) - attr(y,'dimensions') <- d + attr(y,'dimensions') <- attr(x,'dimensions') attr(y,'longitude') <- lon attr(y,'greenwich') <- as.logical(greenwich) attr(y,'history') <- history.stamp(x) diff --git a/R/regrid.R b/R/regrid.R index ca61eeb6..ffc60d1d 100755 --- a/R/regrid.R +++ b/R/regrid.R @@ -332,7 +332,7 @@ regrid.field <- function(x,is=NULL,...,it=NULL,verbose=FALSE,approach="field",cl mlon <- MATCH(attr(x,'longitude'), lon.new) mlat <- MATCH(attr(x,'latitude'), lat.new) #print(mlon); print(mlat) - # If thhose corordinates are the same as the original data, then + # If those coordinates are the same as the original data, then # return with the original data: if ( sum(is.na(c(mlon,mlat))) ==0 ) { if (verbose) print(summary(mlon,mlat)) From 16daa55b8e1ee482f2e16c1b0df1440cebbe9511 Mon Sep 17 00:00:00 2001 From: brasmus Date: Tue, 27 Aug 2024 08:51:23 +0200 Subject: [PATCH 4/4] minor fix in annual.default --- R/as.annual.R | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) 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