R/rtools.R

In metno/esd: Climate analysis and empirical-statistical downscaling (ESD) package for monthly and daily data

#' Simple and handy functions
#'
#' \code{attrcp(x,y)} passes on attributes from \code{x} to \code{y} and
#' returns the \code{y} with new attributes.
#'
#' \code{ensemblemean} returns the ensemble mean for dsensemble objects. The argument \code{FUN} can also be
#' used to estimate other statistics, e.g \code{FUN='q95'} where
#' \code{q95=function(x) apply(x,1,quantile,probs=0.95)}
#'
#' \code{TGW} uses triangulation of pressure measurements to estimate geostrophic wind based on
#' Alexandersson et al. (1998), Glob. Atm. and Oce. Sys.
#'
#' \code{stand} gives a standardised time series.
#'
#' \code{zeros} counts the occurrence of zero values in a vector.
#'
#' \code{nv} returns the number of valid points.
#'
#' \code{missval} computes the percentage of missing data.
#'
#' \code{arec} compares the number of record-breaking events to the number of
#' events for a time series of iid data (\code{sum(1/1:n)})
#' 
#' \code{strstrip} strips off spaces before and after text in a string.
#' 
#' \code{as.decimal} converts between degree-minute-second into decimal value.
#' 
#' \code{cv} computes the coefficient of variation.
#' 
#' \code{nv} count the number of valid data points.
#' 
#' \code{q5}, \code{q95}, \code{q975} and \code{q995} are shortcuts to the 5\%, 95\%, 97.5\% and 99.5\% percentiles.
#' 
#' \code{trend.coef} and \code{trend.pval} return the coefficient and the p-value of the linear trend.
#' 
#' \code{exit} is a handy function for exiting the R session without saving.
#' 
#' \code{figlab} is a handy function for labelling figures (e.g. 'Figure 1')
#' 
#' \code{ndig} estimates the number of digits for round(x,ndig), e.g. in scales for plotting.
#' 
#' \code{factor2numeric} transforms a factor to a numeric object
#' 
#' \code{rmse} and \code{RMSE} calculate the root-mean-square error
#'
#' \code{leapyear} checks if a year is a leap year
#'
#' \code{distance2ocean} estimates the distance from a point to the ocean
#'
#' \code{lastyear} checks the availability of data of a station object and returns the last year(s) of valid data entries
#'
#' \code{difftime.month} returns the number of months between two dates
#' 
#' \code{lastelementrecord} checks last element of the record to see if they are the highest - a record
#' 
#' \code{round2magnitude} applied to a number, rounds up to the nearest order of magnitude
#' 
#' @aliases as.decimal nv cv q5 q95 q975 q995 filt
#' filt.default exit figlab ndig ensemblemean propchange stand rmse RMSE
#' firstyear lastyear eofvar test.num.predictors arec
#' arec.default arec.station lastrains lastdry lastelementrecord strstrip bin
#' factor2numeric zeros missval leapyear distance2ocean difftime.month round2magnitude
#' @seealso attrcp
#'
#' @importFrom stats quantile qgeom qpois dnorm filter 
#'
#' @param x A data.frame or a coredata zoo object.
#' @param na.rm If TRUE, remove NA's from data
#' @param type 'ma' for moving average (box-car), 'gauss' for Gaussian, 'binom'
#' for binomial filter,' parzen' for Parzen filter, 'hanning' for Hanning
#' filter, or 'welch' for Welch filter.
#' @param lowpass True for smoothing, otherwise the highpass results is
#' returned
#' @param triangle a group of three stations with sea-level pressure', e.g.
#' from ECA\&D.
#' @param nbins number of bins/categories
#' @return \item{as.decimal }{Decimal value} \item{trend.coef }{Linear trend
#' per decade}
#' @author A. Mezghani
#' @keywords rtools ~kwd2
#' @examples
#' 
#' ## Monthly mean temperature at 'Oslo - Blindern' station 
#' data(Oslo)
#' ## Compute the linear trend and the p-value on annual aggregated values 
#' tr <- trend.coef(coredata(annual(Oslo)))
#' pval <- trend.pval(coredata(annual(Oslo)))
#' \dontrun{
#' pp <- station(param='slp',cntr='Denmark',src='ecad')
#' wind <- TGW(subset(pp,is=c(1,3,10))
#' plot(wind)
#' ws <- sqrt(wind[,1]^2 + wind[,2]^2)
#' plot(ws)
#' hist(ws)
#' 
#' ## Estimate wind for a larger group of stations
#' wind <- geostrophicwind(pp,nmax=10)
#' u <- subset(wind,param='u')
#' v <- subset(wind,param='u')
#' ws <- sqrt(u^2+v^2)
#' ws <- attrcp(v,ws)
#' class(ws) <- class(v)
#' attr(ws,'variable')='windspeed'
#' attr(ws,'longname')='geostrophic wind speed'
#' map(ws,FUN='quantile',probs=0.98)
#' 
#' ## Test firstyears on HadCRUT4
#' if (!file.exists('~/Downloads/HadCRUT.4.6.0.0.median.nc')) {
#'   print('Download HadCRUT4')
#'   download.file('https://crudata.uea.ac.uk/cru/data/temperature/HadCRUT.4.6.0.0.median.nc',
#'                 dest='~/Downloads/HadCRUT.4.6.0.0.median.nc') 
#' }
#' 
#' Obs <- annual(retrieve('~/Downloads/HadCRUT.4.6.0.0.median.nc',param='temperature_anomaly'))
#' lons <- rep(lon(Obs),length(lat(Obs)))
#' lats <- sort(rep(lat(Obs),length(lon(Obs))))
#' fy <- firstyear(Obs)
#' map(subset(Obs,it=1))
#' points(lons[fy==1850],lats[fy==1850])
#' map(Obs,FUN='firstyear')
#' }
#'
#' ## Check for leap years
#' leapyear(1999)
#' leapyear(2000)
#'
#' ## Calculate the distance to the ocean from the point 10E/60N 
#' distance2ocean(10, 60)
#'
#' @export as.decimal
as.decimal <- function(x=NULL) {
    ## converts from degree min sec format to degrees ...
    ##x is in the form "49 deg 17' 38''"
  if (!is.null(x)) {
        deg <- as.numeric(substr(x,1,2)) 
        min <- as.numeric(substr(x,4,5))
        sec <- as.numeric(substr(x,7,8))     
        x <- deg + min/60 + sec/3600
    }
    return(x)
}

#' @export
nv <- function(x) sum(is.finite(x))

#' @export
ndig <- function(x) {
  i0 <- (x==0) & !is.finite(x)
  if (sum(i0)>0) x[i0] <- 1
  y <- -trunc(log(abs(x))/log(10))
  if (sum(i0)>0) y[i0] <- 0
  return(y)
}

#' @export
strstrip <- function(x) {
  if (is.na(x)) return(NA)
  if (is.factor(x)) x <- as.character(x)
  if (!is.character(x)) return(NA)
  while (substr(x,1,1)==' ') x <- substr(x,2,nchar(x))
  while (substr(x,nchar(x),nchar(x))==' ') x <- substr(x,1,nchar(x)-1)
  return(x)
}

#' @export
eofvar <- function(x) {
  if (inherits(x,c('eof','pca'))) {
    attr(x,'eigenvalues')^2/attr(x,'tot.var')*100
  } else {
    NULL
  }
}

#' @export
firstyear <- function(x,na.rm=FALSE,verbose=FALSE) {
  if (verbose) print('firstyear')
  yrs <- year(x)
  if (verbose) print(range(as.numeric(yrs)))
  if (is.null(dim(x))) y <- min(yrs[is.finite(x)]) else { 
    nv <- apply(x,2,'nv')
    y <- rep(NA,length(nv))
    ok <- (1:length(nv))[nv > 0]
    y[ok] <- apply(x[,ok],2,function(x,yrs=yrs) min(yrs[is.finite(x)]),yrs)
    #for (i in ok) y[i] <- min(yrs[is.finite(x[,i])])
  }
  y[!is.finite(y)] <- NA
  if (verbose) print(table(as.numeric(y)))
  return(y)
}

#' @export
lastyear <- function(x,na.rm=FALSE,verbose=FALSE) {
  if (verbose) print('lastyear')
  yrs <- year(x)
  if (verbose) print(range(as.numeric(yrs)))
  if (is.null(dim(x))) y <- max(yrs[is.finite(x)]) else { 
    nv <- apply(x,2,'nv')
    y <- rep(NA,length(nv))
    ok <- (1:length(nv))[nv > 0]
    y[ok] <- apply(x[,ok],2,function(x,yrs=yrs) max(yrs[is.finite(x)]),yrs)
  }     
  y[!is.finite(y)] <- NA
  if (verbose) print(table(as.numeric(y)))
  return(y)
}

#' @export difftime.month
difftime.month <- function(enddate, startdate, verbose=FALSE) {
  if(verbose) print("difftime.month")
  yend <- as.integer(strftime(enddate,"%Y"))
  ystart <- as.integer(strftime(startdate,"%Y"))
  mend <- as.integer(strftime(enddate,"%m"))
  mstart <- as.integer(strftime(startdate,"%m"))
  return((yend-ystart)*12 + mend-mstart)
}

## Iterate using n number of predictands in the downscaling and retrive the cross-val given the number of predictands
#' @export test.num.predictors
test.num.predictors <- function(x=NA,y=NA,nmax.x=6,nmin.x=3,nmax.y=4,nam.x='NA', nam.y.res='NA', nam.y='NA', nam.x.dom='NA',nam.t='NA',verbose=FALSE) {
  predictor_field <- x
  predictand_field <- y
  max_EOFs_predictor <- nmax.x
  min_EOFs_predictor <- nmin.x
  max_EOFs_predictand <- nmax.y
  min_EOFs_predictand <- 1
  if (verbose) cat('The input data should be a field. \n The field is iteratively downscaling with a different number of predictor EOFs. \n A dataframe with a summary of cross-validation correlation coefficients (Q2) and, if provided, meta-data is returned.\n')
  stopifnot (inherits(predictor_field,'field') & inherits(predictor_field,'zoo'))
  stopifnot (inherits(predictand_field,'field') & inherits(predictand_field,'zoo')) #add not eof not ds
  
  training<-setNames(data.frame(matrix(ncol =max_EOFs_predictand-min_EOFs_predictand+2, nrow = max_EOFs_predictor-min_EOFs_predictor+2), stringsAsFactors = FALSE), c(paste('y.EOF',seq(min_EOFs_predictand,max_EOFs_predictand)),"wQ2_lim0.2"))
  rownames(training)<- c('y R2:', paste(seq(min_EOFs_predictor,max_EOFs_predictor),'x.EOFs Q2:'))

  n_predictand_EOFs <- max_EOFs_predictand 
  if (verbose) cat('\n Downscaling',n_predictand_EOFs,'EOFs from the predictand field. \n',fill = TRUE)
  predictand <- EOF(predictand_field,n=max_EOFs_predictand)
  training[1,1:(ncol(training)-1)]<-round(eofvar(predictand),3)
  
  for (n_predictor_EOFs in rev(seq(min_EOFs_predictor,max_EOFs_predictor))) {
    if (verbose) cat('\n Using',n_predictor_EOFs,'EOFs from the coarse input as predictors.\n',fill = TRUE)
    predictor <- EOF(predictor_field,n=n_predictor_EOFs)
    ds_obj <- DS(predictand,predictor)
  
    #Since ortho one can reduce the DS model later
    xc=cor(attr(ds_obj,'evaluation'),attr(ds_obj,'evaluation'))
    crossvals<-diag(xc[seq(2,nrow(xc),2),seq(1,ncol(xc),2)])
    training[n_predictor_EOFs-1,1:(ncol(training)-1)]<-round(crossvals,3)
    training[n_predictor_EOFs-1,ncol(training)]<-round(sum(eofvar(predictand)[crossvals>0.2]/sum(eofvar(predictand))*crossvals[crossvals>0.2]),3)
  }
  training$y.name=nam.y
  training$x.name=nam.x
  training$y.res=nam.y.res
  training$t.name=nam.t
  training$dom.name=nam.x.dom
  tm<-cbind(setNames(data.frame(rownames(training)),c('names')), data.frame(training, row.names=NULL))
  tmn<-tm[,c(ncol(tm)-3, ncol(tm)-2, ncol(tm)-1,  ncol(tm),1:(ncol(tm)-4))]
  if (verbose) {cat('\n .')  
    print(tmn) }
  invisible(tmn)
}

## compute the percentage of missing data in x
#' @export
missval <- function(x) sum(is.na(coredata(x)))/length(coredata(x))

## compute the quantile 95% of x
#' @export
q95 <- function(x,na.rm=TRUE) quantile(x,probs=.95,na.rm=na.rm)

## compute the quantile 5% of x
#' @export
q5 <- function(x,na.rm=TRUE) quantile(x,probs=.05,na.rm=na.rm)

## compute the quantile 5% of x
#' @export
q995 <- function(x,na.rm=TRUE) quantile(x,probs=.995,na.rm=na.rm)

## compute the quantile 5% of x
#' @export
q975 <- function(x,na.rm=TRUE) quantile(x,probs=.975,na.rm=na.rm)

## count the number of valid data points
#' @export
nv <- function(x,...) sum(is.finite(x))

## Compute the coefficient of variation of x
#' @export
cv <- function(x,na.rm=TRUE) {sd(x,na.rm=na.rm)/mean(x,na.rm=na.rm)}

#' @export
stand <- function(x) (x - mean(x,na.rm=TRUE))/sd(x,na.rm=TRUE)

#' @export
zeros <- function(x) (sum(is.infinite(1/x)) > 0)

## Estimate the root-mean-squared-error
#' @export
rmse <- function(x,y,na.rm=TRUE) {
  z <- (x - y)^2
  z <- sum(z,na.rm=na.rm)/sum(is.finite(z),na.rm=na.rm)
  z <- sqrt(z)
  return(z)             
}


#' @export
RMSE <- function(x,y,...) return(rmse(x,y,...))

#' Wrap-around for lag.zoo to work on station and field objects
#'
#' @aliases lag.field
#'
#' @importFrom stats lag
#'
#' @export lag.station
lag.station <- function(x,...) {
  y <- lag(zoo(x),...)
  y <- attrcp(x,y)
  class(y) <- class(x)
  invisible(y)
}

#' @export lag.field
lag.field <- function(x,...) lag.station(x,...)

#' @export
exit <- function() q(save="no")

#' @export
filt <- function(x,n,type='ma',lowpass=TRUE) UseMethod("filt")

#' @exportS3Method
#' @export filt.default
filt.default <- function(x,n,type='ma',lowpass=TRUE) {
  
  # A number of different filters using different window
  # shapes.
  #
  # R.E. Benestad, July, 2002, met.no.
  #
  # ref: Press et al. (1989), Numerical Recipes in Pascal, pp. 466
  #library(ts)

  # Moving-average (box-car) filter
  ma.filt <- function(x,n) {
    y <- filter(x,rep(1,n)/n)
    y
  }

  # Gaussian filter with cut-off at 0.025 and 0.975 of the area.
  gauss.filt <- function(x,n) {
    i <- seq(0,qnorm(0.975),length=n/2)
    win <- dnorm(c(sort(-i),i))
    win <- win/sum(win)
    y <- filter(x,win)
    y
  }

  # Binomial filter
  binom.filt <- function(x,n) {
    win <- choose(n-1,0:(n-1))
    win <- win/max(win,na.rm=T)
    win[is.na(win)] <- 1
    win <- win/sum(win,na.rm=T)
    y <- filter(x,win)
    y
  }

  # Parzen filter (Press,et al. (1989))
  parzen.filt  <-  function(x,n) {
    j <- 0:(n-1)
    win <- 1 - abs((j - 0.5*(n-1))/(0.5*(n+1)))
    win <- win/sum(win)
    y <- filter(x,win)
    y
  }

  # Hanning filter (Press,et al. (1989))
  hanning.filt  <-  function(x,n) {
    j <- 0:(n-1)
    win <- 0.5*(1-cos(2*pi*j/(n-1)))
    win <- win/sum(win)
    y <- filter(x,win)
    y
  }

  # Welch filter (Press,et al. (1989))
  welch.filt  <-  function(x,n) {
    j <- 0:(n-1)
    win <- 1 - ((j - 0.5*(n-1))/(0.5*(n+1)))^2
    win <- win/sum(win)
    y <- filter(x,win)
    y
  }

  y <- coredata(x)
  z <- eval(parse(text=paste(type,'filt(y,n)',sep='.')))
  hp <- as.numeric(y - coredata(z))
  if (!is.null(dim(x))) dim(hp) <- dim(x)
  if (lowpass) coredata(x) <- coredata(z) else
               coredata(x) <- hp
  attr(x,'history') <- history.stamp(x)
  return(x)
}

#' @export
figlab <- function(x,xpos=0.001,ypos=0.001) {
  par(new=TRUE,fig=c(0,1,0,1),xaxt='n',yaxt='n',bty='n',mar=rep(0,4))
  plot(c(0,1),c(0,1),type='n')
  text(xpos,ypos,x,cex=1.2,pos=4,col='grey30')
}

#' @export
ensemblemean <- function(x,FUN='rowMeans') {
  if (inherits(x,'pca')) z <- as.station(x) else z <- x
  ## Estimate the ensemble mean
  zz <- unlist(lapply(coredata(z),FUN=FUN))
  zm <- matrix(zz,length(index(z[[1]])),length(z))
  zm <- zoo(zm,order.by=index(z[[1]]))
  zm <- as.station(zm,param=varid(z),unit=unit(z),
                   loc=unlist(lapply(z,loc)),lon=unlist(lapply(z,lon)),
                   lat=unlist(lapply(z,lat)),alt=unlist(lapply(z,alt)),
                   longname=attr(x,'longname'),aspect=attr(x,'aspect'),
                   info='Ensemble mean ESD')
  invisible(zm)
}

#' @export
propchange <- function(x,it0=c(1979,2013)) {
  z <- coredata(x)
  if (is.null(dim(z)))
      z <- 100*(z/mean(coredata(subset(x,it=it0)),na.rm=TRUE)) else
      z <- 100*t(t(z)/apply(coredata(subset(x,it=it0)),2,'mean',na.rm=TRUE))
  attributes(z) <- NULL
  z -> coredata(x)  
  x
}

#' @export
arec <- function(x,...) UseMethod("arec")

#' @exportS3Method
#' @export arec.default
arec.default <- function(x,...) {
  y <- length(records(x))/sum(1/(1:nv(x)))
  return(y)
}

#' @exportS3Method
arec.station <- function(x,...) {
  y <- unlist(lapply(records(x),length))/apply(x,2,function(x) sum(1/(1:nv(x))))
  return(y)
}

#' @export
lastrains <- function(x,x0=1,uptodate=TRUE,verbose=FALSE) {
  if (verbose) print('lastrains')
  ## Clean up missing values
  x <- x[is.finite(x)]
  y <- cumsum(rev(coredata(x)))
  z <- sum(y < x0,na.rm=TRUE)
  if (uptodate) if (Sys.Date() - end(x) > 1) z <- NA 
  return(z)
}

#' @export
lastdry <- function(x,x0=1,uptodate=TRUE,verbose=FALSE) {
  if (verbose) print('lastdry')
  ## Clean up missing values
  x <- x[is.finite(x)]
  y <- rev(coredata(x))
  z <- (1:length(y))[y < x0][1] - 1
  if (uptodate) if (Sys.Date() - end(x) > 1) z <- NA 
  return(z)
}

#' @export
lastelementrecord <- function(x,verbose=FALSE) {
  ## Checks last element of the record to see if they are the highest - a record
  if (verbose) print('lastelementrecord')
  ## If minimum, then multiply x with -1
  y <- coredata(x)
  nt <- length(index(x))
  if (length(dim(y)) == 2) {
    z <- rep(0,dim(y)[2])
    validlast <- is.finite(y[nt,])
    if (verbose) print(sum(validlast))
    if (sum(validlast)>1)
      z[validlast] <- apply(y[,validlast],2,function(x) if (x[length(x)] == max(x,na.rm=TRUE)) 1 else 0) else
    if (sum(validlast)>1) if (y[nt,validlast]==max(y[,validlast],na.rm=TRUE)) z <- 1    
  } else {
    z <- 0
    if(is.finite(y[nt])) 
      if (y[nt]==max(y,na.rm=TRUE)) z <- 1
  }
  return(z)
}

#' @export
bin <- function(x,nbins=5,labels=NULL,na.rm=TRUE) {
  if (na.rm) good <- is.finite(x) else good <- rep(TRUE,length(x))
  rank <- order(x[good])
  y <- (rank -1) %/% trunc(length(rank)/nbins) + 1
  z <- rep(NA,length(x))
  z[good] <- y
  if (is.null(labels)) names(z) <- quantile(x,probs = (1:nbins -0.5)/nbins)
  return(z)
}


#' @export
factor2numeric <- function(f) {
  if(!is.null(levels(f))) {return(as.numeric(levels(f))[f])
  } else return(as.numeric(f))
}

#' @export
leapyear <- function(years) {
  is.leap <- (years %% 4 == 0) & 
             ( years %% 100 != 0 | 
              (years %% 100 == 0 & years %% 400 == 0) )
  return(is.leap)
}

#' @export
distance2ocean <- function(lon1, lat1, dlon=c(-10,10), dlat=c(-5,5), 
                           delta=0.4) {
  ## regridding to a coarser resolution to avoid problems with Fjords
  data(etopo5, envir=environment())
  d <- NULL; lon.d <- NULL; lat.d <- NULL 
  for(j in seq_along(lon1)) {
    is <- list(lon=lon1[j]+dlon, lat=lat1[j]+dlat)
    h <- regrid(etopo5, is=list(lon=seq(min(is$lon),max(is$lon),delta), 
                                lat=seq(min(is$lat),max(is$lat),delta)))
    ocean <- which(h < -100, arr.ind=TRUE)
    if(h[which.min(abs(lon1[j]-lon(h))),which.min(abs(lat1[j]-lat(h)))]>0) {
      di <- distAB(lon1[j],lat1[j],lon(h)[ocean[,1]],lat(h)[ocean[,2]])*1E-3
      lon.dj <- lon(h)[ocean[,1]][which.min(di)]
      lat.dj <- lat(h)[ocean[,2]][which.min(di)]
      dj <- di[which.min(di)]
    } else {
      dj <- 0
      lon.dj <- lon1[j]
      lat.dj <- lat1[j]
    }
    d <- c(d, dj)
    lon.d <- c(lon.d, lon.dj)
    lat.d <- c(lat.d, lat.dj)
  }
  return(list("distance"=d, "lon"=lon.d, "lat"=lat.d))
}

#' @export
year2date <- function(y,toy='-01-01') {
  if (inherits(y,'zoo')) {
    index(y) <- as.Date(paste0(index(y),toy))
    return(y)
  } else if(is.numeric(y)) return(as.Date(paste0(y,toy)))
}

#' @export
yyyymmddhh2time <- function(y,option='POSIXct',verbose=FALSE) {
  stopifnot(is.character(y))
  ti <- paste(substr(y,1,4),substr(y,5,6),substr(y,7,8),sep='-')
  if (nchar(y) ==10)
    ti <- paste0(paste(ti,substr(y,9,10),sep=' '),':00:00')
  if (nchar(y) ==12)
    ti <- paste0(paste(ti,substr(y,11,12),sep=':'),':00')
  if (nchar(y) ==14)
    ti <- paste(ti,substr(y,13,14),sep=':')
  if (verbose) print(ti)
  if (option=='POSIXct') y <- as.POSIXct(ti) else
    y <- as.POSIXlt(ti)
  return(y)
}


## Function to check the time period covered
#' @export
getinterval <- function(file,param='tas',verbose=FALSE) {
  if (!file.exists(file)) return(NULL)
  ncid <- nc_open(file)
  time <- check.ncdf4(ncid,param=param,verbose=verbose)$time 
  nc_close(ncid)
  print(range(time$vdate))
  return(time$vdate)
}

## Function to split the CMIP file names into model, ssp and RIPF
#' @export
decipher <- function(x) {
  i <- gregexpr('_',x)[[1]]
  model <- substr(x,i[2]+1,i[3]-1)
  ssp <- substr(x,i[3]+1,i[4]-1)
  ripf <- substr(x,i[4]+1,i[4]+8)
  return(c(model,ssp,ripf))
}

# Function to round a number up to its order of magnitude 
#' @export
round2magnitude <- function(x) {
  10^(ceiling(log10(x)))
}