Ecadevclimdex

.Rbuildignore

@@ -5,3 +5,5 @@ man-roxygen
 comparison
 README.rst
 CONTRIBUTING.rst
+^.*\.Rproj$
+^\.Rproj\.user$

.gitignore

@@ -0,0 +1,6 @@
+.Rproj.user
+.Rhistory
+.RData
+src/*.o
+src/*.so
+src/*.dll

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

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)

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,[email protected]$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, [email protected][[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, [email protected][[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, [email protected][[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), [email protected][[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))
+
+}

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

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))
+
+

tests/testthat.R

@@ -0,0 +1,4 @@
+library(testthat)
+library(climdex.pcic)
+
+test_check("climdex.pcic")