Rewrote weather generator for daily precipitation WG.fwmu.day.precip

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: esd
-Version: 1.10.82
-Date: 2024-05-24
+Version: 1.10.83
+Date: 2024-05-27
 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]>

NAMESPACE

@@ -35,7 +35,7 @@ S3method(PCA,radiosonde)
 S3method(PCA,station)
 S3method(PCA,trajectory)
 S3method(WG,FT.day.t2m)
-S3method(WG,fw.day.precip)
+S3method(WG,fwmu.day.precip)
 S3method(WG,pca.day.t2m.precip)
 S3method(WG,station)
 S3method(aggregate,comb)
@@ -294,7 +294,7 @@ S3method(seasevol,station)
 S3method(season,default)
 S3method(sort,station)
 S3method(spell,default)
-S3method(spell,new)
+S3method(spell,old)
 S3method(spell,station)
 S3method(spell,test)
 S3method(station,default)
@@ -422,7 +422,7 @@ export(TGW)
 export(UTM2LatLon)
 export(WG)
 export(WG.FT.day.t2m)
-export(WG.fw.day.precip)
+export(WG.fwmu.day.precip)
 export(WG.pca.day.t2m.precip)
 export(WG.station)
 export(aggregate.area)
@@ -826,6 +826,7 @@ export(sp2np)
 export(sparseMproduct)
 export(spell)
 export(spell.default)
+export(spell.old)
 export(spell.station)
 export(src)
 export(sst.DNMI)
@@ -861,7 +862,7 @@ export(t2m.NCEP)
 export(t2m.NorESM.M)
 export(t2m.Pr)
 export(t2m.vul)
-export(test.WG.fw.day.precip)
+export(test.WG.fwmu.day.precip)
 export(test.aggregate.grid)
 export(test.ds.field)
 export(test.howsimilar)

R/addstations2nc.R

@@ -20,9 +20,8 @@
 #' @keywords utilities
 #'
 #' @examples 
-#'
-#' data(ferder)
-#' plot(anomaly(ferder))
+#' ncfiles <- list.files(path='~/Downloads/ecad2nc4',pattern='.ecad.nc',full.names=TRUE)
+#' nc4combine.stations(ncfiles,out=paste0('~/Downloads/',variable,'.ecad.nc'))
 #' 
 #' @export
 nc4add.stations <- function(x,fname,verbose=FALSE) {
@@ -136,7 +135,30 @@ nc4combine.stations <- function(fname,out='nc4combine.stations.nc',namelength=24
                           namelength,calendar,doi,namingauthority,processinglevel,creatortype,
                           creatoremail,institution,publishername,publisheremail,
                           publisherurl,project,insufficient)
-
 
   if (verbose) print('nc4combine.stations successfuly finished')
 }
+
+## This function copies information in netCDF files into one assuming the same time
+## coordinates - used for combining station data
+ncrcat <- function(fnames,ncout) {
+  file.copy(fnames[1],ncout)
+  ncid <- nc_open(ncout,write=TRUE)
+  stid <- ncvar_get(ncout,'stid')
+  for (fname in fnames) {
+    ncadd <- nc_open(fname)
+    for (var in names(ncadd$var)) {
+      x <- ncvar_get(ncadd,var)
+      scale <- try(ncatt_get(ncadd,var,'scale_factor'))
+      offset  <- try(ncatt_get(ncadd,var,'add_offset'))
+      if (!inherits(scale,'try-error') & !inherits(offset,'try-error'))
+        x <- (x - offset)/scale
+      ## Locations
+      d <- dim(x); if (is.null(d)) d <- length(d)
+      is <- ncvar_get(ncadd,'stid')
+      ncvar_put(ncout,var,x,start=c(),count=d)
+    }
+
+  }
+  nc_close(ncid)
+}

R/spell.R

@@ -54,8 +54,65 @@ spell <- function(x,threshold,...) UseMethod("spell")
 #' @exportS3Method
 #' @export spell.default
 spell.default <- function(x,threshold,upper=NULL,verbose=FALSE,...) {
+  ## Simpler algorithm for spell duration statistics
+  t0 <- Sys.time()
+  if (verbose) print(paste('spell.default - threshold=',threshold))
+  xgt <- coredata(x >= threshold)
+  nt <- length(x)
+  ## t is time index and y holds the length of spell length
+  t <- rep(NA,nt); y <- matrix(rep(NA,2*nt),nt,2)
+  t[1] <- index(x)[1]
+  ie <- 1; ii <- 1; L <- 0; 
+  if (verbose) print(table(xgt))
+  for (i in 1:length(xgt)) {
+    if (is.finite(x[i])) L <- L + 1 else L <- NA   ## Increment the spell length counter
+    ## Check for changed from above to below or the other way - store the length
+    ## of spell
+    if (verbose) print(paste('i=',i,'#event=',ie,'ii=',ii,year(x[i]),x[i],xgt[i],L))
+    if (xgt[i] != xgt[ii]) {
+      if (verbose) print('...')
+      y[ie,c(!xgt[i],xgt[i])] <- L
+      t[ie] <- index(x)[i]
+      if (is.na(L)) t[ie] <- NA
+      L <- 0       ## Reset spell length counter
+      ie <- ie + 1 ## Increment the event counter
+    }
+    ii <- i
+  }
+  y <- y[1:ie,]; t <- t[1:ie]; t[1] <- NA; t[ie] <- NA
+  good <- is.finite(t) 
+  y <- zoo(y[good,],order.by=as.Date(t)[good])
+
+  if (is.T(x)) {
+    attr(y,'variable') <-  c("warm","cold") 
+    attr(y,'longname') <-  c("duration of warm spells","duration of cold spells") 
+  } else {
+    attr(y,'variable') <-  c("wet","dry")
+    attr(y,'longname') <-  c("duration of wet spells","duration of dry spells") 
+  }
+  attr(y,'location') <- rep(loc(x),2)
+  attr(y,'station_id') <- rep(stid(x),2)
+  attr(y,'longitude') <- rep(lon(x),2)
+  attr(y,'latitude') <- rep(lat(y),2)
+  attr(y,'altitude') <- rep(alt(x),2)
+  attr(y,'unit') <- rep("days",2)
+  attr(y,'threshold') <- rep(threshold,2)
+  attr(y,'threshold.unit') <- rep(attr(x,'unit'),2)
+  # attr(y,'chksum') <- rep(chksum,2)
+  # attr(y,'uncredibly.high') <- t[ignoreh]
+  # attr(y,'uncredibly.low') <- t[ignorel]
+  # attr(y,'p.above') <- rep(sum(above)/length(above),2)
+  # attr(y,'interpolated.missing') <- index(x)[missing]
+  class(y) <- c("spell",class(x))
+  if (verbose) print(paste('spell.default: Time lapsed',Sys.time()-1))
+  invisible(y)
+}
+
+#' @exportS3Method
+#' @export spell.old
+spell.old <- function(x,threshold,upper=NULL,verbose=FALSE,...) {
 
-  if (verbose) print('spell.default')
+  if (verbose) print('spell.old')
 
   ## Deal with missing data
   missing <- (1:length(x))[!is.finite(x)]
@@ -219,7 +276,7 @@ spell.default <- function(x,threshold,upper=NULL,verbose=FALSE,...) {
 spell.station <-  function(x,threshold,upper=150,verbose=FALSE,...) {
   if (verbose) print('spell.station')
   if (!is.null(dim(x))) {
-    if (verbose) print('group of stations')
+    if ( verbose & (dim(x)[2]>1) ) print('group of stations')
     for (is in 1:dim(x)[2]) {
       y1 <- spell.default(subset(x,is=is),threshold=threshold,upper=upper,verbose=verbose,...)
       if (!inherits(y1,'spell')) {
@@ -459,64 +516,19 @@ GDD <- function(x,x0=10,na.rm=TRUE) {
   return(gdd)
 }
 
-## New version of spell that is more 'brute force' 
-## Set first and last estimate to NA as we don't know if the spells continue beyond
-## the data period. Also set estimates adjacent to NAs as NA to reduce potential
-## wrong estimates.
-#' @export
-spell.new <- function(x,threshold,upper=NULL,verbose=FALSE,...) {
-  if (verbose) {
-    print('spell.new')
-    t0 <- Sys.time()
-  }
-
-  if (!is.null(dim(x))) {
-  ia <- x >= threshold ## TRUE if x >= threshold
-  n <- length(x)
-  j <- 1    # number of event
-  A <- rep(NA,n); B <- A  # Length of events: A=above, B=below
-  k <- L
-  for (i in 1:n) {
-    if (ia[i]) {
-      A[j] <- A[j] + 1
-    } else {
-      B[j] <- B[j] + 1
-      ## increase the count every time x goes below threshold
-      j <- j+1
-    }
-    ## Remove the first and last estimates
-    A <- A[2:j-1] 
-    B <- B[2:j-1]
-
-  }
-  } else {
-    y <- apply(x,2,'spell.new')
-  }
-
-  y1 <- zoo(as.matrix(rep(NA,4),2,2),order.by=range(index(x)))
-  if (is.T(x)) {
-    attr(y1,'variable') <-  c("warm","cold") 
-    attr(y1,'longname') <-  c("duration of warm spells","duration of cold spells") 
-  } else {
-    attr(y1,'variable') <-  c("wet","dry")
-    attr(y1,'longname') <-  c("duration of wet spells","duration of dry spells") 
-  }
-  attr(y1,'location') <- loc(x)[is]
-  attr(y1,'station_id') <- stid(x)[is]
-  attr(y1,'longitude') <- lon(x)[is]
-  attr(y1,'latitude') <- lat(y)[is]
-  attr(y1,'altitude') <- alt(x)[is]
-  attr(y1,'unit') <- "days"
-  attr(y1,'threshold') <- rep(threshold,2)
-  attr(y1,'threshold.unit') <- rep(attr(x,'unit'),2)
-  attr(y1,'chksum') <- NA
-  attr(y1,'uncredibly.high') <- NA
-  attr(y1,'uncredibly.low') <- NA
-  attr(y1,'p.above') <- NA
-  attr(y1,'interpolated.missing') <- NA
-  class(y1) <- c("spell",class(x))
-  if (verbose) print(paste('Time taken is',Sys.time() - t0))
-}
-
 #' @export
-spell.test <- function() {}
+spell.test <- function(x,threshold=1,verbose=FALSE) {
+  t0 <- Sys.time()
+  S.old <- spell.old(x,threshold,verbose=verbose)
+  print(paste('spell.old: time lapsed',round(Sys.time()-t0,2),'s'))
+  t0 <- Sys.time()
+  S <- spell.default(x,threshold,verbose=verbose)
+  print(paste('spell.default: time lapsed',round(Sys.time()-t0,2),'s'))
+  par(mfrow=c(2,1),mar=c(2,4,2,1))
+  plot(zoo(subset(S.old,is=1)),lwd=3,col='grey70',
+       main=paste(loc(x),varid(subset(S,is=1))),ylab='Days',xlab='')
+  lines(zoo(subset(S,is=1)),col='red',lty=2)
+  plot(zoo(subset(S.old,is=2)),lwd=3,col='grey70',
+       main=paste(loc(x),varid(subset(S,is=2))),ylab='Days',xlab='')
+  lines(zoo(subset(S,is=2)),col='red',lty=2)
+}