diff --git a/.Rbuildignore b/.Rbuildignore index 1c78466..268525c 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -5,3 +5,5 @@ man-roxygen comparison README.rst CONTRIBUTING.rst +^.*\.Rproj$ +^\.Rproj\.user$ diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..0a74739 --- /dev/null +++ b/.gitignore @@ -0,0 +1,6 @@ +.Rproj.user +.Rhistory +.RData +src/*.o +src/*.so +src/*.dll diff --git a/DESCRIPTION b/DESCRIPTION index d8534b9..48cc36d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -11,13 +11,16 @@ Depends: Imports: caTools (>= 1.0), methods, - Rcpp (>= 0.11.4) + Rcpp (>= 0.11.4), + SPEI Suggests: compiler, - RUnit + RUnit, + testthat Description: PCIC's implementation of Climdex routines for computation of extreme climate indices. License: GPL-3 URL: http://www.r-project.org LazyData: yes BugReports: https://github.com/pacificclimate/climdex.pcic/issues/ +RoxygenNote: 5.0.1 diff --git a/NAMESPACE b/NAMESPACE index 960b6ff..8df3ffa 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,39 +1,87 @@ +# Generated by roxygen2: do not edit by hand + +export(climdex.HI) +export(climdex.cc) +export(climdex.cc2) +export(climdex.cc6) +export(climdex.cd) export(climdex.cdd) +export(climdex.cfd) export(climdex.csdi) +export(climdex.csu) +export(climdex.cw) export(climdex.cwd) +export(climdex.ddeast) +export(climdex.ddnorth) +export(climdex.ddsouth) +export(climdex.ddwest) export(climdex.dtr) export(climdex.fd) +export(climdex.fg) +export(climdex.fg6bft) +export(climdex.fgcalm) +export(climdex.fxstorm) +export(climdex.fxx) export(climdex.get.available.indices) export(climdex.gsl) +export(climdex.hd17) export(climdex.id) +export(climdex.nsd) +export(climdex.nss) +export(climdex.nsx) export(climdex.prcptot) export(climdex.quantile) export(climdex.r10mm) export(climdex.r20mm) +export(climdex.r75p) +export(climdex.r75ptot) +export(climdex.r95p) export(climdex.r95ptot) +export(climdex.r99p) export(climdex.r99ptot) export(climdex.rnnmm) export(climdex.rx1day) export(climdex.rx5day) +export(climdex.sd) +export(climdex.sdd) export(climdex.sdii) +export(climdex.sdx) +export(climdex.spi3) +export(climdex.spi6) +export(climdex.ss) export(climdex.su) +export(climdex.sun_cloudy) +export(climdex.sun_relmean) +export(climdex.sun_sunny) +export(climdex.tmndaymax) +export(climdex.tmndaymin) export(climdex.tn10p) export(climdex.tn90p) export(climdex.tnn) +export(climdex.tnndaymax) +export(climdex.tnndaymin) +export(climdex.tnnp) export(climdex.tnx) export(climdex.tr) export(climdex.tx10p) export(climdex.tx90p) export(climdex.txn) +export(climdex.txndaymax) +export(climdex.txndaymin) +export(climdex.txnp) export(climdex.txx) +export(climdex.wd) export(climdex.wsdi) +export(climdex.ww) export(climdexInput.csv) export(climdexInput.raw) +export(eca.input) export(get.last.monthday.of.year) export(get.outofbase.quantiles) export(get.series.lengths.at.ends) export(growing.season.length) export(nday.consec.prec.max) +export(nday.consec.temp.mean) export(number.days.op.threshold) export(percent.days.op.threshold) export(select.blocks.gt.length) @@ -46,4 +94,5 @@ import(PCICt) import(Rcpp) import(caTools) import(methods) +importFrom(SPEI,spi) useDynLib(climdex.pcic) diff --git a/R/climdex.r b/R/climdex.r index cb022cb..49b00a6 100644 --- a/R/climdex.r +++ b/R/climdex.r @@ -43,7 +43,7 @@ tapply.fast <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) { ## Check that climdexInput data structure is valid. valid.climdexInput <- function(x) { temp.quantiles <- c(10, 90) - prec.quantiles <- c(95, 99) + prec.quantiles <- c(75, 95, 99) errors <- c() separate.base <- c(tmax=T, tmin=T, tavg=T, prec=F) @@ -57,9 +57,11 @@ valid.climdexInput <- function(x) { errors <- c(errors, "Data fields, dates, and date factors must all be of the same length") ## Check that namasks have columns for each of the variables - if(!all(c("annual", "monthly") %in% names(x@namasks)) || !all(present.data.vars %in% names(x@namasks$annual) & present.data.vars %in% names(x@namasks$monthly))) - errors <- c(errors, "NA mask for monthly and annual must contain data for all variables supplied.") - + if(!all(c("annual", "halfyear", "seasonal", "monthly") %in% names(x@namasks)) || + !all(present.data.vars %in% names(x@namasks$annual) & present.data.vars %in% names(x@namasks$halfyear) & + present.data.vars %in% names(x@namasks$seasonal) & present.data.vars %in% names(x@namasks$monthly))) + errors <- c(errors, "NA mask for monthly, seasonal, halfyear and annual must contain data for all variables supplied.") + ## Check that appropriate thresholds are present. need.base.data <- get.num.days.in.range(x@dates, x@base.range) > 0 errors <- do.call(c, c(list(errors), lapply(intersect(present.data.vars, c("tmax", "tmin", "prec")), function(n) { @@ -81,11 +83,13 @@ valid.climdexInput <- function(x) { } ## Class definition declaration -#' climdexInput +#' @title climdexInput #' +#' @description #' The climdexInput class contains all the data necessary to compute the #' climdex indices. #' +#' @details #' The \code{climdexInput} class consists of all the data necessary to compute #' the climdex indices. Users will not need to modify any of the slots in this #' class. That being said, users may want or need to repurpose this data for @@ -136,6 +140,11 @@ valid.climdexInput <- function(x) { #' length is the growing season starting in the beginning of July of the year #' indicated, running to the end of June of the following year. #' +#' The \code{max.missing.days} slot is a vector consisting of 'annual' +#' (the number of days that can be missing in a year) and 'monthly' (the +#' number of days that can be missing in a month. If one month in a year fails +#' the test, the corresponding year will be omitted. +#' #' @name climdexInput #' @aliases climdexInput-class #' @docType class @@ -274,7 +283,14 @@ get.num.days.in.range <- function(x, date.range) { ## Check that arguments to climdexInput.raw et al are complete enough and valid enough. -check.basic.argument.validity <- function(tmax, tmin, prec, tmax.dates, tmin.dates, prec.dates, base.range=c(1961, 1990), n=5, tavg=NULL, tavg.dates=NULL) { +#check.basic.argument.validity <- function(tmax, tmin, prec, tmax.dates, tmin.dates, prec.dates, base.range=c(1961, 1990), + # n=5, tavg=NULL, tavg.dates=NULL) { +check.basic.argument.validity <- function(tmax, tmin, prec, snow, snow_new, wind, wind_gust, wind_dir, cloud, sun, sun_rel, + tmax.dates, tmin.dates, prec.dates, snow.dates, snow_new.dates, wind.dates, + wind_gust.dates, + wind_dir.dates, cloud.dates, sun.dates, sun_rel.dates,base.range=c(1961, 1990), n=5, tavg=NULL, tavg.dates=NULL) { + + check.var <- function(var, var.dates, var.name) { if(is.null(var) != is.null(var.dates)) stop(paste("If passing in", var, ", must pass in", var, "dates too..")) @@ -290,8 +306,21 @@ check.basic.argument.validity <- function(tmax, tmin, prec, tmax.dates, tmin.dat check.var(tmin, tmin.dates, "tmin") check.var(tavg, tavg.dates, "tavg") check.var(prec, prec.dates, "prec") - - if(all(c(is.null(tmax), is.null(tmin), is.null(prec), is.null(tavg)))) + check.var(snow, snow.dates, "snow") + check.var(snow_new, snow_new.dates, "snow_new") + check.var(wind, wind.dates, "wind") + check.var(wind_gust, wind_gust.dates, "wind_gust") + check.var(wind_dir, wind_dir.dates, "wind_dir") + check.var(cloud, cloud.dates, "cloud") + check.var(sun, sun.dates, "sun") + check.var(sun_rel, sun_rel.dates, "sun_rel") + + if(all(c(is.null(tmax), is.null(tmin), is.null(tavg), + is.null(prec), + is.null(snow), is.null(snow_new), + is.null(wind), is.null(wind_gust), is.null(wind_dir), + is.null(cloud), + is.null(sun),is.null(sun_rel)))) stop("Must supply at least one variable to calculate indices upon.") if(!(length(base.range) == 2 && is.numeric(base.range))) @@ -321,8 +350,8 @@ check.quantile.validity <- function(quantiles, present.vars, days.in.base) { if(any(days.in.base > 0) && !all(sapply(quantiles[names(quantiles) %in% intersect(intersect(present.vars, c("tmax", "tmin")), names(days.in.base)[days.in.base > 0])], function(x) { "inbase" %in% names(x) && all(c("q10", "q90") %in% names(x$inbase)) }))) stop("Temperature in-base quantiles must contain 10th and 90th percentiles.\n") - if("prec" %in% names(quantiles) && !all(c("q95", "q99") %in% names(quantiles$prec))) - stop("Precipitation quantiles must contain 95th and 99th percentiles.\n") + if("prec" %in% names(quantiles) && !all(c("q75", "q95", "q99") %in% names(quantiles$prec))) + stop("Precipitation quantiles must contain 75th, 95th and 99th percentiles.\n") } get.temp.var.quantiles <- function(filled.data, date.series, bs.date.series, qtiles, bs.date.range, n, in.base=FALSE, min.base.data.fraction.present=0.1) { @@ -334,7 +363,7 @@ get.temp.var.quantiles <- function(filled.data, date.series, bs.date.series, qti return(list(outbase=zhang.running.qtile(base.data, dates.base=bs.date.series, qtiles=qtiles, bootstrap.range=bs.date.range, n=n, min.fraction.present=min.base.data.fraction.present))) } -get.prec.var.quantiles <- function(filled.prec, date.series, bs.date.range, qtiles=c(0.95, 0.99)) { +get.prec.var.quantiles <- function(filled.prec, date.series, bs.date.range, qtiles=c(0.75, 0.95, 0.99)) { wet.days <- !(is.na(filled.prec) | filled.prec < 1) inset <- date.series >= bs.date.range[1] & date.series <= bs.date.range[2] & !is.na(filled.prec) & wet.days pq <- quantile(filled.prec[inset], qtiles, type=8) @@ -342,16 +371,18 @@ get.prec.var.quantiles <- function(filled.prec, date.series, bs.date.range, qtil return(pq) } -#' Method for getting threshold quantiles for use in computing indices +#' @title Method for getting threshold quantiles for use in computing indices #' +#' @description #' This function creates threshold quantiles for use with climdexInput.raw #' or climdexInput.csv. #' +#' @details #' This function takes input climate data at daily resolution, and produces as #' output a set of threshold quantiles. This data structure can then be passed #' to climdexInput.raw or climdexInput.csv. #' -#' For more details on arguments, see \code{link{climdexInput.raw}}. +#' For more details on arguments, see \code{\link{climdexInput.raw}}. #' #' @seealso \code{\link{climdex.pcic-package}}, \code{\link{climdexInput.raw}}. #' @references \url{http://etccdi.pacificclimate.org/list_27_indices.shtml} @@ -389,9 +420,12 @@ get.prec.var.quantiles <- function(filled.prec, date.series, bs.date.range, qtil #' tmax.dates, tmin.dates, prec.dates, base.range=c(1971, 2000)) #' #' @export -get.outofbase.quantiles <- function(tmax=NULL, tmin=NULL, prec=NULL, tmax.dates=NULL, tmin.dates=NULL, prec.dates=NULL, base.range=c(1961, 1990), n=5, temp.qtiles=c(0.10, 0.90), prec.qtiles=c(0.95, 0.99), min.base.data.fraction.present=0.1) { +get.outofbase.quantiles <- function(tmax=NULL, tmin=NULL, prec=NULL, tmax.dates=NULL, tmin.dates=NULL, prec.dates=NULL, + base.range=c(1961, 1990), n=5, temp.qtiles=c(0.10, 0.90), prec.qtiles=c(0.75, 0.95, 0.99), + min.base.data.fraction.present=0.1) { days.threshold <- 359 check.basic.argument.validity(tmax, tmin, prec, tmax.dates, tmin.dates, prec.dates, base.range, n) + #check.basic.argument.validity(tmax, tmin, tavg, prec, tmax.dates, tmin.dates, tavg.dates, prec.dates, base.range, n) d.list <- list(tmin.dates, tmax.dates, prec.dates) all.dates <- do.call(c, d.list[!sapply(d.list, is.null)]) @@ -418,6 +452,7 @@ get.outofbase.quantiles <- function(tmax=NULL, tmin=NULL, prec=NULL, tmax.dates= filled.tmin <- create.filled.series(tmin, trunc(tmin.dates, "days"), date.series) quantiles$tmin <- get.temp.var.quantiles(filled.tmin, date.series, bs.date.series, temp.qtiles, bs.date.range, n) } + if(!is.null(prec)) { if(get.num.days.in.range(prec.dates, bs.date.range) <= days.threshold) @@ -428,11 +463,13 @@ get.outofbase.quantiles <- function(tmax=NULL, tmin=NULL, prec=NULL, tmax.dates= return(quantiles) } -#' Method for creating climdexInput object from vectors of data +#' @title Method for creating climdexInput object from vectors of data #' +#' @description #' This function creates a climdexInput object from data already ingested into #' R. #' +#' @details #' This function takes input climate data at daily resolution, and produces as #' output a ClimdexInput data structure. This data structure can then be passed #' to any of the routines used to compute the Climdex indices. The indices @@ -527,16 +564,62 @@ get.outofbase.quantiles <- function(tmax=NULL, tmin=NULL, prec=NULL, tmax.dates= #' tmax.dates, tmin.dates, prec.dates, base.range=c(1971, 2000)) #' #' @export -climdexInput.raw <- function(tmax=NULL, tmin=NULL, prec=NULL, tmax.dates=NULL, tmin.dates=NULL, prec.dates=NULL, +# climdexInput.raw <- function(tmax=NULL, tmax.dates=NULL, +# tmin=NULL, tmin.dates=NULL, +# tavg=NULL, tavg.dates=NULL, +# prec=NULL, prec.dates=NULL, +# quantiles=NULL, temp.qtiles=c(0.10, 0.90), prec.qtiles=c(0.75, 0.95, 0.99), +# base.range=c(1961, 1990), n=5, northern.hemisphere=TRUE, +# max.missing.days=c(annual=15, halfyear=10, seasonal=8, monthly=3), +# min.base.data.fraction.present=0.1) { +climdexInput.raw <- function(tmax=NULL, tmax.dates=NULL, + tmin=NULL, tmin.dates=NULL, + tavg=NULL, tavg.dates=NULL, + prec=NULL, prec.dates=NULL, + snow=NULL, snow.dates=NULL, + snow_new=NULL, snow_new.dates=NULL, + wind=NULL, wind.dates=NULL, + wind_gust=NULL, wind_gust.dates=NULL, + wind_dir=NULL, wind_dir.dates=NULL, + cloud=NULL, cloud.dates=NULL, + sun=NULL, sun.dates=NULL, + sun_rel=NULL, sun_rel.dates=NULL, + quantiles=NULL, temp.qtiles=c(0.10, 0.90), prec.qtiles=c(0.75, 0.95, 0.99), base.range=c(1961, 1990), n=5, northern.hemisphere=TRUE, - tavg=NULL, tavg.dates=NULL, quantiles=NULL, temp.qtiles=c(0.10, 0.90), prec.qtiles=c(0.95, 0.99), max.missing.days=c(annual=15, monthly=3), min.base.data.fraction.present=0.1) { + max.missing.days=c(annual=15, halfyear=10, seasonal=8, monthly=3), + min.base.data.fraction.present=0.1) { + ## Make sure all of these arguments are valid... - check.basic.argument.validity(tmax, tmin, prec, tmax.dates, tmin.dates, prec.dates, base.range, n, tavg, tavg.dates) +# check.basic.argument.validity(tmax, tmin, prec, tmax.dates, tmin.dates, prec.dates, base.range, n, tavg, tavg.dates) +# + check.basic.argument.validity(tmax=tmax, tmax.dates=tmax.dates, + tmin=tmin, tmin.dates=tmin.dates, + tavg=tavg, tavg.dates=tavg.dates, + prec=prec, prec.dates=prec.dates, + snow=snow, snow.dates=snow.dates, + snow_new=snow_new, snow_new.dates=snow_new.dates, + wind=wind, wind.dates=wind.dates, + wind_gust=wind_gust, wind_gust.dates=wind_gust.dates, + wind_dir=wind_dir, wind_dir.dates=wind_dir.dates, + cloud=cloud, cloud.dates=cloud.dates, + sun=sun, sun.dates=sun.dates, + sun_rel=sun_rel, sun_rel.dates=sun_rel.dates, + base.range=base.range, + n=n) + + stopifnot(length(max.missing.days) == 4 && all(c("annual", "halfyear", "seasonal", "monthly") %in% names(max.missing.days))) - stopifnot(length(max.missing.days) == 2 && all(c("annual", "monthly") %in% names(max.missing.days))) stopifnot(is.numeric(min.base.data.fraction.present) && length(min.base.data.fraction.present) == 1) - d.list <- list(tmin.dates, tmax.dates, prec.dates, tavg.dates) +# d.list <- list(tmin.dates, tmax.dates, prec.dates, tavg.dates) +# + d.list <- list(tmin.dates, tmax.dates, tavg.dates, + prec.dates, + snow.dates, snow_new.dates, + wind.dates, wind_gust.dates, wind_dir.dates, + cloud.dates, + sun.dates, sun_rel.dates) + all.dates <- do.call(c, d.list[!sapply(d.list, is.null)]) last.day.of.year <- get.last.monthday.of.year(all.dates) cal <- attr(all.dates, "cal") @@ -551,12 +634,38 @@ climdexInput.raw <- function(tmax=NULL, tmin=NULL, prec=NULL, tmax.dates=NULL, t jdays <- get.jdays.replaced.feb29(get.jdays(date.series)) ## Factors for dividing data up - date.factors <- list(annual=factor(format(date.series, format="%Y", tz="GMT")), monthly=factor(format(date.series, format="%Y-%m", tz="GMT"))) - + date.months <- as.numeric(format(date.series, format="%m", tz="GMT")) + date.years <- as.numeric(format(date.series, format="%Y", tz="GMT")) + # get factors for seasons + # Winter month D of prev year and JF of next year belong together, Year belongs to Jan => increase year of prev Dec by 1 + seas.years <- date.years + seas.seas <- date.months %/% 3 + 1 + seas.idx <- which(seas.seas == 5) + seas.years[seas.idx] <- seas.years[seas.idx]+1 + seas.seas[seas.idx] <- 1 + # get factors for half years (winter (ONDJFM) & summer (APJJAS)) + # winter months OND of prev year and JFM of next year belong together, Year belongs to Jan => increase year of prev OND by 1 + half.years <- date.years + half.half <- (date.months+2) %/% 6 + 1 + half.idx <- which(half.half == 3) + half.years[half.idx] <- half.years[half.idx]+1 + half.half[half.idx] <- 1 + + # set up date.factors list + date.factors <- list(annual=factor(format(date.series, format="%Y", tz="GMT")), + halfyear=factor(paste(half.years,half.half,sep="-")), + seasonal=factor(paste(seas.years,seas.seas,sep="-")), + monthly=factor(format(date.series, format="%Y-%m", tz="GMT"))) + ## Filled data... - var.list <- c("tmax", "tmin", "prec", "tavg") +# var.list <- c("tmax", "tmin", "prec", "tavg") + var.list <- c("tmax", "tmin", "tavg", "prec", "snow", "snow_new", "wind", "wind_gust", "wind_dir", + "cloud", "sun", "sun_rel") + present.var.list <- var.list[sapply(var.list, function(x) !is.null(get(x)))] - filled.list <- sapply(present.var.list, function(x) { return(create.filled.series(get(x), trunc(get(paste(x, "dates", sep="."))), date.series)) }, simplify=FALSE) + + filled.list <- sapply(present.var.list, function(x) { + return(create.filled.series(get(x), trunc(get(paste(x, "dates", sep="."))), date.series)) }, simplify=FALSE) if(is.null(tavg) && !is.null(tmin) && !is.null(tmax)) filled.list$tavg <- (filled.list$tmax + filled.list$tmin) / 2 @@ -573,9 +682,11 @@ climdexInput.raw <- function(tmax=NULL, tmin=NULL, prec=NULL, tmax.dates=NULL, t have.quantiles <- all(present.var.list %in% names(quantiles)) ## NA masks - namasks <- list(annual=lapply(filled.list, get.na.mask, date.factors$annual, max.missing.days['annual']), monthly=lapply(filled.list, get.na.mask, date.factors$monthly, max.missing.days['monthly'])) - namasks$annual <- lapply(names(namasks$annual), function(v) { d <- namasks$annual[[v]] * as.numeric(tapply(namasks$monthly[[v]], rep(seq_along(namasks$annual[[v]]), each=12), prod)); dimnames(d) <- dim(d) <- NULL; d }) - names(namasks$annual) <- names(namasks$monthly) + ## NA masks + namasks <- list(annual=lapply(filled.list, get.na.mask, date.factors$annual, max.missing.days['annual']), + halfyear=lapply(filled.list, get.na.mask, date.factors$halfyear, max.missing.days['halfyear']), + seasonal=lapply(filled.list, get.na.mask, date.factors$seasonal, max.missing.days['seasonal']), + monthly=lapply(filled.list, get.na.mask, date.factors$monthly, max.missing.days['monthly'])) ## Pad data passed as base if we're missing endpoints... if(!have.quantiles) { @@ -590,14 +701,16 @@ climdexInput.raw <- function(tmax=NULL, tmin=NULL, prec=NULL, tmax.dates=NULL, t } else { quantiles <- as.environment(quantiles) } - + return(new("climdexInput", data=filled.list, quantiles=quantiles, namasks=namasks, dates=date.series, jdays=jdays, base.range=bs.date.range, date.factors=date.factors, northern.hemisphere=northern.hemisphere, max.missing.days=max.missing.days)) } -#' Method for creating climdexInput object from CSV files +#' @title Method for creating climdexInput object from CSV files #' +#' @description #' This function creates a climdexInput object from data in CSV files. #' +#' @details #' This function takes input climate data in CSV files at daily resolution, #' and produces as output a ClimdexInput data structure. This data structure #' can then be passed to any of the routines used to compute the Climdex @@ -628,7 +741,7 @@ climdexInput.raw <- function(tmax=NULL, tmin=NULL, prec=NULL, tmax.dates=NULL, t #' of names consisting of the columns to be concatenated together with spaces. #' The \code{format} item is a date format as taken by \code{strptime}. #' -#' For more details on arguments, see \code{link{climdexInput.raw}}. +#' For more details on arguments, see \code{\link{climdexInput.raw}}. #' #' @seealso \code{\link{climdex.pcic-package}}, \code{\link{climdexInput.raw}}. #' @references \url{http://etccdi.pacificclimate.org/list_27_indices.shtml} @@ -638,6 +751,14 @@ climdexInput.raw <- function(tmax=NULL, tmin=NULL, prec=NULL, tmax.dates=NULL, t #' @param tmin.file Name of file containing daily minimum temperature data. #' @param prec.file Name of file containing daily total precipitation data. #' @param tavg.file Name of file containing daily mean temperature data. +#' @param snow.file Name of file containing daily mean snow height data. +#' @param snow_new.file Name of file containing daily mean new snow height data. +#' @param wind.file Name of file containing daily mean wind speed data. +#' @param wind_gust.file Name of file containing daily wind gust data. +#' @param wind_dir.file Name of file containing daily mean wind direction data. +#' @param cloud.file Name of file containing daily mean cloud cover data. +#' @param sun.file Name of file containing daily mean sunshine duration data. +#' @param sun_rel.file Name of file containing daily mean relative sunshine duration data. #' @param data.columns Column names for tmin, tmax, and prec data. #' @param date.types Column names for tmin, tmax, and prec data (see notes). #' @param na.strings Strings used for NA values; passed to @@ -659,10 +780,29 @@ climdexInput.raw <- function(tmax=NULL, tmin=NULL, prec=NULL, tmax.dates=NULL, t #' prec.filename, date.types=list(list(fields=c("date"), format="%Y-%m-%d")))} #' #' @export -climdexInput.csv <- function(tmax.file=NULL, tmin.file=NULL, prec.file=NULL, - data.columns=list(tmin="tmin", tmax="tmax", prec="prec"), base.range=c(1961, 1990), - na.strings=NULL, cal="gregorian", date.types=NULL, n=5, northern.hemisphere=TRUE, - tavg.file=NULL, quantiles=NULL, temp.qtiles=c(0.10, 0.90), prec.qtiles=c(0.95, 0.99), max.missing.days=c(annual=15, monthly=3), min.base.data.fraction.present=0.1) { +# climdexInput.csv <- function(tmax.file=NULL, tmin.file=NULL, prec.file=NULL, +# data.columns=list(tmin="tmin", tmax="tmax", prec="prec"), base.range=c(1961, 1990), +# na.strings=NULL, cal="gregorian", date.types=NULL, n=5, northern.hemisphere=TRUE, +# tavg.file=NULL, quantiles=NULL, temp.qtiles=c(0.10, 0.90), prec.qtiles=c(0.75, 0.95, 0.99), +# max.missing.days=c(annual=15, monthly=3), min.base.data.fraction.present=0.1) { +climdexInput.csv <- function(tmax.file=NULL, tmin.file=NULL, tavg.file=NULL, prec.file=NULL, + snow.file=NULL, snow_new.file=NULL, + wind.file=NULL, wind_gust.file=NULL, wind_dir.file=NULL, + cloud.file=NULL, + sun.file=NULL, sun_rel.file=NULL, + data.columns=list(tmin="tmin", tmax="tmax", tavg="tavg", + prec="prec", + snow="snow", snow_new="snow_new", + wind="wind", wind_gust="wind_gust", wind_dir="wind_dir", + cloud="cloud", + sun="sun", sun_rel="sun_rel"), + base.range=c(1961, 1990), + na.strings=NULL, cal="gregorian", + date.types=NULL, n=5, northern.hemisphere=TRUE, + quantiles=NULL, temp.qtiles=c(0.10, 0.90), prec.qtiles=c(0.75, 0.95, 0.99), + max.missing.days=c(annual=15, halfyear=10, seasonal=8, monthly=3), + min.base.data.fraction.present=0.1) { + get.and.check.data <- function(fn, datacol) { if(!is.null(fn)) { dat <- read.csv(fn, na.strings=na.strings) @@ -683,8 +823,34 @@ climdexInput.csv <- function(tmax.file=NULL, tmin.file=NULL, prec.file=NULL, tmax <- get.and.check.data(tmax.file, data.columns$tmax) tavg <- get.and.check.data(tavg.file, data.columns$tavg) prec <- get.and.check.data(prec.file, data.columns$prec) - - return(climdexInput.raw(tmax=tmax$dat, tmin=tmin$dat, prec=prec$dat, tmax.dates=tmax$dates, tmin.dates=tmin$dates, prec.dates=prec$dates, base.range=base.range, n=n, northern.hemisphere=northern.hemisphere, tavg=tavg$dat, tavg.dates=tavg$dates, quantiles=quantiles, temp.qtiles=temp.qtiles, prec.qtiles=prec.qtiles, max.missing.days=max.missing.days, min.base.data.fraction.present=min.base.data.fraction.present)) + snow <- get.and.check.data(snow.file, data.columns$snow) + snow_new <- get.and.check.data(snow_new.file, data.columns$snow_new) + wind <- get.and.check.data(wind.file, data.columns$wind) + wind_gust <- get.and.check.data(wind_gust.file, data.columns$wind_gust) + wind_dir <- get.and.check.data(wind_dir.file, data.columns$wind_dir) + cloud <- get.and.check.data(cloud.file, data.columns$cloud) + sun <- get.and.check.data(sun.file, data.columns$sun) + sun_rel <- get.and.check.data(sun_rel.file, data.columns$sun_rel) +# + # return(climdexInput.raw(tmax=tmax$dat, tmin=tmin$dat, prec=prec$dat, tmax.dates=tmax$dates, tmin.dates=tmin$dates, prec.dates=prec$dates, + # base.range=base.range, n=n, northern.hemisphere=northern.hemisphere, tavg=tavg$dat, tavg.dates=tavg$dates, quantiles=quantiles, + # temp.qtiles=temp.qtiles, prec.qtiles=prec.qtiles, max.missing.days=max.missing.days, min.base.data.fraction.present=min.base.data.fraction.present)) + + return(climdexInput.raw(tmax=tmax$dat, tmax.dates=tmax$dates, + tmin=tmin$dat, tmin.dates=tmin$dates, + tavg=tavg$dat, tavg.dates=tavg$dates, + prec=prec$dat, prec.dates=prec$dates, + snow=snow$dat, snow.dates=snow$dates, + snow_new=snow_new$dat, snow_new.dates=snow_new$dates, + wind=wind$dat, wind.dates=wind$dates, + wind_gust=wind_gust$dat, wind_gust.dates=wind_gust$dates, + wind_dir=wind_dir$dat, wind_dir.dates=wind_dir$dates, + cloud=cloud$dat, cloud.dates=cloud$dates, + sun=sun$dat, sun.dates=sun$dates, + sun_rel=sun_rel$dat, sun_rel.dates=sun_rel$dates, + base.range=base.range, n=n, northern.hemisphere=northern.hemisphere, + quantiles=quantiles, temp.qtiles=temp.qtiles, prec.qtiles=prec.qtiles, + max.missing.days=max.missing.days, min.base.data.fraction.present=min.base.data.fraction.present)) } #' Frost Days @@ -704,7 +870,9 @@ climdexInput.csv <- function(tmax.file=NULL, tmin.file=NULL, prec.file=NULL, #' @template get_generic_example #' #' @export -climdex.fd <- function(ci) { stopifnot(!is.null(ci@data$tmin)); return(number.days.op.threshold(ci@data$tmin, ci@date.factors$annual, 0, "<") * ci@namasks$annual$tmin) } +climdex.fd <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$tmin)); + return(number.days.op.threshold(ci@data$tmin, ci@date.factors[[match.arg(freq)]], 0, "<") * ci@namasks[[match.arg(freq)]]$tmin) } #' Summer Days #' @@ -722,7 +890,9 @@ climdex.fd <- function(ci) { stopifnot(!is.null(ci@data$tmin)); return(number.da #' @template get_generic_example #' #' @export -climdex.su <- function(ci) { stopifnot(!is.null(ci@data$tmax)); return(number.days.op.threshold(ci@data$tmax, ci@date.factors$annual, 25, ">") * ci@namasks$annual$tmax) } +climdex.su <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$tmax)); + return(number.days.op.threshold(ci@data$tmax, ci@date.factors[[match.arg(freq)]], 25, ">") * ci@namasks[[match.arg(freq)]]$tmax) } #' Icing Days #' @@ -740,7 +910,9 @@ climdex.su <- function(ci) { stopifnot(!is.null(ci@data$tmax)); return(number.da #' @template get_generic_example #' #' @export -climdex.id <- function(ci) { stopifnot(!is.null(ci@data$tmax)); return(number.days.op.threshold(ci@data$tmax, ci@date.factors$annual, 0, "<") * ci@namasks$annual$tmax) } +climdex.id <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$tmax)); + return(number.days.op.threshold(ci@data$tmax, ci@date.factors[[match.arg(freq)]], 0, "<") * ci@namasks[[match.arg(freq)]]$tmax) } #' Tropical Nights #' @@ -758,12 +930,16 @@ climdex.id <- function(ci) { stopifnot(!is.null(ci@data$tmax)); return(number.da #' @template get_generic_example #' #' @export -climdex.tr <- function(ci) { stopifnot(!is.null(ci@data$tmin)); return(number.days.op.threshold(ci@data$tmin, ci@date.factors$annual, 20, ">") * ci@namasks$annual$tmin) } +climdex.tr <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$tmin)); + return(number.days.op.threshold(ci@data$tmin, ci@date.factors[[match.arg(freq)]], 20, ">") * ci@namasks[[match.arg(freq)]]$tmin) } -#' Growing Season Length +#' @title Growing Season Length #' +#' @description #' This function computes the growing season length (GSL) given the input. #' +#' @details #' This function takes a climdexInput object as input and computes the growing #' season length based on this data. #' @@ -847,7 +1023,9 @@ climdex.gsl <- function(ci, gsl.mode=c("GSL", "GSL_first", "GSL_max", "GSL_sum") #' @template get_generic_example #' #' @export -climdex.txx <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null(ci@data$tmax)); return(suppressWarnings(tapply.fast(ci@data$tmax, ci@date.factors[[match.arg(freq)]], max, na.rm=TRUE)) * ci@namasks[[match.arg(freq)]]$tmax) } +climdex.txx <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$tmax)); + return(suppressWarnings(tapply.fast(ci@data$tmax, ci@date.factors[[match.arg(freq)]], max, na.rm=TRUE)) * ci@namasks[[match.arg(freq)]]$tmax) } #' Monthly Maximum of Daily Minimum Temperature #' @@ -865,7 +1043,9 @@ climdex.txx <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null(ci #' @template get_generic_example #' #' @export -climdex.tnx <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null(ci@data$tmin)); return(suppressWarnings(tapply.fast(ci@data$tmin, ci@date.factors[[match.arg(freq)]], max, na.rm=TRUE)) * ci@namasks[[match.arg(freq)]]$tmin) } +climdex.tnx <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$tmin)); + return(suppressWarnings(tapply.fast(ci@data$tmin, ci@date.factors[[match.arg(freq)]], max, na.rm=TRUE)) * ci@namasks[[match.arg(freq)]]$tmin) } #' Monthly Minimum of Daily Maximum Temperature #' @@ -883,7 +1063,9 @@ climdex.tnx <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null(ci #' @template get_generic_example #' #' @export -climdex.txn <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null(ci@data$tmax)); return(suppressWarnings(tapply.fast(ci@data$tmax, ci@date.factors[[match.arg(freq)]], min, na.rm=TRUE)) * ci@namasks[[match.arg(freq)]]$tmax) } +climdex.txn <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$tmax)); + return(suppressWarnings(tapply.fast(ci@data$tmax, ci@date.factors[[match.arg(freq)]], min, na.rm=TRUE)) * ci@namasks[[match.arg(freq)]]$tmax) } #' Monthly Minimum of Daily Minimum Temperature #' @@ -901,7 +1083,9 @@ climdex.txn <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null(ci #' @template get_generic_example #' #' @export -climdex.tnn <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null(ci@data$tmin)); return(suppressWarnings(tapply.fast(ci@data$tmin, ci@date.factors[[match.arg(freq)]], min, na.rm=TRUE)) * ci@namasks[[match.arg(freq)]]$tmin) } +climdex.tnn <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$tmin)); + return(suppressWarnings(tapply.fast(ci@data$tmin, ci@date.factors[[match.arg(freq)]], min, na.rm=TRUE)) * ci@namasks[[match.arg(freq)]]$tmin) } ## Our implementation currently follows the example set by fclimdex for dealing with missing values, which is wrong; it biases results upwards when missing values are present. @@ -923,7 +1107,9 @@ climdex.tnn <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null(ci #' @template get_generic_example #' #' @export -climdex.tn10p <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null(ci@data$tmin) && !is.null(ci@quantiles$tmin)); return(percent.days.op.threshold(ci@data$tmin, ci@dates, ci@jdays, ci@date.factors[[match.arg(freq)]], ci@quantiles$tmin$outbase$q10, ci@quantiles$tmin$inbase$q10, ci@base.range, "<", ci@max.missing.days[match.arg(freq)]) * ci@namasks[[match.arg(freq)]]$tmin) } +climdex.tn10p <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$tmin) && !is.null(ci@quantiles$tmin)); + return(percent.days.op.threshold(ci@data$tmin, ci@dates, ci@jdays, ci@date.factors[[match.arg(freq)]], ci@quantiles$tmin$outbase$q10, ci@quantiles$tmin$inbase$q10, ci@base.range, "<", ci@max.missing.days[match.arg(freq)]) * ci@namasks[[match.arg(freq)]]$tmin) } #' Percent of Values Below 10th Percentile Daily Maximum Temperature #' @@ -943,7 +1129,9 @@ climdex.tn10p <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null( #' @template get_generic_example #' #' @export -climdex.tx10p <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null(ci@data$tmax) && !is.null(ci@quantiles$tmax)); return(percent.days.op.threshold(ci@data$tmax, ci@dates, ci@jdays, ci@date.factors[[match.arg(freq)]], ci@quantiles$tmax$outbase$q10, ci@quantiles$tmax$inbase$q10, ci@base.range, "<", ci@max.missing.days[match.arg(freq)]) * ci@namasks[[match.arg(freq)]]$tmax) } +climdex.tx10p <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$tmax) && !is.null(ci@quantiles$tmax)); + return(percent.days.op.threshold(ci@data$tmax, ci@dates, ci@jdays, ci@date.factors[[match.arg(freq)]], ci@quantiles$tmax$outbase$q10, ci@quantiles$tmax$inbase$q10, ci@base.range, "<", ci@max.missing.days[match.arg(freq)]) * ci@namasks[[match.arg(freq)]]$tmax) } #' Percent of Values Above 90th Percentile Daily Minimum Temperature #' @@ -963,7 +1151,9 @@ climdex.tx10p <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null( #' @template get_generic_example #' #' @export -climdex.tn90p <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null(ci@data$tmin) && !is.null(ci@quantiles$tmin)); return(percent.days.op.threshold(ci@data$tmin, ci@dates, ci@jdays, ci@date.factors[[match.arg(freq)]], ci@quantiles$tmin$outbase$q90, ci@quantiles$tmin$inbase$q90, ci@base.range, ">", ci@max.missing.days[match.arg(freq)]) * ci@namasks[[match.arg(freq)]]$tmin) } +climdex.tn90p <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$tmin) && !is.null(ci@quantiles$tmin)); + return(percent.days.op.threshold(ci@data$tmin, ci@dates, ci@jdays, ci@date.factors[[match.arg(freq)]], ci@quantiles$tmin$outbase$q90, ci@quantiles$tmin$inbase$q90, ci@base.range, ">", ci@max.missing.days[match.arg(freq)]) * ci@namasks[[match.arg(freq)]]$tmin) } #' Percent of Values Above 90th Percentile Daily Maximum Temperature #' @@ -983,12 +1173,18 @@ climdex.tn90p <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null( #' @template get_generic_example #' #' @export -climdex.tx90p <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null(ci@data$tmax) && !is.null(ci@quantiles$tmax)); return(percent.days.op.threshold(ci@data$tmax, ci@dates, ci@jdays, ci@date.factors[[match.arg(freq)]], ci@quantiles$tmax$outbase$q90, ci@quantiles$tmax$inbase$q90, ci@base.range, ">", ci@max.missing.days[match.arg(freq)]) * ci@namasks[[match.arg(freq)]]$tmax) } +climdex.tx90p <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$tmax) && !is.null(ci@quantiles$tmax)); + return(percent.days.op.threshold(ci@data$tmax, ci@dates, ci@jdays, ci@date.factors[[match.arg(freq)]], + ci@quantiles$tmax$outbase$q90, ci@quantiles$tmax$inbase$q90, ci@base.range, ">", + ci@max.missing.days[match.arg(freq)]) * ci@namasks[[match.arg(freq)]]$tmax) } -#' Warm Spell Duration Index +#' @title Warm Spell Duration Index #' +#' @description #' This function computes the climdex index WSDI. #' +#' @details #' This function takes a climdexInput object as input and computes the climdex #' index WSDI (Warm Spell Duration Index). #' @@ -1008,12 +1204,19 @@ climdex.tx90p <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null( #' @template get_generic_example #' #' @export -climdex.wsdi <- function(ci, spells.can.span.years=FALSE) { stopifnot(!is.null(ci@data$tmax) && !is.null(ci@quantiles$tmax)); return(threshold.exceedance.duration.index(ci@data$tmax, ci@date.factors$annual, ci@jdays, ci@quantiles$tmax$outbase$q90, ">", spells.can.span.years=spells.can.span.years, max.missing.days=ci@max.missing.days['annual']) * ci@namasks$annual$tmax) } +climdex.wsdi <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal"), spells.can.span.years=FALSE) { + stopifnot(!is.null(ci@data$tmax) && !is.null(ci@quantiles$tmax)) + max_missing_days = ci@max.missing.days + stopifnot(sum(match.arg(freq) == names(max_missing_days)) == 1) # Check if we only get one TRUE, behavior is undefined if not + current_max_missing_days = max_missing_days[match.arg(freq) == names(max_missing_days)] + return(threshold.exceedance.duration.index(ci@data$tmax, ci@date.factors[[match.arg(freq)]], ci@jdays, ci@quantiles$tmax$outbase$q90, ">", spells.can.span.years=spells.can.span.years, max.missing.days=current_max_missing_days) * ci@namasks[[match.arg(freq)]]$tmax) } -#' Cold Spell Duration Index +#' @title Cold Spell Duration Index #' +#' @description #' This function computes the climdex index CSDI. #' +#' @details #' This function takes a climdexInput object as input and computes the climdex #' index CSDI (Cold Spell Duration Index). #' @@ -1033,7 +1236,13 @@ climdex.wsdi <- function(ci, spells.can.span.years=FALSE) { stopifnot(!is.null(c #' @template get_generic_example #' #' @export -climdex.csdi <- function(ci, spells.can.span.years=FALSE) { stopifnot(!is.null(ci@data$tmin) && !is.null(ci@quantiles$tmin)); return(threshold.exceedance.duration.index(ci@data$tmin, ci@date.factors$annual, ci@jdays, ci@quantiles$tmin$outbase$q10, "<", spells.can.span.years=spells.can.span.years, max.missing.days=ci@max.missing.days['annual']) * ci@namasks$annual$tmin) } +climdex.csdi <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal"), spells.can.span.years=FALSE) { + stopifnot(!is.null(ci@data$tmin) && !is.null(ci@quantiles$tmin)); + max_missing_days = ci@max.missing.days + stopifnot(sum(match.arg(freq) == names(max_missing_days)) == 1) # Check if we only get one TRUE, behavior is undefined if not + current_max_missing_days = max_missing_days[match.arg(freq) == names(max_missing_days)] + return(threshold.exceedance.duration.index(ci@data$tmin, ci@date.factors[[match.arg(freq)]], ci@jdays, ci@quantiles$tmin$outbase$q10, "<", + spells.can.span.years=spells.can.span.years, max.missing.days=current_max_missing_days) * ci@namasks[[match.arg(freq)]]$tmin) } #' Mean Diurnal Temperature Range #' @@ -1054,7 +1263,9 @@ climdex.csdi <- function(ci, spells.can.span.years=FALSE) { stopifnot(!is.null(c #' @template get_generic_example #' #' @export -climdex.dtr <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null(ci@data$tmin) && !is.null(ci@data$tmax) && !is.null(ci@data$tavg)); return(mean.daily.temp.range(ci@data$tmax, ci@data$tmin, ci@date.factors[[match.arg(freq)]]) * ci@namasks[[match.arg(freq)]]$tavg) } +climdex.dtr <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$tmin) && !is.null(ci@data$tmax) && !is.null(ci@data$tavg)); + return(mean.daily.temp.range(ci@data$tmax, ci@data$tmin, ci@date.factors[[match.arg(freq)]]) * ci@namasks[[match.arg(freq)]]$tavg) } #' Monthly Maximum 1-day Precipitation #' @@ -1072,7 +1283,9 @@ climdex.dtr <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null(ci #' @template get_generic_example #' #' @export -climdex.rx1day <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null(ci@data$prec)); return(nday.consec.prec.max(ci@data$prec, ci@date.factors[[match.arg(freq)]], 1) * ci@namasks[[match.arg(freq)]]$prec) } +climdex.rx1day <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$prec)); + return(nday.consec.prec.max(ci@data$prec, ci@date.factors[[match.arg(freq)]], 1) * ci@namasks[[match.arg(freq)]]$prec) } #' Monthly Maximum 5-day Consecutive Precipitation #' @@ -1092,7 +1305,9 @@ climdex.rx1day <- function(ci, freq=c("monthly", "annual")) { stopifnot(!is.null #' @template get_generic_example #' #' @export -climdex.rx5day <- function(ci, freq=c("monthly", "annual"), center.mean.on.last.day=FALSE) { stopifnot(!is.null(ci@data$prec)); return(nday.consec.prec.max(ci@data$prec, ci@date.factors[[match.arg(freq)]], 5, center.mean.on.last.day) * ci@namasks[[match.arg(freq)]]$prec) } +climdex.rx5day <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal"), center.mean.on.last.day=FALSE) { + stopifnot(!is.null(ci@data$prec)); + return(nday.consec.prec.max(ci@data$prec, ci@date.factors[[match.arg(freq)]], 5, center.mean.on.last.day) * ci@namasks[[match.arg(freq)]]$prec) } #' Simple Precpitation Intensity Index #' @@ -1113,7 +1328,9 @@ climdex.rx5day <- function(ci, freq=c("monthly", "annual"), center.mean.on.last. #' @template get_generic_example #' #' @export -climdex.sdii <- function(ci) { stopifnot(!is.null(ci@data$prec)); return(simple.precipitation.intensity.index(ci@data$prec, ci@date.factors$annual) * ci@namasks$annual$prec) } +climdex.sdii <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$prec)); + return(simple.precipitation.intensity.index(ci@data$prec, ci@date.factors[[match.arg(freq)]]) * ci@namasks[[match.arg(freq)]]$prec) } #' Precipitation Exceeding 10mm Per Day #' @@ -1130,7 +1347,9 @@ climdex.sdii <- function(ci) { stopifnot(!is.null(ci@data$prec)); return(simple. #' @template get_generic_example #' #' @export -climdex.r10mm <- function(ci) { stopifnot(!is.null(ci@data$prec)); return(number.days.op.threshold(ci@data$prec, ci@date.factors$annual, 10, ">=") * ci@namasks$annual$prec) } +climdex.r10mm <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$prec)); + return(number.days.op.threshold(ci@data$prec, ci@date.factors[[match.arg(freq)]], 10, ">=") * ci@namasks[[match.arg(freq)]]$prec) } #' Precipitation Exceeding 20mm Per Day #' @@ -1147,7 +1366,9 @@ climdex.r10mm <- function(ci) { stopifnot(!is.null(ci@data$prec)); return(number #' @template get_generic_example #' #' @export -climdex.r20mm <- function(ci) { stopifnot(!is.null(ci@data$prec)); return(number.days.op.threshold(ci@data$prec, ci@date.factors$annual, 20, ">=") * ci@namasks$annual$prec) } +climdex.r20mm <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$prec)); + return(number.days.op.threshold(ci@data$prec, ci@date.factors[[match.arg(freq)]], 20, ">=") * ci@namasks[[match.arg(freq)]]$prec) } #' Precipitation Exceeding A Specified Amount Per Day #' @@ -1165,11 +1386,11 @@ climdex.r20mm <- function(ci) { stopifnot(!is.null(ci@data$prec)); return(number #' @template get_generic_example #' #' @export -climdex.rnnmm <- function(ci, threshold=1) { +climdex.rnnmm <- function(ci, threshold=1, freq=c("monthly", "annual", "halfyear", "seasonal")) { stopifnot(!is.null(ci@data$prec)); if(!is.numeric(threshold) || length(threshold) != 1) stop("Please specify a single numeric threshold value."); - return(number.days.op.threshold(ci@data$prec, ci@date.factors$annual, threshold, ">=") * ci@namasks$annual$prec) + return(number.days.op.threshold(ci@data$prec, ci@date.factors[[match.arg(freq)]], threshold, ">=") * ci@namasks[[match.arg(freq)]]$prec) } #' Maximum Consecutive Dry Days @@ -1186,7 +1407,9 @@ climdex.rnnmm <- function(ci, threshold=1) { #' @template get_generic_example #' #' @export -climdex.cdd <- function(ci, spells.can.span.years=TRUE) { stopifnot(!is.null(ci@data$prec)); return(spell.length.max(ci@data$prec, ci@date.factors$annual, 1, "<", spells.can.span.years) * ci@namasks$annual$prec) } +climdex.cdd <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal"), spells.can.span.years=TRUE) { + stopifnot(!is.null(ci@data$prec)); + return(spell.length.max(ci@data$prec, ci@date.factors[[match.arg(freq)]], 1, "<", spells.can.span.years) * ci@namasks[[match.arg(freq)]]$prec) } #' Maximum Consecutive Wet Days #' @@ -1203,14 +1426,105 @@ climdex.cdd <- function(ci, spells.can.span.years=TRUE) { stopifnot(!is.null(ci@ #' @template get_generic_example #' #' @export -climdex.cwd <- function(ci, spells.can.span.years=TRUE) { stopifnot(!is.null(ci@data$prec)); return(spell.length.max(ci@data$prec, ci@date.factors$annual, 1, ">=", spells.can.span.years) * ci@namasks$annual$prec) } +climdex.cwd <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal"), spells.can.span.years=TRUE) { + stopifnot(!is.null(ci@data$prec)); + return(spell.length.max(ci@data$prec, ci@date.factors[[match.arg(freq)]], 1, ">=", spells.can.span.years) * ci@namasks[[match.arg(freq)]]$prec) } + + +#' Total Daily Precipitation Exceeding 75\%ile Threshold +#' +#' This function computes the climdex index R75p. +#' +#' This function takes a climdexInput object as input and computes the climdex +#' index R75p: the annual sum of precipitation in days where daily precipitation exceeds the +#' 75th percentile of daily precipitation in the base period. +#' +#' @param ci Object of type climdexInput. +#' @return A vector containing an annual timeseries of precipitation exceeding +#' the threshold. +#' @template generic_seealso_references +#' @templateVar cdxvar r75p +#' @templateVar cdxdescription an annual timeseries of the R75p index. +#' @template get_generic_example +#' +#' @export +climdex.r75p <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$prec) && !is.null(ci@quantiles$prec)); + return(total.precip.op.threshold(ci@data$prec, ci@date.factors[[match.arg(freq)]], ci@quantiles$prec['q75'], ">") * ci@namasks[[match.arg(freq)]]$prec) } + #' Total Daily Precipitation Exceeding 95\%ile Threshold #' +#' This function computes the climdex index R95p. +#' +#' This function takes a climdexInput object as input and computes the climdex +#' index R95p: the annual sum of precipitation in days where daily precipitation exceeds the +#' 95th percentile of daily precipitation in the base period. +#' +#' @param ci Object of type climdexInput. +#' @return A vector containing an annual timeseries of precipitation exceeding +#' the threshold. +#' @template generic_seealso_references +#' @templateVar cdxvar r95p +#' @templateVar cdxdescription an annual timeseries of the R95p index. +#' @template get_generic_example +#' +#' @export +climdex.r95p <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$prec) && !is.null(ci@quantiles$prec)); + return(total.precip.op.threshold(ci@data$prec, ci@date.factors[[match.arg(freq)]], ci@quantiles$prec['q95'], ">") * ci@namasks[[match.arg(freq)]]$prec) } + +#' Total Daily Precipitation Exceeding 99\%ile Threshold +#' +#' This function computes the climdex index R99p. +#' +#' This function takes a climdexInput object as input and computes the climdex +#' index R99p: the annual sum of precipitation in days where daily precipitation exceeds the +#' 99th percentile of daily precipitation in the base period. +#' +#' @param ci Object of type climdexInput. +#' @return A vector containing an annual timeseries of precipitation exceeding +#' the threshold. +#' @template generic_seealso_references +#' @templateVar cdxvar r99p +#' @templateVar cdxdescription an annual timeseries of the R99p index. +#' @template get_generic_example +#' +#' @export +climdex.r99p <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$prec) && !is.null(ci@quantiles$prec)); + return(total.precip.op.threshold(ci@data$prec, ci@date.factors[[match.arg(freq)]], ci@quantiles$prec['q99'], ">") * ci@namasks[[match.arg(freq)]]$prec) } + +#' Precipitation fraction due to moderate wet days (exceeding 75\%ile Threshold) +#' +#' This function computes the climdex index R75pTOT. +#' +#' This function takes a climdexInput object as input and computes the climdex +#' index R75pTOT: the fraction of the sum of precipitation in days where daily precipitation exceeds the +#' 75th percentile of daily precipitation in the base period. +#' +#' @param ci Object of type climdexInput. +#' @return A vector containing an annual timeseries of precipitation exceeding +#' the threshold. +#' @template generic_seealso_references +#' @templateVar cdxvar r75ptot +#' @templateVar cdxdescription an annual timeseries of the R75pTOT index. +#' @template get_generic_example +#' +#' @export +climdex.r75ptot <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + + stopifnot(!is.null(ci@data$prec),!is.null(ci@quantiles$prec)); + prcptot <- total.precip.op.threshold(ci@data$prec, ci@date.factors[[match.arg(freq)]], 1, ">=") * ci@namasks[[match.arg(freq)]]$prec + r75p <- total.precip.op.threshold(ci@data$prec, ci@date.factors[[match.arg(freq)]], ci@quantiles$prec['q75'], ">") * ci@namasks[[match.arg(freq)]]$prec + return(100*r75p/prcptot) } + +#' Precipitation fraction due to very wet days (exceeding 95\%ile Threshold) +#' #' This function computes the climdex index R95pTOT. #' #' This function takes a climdexInput object as input and computes the climdex -#' index R95pTOT: the annual sum of precipitation in days where daily precipitation exceeds the +#' index R95pTOT: the fraction of the sum of precipitation in days where daily precipitation exceeds the #' 95th percentile of daily precipitation in the base period. #' #' @param ci Object of type climdexInput. @@ -1222,14 +1536,19 @@ climdex.cwd <- function(ci, spells.can.span.years=TRUE) { stopifnot(!is.null(ci@ #' @template get_generic_example #' #' @export -climdex.r95ptot <- function(ci) { stopifnot(!is.null(ci@data$prec) && !is.null(ci@quantiles$prec)); return(total.precip.op.threshold(ci@data$prec, ci@date.factors$annual, ci@quantiles$prec['q95'], ">") * ci@namasks$annual$prec) } +climdex.r95ptot <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + + stopifnot(!is.null(ci@data$prec),!is.null(ci@quantiles$prec)); + prcptot <- total.precip.op.threshold(ci@data$prec, ci@date.factors[[match.arg(freq)]], 1, ">=") * ci@namasks[[match.arg(freq)]]$prec + r95p <- total.precip.op.threshold(ci@data$prec, ci@date.factors[[match.arg(freq)]], ci@quantiles$prec['q95'], ">") * ci@namasks[[match.arg(freq)]]$prec + return(100*r95p/prcptot) } -#' Total Daily Precipitation Exceeding 99\%ile Threshold +#' Precipitation fraction due to extremly wet days (exceeding 99\%ile Threshold) #' #' This function computes the climdex index R99pTOT. #' #' This function takes a climdexInput object as input and computes the climdex -#' index R99pTOT: the annual sum of precipitation in days where daily precipitation exceeds the +#' index R99pTOT: the fraction of the sum of precipitation in days where daily precipitation exceeds the #' 99th percentile of daily precipitation in the base period. #' #' @param ci Object of type climdexInput. @@ -1241,7 +1560,12 @@ climdex.r95ptot <- function(ci) { stopifnot(!is.null(ci@data$prec) && !is.null(c #' @template get_generic_example #' #' @export -climdex.r99ptot <- function(ci) { stopifnot(!is.null(ci@data$prec) && !is.null(ci@quantiles$prec)); return(total.precip.op.threshold(ci@data$prec, ci@date.factors$annual, ci@quantiles$prec['q99'], ">") * ci@namasks$annual$prec) } +climdex.r99ptot <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + + stopifnot(!is.null(ci@data$prec),!is.null(ci@quantiles$prec)); + prcptot <- total.precip.op.threshold(ci@data$prec, ci@date.factors[[match.arg(freq)]], 1, ">=") * ci@namasks[[match.arg(freq)]]$prec + r99p <- total.precip.op.threshold(ci@data$prec, ci@date.factors[[match.arg(freq)]], ci@quantiles$prec['q99'], ">") * ci@namasks[[match.arg(freq)]]$prec + return(100*r99p/prcptot) } #' Total Daily Precipitation #' @@ -1259,7 +1583,1020 @@ climdex.r99ptot <- function(ci) { stopifnot(!is.null(ci@data$prec) && !is.null(c #' @template get_generic_example #' #' @export -climdex.prcptot <- function(ci) { stopifnot(!is.null(ci@data$prec)); return(total.precip.op.threshold(ci@data$prec, ci@date.factors$annual, 1, ">=") * ci@namasks$annual$prec) } +climdex.prcptot <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$prec)); + return(total.precip.op.threshold(ci@data$prec, ci@date.factors[[match.arg(freq)]], 1, ">=") * ci@namasks[[match.arg(freq)]]$prec) } + +#' Standardized Precipitation Index 3-mon +#' @description Adapted by ClimPACT2 and SPEI (arguments taken from here) libraries. +#' Given a timeseries of daily precipitation amounts [mm], the function calculates the Standardized Precipitation Index (SPI) as it is created and defined by the SPEI library. +#' This function computes the climdex index spi at multiple scales. +#' +#' @param ci Object of type climdexInput (representing the daily precipitation [mm]) +#' @param freq Time frequency to aggregate to. Allowed only monthly +#' @param scale an integer, representing the time scale at which the SPI will be computed. Default is 3 +#' @param distribution name of the distribution function to be used for computing the SPI (default ’Gamma’) +#' @param fit name of the method used for computing the distribution function parameters (default 'ub-pwm') +#' @param kernel optional, a list defining the type of kernel used for computing the SPI atscales higher than one. Defaults to unshifted rectangular kernel. +#' @param ref.start optional, starting point of the reference period used for computing the index. Defaults to NULL, indicating that the first value in data will be used as starting point. +#' @param ref.end optional, ending point of the reference period used for computing the index. Defaults to NULL, indicating that the last value in data will be used as ending point. +#' @param ... For more details please refer to the SPEI +#' @return A vector containing a monthly SPI at the selected scale +#' @author Christiana Photiadou (KNMI) +#' @references \url{http://www.ecad.eu/documents/atbd.pdf} +#' @references \url{https://cran.r-project.org/web/packages/SPEI/SPEI.pdf} +#' @importFrom SPEI spi +#' +#' @export +climdex.spi3 <- function(ci, freq=c("monthly"), scale=3, distribution="Gamma", fit="ub-pwm", kernel=list(type="rectangular",shift=1), ref.start=NULL, ref.end=NULL){ + + spiprec <- ci@data$prec + spifactor <- ci@date.factors$monthly + prec_sum <- as.numeric(tapply.fast(spiprec, spifactor, sum, na.rm=TRUE)) + + ts.start <- c(as.numeric(format(ci@dates[1],format="%Y")),1) + ts.end <- c(as.numeric(format(ci@dates[length(ci@dates)],format="%Y")),12) + + data.spi <- ts(prec_sum,freq=12,start=ts.start,end=ts.end) + + spi_col <- spi(data.spi, scale=scale,ref.start=ref.start,ref.end=ref.end, + distribution=distribution,fit=fit,kernal=kernal,na.rm=TRUE) + + tmpvar <- (spi_col$fitted) + tmpvar <- ifelse(tmpvar=="-Inf",NA,tmpvar) + tmpvar <- ifelse(tmpvar=="Inf",NA,tmpvar) + + tmpvar <- ifelse(tmpvar=="NaNf",NA,tmpvar) + tmpvar <- ifelse(tmpvar=="NaN",NA,tmpvar) + + x <- as.numeric(tmpvar) + names(x) <- unique(spifactor) + return(x) +} + +#' Standardized Precipitation Index 6-mon +#' @description Adapted by ClimPACT2 and SPEI (arguments taken from here) libraries. +#' Given a timeseries of daily precipitation amounts [mm], the function calculates the Standardized Precipitation Index (SPI) as it is created and defined by the SPEI library. +#' This function computes the climdex index spi at multiple scales. +#' +#' @param ci Object of type climdexInput (representing the daily precipitation [mm]) +#' @param freq Time frequency to aggregate to. Allowed only monthly +#' @param scale an integer, representing the time scale at which the SPI will be computed. Default is 3 +#' @param distribution name of the distribution function to be used for computing the SPI (default ’Gamma’) +#' @param fit name of the method used for computing the distribution function parameters (default 'ub-pwm') +#' @param kernel optional, a list defining the type of kernel used for computing the SPI atscales higher than one. Defaults to unshifted rectangular kernel. +#' @param ref.start optional, starting point of the reference period used for computing the index. Defaults to NULL, indicating that the first value in data will be used as starting point. +#' @param ref.end optional, ending point of the reference period used for computing the index. Defaults to NULL, indicating that the last value in data will be used as ending point. +#' @param ... For more details please refer to the SPEI +#' @return A vector containing a monthly SPI at the selected scale +#' @author Christiana Photiadou (KNMI) +#' @references \url{http://www.ecad.eu/documents/atbd.pdf} +#' @references \url{https://cran.r-project.org/web/packages/SPEI/SPEI.pdf} +#' @importFrom SPEI spi +#' +#' @export +climdex.spi6 <- function(ci, freq=c("monthly"), scale=6, distribution="Gamma", fit="ub-pwm", kernel=list(type="rectangular",shift=1), ref.start=NULL, ref.end=NULL){ + + spiprec <- ci@data$prec + spifactor <- ci@date.factors$monthly + prec_sum <- as.numeric(tapply.fast(spiprec, spifactor, sum, na.rm=TRUE)) + + ts.start <- c(as.numeric(format(ci@dates[1],format="%Y")),1) + ts.end <- c(as.numeric(format(ci@dates[length(ci@dates)],format="%Y")),12) + + data.spi <- ts(prec_sum,freq=12,start=ts.start,end=ts.end) + + spi_col <- spi(data.spi, scale=scale,ref.start=ref.start,ref.end=ref.end, + distribution=distribution,fit=fit,kernal=kernal,na.rm=TRUE) + + tmpvar <- (spi_col$fitted) + tmpvar <- ifelse(tmpvar=="-Inf",NA,tmpvar) + tmpvar <- ifelse(tmpvar=="Inf",NA,tmpvar) + + tmpvar <- ifelse(tmpvar=="NaNf",NA,tmpvar) + tmpvar <- ifelse(tmpvar=="NaN",NA,tmpvar) + + x <- as.numeric(tmpvar) + names(x) <- unique(spifactor) + return(x) +} + + +#' Consecutive Summer Days +#' @description +#' This function takes a climdexInput object as input and computes the climdex +#' index csu: the annual (or at different periods) count of consecutive summer days (TX >25C) +#' +#' @param ci Object of type climdexInput. Here the daily maximum temperature. +#' @param freq Time frequency to aggregate to. Allowed only monthly, annual, halfyear, seasonal. +#' @param spells.can.span.years Default FALSE +#' @return A vector containing a timeseries of the number of consecutive summer days in a given period (freq). +#' @author Christiana Photiadou (KNMI) +#' @references \url{http://www.ecad.eu/documents/atbd.pdf} +#' +#' @export +climdex.csu <- function(ci,freq=c("monthly","annual", "halfyear", "seasonal"),spells.can.span.years=FALSE) { + stopifnot(!is.null(ci@data$tmax)); + return(spell.length.max(ci@data$tmax, ci@date.factors[[match.arg(freq)]], 25, ">", spells.can.span.years) * ci@namasks[[match.arg(freq)]]$tmax) +} + + +#' Consecutive Frost Days +#' @description +#' This function takes a climdexInput object as input and computes the climdex +#' index cfd: the annual (or at different periods) count of consecutive frost days (TN < 0C) +#' +#' @param ci Object of type climdexInput. Here the daily maximum temperature. +#' @param freq Time frequency to aggregate to. Allowed only monthly, annual, halfyear, seasonal. +#' @param spells.can.span.years Default FALSE +#' @return A vector containing a timeseries of the number of consecutive summer days in a given period (freq). +#' @author Christiana Photiadou (KNMI) +#' @references \url{http://www.ecad.eu/indicesextremes/indicesdictionary.php} +#' +#' @export +climdex.cfd <- function(ci,freq=c("monthly","annual", "halfyear", "seasonal"),spells.can.span.years=FALSE) { + stopifnot(!is.null(ci@data$tmin)); + return(spell.length.max(ci@data$tmin, ci@date.factors[[match.arg(freq)]], 0, "<", spells.can.span.years) * ci@namasks[[match.arg(freq)]]$tmin) } + + + +#' Heating Degree Days +#' @description +#' This function takes a climdexInput object as input and computes the climdex +#' index hd17: the annual (or at different periods) sum of heating degree days (17-tavg) +#' @param ci Object of type climdexInput. Here the daily maximum temperature. +#' @param freq Time frequency to aggregate to. Allowed only monthly, annual, halfyear, seasonal. +#' @param spells.can.span.years Default FALSE +#' @return A vector containing a timeseries of the number of consecutive summer days in a given period (freq). +#' @author Christiana Photiadou (KNMI) +#' @references \url{http://www.ecad.eu/indicesextremes/indicesdictionary.php} +#' +#' @export +climdex.hd17 <- function(ci, freq=c("monthly","annual", "halfyear", "seasonal")) { + stopifnot(!is.null(ci@data$tavg)); + return(tapply((17 - ci@data$tavg), ci@date.factors[[match.arg(freq)]], sum)* ci@namasks[[match.arg(freq)]]$tavg) +} + +#' Lengths of strings of TRUE (1 & 0) values +## -- introduced by C. Photiadou (KNMI), September 2015 +#' Computes which days are above or below the baseline threshold. +#' +#' This function computes which days are above or below baseline thresholds. +#' It is used to implement the compound indices. +#' It is based on the "percent.days.op.threshold" +#' +#' @param temp Sequence of temperature values. +#' @param dates Sequence of associated dates. +#' @param jdays Sequence of associated days of year. +#' @param date.factor Factor to aggregate data using. +#' @param threshold.outside.base Sequence of thresholds to be used for data outside the base period. +#' @param base.thresholds Data structure containing sets of thresholds to be used inside the base period; see \link{climdexInput-class}. +#' @param base.range Date range (type PCICt) of the baseline period. +#' @param op Comparison operator to use. +#' @param max.missing.days Maximum number of NA values per time period. +#' @return A vector consisting of the mean fraction of days above or below the supplied set of thresholds. +#' @note If date.factor is omitted, daily series will be returned. +#' @seealso \link{climdexInput-class}. +#' @keywords ts climate +days.op.threshold <- function(temp, dates, jdays, date.factor, threshold.outside.base, base.thresholds, base.range, op='<') { + f <- match.fun(op) + dat <- f(temp, threshold.outside.base[jdays]) + + inset <- dates >= base.range[1] & dates <= base.range[2] + ## Don't use in-base thresholds with data shorter than two years; no years to replace with. + if(sum(inset) > 0 && length(dates) >= 360 * 2) { + jdays.base <- jdays[inset] + years.base <- climdex.pcic:::get.years(dates[inset]) + + ## Get number of base years, subset temp data to base period only. + temp.base <- temp[inset] + years.base.range <- range(years.base) + byrs <- (years.base.range[2] - years.base.range[1] + 1) + + ## Linearize thresholds, then compare them to the temperatures + bdim <- dim(base.thresholds) + dim(base.thresholds) <- c(bdim[1] * bdim[2], bdim[3]) + yday.byr.indices <- jdays.base + (years.base - climdex.pcic:::get.years(base.range)[1]) * bdim[1] + f.result <- f(rep(temp.base, byrs - 1), base.thresholds[yday.byr.indices,]) + dim(f.result) <- c(length(yday.byr.indices), bdim[3]) + + ## Chop up data along the 2nd dim into a list; sum elements of the list + dat[inset] <- rowSums(f.result, na.rm=TRUE) / (byrs - 1) + } + return(dat) +} + +## Climate Compound Indices +## -- introduced by C. Photiadou (KNMI), September 2015 +#' Compound Indices +#' +#' Cold-dry days +#' +#' This function computes the climdex index CD. +#' +#' This function takes a climdexInput object as input and computes the climdex +#' index CD: the number of days where TG<25 & RR<25 (for wet days: days where precipitation is at least 1mm). +#' Note: this function doesnt use directly tavg but derives it from tmax & tmin (concistent with .ncdf) +#' @param ci Object of type climdexInput (representing the daily precipitation [mm] and the averaged daily temperature [C]) +#' @return A vector containing an annual timeseries of precipitation in wet days. +#' @template generic_seealso_references +#' +#' @export +climdex.cd <- function(ci, freq=c("annual", "monthly", "halfyear", "seasonal"), precip.thresh="q25", precip.op="<", temp.thresh="q25", temp.op="<") { + + stopifnot(!is.null(ci@data$prec) && !is.null(ci@quantiles$prec) && !is.null(ci@data$tavg) && !is.null(ci@quantiles$tavg)) + daily.prec <- ci@data$prec + daily.temp <- ci@data$tavg + q.precip <- ci@quantiles$prec[[precip.thresh]] + f.prec <- match.fun(precip.op) + f.temp<- days.op.threshold(daily.temp, ci@dates, ci@jdays, ci@date.factors[[match.arg(freq)]], + ci@quantiles$tavg$outbase[[temp.thresh]], + ci@quantiles$tavg$inbase[[temp.thresh]], ci@base.range, temp.op) + #Convert from 1/0 to TRUE/FALSE + logic.f.temp <- f.temp==1 + df.for.data <- data.frame(daily.prec, logic.f.temp, date.factor = ci@date.factors[[match.arg(freq)]]) + result <- lapply(split(df.for.data, df.for.data$date.factor), function(chunk) { + sum(f.prec(chunk$daily.prec, q.precip) & chunk$logic.f.temp, na.rm=FALSE) + }) + return(unlist(result)) +} + +#' Cold-wet days +#' +#' This function computes the climdex index CW. +#' +#' This function takes a climdexInput object as input and computes the climdex +#' index CW: the number of days where TG<25 & RR>75 (for wet days: days where precipitation is at least 1mm). +#' +#' @param ci Object of type climdexInput (representing the daily precipitation [mm] and the averaged daily temperature [C]) +#' @return A vector containing an annual timeseries of precipitation in wet days. +#' @template generic_seealso_references +#' +#' @export +climdex.cw <- function(ci, freq=c("annual", "monthly", "halfyear", "seasonal"), precip.thresh="q75", precip.op=">", temp.thresh="q25", temp.op="<") { + + stopifnot(!is.null(ci@data$prec) && !is.null(ci@quantiles$prec) && !is.null(ci@data$tavg) && !is.null(ci@quantiles$tavg)) + daily.prec <- ci@data$prec + daily.temp <- ci@data$tavg + q.precip <- ci@quantiles$prec[[precip.thresh]] + f.prec <- match.fun(precip.op) + f.temp<- days.op.threshold(daily.temp, ci@dates, ci@jdays, ci@date.factors[[match.arg(freq)]], + ci@quantiles$tavg$outbase[[temp.thresh]], + ci@quantiles$tavg$inbase[[temp.thresh]], ci@base.range, temp.op) + #Convert from 1/0 to TRUE/FALSE + logic.f.temp <- f.temp==1 + df.for.data <- data.frame(daily.prec, logic.f.temp, date.factor = ci@date.factors[[match.arg(freq)]]) + result <- lapply(split(df.for.data, df.for.data$date.factor), function(chunk) { + sum(f.prec(chunk$daily.prec, q.precip) & chunk$logic.f.temp, na.rm=FALSE) + }) + return(unlist(result)) +} + +#' Warm-dry days +#' +#' This function computes the climdex index WD. +#' +#' This function takes a climdexInput object as input and computes the climdex +#' index WD: the number of days where TG>75 & RR<25 (for wet days: days where precipitation is at least 1mm). +#' +#' @param ci Object of type climdexInput (representing the daily precipitation [mm] and the averaged daily temperature [C]) +#' @return A vector containing an annual timeseries of precipitation in wet days. +#' @template generic_seealso_references +#' +#' @export +climdex.wd <- function(ci, freq=c("annual", "monthly", "halfyear", "seasonal"), precip.thresh="q25", precip.op="<", temp.thresh="q75", temp.op=">") { + + stopifnot(!is.null(ci@data$prec) && !is.null(ci@quantiles$prec) && !is.null(ci@data$tavg) && !is.null(ci@quantiles$tavg)) + daily.prec <- ci@data$prec + daily.temp <- ci@data$tavg + q.precip <- ci@quantiles$prec[[precip.thresh]] + f.prec <- match.fun(precip.op) + f.temp<- days.op.threshold(daily.temp, ci@dates, ci@jdays, ci@date.factors[[match.arg(freq)]], + ci@quantiles$tavg$outbase[[temp.thresh]], + ci@quantiles$tavg$inbase[[temp.thresh]], ci@base.range, temp.op) + #Convert from 1/0 to TRUE/FALSE + logic.f.temp <- f.temp==1 + df.for.data <- data.frame(daily.prec, logic.f.temp, date.factor = ci@date.factors[[match.arg(freq)]]) + result <- lapply(split(df.for.data, df.for.data$date.factor), function(chunk) { + sum(f.prec(chunk$daily.prec, q.precip) & chunk$logic.f.temp, na.rm=FALSE) + }) + return(unlist(result)) +} + +#' Warm-wet days +#' +#' This function computes the climdex index WW +#' +#' This function takes a climdexInput object as input and computes the climdex +#' index WW: the number of days where TG>75 & RR>75 (for wet days: days where precipitation is at least 1mm). +#' +#' @param ci Object of type climdexInput (representing the daily precipitation [mm] and the averaged daily temperature [C]) +#' @return A vector containing an annual timeseries of precipitation in wet days. +#' @template generic_seealso_references +#' +#' @export +climdex.ww <- function(ci, freq=c("annual", "monthly", "halfyear", "seasonal"), precip.thresh="q75", precip.op=">", temp.thresh="q75", temp.op=">") { + + stopifnot(!is.null(ci@data$prec) && !is.null(ci@quantiles$prec) && !is.null(ci@data$tavg) && !is.null(ci@quantiles$tavg)) + daily.prec <- ci@data$prec + daily.temp <- ci@data$tavg + q.precip <- ci@quantiles$prec[[precip.thresh]] + f.prec <- match.fun(precip.op) + f.temp<- days.op.threshold(daily.temp, ci@dates, ci@jdays, ci@date.factors[[match.arg(freq)]], + ci@quantiles$tavg$outbase[[temp.thresh]], + ci@quantiles$tavg$inbase[[temp.thresh]], ci@base.range, temp.op) + #Convert from 1/0 to TRUE/FALSE + logic.f.temp <- f.temp==1 + df.for.data <- data.frame(daily.prec, logic.f.temp, date.factor = ci@date.factors[[match.arg(freq)]]) + result <- lapply(split(df.for.data, df.for.data$date.factor), function(chunk) { + sum(f.prec(chunk$daily.prec, q.precip) & chunk$logic.f.temp, na.rm=FALSE) + }) + return(unlist(result)) +} + +## Climate indices involving wind timeseries +## -- introduced by R. Posselt (MeteoSwiss), July 2015 + +#' Mean wind +#' @description +#' This function computes the climdex index FG: the mean wind speed measured within a period. +#' +#' @param ci Object of type climdexInput (representing the daily mean wind speed in [m/s]) +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear", "seasonal" or "monthly". Default: "annual". +#' @return A vector containing the time series of mean wind speed. +#' @author Rebekka Posselt (MeteoSwiss) +#' @references \url{http://www.ecad.eu/documents/atbd.pdf} +#' +#' @export +climdex.fg <- function(ci,freq=c("annual","halfyear","seasonal","monthly")) { + stopifnot(!is.null(ci@data$wind)) + return(tapply.fast(ci@data$wind, ci@date.factors[[match.arg(freq)]], mean, na.rm=TRUE) * ci@namasks[[match.arg(freq)]]$wind) +} + +#' Calm days +#' @description +#' This function computes the climdex index FGcalm days: The number of days with mean wind lower than or equal 2 m/s. +#' +#' @param ci Object of type climdexInput (representing the daily mean wind speed in [m/s]) +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear", "seasonal" or "monthly". Default: "annual". +#' @return A vector containing the number of calm days. +#' @author Rebekka Posselt (MeteoSwiss) +#' @references \url{http://www.ecad.eu/indicesextremes/indicesdictionary.php} +#' +#' @export +#' +climdex.fgcalm <- function(ci,freq=c("annual","halfyear","seasonal","monthly")) { + stopifnot(!is.null(ci@data$wind)) + return(number.days.op.threshold(temp=ci@data$wind, date.factor=ci@date.factors[[match.arg(freq)]], + threshold=2, op="<=") * ci@namasks[[match.arg(freq)]]$wind) +} + +#' Windy days +#' @description +#' This function computes the climdex index FG6bft: The number of days with mean wind greater than or equal 10.8 m/s (~6 Bft). +#' +#' @param ci Object of type climdexInput (representing the daily mean wind speed in [m/s]) +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear", "seasonal" or "monthly". Default: "annual". +#' @return A vector containing the number of calm days. +#' @author Rebekka Posselt (MeteoSwiss) +#' @references \url{http://www.ecad.eu/indicesextremes/indicesdictionary.php} +#' +#' @export +#' +climdex.fg6bft <- function(ci,freq=c("annual","halfyear","seasonal","monthly")) { + stopifnot(!is.null(ci@data$wind)) + return(number.days.op.threshold(temp=ci@data$wind, date.factor=ci@date.factors[[match.arg(freq)]], + threshold=10.8, op=">=") * ci@namasks[[match.arg(freq)]]$wind) +} + +#' Storm days +#' +#' This function computes the climdex index FXstorm: The number of days with wind gusts greater than or equal 20.8 m/s (75 km/h). +#' +#' @param ci Object of type climdexInput (representing the daily maximum wind gust in [m/s]) +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear", "seasonal" or "monthly". Default: "annual". +#' @return A vector containing the number of storm days. +#' @author Rebekka Posselt (MeteoSwiss) +#' @references \url{http://www.ecad.eu/indicesextremes/indicesdictionary.php} +#' +#' @export +#' +climdex.fxstorm <- function(ci,freq=c("annual","halfyear","seasonal","monthly")) { + stopifnot(!is.null(ci@data$wind_gust)) + return(number.days.op.threshold(temp=ci@data$wind_gust, date.factor=ci@date.factors[[match.arg(freq)]], + threshold=20.8, op=">=") * ci@namasks[[match.arg(freq)]]$wind_gust) +} + +#' Maximum wind gust +#' +#' This function computes the climdex index FXx: The maximum of the maximum daily wind gust. +#' +#' @param ci Object of type climdexInput (representing the daily maximum wind gust in [m/s]) +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear", "seasonal" or "monthly". Default: "annual". +#' @return A vector containing the maximum of the maximum daily wind gust. +#' @author Rebekka Posselt (MeteoSwiss) +#' @references \url{http://www.ecad.eu/indicesextremes/indicesdictionary.php} +#' +#' @export +#' +climdex.fxx <- function(ci,freq=c("annual","halfyear","seasonal","monthly")) { + stopifnot(!is.null(ci@data$wind_gust)) + return(tapply.fast(ci@data$wind_gust, ci@date.factors[[match.arg(freq)]], max, na.rm=TRUE) * + ci@namasks[[match.arg(freq)]]$wind_gust) +} + +#' Northerly winds +#' @description +#' This function computes the climdex index DDnorth: Days with northerly wind (-45deg (315deg) < wind_dir <= 45deg). +#' +#' @param ci Object of type climdexInput (representing the daily mean wind direction in [deg] with 0deg being wind from north) +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear", "seasonal" or "monthly". Default: "annual". +#' @return A vector containing the time series of days with northerly winds. +#' @author Rebekka Posselt (MeteoSwiss) +#' @references \url{http://www.ecad.eu/indicesextremes/indicesdictionary.php} +#' +#' @export +#' +climdex.ddnorth <- function(ci,freq=c("annual","halfyear","seasonal","monthly")) { + stopifnot(!is.null(ci@data$wind_dir)) + wind_northerly <- (number.days.op.threshold(temp=ci@data$wind_dir, date.factor=ci@date.factors[[match.arg(freq)]], + threshold=315, op=">") + + number.days.op.threshold(temp=ci@data$wind_dir, date.factor=ci@date.factors[[match.arg(freq)]], + threshold=45, op="<=")) * ci@namasks[[match.arg(freq)]]$wind_dir + return(wind_northerly) +} + +#' Easterly winds +#' @description +#' This function computes the climdex index DDeast: Days with easterly wind (45deg < wind_dir <= 135deg). +#' +#' @param ci Object of type climdexInput (representing the daily mean wind direction in [deg] with 90deg being wind from east) +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear", "seasonal" or "monthly". Default: "annual". +#' @return A vector containing the time series of days with easterly winds. +#' @author Rebekka Posselt (MeteoSwiss) +#' @references \url{http://www.ecad.eu/indicesextremes/indicesdictionary.php} +#' +#' @export +#' +climdex.ddeast <- function(ci,freq=c("annual","halfyear","seasonal","monthly")) { + stopifnot(!is.null(ci@data$wind_dir)) + wind_easterly <- (number.days.op.threshold(temp=ci@data$wind_dir, date.factor=ci@date.factors[[match.arg(freq)]], + threshold=45, op=">") - + number.days.op.threshold(temp=ci@data$wind_dir, date.factor=ci@date.factors[[match.arg(freq)]], + threshold=135, op=">")) * ci@namasks[[match.arg(freq)]]$wind_dir + return(wind_easterly) +} + +#' Southerly winds +#' @description +#' This function computes the climdex index DDsouth: Days with southerly wind (135deg < wind_dir <= 225deg). +#' +#' @param ci Object of type climdexInput (representing the daily mean wind direction in [deg] with 180deg being wind from south) +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear", "seasonal" or "monthly". Default: "annual". +#' @return A vector containing the time series of days with southerly winds. +#' @author Rebekka Posselt (MeteoSwiss) +#' @references \url{http://www.ecad.eu/indicesextremes/indicesdictionary.php} +#' +#' @export +#' +climdex.ddsouth <- function(ci,freq=c("annual","halfyear","seasonal","monthly")) { + stopifnot(!is.null(ci@data$wind_dir)) + wind_southerly <- (number.days.op.threshold(temp=ci@data$wind_dir, date.factor=ci@date.factors[[match.arg(freq)]], + threshold=135, op=">") - + number.days.op.threshold(temp=ci@data$wind_dir, date.factor=ci@date.factors[[match.arg(freq)]], + threshold=225, op=">")) * ci@namasks[[match.arg(freq)]]$wind_dir + return(wind_southerly) +} + +#' Westerly winds +#' @description +#' This function computes the climdex index DDwest: Days with westerly wind (225deg < wind_dir <= 315deg). +#' +#' @param ci Object of type climdexInput (representing the daily mean wind direction in [deg] with 270deg being wind from west) +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear", "seasonal" or "monthly". Default: "annual". +#' @return A vector containing the time series of days with westerly winds. +#' @author Rebekka Posselt (MeteoSwiss) +#' @references \url{http://www.ecad.eu/indicesextremes/indicesdictionary.php} +#' +#' @export +#' +climdex.ddwest <- function(ci,freq=c("annual","halfyear","seasonal","monthly")) { + stopifnot(!is.null(ci@data$wind_dir)) + wind_westerly <- (number.days.op.threshold(temp=ci@data$wind_dir, date.factor=ci@date.factors[[match.arg(freq)]], + threshold=225, op=">") - + number.days.op.threshold(temp=ci@data$wind_dir, date.factor=ci@date.factors[[match.arg(freq)]], + threshold=315, op=">")) * ci@namasks[[match.arg(freq)]]$wind_dir + return(wind_westerly) +} + + +## Climate indices involving temperature timeseries +#' Min/Max of mean n-day Min/Avg/Max.-Temperature +#' +#' This functions compute the climdex indeces Txndaymax, Txndaymin, +#' Tnndaymax, Tnndaymin, +#' Tmndaymax, Tmndaymin. +#' +#' This function takes a climdexInput object as input and computes the climdex +#' index T[x|m|n]nday[min|max]: min/max of mean n-day [max|avg|min] temperature. +#' +#' @param ci Object of type climdexInput (representing daily Tmax/Tavg/Tmin). +#' @param ndays number of days to consider (default=5). +#' @param freq Time frequency to aggregate to (default="monthly"). +#' @param center.mean.on.last.day Whether to center the n-day running mean on +#' the last day of the window, instead of the center day. +#' +#' @template generic_seealso_references +#' +#' @templateVar cdxvar txndaymax +#' @templateVar cdxdescription a time series of the max of average maximum n-day temperature. +#' @template get_generic_example +#' @author Rebekka Posselt (MeteoSwiss) +#' +#' @name climdex.tnday +NULL + +#' +#' @rdname climdex.tnday +#' @export +#' +climdex.txndaymax <- function(ci, ndays=5, freq=c("monthly", "annual", "halfyear", "seasonal"), center.mean.on.last.day=FALSE) { + stopifnot(!is.null(ci@data$tmax)) + return(nday.consec.temp.mean(ci@data$tmax, ci@date.factors[[match.arg(freq)]], + ndays=ndays, freq.fun="max", center.mean.on.last.day) * + ci@namasks[[match.arg(freq)]]$tmax) +} + +#' +#' @rdname climdex.tnday +#' @export +#' +climdex.txndaymin <- function(ci, ndays=5, freq=c("monthly", "annual", "halfyear", "seasonal"), center.mean.on.last.day=FALSE) { + stopifnot(!is.null(ci@data$tmax)) + return(nday.consec.temp.mean(ci@data$tmax, ci@date.factors[[match.arg(freq)]], + ndays=ndays, freq.fun="min", center.mean.on.last.day) * + ci@namasks[[match.arg(freq)]]$tmax) +} + +#' +#' @rdname climdex.tnday +#' @export +#' +climdex.tnndaymax <- function(ci, ndays=5, freq=c("monthly", "annual", "halfyear", "seasonal"), center.mean.on.last.day=FALSE) { + stopifnot(!is.null(ci@data$tmin)) + return(nday.consec.temp.mean(ci@data$tmin, ci@date.factors[[match.arg(freq)]], + ndays=ndays, freq.fun="max", center.mean.on.last.day) * + ci@namasks[[match.arg(freq)]]$tmin) +} + +#' +#' @rdname climdex.tnday +#' @export +#' +climdex.tnndaymin <- function(ci, ndays=5, freq=c("monthly", "annual", "halfyear", "seasonal"), center.mean.on.last.day=FALSE) { + stopifnot(!is.null(ci@data$tmin)) + return(nday.consec.temp.mean(ci@data$tmin, ci@date.factors[[match.arg(freq)]], + ndays=ndays, freq.fun="min", center.mean.on.last.day) * + ci@namasks[[match.arg(freq)]]$tmin) +} + +#' +#' @rdname climdex.tnday +#' @export +#' +climdex.tmndaymax <- function(ci, ndays=5, freq=c("monthly", "annual", "halfyear", "seasonal"), center.mean.on.last.day=FALSE) { + stopifnot(!is.null(ci@data$tavg)) + return(nday.consec.temp.mean(ci@data$tavg, ci@date.factors[[match.arg(freq)]], + ndays=ndays, freq.fun="max", center.mean.on.last.day) * + ci@namasks[[match.arg(freq)]]$tavg) +} + +#' +#' @rdname climdex.tnday +#' @export +#' +climdex.tmndaymin <- function(ci, ndays=5, freq=c("monthly", "annual", "halfyear", "seasonal"), center.mean.on.last.day=FALSE) { + stopifnot(!is.null(ci@data$tavg)) + return(nday.consec.temp.mean(ci@data$tavg, ci@date.factors[[match.arg(freq)]], + ndays=ndays, freq.fun="min", center.mean.on.last.day) * + ci@namasks[[match.arg(freq)]]$tavg) +} + +#' Percent of Values Above nth Percentile Daily Maximum Temperature +#' +#' This function computes the climdex index TXnp. +#' +#' This function takes a climdexInput object as input and computes the +#' monthly or annual percent of values above/below the nth percentile of baseline +#' daily maximum temperature. +#' +#' @template threshold_indices_block +#' @template threshold_indices_args +#' @param quant Quantile to used (must be stated and calculated within \code{\link{climdexInput.raw}}). Default: 0.9 +#' @param op Operator to use for comparison. Default: ">" +#' @template missing_values_caveat +#' +#' @template generic_seealso_references +#' +#' @export +#' +climdex.txnp <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal"), quant=0.9, op=">") { + quant.str <- paste0("q",as.integer(quant*100)) + stopifnot(!is.null(ci@data$tmax) && !is.null(ci@quantiles$tmax) && !is.null(ci@quantiles$tmax$outbase[[quant.str]])) + return(percent.days.op.threshold(ci@data$tmax, ci@dates, ci@jdays, ci@date.factors[[match.arg(freq)]], + ci@quantiles$tmax$outbase[[quant.str]], ci@quantiles$tmax$inbase[[quant.str]], ci@base.range, op, + ci@max.missing.days[match.arg(freq)]) * ci@namasks[[match.arg(freq)]]$tmax) +} + +#' Percent of Values Above nth Percentile Daily Minimum Temperature +#' +#' This function computes the climdex index TNnp. +#' +#' This function takes a climdexInput object as input and computes the +#' monthly or annual percent of values above the 90th percentile of baseline +#' daily minimum temperature. +#' +#' @template threshold_indices_block +#' @template threshold_indices_args +#' @param quant Quantile to used (must be stated and calculated within \code{\link{climdexInput.raw}}). Default: 0.9 +#' @param op Operator to use for comparison. Default: ">" +#' @template missing_values_caveat +#' +#' @template generic_seealso_references +#' +#' @export +#' +climdex.tnnp <- function(ci, freq=c("monthly", "annual", "halfyear", "seasonal"), quant=0.9, op=">") { + quant.str <- paste0("q",as.integer(quant*100)) + stopifnot(!is.null(ci@data$tmin) && !is.null(ci@quantiles$tmin) && !is.null(ci@quantiles$tmin$outbase[[quant.str]])) + return(percent.days.op.threshold(ci@data$tmin, ci@dates, ci@jdays, ci@date.factors[[match.arg(freq)]], + ci@quantiles$tmin$outbase[[quant.str]], ci@quantiles$tmin$inbase[[quant.str]], ci@base.range, op, + ci@max.missing.days[match.arg(freq)]) * ci@namasks[[match.arg(freq)]]$tmin) +} + +## Climate indices involving sun timeseries +#' Sunshine duration +#' +#' This function computes the climdex index SS: The sum of the sunshine duration hours within a period. +#' +#' @param ci Object of type climdexInput (representing the daily sunshine duration in [h]) +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear", "seasonal" or "monthly". Default: "annual". +#' @return A vector containing the time series of sunshine duration. +#' +#' @template generic_seealso_references +#' @author Rebekka Posselt (MeteoSwiss) +#' +#' @export +climdex.ss <- function(ci,freq=c("annual","halfyear","seasonal","monthly")) { + stopifnot(!is.null(ci@data$sun)) + return(tapply.fast(ci@data$sun, ci@date.factors[[match.arg(freq)]], sum, na.rm=TRUE) * ci@namasks[[match.arg(freq)]]$sun) +} + +#' Relative sunshine duration +#' +#' This function computes the climdex index SUN_RELMEAN. +#' +#' This function takes a climdexInput object as input and computes the +#' SUN_RELMEAN index: The mean of the relative sunshine duration within a period. +#' +#' @param ci Object of type climdexInput (representing the daily relative sunshine duration in [\%]) +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear", "seasonal" or "monthly". Default: "annual". +#' @return A vector containing the time series of the relative sunshine duration. +#' +#' @template generic_seealso_references +#' +#' @export +climdex.sun_relmean <- function(ci,freq=c("annual","halfyear","seasonal","monthly")) { + stopifnot(!is.null(ci@data$sun_rel)) + return(tapply.fast(ci@data$sun_rel, ci@date.factors[[match.arg(freq)]], mean, na.rm=TRUE) * ci@namasks[[match.arg(freq)]]$sun_rel) +} + +#' Mostly cloudy days +#' +#' This function computes the climdex index SUN_CLOUDY. +#' +#' This function takes a climdexInput object as input and computes the +#' SUN_CLOUDY index: The number of mostly cloudy days (relative sunshine duration < 20 \%) within a period. +#' +#' @param ci Object of type climdexInput (representing the daily relative sunshine duration in [\%]) +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear", "seasonal" or "monthly". Default: "annual". +#' @return A vector containing the number of mostly cloudy days. +#' +#' @template generic_seealso_references +#' +#' @export +climdex.sun_cloudy <- function(ci,freq=c("annual","halfyear","seasonal","monthly")) { + stopifnot(!is.null(ci@data$sun_rel)) + return(number.days.op.threshold(temp=ci@data$sun_rel, date.factor=ci@date.factors[[match.arg(freq)]], + threshold=20, op="<") * ci@namasks[[match.arg(freq)]]$sun_rel) +} + +#' Mostly sunny days +#' +#' This function computes the climdex index SUN_SUNNY. +#' +#' This function takes a climdexInput object as input and computes the +#' SUN_SUNNY index: The number of mostly sunny days (relative sunshine duration > 80 \%) within a period. +#' +#' @param ci Object of type climdexInput (representing the daily relative sunshine duration in [\%]) +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear", "seasonal" or "monthly". Default: "annual". +#' @return A vector containing the number of mostly sunny days. +#' +#' @template generic_seealso_references +#' +#' @export +climdex.sun_sunny <- function(ci,freq=c("annual","halfyear","seasonal","monthly")) { + stopifnot(!is.null(ci@data$sun)) + return(number.days.op.threshold(temp=ci@data$sun_rel, date.factor=ci@date.factors[[match.arg(freq)]], + threshold=80, op="<") * ci@namasks[[match.arg(freq)]]$sun) +} + +## Climate indices involving snow timeseries +#' Number of snow days +#' +#' This function computes the climdex index SDD: The number of days with a snow depth > Y cm. +#' +#' @param ci Object of type climdexInput (representing the daily snow depth timeseries in [cm]). +#' @param threshold Snow depth threshhold [in cm]. (Default: 1) +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear" or "seasonal". Default: "annual". +#' @return A vector containing an annual timeseries of the number of snow days. +#' @author Rebekka Posselt (MeteoSwiss) +#' @references \url{http://www.ecad.eu/indicesextremes/indicesdictionary.php} +#' +#' @export +#' +climdex.sdd <- function(ci,threshold=1,freq=c("annual", "monthly", "halfyear","seasonal")) { + stopifnot(!is.null(ci@data$snow)) + return(number.days.op.threshold(temp=ci@data$snow, date.factor=ci@date.factors[[match.arg(freq)]], + threshold=threshold, op=">=") * ci@namasks[[match.arg(freq)]]$snow) +} + +#' Maximum snow depth +#' +#' This function computes the climdex index SDx: The maximum snow depth measured within a year. +#' +#' @param ci Object of type climdexInput (representing the daily snow depth timeseries in [cm]). +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear" or "seasonal". Default: "annual". +#' @return A vector containing an annual timeseries of the maximum snow depth. +#' @author Rebekka Posselt (MeteoSwiss) +#' +#' @export +#' +climdex.sdx <- function(ci,freq=c("annual","monthly", "halfyear","seasonal")) { + stopifnot(!is.null(ci@data$snow)) + return(tapply.fast(ci@data$snow, ci@date.factors[[match.arg(freq)]], max, na.rm=TRUE) * + ci@namasks[[match.arg(freq)]]$snow) +} + +#' Mean snow depth +#' +#' This function computes the climdex index SD: The mean snow depth measured within a period. +#' +#' @param ci Object of type climdexInput (representing the daily snow depth timeseries in [cm]). +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear" or "seasonal". Default: "annual". +#' @return A vector containing the timeseries of the mean snow depth. +#' @author Rebekka Posselt (MeteoSwiss) +#' @references \url{http://www.ecad.eu/indicesextremes/indicesdictionary.php} +#' +#' @export +#' +climdex.sd <- function(ci,freq=c("annual","monthly", "halfyear","seasonal")) { + stopifnot(!is.null(ci@data$snow)) + return(tapply.fast(ci@data$snow, ci@date.factors[[match.arg(freq)]], mean, na.rm=TRUE) * + ci@namasks[[match.arg(freq)]]$snow) +} + +#' Number of new snow days +#' +#' This function computes the climdex index NSD: The number of days with a (new) snowfall > Y cm. +#' +#' @param ci Object of type climdexInput (representing the daily snowfall timeseries in [cm]). +#' @param threshold Snow depth threshhold [in cm]. (Default: 1) +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear" or "seasonal". Default: "annual". +#' @return A vector containing the timeseries of the number of new snow days. +#' @author Rebekka Posselt (MeteoSwiss) +#' +#' @export +#' +climdex.nsd <- function(ci,threshold=1,freq=c("annual","monthly", "halfyear","seasonal")) { + stopifnot(!is.null(ci@data$snow_new)) + return(number.days.op.threshold(temp=ci@data$snow_new, date.factor=ci@date.factors[[match.arg(freq)]], + threshold=threshold, op=">=") * ci@namasks[[match.arg(freq)]]$snow_new) +} + +#' Maximum new snow +#' +#' This function computes the climdex index NSX: The maximum snowfall measured within a year. +#' +#' @param ci Object of type climdexInput (representing the daily snowfall timeseries in [cm]). +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear" or "seasonal". Default: "annual". +#' @return A vector containing the timeseries of the maximum snowfall. +#' @author Rebekka Posselt (MeteoSwiss) +#' +#' @export +#' +climdex.nsx <- function(ci,freq=c("annual","monthly", "halfyear","seasonal")) { + stopifnot(!is.null(ci@data$snow_new)) + return(tapply.fast(ci@data$snow_new, ci@date.factors[[match.arg(freq)]], max, na.rm=TRUE) * + ci@namasks[[match.arg(freq)]]$snow_new) +} + +#' New snow sum +#' +#' This function computes the climdex index NSS: The sum of all snowfall measured within a year. +#' +#' @param ci Object of type climdexInput (representing the daily snowfall timeseries in [cm]). +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear" or "seasonal". Default: "annual". +#' @return A vector containing the timeseries of the snowfall sum. +#' @author Rebekka Posselt (MeteoSwiss) +#' +#' @export +#' +climdex.nss <- function(ci,freq=c("annual","monthly", "halfyear","seasonal")) { + stopifnot(!is.null(ci@data$snow_new)) + return(tapply.fast(ci@data$snow_new, ci@date.factors[[match.arg(freq)]], sum, na.rm=TRUE) * + ci@namasks[[match.arg(freq)]]$snow_new) +} + +## Climate indices involving cloud timeseries +#' Mean cloud cover +#' +#' This function computes the climdex index CC +#' +#' This function takes a climdexInput object as input and computes the +#' CLOUD_MEAN index: The mean cloud cover measured within a period. +#' +#' @param ci Object of type climdexInput (representing the daily mean cloud cover in octa or percent) +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear", "seasonal" or "monthly". Default: "annual". +#' @return A vector containing the time series of mean cloud cover. +#' @author Rebekka Posselt (MeteoSwiss) +#' @references \url{http://www.ecad.eu/indicesextremes/indicesdictionary.php} +#' +#' @export +climdex.cc <- function(ci,freq=c("annual","halfyear","seasonal","monthly")) { + stopifnot(!is.null(ci@data$cloud)) + return(tapply.fast(ci@data$cloud, ci@date.factors[[match.arg(freq)]], mean, na.rm=TRUE) * ci@namasks[[match.arg(freq)]]$cloud) +} + +#' Mostly cloudy days +#' +#' This function computes the climdex index CC6: The number of mostly cloudy days (Cloud cover >= 6octa / 80%) within a period. +#' +#' @param ci Object of type climdexInput (representing the daily mean cloud cover in octa or percent) +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear", "seasonal" or "monthly". Default: "annual". +#' @param unit unit of the cloud cover. Allowed are: "octa" or "percent". Default: "octa". +#' @return A vector containing the number of mostly cloudy days. +#' @author Rebekka Posselt (MeteoSwiss) +#' @references \url{http://www.ecad.eu/indicesextremes/indicesdictionary.php} +#' +#' @export +climdex.cc6 <- function(ci,freq=c("annual","halfyear","seasonal","monthly"),unit="octa") { + stopifnot(!is.null(ci@data$cloud)) + if (unit=="octa"){ + threshold=6 + }else if (unit=="percent"){ + threshold=80 + } else{ + stop("unit should be either 'octa' or 'percent'") + } + return(number.days.op.threshold(temp=ci@data$cloud, date.factor=ci@date.factors[[match.arg(freq)]], + threshold=threshold, op=">=") * ci@namasks[[match.arg(freq)]]$cloud) +} + +#' Mostly sunny days +#' +#' This function computes the climdex index CC2: The number of mostly sunny days (Cloud cover <= 2octa / 20%) within a period. +#' +#' @param ci Object of type climdexInput (representing the daily mean cloud cover in octa or percent) +#' @param freq Time frequency to aggregate to. Allowed are: "annual","halfyear", "seasonal" or "monthly". Default: "annual". +#' @param unit unit of the cloud cover. Allowed are: "octa" or "percent". Default: "octa". +#' @return A vector containing the number of mostly sunny days. +#' @author Rebekka Posselt (MeteoSwiss) +#' @references \url{http://www.ecad.eu/indicesextremes/indicesdictionary.php} +#' +#' @export +climdex.cc2 <- function(ci,freq=c("annual","halfyear","seasonal","monthly"),unit="octa") { + stopifnot(!is.null(ci@data$cloud)) + if (unit=="octa"){ + threshold=2 + }else if (unit=="percent"){ + threshold=20 + } else{ + stop("unit should be either 'octa' or 'percent'") + } + return(number.days.op.threshold(temp=ci@data$cloud, date.factor=ci@date.factors[[match.arg(freq)]], + threshold=threshold, op="<=") * ci@namasks[[match.arg(freq)]]$cloud) +} + +# ###################################################################################################################### +# #' HUGLIN INDEX (only for climdex.pcic.ncdf) +#' Introduces by C.Photiadou (KNMI) +# #' This function is not ready yet. It is uses coefficient based on the latitude +# #' For this index I had to curry the cdx.funcs to be able to include the subset. Later I realised +# #' I need also to include the latitude. This would be used in compute.indices.for.stripe together with get.lat +# #' to retrieve subset & latitude +# #' I didn't proceed with finisheing the eca.HI function. I thought to ask you first if its possible +# #' to adapt compute.indices.for.stripe so it can include the currying and the latitude. Or if you had a better idea on this +# #' please let me know. +# +# #' Function to Curry a cxd.funcs for subset (now at cur_sub) +# #' used only for Huglin Index +# curry_in_subset_for_huglin <- function(cdx.funcs, cur_sub){ +# cdx.names = names(cdx.funcs) +# cdx.funcs <- lapply(cdx.names, function(function_name) { +# f = cdx.funcs[[function_name]] +# if(grepl('^hi', function_name)) { +# return(functional::Curry(f, cur_sub = cur_sub)) +# } else { +# return(f) +# } +# }) +# names(cdx.funcs) = cdx.names +# return(cdx.funcs) +# } + +# ### Get latitude function +# get.lat <- function(open_file_list, variable.name.map) { +# #var.name <- variable.name.map[[names(v.f.idx)[1]]] +# y.dim <- ncdf4.helpers::nc.get.dim.for.axis(open_file_list[[1]], variable.name.map, "Y") +# return(y.dim$vals) +# } + +#' Huglin Index +#' +#' This function computes the climdex index HI: +#' @param ci Object of type climdexInput (representing the daily mean and daily max temperature) +#' @param freq Time frequency to aggregate to. Allowed are only "annual" +#' @param unit. Allowed are only deg C. +#' @return A vector containing the HI +#' @author Rebekka Posselt (MeteoSwiss) +#' @references \url{http://www.ecad.eu/indicesextremes/indicesdictionary.php} +#' +#' @export +climdex.HI <- function(ci,freq=c("annual"),cur_sub){ + + tempavg <- ci@data$tavg + tempmax <- ci@data$tmax + + month.series <- get.months(ci@dates) + year.series <- get.years(ci@dates) + valid.months <- month.series >=4 & month.series <=9 + hi_coef <- if (cur_sub <=40) {hi_coeff <- 1 + }else if(cur_sub >40 & cur_sub <42) {hi_coef <- 1.02 + }else if(cur_sub >42 & cur_sub <44) {hi_coef <- 1.03 + }else if(cur_sub >44 & cur_sub <48) {hi_coef <- 1.04 + }else if(cur_sub >46 & cur_sub <48) {hi_coef <- 1.05 + }else if(cur_sub >48 & cur_sub <50) {hi_coef <- 1.06 + }else if(cur_sub >=50){hi_coef <- 1} + valid.sel<- year.series[valid.months] + tempdata <- ((((tempavg -10) + (tempmax -10)) /2) * hi_coef) + dat_final <- tempdata[valid.months] + + return(tapply(dat_final,valid.sel,sum)) +} +####################################################################33 +## Functions for eca\&d station files +## -- introduced by C. Photiadou (KNMI), November 2015 +# Secondary function for eca.input +# +# Typical ECA\&D files contain a header and this functions finds the line number where the data starts. +# +find.start.of.data.index = function(fname) { + matched.indices = which(grepl('^SOUID', gsub(" ", "", readLines(fname, n = 50)))) + if (length(matched.indices) > 1) stop('ECA fileformat error: cannot determine start of data, multiple header lines') + if (length(matched.indices) == 0) stop('ECA fileformat error: cannot find start of data: cannot find header line') + return(matched.indices - 1) +} + +## Functions for eca\&d station files +#' ECA &D station data files +#' +#' This function reads and prepares the station data files from the European Climate Assessment and Dataset (ECA&D) for use in climdex +#' +#' @param filename File name and path of station data file. +#' @param var.name A varialbe name from the ECA\&D variavle list: prec (RR), tavg (TG), tmax (TX), tmin(TN), sun (SS), +#' wind_gust (FX), wind (FG). +#' @param data.name Always DATE +#' @return A data frame containing two columns: DATE with dates in PCICt format and "var.name" the varialbe name. +#' @template get_generic_example +#' @author Christiana Photiadou (KNMI) +#' @references \url{http://www.ecad.eu/} +#' @export +eca.input <- function(filename, var.name, date.name){ + + ifile.eca <- read.table(filename, skip=find.start.of.data.index(filename), sep=",", header=T)[,c(date.name,var.name)] + + ifile.eca[[date.name]] <- as.PCICt(strptime(as.character(ifile.eca[[date.name]]),"%Y%m%d"), cal="gregorian") + + ifile.eca[[var.name]][ifile.eca[[var.name]]==-9999]<- NA + + if(!(var.name %in% c("TG", "TX" ,"TN", "RR", "SS", "FX", "FG"))){ + return(ifile.eca) + } else + ifile.eca[[var.name]] <- ifile.eca[[var.name]]*0.1 + + return(ifile.eca) +} #' Get available indices by name #' @@ -1296,10 +2633,19 @@ climdex.prcptot <- function(ci) { stopifnot(!is.null(ci@data$prec)); return(tota #' func.names <- climdex.get.available.indices(ci) #' @export climdex.get.available.indices <- function(ci, function.names=TRUE) { - available.indices <- list(tmax=c('su', 'id', 'txx', 'txn', 'tx10p', 'tx90p', 'wsdi'), - tmin=c('fd', 'tr', 'tnx', 'tnn', 'tn10p', 'tn90p', 'csdi'), - tavg=c('gsl', 'dtr'), - prec=c('rx1day', 'rx5day', 'sdii', 'r10mm', 'r20mm', 'rnnmm', 'cdd', 'cwd', 'r95ptot', 'r99ptot', 'prcptot')) + available.indices <- list(tmax=c('su', 'id', 'txx', 'txn', 'tx10p', 'tx90p', 'wsdi', 'csu', 'txndaymin','txndaymax'), + tmin=c('fd', 'tr', 'tnx', 'tnn', 'tn10p', 'tn90p', 'csdi', 'cfd', 'tnndaymin','tnndaymax'), + tavg=c('gsl', 'dtr', 'hd17', 'tmndaymin','tmndaymax', 'cd', 'cw', 'wd', 'ww'), + prec=c('rx1day', 'rx5day', 'sdii', 'r10mm', 'r20mm', 'rnnmm', 'cdd', 'cwd', + 'r75p', 'r95p', 'r99p', 'r75ptot', 'r95ptot', 'r99ptot', 'prcptot', 'spi3', 'spi6'), + snow=c('sdd','sdx','sd'), + snow_new=c('nsd','nsx','nss'), + wind=c("fg","fgcalm","fg6bft"), + wind_gust=c('fxstorm','fxx'), + wind_dir=c('ddnorth','ddeast','ddsouth','ddwest'), + cloud=c('cc','cc6','cc2'), + sun=c("ss"), + sun_rel=c("sun_cloudy","sun_sunny","sun_relmean")) if(function.names) { return(paste("climdex", unlist(available.indices[names(ci@data)]), sep=".")) } else { @@ -1312,7 +2658,7 @@ climdex.get.available.indices <- function(ci, function.names=TRUE) { ## #' Get series length at ends -#' +#' #' This function takes a series of boolean values and returns a list of #' integers of the same length corresponding to the lengths at the ends of #' sequences of TRUE values. @@ -1391,11 +2737,13 @@ number.days.op.threshold <- function(temp, date.factor, threshold, op="<") { return(tapply.fast(match.fun(op)(temp, threshold), date.factor, sum, na.rm=TRUE)) } -#' Flexible GSL function +#' @title Flexible GSL function #' +#' @description #' This function computes the growing season length (GSL) given the input, #' which is allowed to vary considerably from the ETCCDI definitions. #' +#' @details #' This function is the function used to implement \code{\link{climdex.gsl}}. #' It's designed to be flexible to allow for experimentation and testing of new #' thresholds and methods. @@ -1567,11 +2915,13 @@ percent.days.op.threshold <- function(temp, dates, jdays, date.factor, threshold return(ret) } -#' Sum of spell lengths exceeding daily threshold +#' @title Sum of spell lengths exceeding daily threshold #' +#' @description #' This function returns the number of spells of more than \code{min.length} #' days which exceed or are below the given threshold. #' +#' @details #' This routine compares data to the thresholds using the given operator, #' generating a series of TRUE or FALSE values; these values are then filtered #' to remove any sequences of less than \code{min.length} days of TRUE values. @@ -1633,10 +2983,12 @@ mean.daily.temp.range <- function(daily.max.temp, daily.min.temp, date.factor) { dat } -#' Number of days (less than, greater than, etc) a threshold +#' @title Number of days (less than, greater than, etc) a threshold #' +#' @description #' Produces sums of values that exceed (or are below) the specified threshold. #' +#' @details #' This function takes a data series, the number of days in the running window, #' a date factor to aggregate by, and an optional modifier parameter #' (center.mean.on.last.day). It computes the n-day running sum of @@ -1674,15 +3026,74 @@ mean.daily.temp.range <- function(daily.max.temp, daily.min.temp, date.factor) { nday.consec.prec.max <- function(daily.prec, date.factor, ndays, center.mean.on.last.day=FALSE) { if(ndays == 1) { return(suppressWarnings(tapply.fast(daily.prec, date.factor, max, na.rm=TRUE))) - } else { - ## Ends of the data will be de-emphasized (padded with zero precip data); NAs replaced with 0 - new.series <- c(rep(0, floor(ndays / 2)), daily.prec, rep(0, floor(ndays / 2))) - new.series[is.na(new.series)] <- 0 - prec.runsum <- runmean(new.series, k=ndays, endrule="trim") * ndays - if(center.mean.on.last.day) - prec.runsum <- c(rep(0, floor(ndays / 2)), prec.runsum[1:(length(prec.runsum) - floor(ndays / 2))]) - return(tapply.fast(prec.runsum, date.factor, max)) } + ## Ends of the data will be de-emphasized (padded with zero precip data); NAs replaced with 0 + daily.prec[is.na(daily.prec)] <- 0 + prec.runsum <- runmean(daily.prec, k=ndays, endrule="NA") + prec.runsum[is.na(prec.runsum)] <- 0 + if(center.mean.on.last.day) { + k2 = ndays %/% 2 + prec.runsum <- c(rep(0, k2), prec.runsum[1:(length(prec.runsum) - k2)]) + } + return(tapply.fast(prec.runsum, date.factor, max) * ndays) +} +#' @title Number of days (less than, greater than, etc) a threshold +#' +#' @description +#' Produces means of values that exceed (or are below) the specified threshold. +#' +#' @details +#' This function takes a data series, the number of days in the running window, +#' a date factor to aggregate by, and an optional modifier parameter +#' (center.mean.on.last.day). It computes the n-day running mean of +#' temperature and returns the mean n-day temperature per unit +#' time, as defined by \code{date.factor}. +#' +#' @param daily.temp Daily timeseries of temperature (tmax, tavg or tmin). +#' @param date.factor Factor to aggregate by. +#' @param ndays Number of days in the running window. +#' @param freq.fun Function to apply on the choosen period (e.g., max, min, ...) +#' @param center.mean.on.last.day Whether to center the n-day running mean on +#' the last day of the series, instead of the middle day. +#' @return A vector consisting of the mean n-day temp per +#' time interval. +#' @keywords ts climate +#' +#' @examples +#' library(PCICt) +#' +#' ## Parse the dates into PCICt. +#' tmax.dates <- as.PCICt(do.call(paste, ec.1018935.tmax[,c("year","jday")]), format="%Y %j", cal="gregorian") +#' tmin.dates <- as.PCICt(do.call(paste, ec.1018935.tmin[,c("year","jday")]), format="%Y %j", cal="gregorian") +#' prec.dates <- as.PCICt(do.call(paste, ec.1018935.prec[,c("year","jday")]), format="%Y %j", cal="gregorian") +#' +#' ## Load the data in. +#' ci <- climdexInput.raw(tmax=ec.1018935.tmax$MAX_TEMP, +#' tmin=ec.1018935.tmin$MIN_TEMP, +#' prec=ec.1018935.prec$ONE_DAY_PRECIPITATION, +#' tmax.dates=tmax.dates, +#' tmin.dates=tmin.dates, +#' prec.dates=prec.dates, +#' base.range=c(1971, 2000)) +#' +#' ## Compute txnday on a monthly basis. +#' tx5daymax <- nday.consec.temp.mean(ci@@data$tmax, ci@@date.factors$monthly, freq.fun="max", ndays=5) +#' +#' @export +#' +nday.consec.temp.mean <- function(daily.temp, date.factor, ndays, freq.fun="max", center.mean.on.last.day=FALSE) { + if(ndays == 1) { + return(suppressWarnings(tapply.fast(daily.temp, date.factor, mean, na.rm=TRUE))) + } + ## Ends of the data will be de-emphasized (padded with zero precip data); NAs replaced with 0 + daily.temp[is.na(daily.temp)] <- 0 + temp.runmean <- runmean(daily.temp, k=ndays, endrule="NA") + temp.runmean[is.na(temp.runmean)] <- 0 + if(center.mean.on.last.day) { + k2 = ndays %/% 2 + temp.runmean <- c(rep(0, k2), temp.runmean[1:(length(temp.runmean) - k2)]) + } + return(tapply.fast(temp.runmean, date.factor, match.fun(freq.fun))) } #' Simple Precipitation Intensity Index @@ -1710,11 +3121,13 @@ simple.precipitation.intensity.index <- function(daily.prec, date.factor) { return(tapply.fast(daily.prec, date.factor, function(prec) { idx <- prec >= 1 & !is.na(prec); if(sum(idx) == 0) { return(0); } else { return(sum(prec[idx], na.rm=TRUE) / sum(idx)) } } )) } -#' Maximum spell length +#' @title Maximum spell length #' +#' @description #' This function returns the longest string of days which exceed or are below #' the given threshold. #' +#' @details #' This routine compares data to the threshold using the given operator, #' generating a series of TRUE or FALSE values. It then computes the lengths of #' sequences of TRUE values (spells) and chooses the longest spell in each @@ -1789,6 +3202,7 @@ total.precip.op.threshold <- function(daily.prec, date.factor, threshold, op) { return(tapply.fast(daily.prec, date.factor, function(pr) { return(sum(pr[f(pr, threshold)], na.rm=TRUE)) } )) } + ## Returns an n-day running quantile for each day of data (dimensions c(dpy, q)) running.quantile <- function(data, n, q, dpy, min.fraction) { ret <- .Call("running_quantile_windowed", data, n, q, dpy, min.fraction, PACKAGE='climdex.pcic') diff --git a/R/get_climdex_function_eobs.R b/R/get_climdex_function_eobs.R new file mode 100644 index 0000000..932b2ab --- /dev/null +++ b/R/get_climdex_function_eobs.R @@ -0,0 +1,277 @@ + +#' Returns a list of Climdex variables given constraints +#' Added: monthly only for SPI indices +#' Changed annual only is changed +#' Added huglin index in tavg although it takes two input variables (tmax and tavg) +#' Changed fre.list to accept also only months and added helper_fun for same reason +get.climdex.variable.list.eobs <- function(source.data.present, time.resolution=c("all", "annual", "monthly"), climdex.vars.subset=NULL) { + time.res <- match.arg(time.resolution) + #annual.only <- c("fdETCCDI", "suETCCDI", "idETCCDI", "trETCCDI", "gslETCCDI", "wsdiETCCDI", "csdiETCCDI", "sdiiETCCDI", "r10mmETCCDI", "r20mmETCCDI", "r1mmETCCDI", "cddETCCDI", "cwdETCCDI", "r95pETCCDI", "r99pETCCDI", "prcptotETCCDI", "altcddETCCDI", "altcwdETCCDI", "altcsdiETCCDI", "altwsdiETCCDI") + annual.only <- c("hiETCCDI","gslETCCDI", "wsdiETCCDI", "csdiETCCDI", "cddETCCDI", "cwdETCCDI", "altcddETCCDI", "altcwdETCCDI", "altcsdiETCCDI", "altwsdiETCCDI") + monthly.only <- c("spi3ETCCDI", "spi6ETCCDI") + vars.by.src.data.reqd <- list(tmax=c("csuETCCDI", "suETCCDI", "idETCCDI", "txxETCCDI", "txnETCCDI", "tx10pETCCDI", "tx90pETCCDI", "wsdiETCCDI", "altwsdiETCCDI"), + tmin=c("cfdETCCDI", "fdETCCDI", "trETCCDI", "tnxETCCDI", "tnnETCCDI", "tn10pETCCDI", "tn90pETCCDI", "csdiETCCDI", "altcsdiETCCDI"), + prec=c("spi3ETCCDI", "spi6ETCCDI", "rx1dayETCCDI", "rx5dayETCCDI", "sdiiETCCDI", "r10mmETCCDI", "r20mmETCCDI", "r1mmETCCDI", + "cddETCCDI", "cwdETCCDI", "r95pETCCDI", "r99pETCCDI", "prcptotETCCDI", "altcddETCCDI", "altcwdETCCDI"), + tavg=c("hd17ETCCDI", "hiETCCDI", "gslETCCDI", "dtrETCCDI") ) + + if(any(!(source.data.present %in% c("tmin", "tmax", "tavg", "prec")))) + stop("Invalid variable listed in source.data.present.") + + if(all(c("tmax", "tmin") %in% source.data.present) && !("tavg" %in% source.data.present)) + source.data.present <- c(source.data.present, "tavg") + + climdex.vars <- unlist(vars.by.src.data.reqd[source.data.present]) + if(!is.null(climdex.vars.subset)) + climdex.vars <- climdex.vars[climdex.vars %in% paste(climdex.vars.subset, "ETCCDI", sep="")] + + #freq.lists <- list(c("mon", "yr"), c("yr")) + freq.lists <- list(c("mon", "yr"), c("yr"),c("mon")) + + helper_fun <- function(climdex.vars, annual.only, month.only) { + if (climdex.vars %in% annual.only) { + return(paste(climdex.vars, 'yr', sep = '_')) + } else if (climdex.vars %in% monthly.only) { + return(paste(climdex.vars, 'mon', sep = '_')) + } else { + return(paste(climdex.vars, c('mon', 'yr'), sep = '_')) + } + } + +# dat <- switch(time.res, +# all=unlist(lapply(climdex.vars, function(x) { paste(x, freq.lists[[(x %in% annual.only) + 1]], sep="_") })), +# annual=paste(climdex.vars, "yr", sep="_"), +# monthly=paste(climdex.vars[!(climdex.vars %in% annual.only)], "mon", sep="_")) + + dat <- switch(time.res, + all=unlist(lapply(climdex.vars, helper_fun, annual.only = annual.only, month.only = month.only)), + annual=paste(climdex.vars[!(climdex.vars %in% monthly.only)], "yr", sep="_"), + monthly=paste(climdex.vars[!(climdex.vars %in% annual.only)], "mon", sep="_")) + + names(dat) <- NULL + + return(dat) +} + + +#' Returns a list of Climdex functions, with parameters curried in. +#' Create func.names.eca, options.eca and func.eca +#' +get.climdex.functions.eobs <- function(vars.list, fclimdex.compatible=TRUE) { + func.names <- c("climdex.fd", "climdex.su", "climdex.id", "climdex.tr", "climdex.gsl", + "climdex.fd", "climdex.su", "climdex.id", "climdex.tr", + + "climdex.txx", "climdex.tnx", "climdex.txn", "climdex.tnn", "climdex.tn10p", "climdex.tx10p", "climdex.tn90p", "climdex.tx90p", + "climdex.txx", "climdex.tnx", "climdex.txn", "climdex.tnn", "climdex.tn10p", "climdex.tx10p", "climdex.tn90p", "climdex.tx90p", + + "climdex.wsdi", "climdex.csdi", + + "climdex.dtr", "climdex.rx1day", "climdex.rx5day", + "climdex.dtr", "climdex.rx1day", "climdex.rx5day", + + "climdex.sdii", "climdex.r10mm", "climdex.r20mm", "climdex.rnnmm", "climdex.cdd", "climdex.cwd", "climdex.r95ptot", "climdex.r99ptot", "climdex.prcptot", + "climdex.sdii", "climdex.r10mm", "climdex.r20mm", "climdex.rnnmm", "climdex.r95ptot", "climdex.r99ptot", "climdex.prcptot", + + "climdex.cdd", "climdex.cwd", "climdex.csdi", "climdex.wsdi") + + func.names.eca <- c("eca.csu", "eca.csu", "eca.cfd", "eca.cfd", "eca.hd17", "eca.hd17", "eca.HI", "eca.SPI3", "eca.SPI6") + + el <- list() + af <- list(freq="annual") + mf <- list(freq="monthly") + cwdd.opts <- list(spells.can.span.years=TRUE) + altcwdd.opts <- list(spells.can.span.years=FALSE) + wcsdi.opts <- list(spells.can.span.years=FALSE) + altwcsdi.opts <- list(spells.can.span.years=TRUE) + rx5day.opts <- list(center.mean.on.last.day=fclimdex.compatible) + r1mm.opts <- list(threshold=1) + + options.climdex <- list( af, af, af, af, el, + mf, mf, mf, mf, + + mf, mf, mf, mf, mf, mf, mf, mf, + af, af, af, af, af, af, af, af, + + wcsdi.opts, wcsdi.opts, + + mf, mf, c(mf, rx5day.opts), + af, af, c(af, rx5day.opts), + + af, af, af, c(af, r1mm.opts), cwdd.opts, cwdd.opts, af, af, af, + mf, mf, mf, c(mf, r1mm.opts), mf, mf, mf, + altcwdd.opts, altcwdd.opts, altwcsdi.opts, altwcsdi.opts) + + + options.eca <- list(af,mf, af, mf, af, mf, af, mf, mf) + + func.eca <- lapply(1:(length(func.names.eca)), function(n) do.call(functional::Curry, c(list(getFromNamespace(func.names.eca[n], 'ecaclimdex')), options.eca[[n]]))) + + func.climdex <- lapply(1:(length(func.names)), function(n) do.call(functional::Curry, c(list(getFromNamespace(func.names[n], 'climdex.pcic')), options.climdex[[n]]))) + + func <- c(func.eca,func.climdex) + + names(func) <- c("csuETCCDI_yr", "csuETCCDI_mon", "cfdETCCDI_yr", "cfdETCCDI_mon", "hd17ETCCDI_yr", "hd17ETCCDI_mon", "hiETCCDI_yr", "spi3ETCCDI_mon", "spi6ETCCDI_mon", "fdETCCDI_yr", "suETCCDI_yr", "idETCCDI_yr","trETCCDI_yr", "gslETCCDI_yr", + "fdETCCDI_mon", "suETCCDI_mon", "idETCCDI_mon","trETCCDI_mon", + + "txxETCCDI_mon", "tnxETCCDI_mon", "txnETCCDI_mon", "tnnETCCDI_mon", "tn10pETCCDI_mon", "tx10pETCCDI_mon", "tn90pETCCDI_mon", "tx90pETCCDI_mon", + "txxETCCDI_yr", "tnxETCCDI_yr", "txnETCCDI_yr", "tnnETCCDI_yr", "tn10pETCCDI_yr", "tx10pETCCDI_yr", "tn90pETCCDI_yr", "tx90pETCCDI_yr", + + "wsdiETCCDI_yr", "csdiETCCDI_yr", + + "dtrETCCDI_mon", "rx1dayETCCDI_mon", "rx5dayETCCDI_mon", + "dtrETCCDI_yr", "rx1dayETCCDI_yr", "rx5dayETCCDI_yr", + + "sdiiETCCDI_yr", "r10mmETCCDI_yr", "r20mmETCCDI_yr", "r1mmETCCDI_yr", "cddETCCDI_yr", "cwdETCCDI_yr", "r95pETCCDI_yr", "r99pETCCDI_yr", "prcptotETCCDI_yr", + "sdiiETCCDI_mon", "r10mmETCCDI_mon", "r20mmETCCDI_mon", "r1mmETCCDI_mon", "r95pETCCDI_mon", "r99pETCCDI_mon", "prcptotETCCDI_mon", + + "altcddETCCDI_yr", "altcwdETCCDI_yr", "altcsdiETCCDI_yr", "altwsdiETCCDI_yr") + + return(func[vars.list]) +} + + +#' Returns metadata for specified Climdex variables +#' Add details for additional indices +get.climdex.variable.metadata.eobs <- function(vars.list, template.filename) { + all.data <- data.frame(long.name=c("Annual Consecutive Summer Days", "Monthly Consecutive Summer Days", + "Annual Consecutive Frost Days", "Monthly Consecutive Frost Days", + "Annual Heating Degree days", "Monthly Heating Degree days", "Huglin Index", + "Standardized Precipitation Index 3mon", + "Standardized Precipitation Index 6mon", + "Annual Number of Frost Days", "Annual Number of Summer Days", "Annual Number of Icing Days", "Annual Number of Tropical Nights", "Growing Season Length", + "Monthly Number of Frost Days", "Monthly Number of Summer Days", "Monthly Number of Icing Days", "Monthly Number of Tropical Nights", + + "Monthly Maximum of Daily Maximum Temperature", "Monthly Maximum of Daily Minimum Temperature", + "Monthly Minimum of Daily Maximum Temperature", "Monthly Minimum of Daily Minimum Temperature", + + "Percentage of Days when Daily Minimum Temperature is Below the 10th Percentile", "Percentage of Days when Daily Maximum Temperature is Below the 10th Percentile", + "Percentage of Days when Daily Minimum Temperature is Above the 90th Percentile", "Percentage of Days when Daily Maximum Temperature is Above the 90th Percentile", + "Annual Maximum of Daily Maximum Temperature", "Annual Maximum of Daily Minimum Temperature", + "Annual Minimum of Daily Maximum Temperature", "Annual Minimum of Daily Minimum Temperature", + "Percentage of Days when Daily Minimum Temperature is Below the 10th Percentile", "Percentage of Days when Daily Maximum Temperature is Below the 10th Percentile", + "Percentage of Days when Daily Minimum Temperature is Above the 90th Percentile", "Percentage of Days when Daily Maximum Temperature is Above the 90th Percentile", + + "Warm Spell Duration Index", "Cold Spell Duration Index", + + "Mean Diurnal Temperature Range", "Monthly Maximum 1-day Precipitation", "Monthly Maximum Consecutive 5-day Precipitation", + "Mean Diurnal Temperature Range", "Annual Maximum 1-day Precipitation", "Annual Maximum Consecutive 5-day Precipitation", + + "Annual Simple Precipitation Intensity Index", "Annual Count of Days with At Least 10mm of Precipitation", + "Annual Count of Days with At Least 20mm of Precipitation", "Annual Count of Days with At Least 1mm of Precipitation", + "Maximum Number of Consecutive Days with Less Than 1mm of Precipitation", "Maximum Number of Consecutive Days with At Least 1mm of Precipitation", + "Annual Total Precipitation when Daily Precipitation Exceeds the 95th Percentile of Wet Day Precipitation", + "Annual Total Precipitation when Daily Precipitation Exceeds the 99th Percentile of Wet Day Precipitation", "Annual Total Precipitation in Wet Days", + + "Monthly Simple Precipitation Intensity Index", "Monthly Count of Days with At Least 10mm of Precipitation", + "Monthly Count of Days with At Least 20mm of Precipitation", "Monthly Count of Days with At Least 1mm of Precipitation", + "Monthly Total Precipitation when Daily Precipitation Exceeds the 95th Percentile of Wet Day Precipitation", + "Monthly Total Precipitation when Daily Precipitation Exceeds the 99th Percentile of Wet Day Precipitation", "Monthly Total Precipitation in Wet Days", + + "Maximum Number of Consecutive Days Per Year with Less Than 1mm of Precipitation", "Maximum Number of Consecutive Days Per Year with At Least 1mm of Precipitation", + "Cold Spell Duration Index Spanning Years", "Warm Spell Duration Index Spanning Years"), + + var.name=c("csuETCCDI", "csuETCCDI", "cfdETCCDI", "cfdETCCDI", "hd17ETCCDI", "hd17ETCCDI", "hiETCCDI", "spi3ETCCDI", "spi6ETCCDI", + + "fdETCCDI", "suETCCDI","idETCCDI", "trETCCDI", "gslETCCDI", + "fdETCCDI", "suETCCDI","idETCCDI", "trETCCDI", + + "txxETCCDI", "tnxETCCDI", "txnETCCDI", "tnnETCCDI", "tn10pETCCDI", "tx10pETCCDI", "tn90pETCCDI", "tx90pETCCDI", + "txxETCCDI", "tnxETCCDI", "txnETCCDI", "tnnETCCDI", "tn10pETCCDI", "tx10pETCCDI", "tn90pETCCDI", "tx90pETCCDI", + + "wsdiETCCDI", "csdiETCCDI", + + "dtrETCCDI", "rx1dayETCCDI", "rx5dayETCCDI", + "dtrETCCDI", "rx1dayETCCDI", "rx5dayETCCDI", + + "sdiiETCCDI", "r10mmETCCDI", "r20mmETCCDI", "r1mmETCCDI", "cddETCCDI", "cwdETCCDI", "r95pETCCDI", "r99pETCCDI", "prcptotETCCDI", + "sdiiETCCDI", "r10mmETCCDI", "r20mmETCCDI", "r1mmETCCDI", "r95pETCCDI", "r99pETCCDI", "prcptotETCCDI", + + "altcddETCCDI", "altcwdETCCDI", "altcsdiETCCDI", "altwsdiETCCDI"), + + units=c("days", "days","days","days", "degrees_C", "degrees_C", "degrees_C", "", "", "days", "days", "days","days", "days", + "days", "days", "days","days", + + "degrees_C", "degrees_C", "degrees_C", "degrees_C", "%", "%", "%", "%", + "degrees_C", "degrees_C", "degrees_C", "degrees_C", "%", "%", "%", "%", + + "days", "days", + + "degrees_C", "mm", "mm", + "degrees_C", "mm", "mm", + + "mm d-1", "days", "days", "days", "days", "days", "mm", "mm", "mm", + "mm d-1", "days", "days", "days", "mm", "mm", "mm", + + "days", "days", "days", "days"), + + annual=c(T, F, T, F, T, F, T, F, F, T, T, T, T, T, + F, F, F, F, + + F, F, F, F, F, F, F, F, + T, T, T, T, T, T, T, T, + + T, T, + + F, F, F, + T, T, T, + + T, T, T, T, T, T, T, T, T, + F, F, F, F, F, F, F, + + T, T, T, T), + + base.period.attr=c(F, F, F, F, F, F, F, F, F, F, F, F, F, F, + F, F, F, F, + + F, F, F, F, T, T, T, T, + F, F, F, F, T, T, T, T, + + T, T, + + F, F, F, + F, F, F, + + F, F, F, F, F, F, T, T, F, + F, F, F, F, T, T, F, + + F, F, T, T), + row.names=c("csuETCCDI_yr", "csuETCCDI_mon", "cfdETCCDI_yr", "cfdETCCDI_mon", "hd17ETCCDI_yr", "hd17ETCCDI_mon", "hiETCCDI_yr", + "spi3ETCCDI_mon", "spi6ETCCDI_mon", + "fdETCCDI_yr", "suETCCDI_yr", "idETCCDI_yr", "trETCCDI_yr", "gslETCCDI_yr", + "fdETCCDI_mon", "suETCCDI_mon", "idETCCDI_mon", "trETCCDI_mon", + + "txxETCCDI_mon", "tnxETCCDI_mon", "txnETCCDI_mon", "tnnETCCDI_mon", "tn10pETCCDI_mon", "tx10pETCCDI_mon", "tn90pETCCDI_mon", "tx90pETCCDI_mon", + "txxETCCDI_yr", "tnxETCCDI_yr", "txnETCCDI_yr", "tnnETCCDI_yr", "tn10pETCCDI_yr", "tx10pETCCDI_yr", "tn90pETCCDI_yr", "tx90pETCCDI_yr", + + "wsdiETCCDI_yr", "csdiETCCDI_yr", + + "dtrETCCDI_mon", "rx1dayETCCDI_mon", "rx5dayETCCDI_mon", + "dtrETCCDI_yr", "rx1dayETCCDI_yr", "rx5dayETCCDI_yr", + + "sdiiETCCDI_yr", "r10mmETCCDI_yr", "r20mmETCCDI_yr", "r1mmETCCDI_yr", "cddETCCDI_yr", "cwdETCCDI_yr", "r95pETCCDI_yr", "r99pETCCDI_yr", "prcptotETCCDI_yr", + "sdiiETCCDI_mon", "r10mmETCCDI_mon", "r20mmETCCDI_mon", "r1mmETCCDI_mon", "r95pETCCDI_mon", "r99pETCCDI_mon", "prcptotETCCDI_mon", + + "altcddETCCDI_yr", "altcwdETCCDI_yr", "altcsdiETCCDI_yr", "altwsdiETCCDI_yr"), + + stringsAsFactors=FALSE) + + standard.name.lookup <- c(csuETCCDI="consecutive summer days", cfdETCCDI="consecutive frost days", hd17ETCCDI="heating degree days", hiETCCDI="huglin index", + spi3ETCCDI="standardized precipitation index",spi6ETCCDI="standardized precipitation index", + fdETCCDI="number_frost_days", suETCCDI="number_summer_days", idETCCDI="number_icing_days", trETCCDI="number_tropical_nights", gslETCCDI="growing_season_length", + txxETCCDI="maximum_daily_maximum_temperature", tnxETCCDI="maximum_daily_minimum_temperature", txnETCCDI="minimum_daily_maximum_temperature", + tnnETCCDI="minimum_daily_minimum_temperature", + tn10pETCCDI="percent_days_when_daily_minimum_temperature_below_10p", tx10pETCCDI="percent_days_when_daily_maximum_temperature_below_10p", + tn90pETCCDI="percent_days_when_daily_minimum_temperature_above_90p", tx90pETCCDI="percent_days_when_daily_maximum_temperature_above_90p", + wsdiETCCDI="warm_spell_duration_index", csdiETCCDI="cold_spell_duration_index", dtrETCCDI="diurnal_temperature_range", + altwsdiETCCDI="warm_spell_duration_index", altcsdiETCCDI="cold_spell_duration_index", + rx1dayETCCDI="maximum_1day_precipitation", rx5dayETCCDI="maximum_5day_precipitation", sdiiETCCDI="simple_precipitation_intensity_index", + r10mmETCCDI="count_days_more_than_10mm_precipitation", r20mmETCCDI="count_days_more_than_20mm_precipitation", r1mmETCCDI="count_days_more_than_1mm_precipitation", + cddETCCDI="maximum_number_consecutive_dry_days", cwdETCCDI="maximum_number_consecutive_wet_days", + altcddETCCDI="maximum_number_consecutive_dry_days", altcwdETCCDI="maximum_number_consecutive_wet_days", + r95pETCCDI="total_precipitation_exceeding_95th_percentile", r99pETCCDI="total_precipitation_exceeding_99th_percentile", prcptotETCCDI="total_wet_day_precipitation") + + all.data$standard.name <- standard.name.lookup[all.data$var.name] + + all.data$filename <- create.climdex.eobs.filenames(get.split.filename.eobs(template.filename), rownames(all.data)) + return(all.data[vars.list,]) +} diff --git a/R/indices_eca.R b/R/indices_eca.R new file mode 100644 index 0000000..1683e21 --- /dev/null +++ b/R/indices_eca.R @@ -0,0 +1,156 @@ +#' STANDARDIZE PRECIPITATION INDEX +#' based on the formula in the spi package +#' SPI is only has monthly outputs (see changes in get.climdex.variable.list.eobs) +#' make monthly sums +monthly_sums_spi <- function(temp,date.factor) { + stopifnot(is.numeric(temp) && is.factor(date.factor)) + return(tapply(temp,date.factor,sum,na.rm=TRUE)) +} +#' calcualte eca monthly sums +eca_sums_mon <- function(ci,freq=c("monthly")){ + stopifnot(!is.null(ci@data$prec)) + return(monthly_sums_spi(ci@data$prec,ci@date.factors$monthly) * ci@namasks$monthly$prec) +} +## Makes multi-month averages depending on k (here k=3,6). For spi3 k=3 +# precipitation is a vector of monthly precipitation sums +# returns monthly precipation averaged over current month and prior k-1 months +getPrecOnTimescale <- function(precipitation,k){ + Nt <- length(precipitation) + prec.k <- as.vector(sapply(seq(from=1, to=Nt),function(t) {tm <- max(t-k+1,1); sum(as.vector(precipitation[tm:t]))})) + return(prec.k) +} + +#'In ECA&D if 85% of data are present then we consider them as valid. +#'This was not used in SPI or elsewhere (yet) +# dataLeadsToNAindex <- function(dat) { +# percentage_should_be_present = 0.85 +# number_of_nan = sum(is.na(dat)) +# total_length = length(dat) +# return((number_of_nan / total_length) < percentage_should_be_present) +# } + +#' empirical Standard Precipitation Index 3 (k in getPrecOnTimescale is 3) +eca.SPI3 <- function(ci,freq=c("monthly")){ + dat_mon <- eca_sums_mon(ci, freq=c("monthly")) + dat <- getPrecOnTimescale(dat_mon,3) + #if(dataLeadsToNAindex(dat)) return(NA) + if(all(is.na(dat))) return(NA) + + fit.cdf <- ecdf(dat) + cdfs <- fit.cdf(dat) + spi.t <- qnorm(cdfs) + spi.sym <- spi.t + # drop Inf + spi.sym[which(spi.t == Inf)] <- NA + spi.sym[which(spi.t == -Inf)] <-NA + return(spi.sym) +} + +#' empirical Standard Precipitation Index 6 (k in getPrecOnTimescale is 6) +eca.SPI6 <- function(ci,freq=c("monthly")){ + dat_mon <- eca_sums_mon(ci, freq=c("monthly")) + dat <- getPrecOnTimescale(dat_mon,6) + #if(dataLeadsToNAindex(dat)) return(NA) + if(all(is.na(dat))) return(NA) + + fit.cdf <- ecdf(dat) + cdfs <- fit.cdf(dat) + spi.t <- qnorm(cdfs) + spi.sym <- spi.t + # drop Inf (in the spi package Inf is set to 2.2 and -Inf to -2.2) + spi.sym[which(spi.t == Inf)] <- NA + spi.sym[which(spi.t == -Inf)] <-NA + return(spi.sym) +} + + +#' CONSECUTIVE SUMMER DAYS +#' maximum number of concsecutive summer days (tx > 25C) +#' spell.can.span.years for us is FALSE +eca.csu <- function(ci,freq=c("monthly","annual"),spells.can.span.years=FALSE) { + stopifnot(!is.null(ci@data$tmax)); + return(climdex.pcic:::spell.length.max(ci@data$tmax, ci@date.factors[[match.arg(freq)]], 25, ">", spells.can.span.years) * ci@namasks[[match.arg(freq)]]$tmax) +} + +#' CONSECUTIVE FROST DAYS +#' maximum number of concsecutive frost days (tn < 0C) +#' spell.can.span.years for us is FALSE +eca.cfd <- function(ci,freq=c("monthly","annual"),spells.can.span.years=FALSE) { + stopifnot(!is.null(ci@data$tmin)); + return(climdex.pcic:::spell.length.max(ci@data$tmin, ci@date.factors[[match.arg(freq)]], 0, "<", spells.can.span.years) * ci@namasks[[match.arg(freq)]]$tmin) } + +#' CONSECUTIVE WET DAYS +#' We would This is a trial to have cwd for months also (not sure about this) +#' spells.can.span.years=FALSE +climdex.cwd.eobs <- function(ci, freq=c("monthly","annual"),spells.can.span.years=FALSE) { + stopifnot(!is.null(ci@data$prec)); + return(climdex.pcic:::spell.length.max(ci@data$prec, ci@date.factors[[match.arg(freq)]], 1, ">=", spells.can.span.years) * ci@namasks[[match.arg(freq)]]$prec) +} + +#' HEATING DEGREE DAYS +#' sum of (17C-TG) (C) +#' +eca.hd17 <- function(ci, freq=c("monthly","annual")) { + stopifnot(!is.null(ci@data$tavg)); + return(tapply((17 - ci@data$tavg), ci@date.factors[[match.arg(freq)]], sum)* ci@namasks[[match.arg(freq)]]$tavg) +} + +###################################################################################################################### +#' HUGLIN INDEX +#' This function is not ready yet. It is uses coeeficient based on the latitude +#' For this index I had to curry the cdx.funcs to be able to include the subset. Later I realised +#' I need also to include the latitude. This would be used in compute.indices.for.stripe together with get.lat +#' to retrieve subset & latitude +#' I didn't proceed with finisheing the eca.HI function. I thought to ask you first if its possible +#' to adapt compute.indices.for.stripe so it can include the currying and the latitude. Or if you had a better idea on this +#' please let me know. + +#' Function to Curry a cxd.funcs for subset +#' used only for Huglin Index +curry_in_subset_for_huglin <- function(cdx.funcs, subset){ + cdx.names = names(cdx.funcs) + cdx.funcs <- lapply(cdx.names, function(function_name) { + if(grepl('^hi', function_name)) { + f = cdx.funcs[[function_name]] + return(functional::Curry(f, subset = subset[['Y']])) + } else { + return(f) + } + }) + names(cdx.funcs) = cdx.names + return(cdx.funcs) +} + +#' Function to retrieve the latitude +#' we have different variable names for temp & precip; that's why I've changed the variable.name.map +get.lat <- function(open_file_list, v.f.idx, variable.name.map=c(tmax="tx", tmin="tn", prec="rr", tavg="tg")) { + var.name <- variable.name.map[[names(v.f.idx)[1]]] + y.dim <- ncdf4.helpers::nc.get.dim.for.axis(open_file_list[[1]], var.name, "Y") + return(y.dim$vals) +} + +#' Huglin index +#' It is annual only and valid for months (April-Sep) +eca.HI <- function(ci,freq=c("annual"),subset=-1){ + + tempavg <- ci@data$tavg + tempmax <- ci@data$tmax + + #browser() + month.series <- climdex.pcic:::get.months(ci@dates) + year.series <- climdex.pcic:::get.years(ci@dates) + valid.months <- month.series >=4 & month.series <=9 + #browser() + hi_coef <- if (subset <=40) {hi_coeff <- 1 + }else if(subset >40 & subset <42) {hi_coef <- 1.02 + }else if(subset >42 & subset <44) {hi_coef <- 1.03 + }else if(subset >44 & subset <48) {hi_coef <- 1.04 + }else if(subset >46 & subset <48) {hi_coef <- 1.05 + }else if(subset >48 & subset <50) {hi_coef <- 1.06 + }else if(subset >=50){hi_coef <- 1} + valid.sel<- year.series[valid.months] + tempdata <- ((((tempavg -10) + (tempmax -10)) /2) * hi_coef) + dat_final <- tempdata[valid.months] + return(tapply(dat_final,valid.sel,sum,na.rm=TRUE)) + +} diff --git a/climdex.pcic.Rproj b/climdex.pcic.Rproj new file mode 100644 index 0000000..21a4da0 --- /dev/null +++ b/climdex.pcic.Rproj @@ -0,0 +1,17 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX + +BuildType: Package +PackageUseDevtools: Yes +PackageInstallArgs: --no-multiarch --with-keep.source diff --git a/tests/test_for_eca_input.r b/tests/test_for_eca_input.r new file mode 100644 index 0000000..aa54caa --- /dev/null +++ b/tests/test_for_eca_input.r @@ -0,0 +1,41 @@ +fname = '~/Downloads/TX_STAID000162.txt' + +find_start_of_data_index = function(fname) { + matched_indices = which(grepl('^SOUID', gsub(" ", "", readLines(fname, n = 50)))) + if (length(matched_indices) > 1) stop('ECA fileformat error: cannot determine start of data, multiple header lines') + if (length(matched_indices) == 0) stop('ECA fileformat error: cannot find start of data: cannot find header line') + return(matched_indices - 1) +} + + +tg <- eca.input('~/Documents/deBilt_tavg/TG_STAID000162.txt', 'TG', 'DATE') +tx <- eca.input('~/Documents/deBilt_tmax/TX_STAID000162.txt', 'TX', 'DATE') +tn <- eca.input('~/Documents/deBilt_tmin/TN_STAID000162.txt', 'TN', 'DATE') +prec <- eca.input('~/Documents/deBilt_prec/RR_STAID000162.txt', 'RR', 'DATE') +fx <- eca.input('~/Documents/deBilt_windgust/FX_STAID000162.txt', 'FX', 'DATE') + +ci <- climdexInput.raw(tmax= tx$TX, tmin=tn$TN, tavg=tg$TG, prec=prec$RR, wind_gust = fx$FX, + tmax.dates = tx$DATE, tmin.dates = tn$DATE, tavg.dates = tg$DATE, prec.dates = prec$DATE, wind_gust.dates = fx$DATE, + base.range=c(1961, 1991)) + + +ci <- climdexInput.raw(tmax=tx$TX, tavg=tg$TG, prec=prec$RR, wind_gust = fx$FX, + tmax.dates = tx$DATE, tavg.dates = tg$DATE, prec.dates = prec$DATE, wind_gust.dates = fx$DATE) + + + +ci <- climdexInput.raw(tavg=tg$TG, prec=prec$RR, wind_gust = fx$FX, + tavg.dates = tg$DATE, prec.dates = prec$DATE, wind_gust.dates = fx$DATE) + + +test_max_fx <- climdex.wind_maxgust(ci, "halfyear") + +test_rx1day <- climdex.rx1day(ci, "halfyear") + +test_cd <- climdex.cd(ci, 'annual') + + +rm(list = ls()) +ci <- climdexInput.raw(prec=prec$RR, prec.dates = prec$DATE, base.range=c(1961,1990)) + + diff --git a/tests/testthat.R b/tests/testthat.R new file mode 100644 index 0000000..3754686 --- /dev/null +++ b/tests/testthat.R @@ -0,0 +1,4 @@ +library(testthat) +library(climdex.pcic) + +test_check("climdex.pcic") diff --git a/tests/testthat/outputTests/cc2_deBilt.rds b/tests/testthat/outputTests/cc2_deBilt.rds new file mode 100644 index 0000000..0f399ad Binary files /dev/null and b/tests/testthat/outputTests/cc2_deBilt.rds differ diff --git a/tests/testthat/outputTests/cc6_deBilt.rds b/tests/testthat/outputTests/cc6_deBilt.rds new file mode 100644 index 0000000..8da3187 Binary files /dev/null and b/tests/testthat/outputTests/cc6_deBilt.rds differ diff --git a/tests/testthat/outputTests/cc_deBilt.rds b/tests/testthat/outputTests/cc_deBilt.rds new file mode 100644 index 0000000..e8a23e0 Binary files /dev/null and b/tests/testthat/outputTests/cc_deBilt.rds differ diff --git a/tests/testthat/outputTests/cd_deBilt.rds b/tests/testthat/outputTests/cd_deBilt.rds new file mode 100644 index 0000000..8ae861c Binary files /dev/null and b/tests/testthat/outputTests/cd_deBilt.rds differ diff --git a/tests/testthat/outputTests/cfd_deBilt.rds b/tests/testthat/outputTests/cfd_deBilt.rds new file mode 100644 index 0000000..9749cea Binary files /dev/null and b/tests/testthat/outputTests/cfd_deBilt.rds differ diff --git a/tests/testthat/outputTests/ci_cloud_deBilt.rds b/tests/testthat/outputTests/ci_cloud_deBilt.rds new file mode 100644 index 0000000..72ea664 Binary files /dev/null and b/tests/testthat/outputTests/ci_cloud_deBilt.rds differ diff --git a/tests/testthat/outputTests/ci_rr_deBilt.rds b/tests/testthat/outputTests/ci_rr_deBilt.rds new file mode 100644 index 0000000..576535b Binary files /dev/null and b/tests/testthat/outputTests/ci_rr_deBilt.rds differ diff --git a/tests/testthat/outputTests/ci_snow_Lugano.rds b/tests/testthat/outputTests/ci_snow_Lugano.rds new file mode 100644 index 0000000..74cb2c1 Binary files /dev/null and b/tests/testthat/outputTests/ci_snow_Lugano.rds differ diff --git a/tests/testthat/outputTests/ci_sun_deBilt.rds b/tests/testthat/outputTests/ci_sun_deBilt.rds new file mode 100644 index 0000000..5e32f0f Binary files /dev/null and b/tests/testthat/outputTests/ci_sun_deBilt.rds differ diff --git a/tests/testthat/outputTests/ci_temp_deBilt.rds b/tests/testthat/outputTests/ci_temp_deBilt.rds new file mode 100644 index 0000000..0970db3 Binary files /dev/null and b/tests/testthat/outputTests/ci_temp_deBilt.rds differ diff --git a/tests/testthat/outputTests/ci_wind_deBilt.rds b/tests/testthat/outputTests/ci_wind_deBilt.rds new file mode 100644 index 0000000..a6dae41 Binary files /dev/null and b/tests/testthat/outputTests/ci_wind_deBilt.rds differ diff --git a/tests/testthat/outputTests/cloud_deBilt.rds b/tests/testthat/outputTests/cloud_deBilt.rds new file mode 100644 index 0000000..b23d8a7 Binary files /dev/null and b/tests/testthat/outputTests/cloud_deBilt.rds differ diff --git a/tests/testthat/outputTests/csu_deBilt.rds b/tests/testthat/outputTests/csu_deBilt.rds new file mode 100644 index 0000000..692cf31 Binary files /dev/null and b/tests/testthat/outputTests/csu_deBilt.rds differ diff --git a/tests/testthat/outputTests/cw_deBilt.rds b/tests/testthat/outputTests/cw_deBilt.rds new file mode 100644 index 0000000..8239604 Binary files /dev/null and b/tests/testthat/outputTests/cw_deBilt.rds differ diff --git a/tests/testthat/outputTests/dd_deBilt.rds b/tests/testthat/outputTests/dd_deBilt.rds new file mode 100644 index 0000000..bb2c1eb Binary files /dev/null and b/tests/testthat/outputTests/dd_deBilt.rds differ diff --git a/tests/testthat/outputTests/ddeast_deBilt.rds b/tests/testthat/outputTests/ddeast_deBilt.rds new file mode 100644 index 0000000..6d27d39 Binary files /dev/null and b/tests/testthat/outputTests/ddeast_deBilt.rds differ diff --git a/tests/testthat/outputTests/ddnorth_deBilt.rds b/tests/testthat/outputTests/ddnorth_deBilt.rds new file mode 100644 index 0000000..e38c96e Binary files /dev/null and b/tests/testthat/outputTests/ddnorth_deBilt.rds differ diff --git a/tests/testthat/outputTests/ddsouth_deBilt.rds b/tests/testthat/outputTests/ddsouth_deBilt.rds new file mode 100644 index 0000000..8b9ccb6 Binary files /dev/null and b/tests/testthat/outputTests/ddsouth_deBilt.rds differ diff --git a/tests/testthat/outputTests/ddwest_deBilt.rds b/tests/testthat/outputTests/ddwest_deBilt.rds new file mode 100644 index 0000000..25a57ed Binary files /dev/null and b/tests/testthat/outputTests/ddwest_deBilt.rds differ diff --git a/tests/testthat/outputTests/fg6bft_deBilt.rds b/tests/testthat/outputTests/fg6bft_deBilt.rds new file mode 100644 index 0000000..4cc5032 Binary files /dev/null and b/tests/testthat/outputTests/fg6bft_deBilt.rds differ diff --git a/tests/testthat/outputTests/fg_ci_deBilt.rds b/tests/testthat/outputTests/fg_ci_deBilt.rds new file mode 100644 index 0000000..d11ab8a Binary files /dev/null and b/tests/testthat/outputTests/fg_ci_deBilt.rds differ diff --git a/tests/testthat/outputTests/fg_deBilt.rds b/tests/testthat/outputTests/fg_deBilt.rds new file mode 100644 index 0000000..40caf25 Binary files /dev/null and b/tests/testthat/outputTests/fg_deBilt.rds differ diff --git a/tests/testthat/outputTests/fgcalm_deBilt.rds b/tests/testthat/outputTests/fgcalm_deBilt.rds new file mode 100644 index 0000000..987d6a7 Binary files /dev/null and b/tests/testthat/outputTests/fgcalm_deBilt.rds differ diff --git a/tests/testthat/outputTests/fx_deBilt.rds b/tests/testthat/outputTests/fx_deBilt.rds new file mode 100644 index 0000000..a8c905c Binary files /dev/null and b/tests/testthat/outputTests/fx_deBilt.rds differ diff --git a/tests/testthat/outputTests/fxstorm_deBilt.rds b/tests/testthat/outputTests/fxstorm_deBilt.rds new file mode 100644 index 0000000..ee96edc Binary files /dev/null and b/tests/testthat/outputTests/fxstorm_deBilt.rds differ diff --git a/tests/testthat/outputTests/fxx_deBilt.rds b/tests/testthat/outputTests/fxx_deBilt.rds new file mode 100644 index 0000000..bb89667 Binary files /dev/null and b/tests/testthat/outputTests/fxx_deBilt.rds differ diff --git a/tests/testthat/outputTests/hd17_deBilt.rds b/tests/testthat/outputTests/hd17_deBilt.rds new file mode 100644 index 0000000..e38c3e6 Binary files /dev/null and b/tests/testthat/outputTests/hd17_deBilt.rds differ diff --git a/tests/testthat/outputTests/r99ptot_deBilt.rds b/tests/testthat/outputTests/r99ptot_deBilt.rds new file mode 100644 index 0000000..04e5490 Binary files /dev/null and b/tests/testthat/outputTests/r99ptot_deBilt.rds differ diff --git a/tests/testthat/outputTests/rr_deBilt.rds b/tests/testthat/outputTests/rr_deBilt.rds new file mode 100644 index 0000000..1ba67ca Binary files /dev/null and b/tests/testthat/outputTests/rr_deBilt.rds differ diff --git a/tests/testthat/outputTests/sd_deBilt.rds b/tests/testthat/outputTests/sd_deBilt.rds new file mode 100644 index 0000000..9e6444a Binary files /dev/null and b/tests/testthat/outputTests/sd_deBilt.rds differ diff --git a/tests/testthat/outputTests/sdd_deBilt.rds b/tests/testthat/outputTests/sdd_deBilt.rds new file mode 100644 index 0000000..79f65f8 Binary files /dev/null and b/tests/testthat/outputTests/sdd_deBilt.rds differ diff --git a/tests/testthat/outputTests/sdx_deBilt.rds b/tests/testthat/outputTests/sdx_deBilt.rds new file mode 100644 index 0000000..4377569 Binary files /dev/null and b/tests/testthat/outputTests/sdx_deBilt.rds differ diff --git a/tests/testthat/outputTests/snowD_deBilt.rds b/tests/testthat/outputTests/snowD_deBilt.rds new file mode 100644 index 0000000..84925a4 Binary files /dev/null and b/tests/testthat/outputTests/snowD_deBilt.rds differ diff --git a/tests/testthat/outputTests/spi3_deBilt.rds b/tests/testthat/outputTests/spi3_deBilt.rds new file mode 100644 index 0000000..31d2a1b Binary files /dev/null and b/tests/testthat/outputTests/spi3_deBilt.rds differ diff --git a/tests/testthat/outputTests/spi6_deBilt.rds b/tests/testthat/outputTests/spi6_deBilt.rds new file mode 100644 index 0000000..fafc8bd Binary files /dev/null and b/tests/testthat/outputTests/spi6_deBilt.rds differ diff --git a/tests/testthat/outputTests/ss_deBilt.rds b/tests/testthat/outputTests/ss_deBilt.rds new file mode 100644 index 0000000..aaaf100 Binary files /dev/null and b/tests/testthat/outputTests/ss_deBilt.rds differ diff --git a/tests/testthat/outputTests/sunshine_deBilt.rds b/tests/testthat/outputTests/sunshine_deBilt.rds new file mode 100644 index 0000000..0e2bcb1 Binary files /dev/null and b/tests/testthat/outputTests/sunshine_deBilt.rds differ diff --git a/tests/testthat/outputTests/tg_deBilt.rds b/tests/testthat/outputTests/tg_deBilt.rds new file mode 100644 index 0000000..30db5c4 Binary files /dev/null and b/tests/testthat/outputTests/tg_deBilt.rds differ diff --git a/tests/testthat/outputTests/tmndaymax_deBilt.rds b/tests/testthat/outputTests/tmndaymax_deBilt.rds new file mode 100644 index 0000000..e20fd95 Binary files /dev/null and b/tests/testthat/outputTests/tmndaymax_deBilt.rds differ diff --git a/tests/testthat/outputTests/tmndaymin_deBilt.rds b/tests/testthat/outputTests/tmndaymin_deBilt.rds new file mode 100644 index 0000000..e30893a Binary files /dev/null and b/tests/testthat/outputTests/tmndaymin_deBilt.rds differ diff --git a/tests/testthat/outputTests/tn_deBilt.rds b/tests/testthat/outputTests/tn_deBilt.rds new file mode 100644 index 0000000..4a8e953 Binary files /dev/null and b/tests/testthat/outputTests/tn_deBilt.rds differ diff --git a/tests/testthat/outputTests/tnndaymax_deBilt.rds b/tests/testthat/outputTests/tnndaymax_deBilt.rds new file mode 100644 index 0000000..0fbf76b Binary files /dev/null and b/tests/testthat/outputTests/tnndaymax_deBilt.rds differ diff --git a/tests/testthat/outputTests/tnndaymin_deBilt.rds b/tests/testthat/outputTests/tnndaymin_deBilt.rds new file mode 100644 index 0000000..3861811 Binary files /dev/null and b/tests/testthat/outputTests/tnndaymin_deBilt.rds differ diff --git a/tests/testthat/outputTests/tx_deBilt.rds b/tests/testthat/outputTests/tx_deBilt.rds new file mode 100644 index 0000000..1606137 Binary files /dev/null and b/tests/testthat/outputTests/tx_deBilt.rds differ diff --git a/tests/testthat/outputTests/txndaymax_deBilt.rds b/tests/testthat/outputTests/txndaymax_deBilt.rds new file mode 100644 index 0000000..04ff2d6 Binary files /dev/null and b/tests/testthat/outputTests/txndaymax_deBilt.rds differ diff --git a/tests/testthat/outputTests/txndaymin_deBilt.rds b/tests/testthat/outputTests/txndaymin_deBilt.rds new file mode 100644 index 0000000..04ff2d6 Binary files /dev/null and b/tests/testthat/outputTests/txndaymin_deBilt.rds differ diff --git a/tests/testthat/outputTests/wd_deBilt.rds b/tests/testthat/outputTests/wd_deBilt.rds new file mode 100644 index 0000000..d7d56b5 Binary files /dev/null and b/tests/testthat/outputTests/wd_deBilt.rds differ diff --git a/tests/testthat/outputTests/ww_deBilt.rds b/tests/testthat/outputTests/ww_deBilt.rds new file mode 100644 index 0000000..8de937a Binary files /dev/null and b/tests/testthat/outputTests/ww_deBilt.rds differ diff --git a/tests/testthat/test_eca_indices.r b/tests/testthat/test_eca_indices.r new file mode 100644 index 0000000..b911247 --- /dev/null +++ b/tests/testthat/test_eca_indices.r @@ -0,0 +1,265 @@ +#' Description +#' This testthat tests the new indices added in climdex.pcic +#' 1. Test the eca.input function which can be used for station files from ECA&D +#' 2. Test the climdexInput.raw for the added parameters, added frequencies (halfyear, seasons), and new quantiles +#' 3. Test new indices for precipitation and temperature +#' 4. test new indices for new parameters: cloud, sun, snow, wind. +#' 5. We don't test yet the sun_rel and the snow_new (related indices) + +rm(list = ls()) +library(climdex.pcic) +library(testthat) + +#################################### +### Precipitation indices +#################################### +context("Precipitation indices") + +test_that("Precipitation indices annual", { + prec <- eca.input('~/Documents/deBilt_prec/RR_STAID000162.txt', 'RR', 'DATE') + + expect_equal_to_reference(prec, "./tests/testthat/outputTests/rr_deBilt.rds") + + ## Here we test the additional frequencies (halfyear & seasons) + ci <- climdexInput.raw(prec=prec$RR, prec.dates = prec$DATE, base.range=c(1961,1990)) + expect_equal_to_reference(ci, "./tests/testthat/outputTests/ci_rr_deBilt.rds") + + ## Test new indices for precipitation + ci_spi3 <- climdex.spi3(ci, freq=c("monthly"), scale=3) + ci_spi6 <- climdex.spi6(ci, freq=c("monthly"), scale=6) + ci_r99p <- climdex.r99p(ci, freq=c("monthly")) + ci_r99ptot <- climdex.r99ptot(ci, freq=c("monthly")) + + expect_equal_to_reference(ci_spi3, "./tests/testthat/outputTests/spi3_deBilt.rds") + expect_equal_to_reference(ci_spi6, "./tests/testthat/outputTests/spi6_deBilt.rds") + expect_equal_to_reference(ci_r99p, "./tests/testthat/outputTests/r99p_deBilt.rds") + expect_equal_to_reference(ci_r99ptot, "./tests/testthat/outputTests/r99ptot_deBilt.rds") + + +}) + + +#################################### +### Temperature indices +#################################### + +context("tmax indices") + +test_that("tmax indices annual", { + tg <- eca.input('~/Documents/deBilt_tavg/TG_STAID000162.txt', 'TG', 'DATE') + tx <- eca.input('~/Documents/deBilt_tmax/TX_STAID000162.txt', 'TX', 'DATE') + tn <- eca.input('~/Documents/deBilt_tmin/TN_STAID000162.txt', 'TN', 'DATE') + prec <- eca.input('~/Documents/deBilt_prec/RR_STAID000162.txt', 'RR', 'DATE') + + expect_equal_to_reference(tg, "./tests/testthat/outputTests/tg_deBilt.rds") + expect_equal_to_reference(tx, "./tests/testthat/outputTests/tx_deBilt.rds") + expect_equal_to_reference(tn, "./tests/testthat/outputTests/tn_deBilt.rds") + expect_equal_to_reference(prec, "./tests/testthat/outputTests/rr_deBilt.rds") + + + ## Here we test the additional frequencies (halfyear & seasons) and the additional quantiles for temperature (q25, q75 ) + ci_temp <- climdexInput.raw(tmax= tx$TX, tmin=tn$TN, tavg=tg$TG, tmax.dates = tx$DATE, tmin.dates = tn$DATE, + prec=prec$RR, prec.dates = prec$DATE, + tavg.dates = tg$DATE,base.range=c(1961, 1991)) + + expect_equal_to_reference(ci_temp, "./tests/testthat/outputTests/ci_temp_deBilt.rds") + + ## tmax + ci_csu <- climdex.csu(ci_temp, freq=c("annual")) + ci_txndaymin <- climdex.txndaymin(ci_temp, freq=c("annual")) + ci_txndaymax <- climdex.txndaymin(ci_temp, freq=c("annual")) + + expect_equal_to_reference(ci_csu, "./tests/testthat/outputTests/csu_deBilt.rds") + expect_equal_to_reference(ci_txndaymin, "./tests/testthat/outputTests/txndaymin_deBilt.rds") + expect_equal_to_reference(ci_txndaymax, "./tests/testthat/outputTests/txndaymax_deBilt.rds") +}) + +context("tmin indices") + +test_that("tmin indices annual & monthly", { + tn <- eca.input('~/Documents/deBilt_tmin/TN_STAID000162.txt', 'TN', 'DATE') + + ci_temp <- climdexInput.raw(tmin=tn$TN, tmin.dates = tn$DATE,base.range=c(1961, 1991)) + + ci_cfd <- climdex.cfd(ci_temp, freq=c("monthly")) + ci_tnndaymin <- climdex.tnndaymin(ci_temp, freq=c("monthly")) + ci_tnndaymax <- climdex.tnndaymax(ci_temp, freq=c("annual")) + + expect_equal_to_reference(ci_cfd, "./tests/testthat/outputTests/cfd_deBilt.rds") + expect_equal_to_reference(ci_tnndaymin, "./tests/testthat/outputTests/tnndaymin_deBilt.rds") + expect_equal_to_reference(ci_tnndaymax, "./tests/testthat/outputTests/tnndaymax_deBilt.rds") +}) + +context("tavg indices") + +test_that("tavg indices annual & monthly", { + tn <- eca.input('~/Documents/deBilt_tmin/TN_STAID000162.txt', 'TN', 'DATE') + tg <- eca.input('~/Documents/deBilt_tavg/TG_STAID000162.txt', 'TG', 'DATE') + prec <- eca.input('~/Documents/deBilt_prec/RR_STAID000162.txt', 'RR', 'DATE') + + ci_temp <- climdexInput.raw(tmin = tn$TN, tavg = tg$TG, prec = prec$RR, + tmin.dates = tn$DATE, prec.dates = prec$DATE, tavg.dates = tg$DATE, + base.range=c(1961, 1991)) + + ##tavg + ci_hd17 <- climdex.hd17(ci_temp, freq=c("annual")) + ci_tmndaymin <- climdex.tmndaymin(ci_temp, freq=c("monthly")) + ci_tmndaymax <- climdex.tmndaymax(ci_temp, freq=c("monthly")) + ci_cd <- climdex.cd(ci_temp, freq=c("annual")) + ci_cw <- climdex.cw(ci_temp, freq=c("monthly")) + ci_wd <- climdex.wd(ci_temp, freq=c("monthly")) + ci_ww <- climdex.ww(ci_temp, freq=c("monthly")) + + expect_equal_to_reference(ci_hd17, "./tests/testthat/outputTests/hd17_deBilt.rds") + expect_equal_to_reference(ci_tmndaymin, "./tests/testthat/outputTests/tmndaymin_deBilt.rds") + expect_equal_to_reference(ci_tmndaymax, "./tests/testthat/outputTests/tmndaymax_deBilt.rds") + expect_equal_to_reference(ci_cd, "./tests/testthat/outputTests/cd_deBilt.rds") + expect_equal_to_reference(ci_cw, "./tests/testthat/outputTests/cw_deBilt.rds") + expect_equal_to_reference(ci_wd, "./tests/testthat/outputTests/wd_deBilt.rds") + expect_equal_to_reference(ci_ww, "./tests/testthat/outputTests/ww_deBilt.rds") + +}) + +#################################### +### Wind indices +#################################### + +context("Wind speed") + +test_that("Wind speed annual & monthly", { + fg <- eca.input('~/Documents/deBilt_windspeed/FG_STAID000162.txt', 'FG', 'DATE') + fx <- eca.input('~/Documents/deBilt_windgust/FX_STAID000162.txt', 'FX', 'DATE') + dd <- eca.input('~/Documents/deBilt_winddir/DD_STAID000162.txt', 'DD', 'DATE') + + expect_equal_to_reference(fg, "./tests/testthat/outputTests/fg_deBilt.rds") + expect_equal_to_reference(fx, "./tests/testthat/outputTests/fx_deBilt.rds") + expect_equal_to_reference(dd, "./tests/testthat/outputTests/dd_deBilt.rds") + + + ## Here we test the additional frequencies (halfyear & seasons) and the additional quantiles for temperature (q25, q75 ) + ci_wind <- climdexInput.raw(wind= fg$FG, wind_gust=fx$FX, wind_dir=dd$DD, wind.dates = fg$DATE, wind_gust.dates = fx$DATE, + wind_dir.dates = dd$DATE, base.range=c(1961, 1991)) + + expect_equal_to_reference(ci_wind, "./tests/testthat/outputTests/ci_wind_deBilt.rds") + + ## wind speed + ci_fg <- climdex.fg(ci_wind, freq=c("monthly")) + ci_fgcalm <- climdex.fgcalm(ci_wind, freq=c("monthly")) + ci_fg6bft <- climdex.fg6bft(ci_wind, freq=c("annual")) + + expect_equal_to_reference(ci_fg, "./tests/testthat/outputTests/fg_ci_deBilt.rds") + expect_equal_to_reference(ci_fgcalm, "./tests/testthat/outputTests/fgcalm_deBilt.rds") + expect_equal_to_reference(ci_fg6bft, "./tests/testthat/outputTests/fg6bft_deBilt.rds") + +}) +# +# context("Wind gust") +# +# test_that("Wind gust monthly", { +# fx <- eca.input('~/Documents/deBilt_windgust/FX_STAID000162.txt', 'FX', 'DATE') +# +# ci_wind <- climdexInput.raw(wind_gust=fx$FX, wind_gust.dates = fx$DATE, base.range=c(1961, 1991)) +# +# ## wind gust +# ci_fxstorm <- climdex.fxstorm(ci_wind, freq=c("monthly")) +# ci_fxx <- climdex.fxx(ci_wind, freq=c("monthly")) +# +# expect_equal_to_reference(ci_fxstorm, "./tests/testthat/outputTests/fxstorm_deBilt.rds") +# expect_equal_to_reference(ci_fxx, "./tests/testthat/outputTests/fxx_deBilt.rds") +# }) +# +# context("Wind direction") +# +# test_that("Wind direction annual and monthly", { +# dd <- eca.input('~/Documents/deBilt_winddir/DD_STAID000162.txt', 'DD', 'DATE') +# +# ci_wind <- climdexInput.raw(wind_dir=dd$DD, wind_dir.dates = dd$DATE, base.range=c(1961, 1991)) +# +# ## wind direction +# ci_ddnorth <- climdex.ddnorth(ci_wind, freq=c("monthly")) +# ci_ddeast <- climdex.ddeast(ci_wind, freq=c("monthly")) +# ci_ddsouth <- climdex.ddsouth(ci_wind, freq=c("annual")) +# ci_ddwest <- climdex.ddwest(ci_wind, freq=c("annual")) +# +# expect_equal_to_reference(ci_ddnorth, "./tests/testthat/outputTests/ddnorth_deBilt.rds") +# expect_equal_to_reference(ci_ddeast, "./tests/testthat/outputTests/ddeast_deBilt.rds") +# expect_equal_to_reference(ci_ddsouth, "./tests/testthat/outputTests/ddsouth_deBilt.rds") +# expect_equal_to_reference(ci_ddwest, "./tests/testthat/outputTests/ddwest_deBilt.rds") +# +# }) + +#################################### +### Snow indices +#################################### +# +# context("Snow depth") +# +# test_that("Snow annual & monthly", { +# snowD <- eca.input('~/Documents/Lugano_snowdepth/SD_STAID000242.txt', 'SD', 'DATE') +# +# expect_equal_to_reference(snowD, "./tests/testthat/outputTests/snowD_deBilt.rds") +# +# ci_snow <- climdexInput.raw(snow = snowD$SD, snow.dates = snowD$DATE, base.range=c(1961, 1991)) +# +# expect_equal_to_reference(ci_snow, "./tests/testthat/outputTests/ci_snow_Lugano.rds") +# +# ## snow depth +# ci_sdd <- climdex.sdd(ci_snow, freq=c("monthly")) +# ci_sdx <- climdex.sdx(ci_snow, freq=c("monthly")) +# ci_sd <- climdex.sd(ci_snow, freq=c("annual")) +# +# expect_equal_to_reference(ci_sdd, "./tests/testthat/outputTests/sdd_deBilt.rds") +# expect_equal_to_reference(ci_sdx, "./tests/testthat/outputTests/sdx_deBilt.rds") +# expect_equal_to_reference(ci_sd, "./tests/testthat/outputTests/sd_deBilt.rds") +# +# }) + +#################################### +### Cloud indices +#################################### + +context("Cloud indices") + +test_that("Cloud indices annual & monthly", { + cloud <- eca.input('~/Documents/deBilt_cloud/CC_STAID000162.txt', 'CC', 'DATE') + + expect_equal_to_reference(cloud, "./tests/testthat/outputTests/cloud_deBilt.rds") + + ci_cloud <- climdexInput.raw(cloud = cloud$CC, cloud.dates = cloud$DATE, base.range=c(1961, 1991)) + + expect_equal_to_reference(ci_cloud, "./tests/testthat/outputTests/ci_cloud_deBilt.rds") + + ## snow depth + ci_cc <- climdex.cc(ci_cloud, freq=c("monthly")) + ci_cc6 <- climdex.cc6(ci_cloud, freq=c("monthly")) + ci_cc2 <- climdex.cc2(ci_cloud, freq=c("annual")) + + expect_equal_to_reference(ci_cc, "./tests/testthat/outputTests/cc_deBilt.rds") + expect_equal_to_reference(ci_cc6, "./tests/testthat/outputTests/cc6_deBilt.rds") + expect_equal_to_reference(ci_cc2, "./tests/testthat/outputTests/cc2_deBilt.rds") + +}) + +#################################### +### Sun indices +#################################### + +context("Sun indices") + +test_that("Sun index monthly", { + sunshine <- eca.input('~/Documents/deBilt_sunshine/SS_STAID000162.txt', 'SS', 'DATE') + + expect_equal_to_reference(sunshine, "./tests/testthat/outputTests/sunshine_deBilt.rds") + + ci_sun <- climdexInput.raw(sun = sunshine$SS, sun.dates = sunshine$DATE, base.range=c(1961, 1991)) + + expect_equal_to_reference(ci_sun, "./tests/testthat/outputTests/ci_sun_deBilt.rds") + + ## snow depth + ci_ss <- climdex.ss(ci_sun, freq=c("monthly")) + + expect_equal_to_reference(ci_ss, "./tests/testthat/outputTests/ss_deBilt.rds") + +}) + +