From 26f2f04dbc3e90c23ec5020167c9c6244cf0dd89 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Wed, 9 Oct 2019 17:08:36 +0200
Subject: [PATCH 01/31] corrected varname to dates.tmax

---
 R/DTRRef.R | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/R/DTRRef.R b/R/DTRRef.R
index 008a993..f275ddc 100644
--- a/R/DTRRef.R
+++ b/R/DTRRef.R
@@ -99,12 +99,12 @@ DTRRef <- function(tmax, tmin, by.seasons = TRUE, dates = NULL, timedim = NULL,
     dim_names <- names(dim(tmax))
   }
   if (is.null(dates)) {
-    dates.max <- attr(tmax, 'Variables')$common$time
-    if (is.null(dates.max)) {
+    dates.tmax <- attr(tmax, 'Variables')$common$time
+    if (is.null(dates.tmax)) {
       dates.tmax <- attr(tmax, 'Variables')$dat1$time
     }
-    dates.min <- attr(tmin, 'Variables')$common$time
-    if (is.null(dates.min)) {
+    dates.tmin <- attr(tmin, 'Variables')$common$time
+    if (is.null(dates.tmin)) {
         dates.tmin <- attr(tmin, 'Variables')$dat1$time
     }
     if (length(dates.tmax) != length(dates.tmin)) {
-- 
GitLab


From c9e150388a9cedf9cfacf7c74f8860f7b9adfdef Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Fri, 22 Nov 2019 16:53:05 +0100
Subject: [PATCH 02/31] accept positive values for west longitudes

---
 R/WeightedMean.R | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/R/WeightedMean.R b/R/WeightedMean.R
index 70c33b4..2ab1a1a 100644
--- a/R/WeightedMean.R
+++ b/R/WeightedMean.R
@@ -120,7 +120,8 @@ WeightedMean <- function(data, lon, lat, region = NULL, mask = NULL, londim = NU
   cosphi <- t(array(cos(lat * pi / 180), dim = c(length(lat), length(lon))))
   nblat <- length(lat)
   nblon <- length(lon)
-  dlon <- abs(c(lon[2 : nblon] - lon[1 : nblon - 1])) * pi / 180
+  lon[lon > 180] = lon[lon > 180] - 360
+  dlon <- abs(c(abs(lon[2 : nblon]) - abs(lon[1 : nblon - 1]))) * pi / 180
   dlon <- c(dlon, dlon[1])
   dlon <- array(dlon, dim = c(nblon, nblat))
   dlat <- abs(c(lat[2 : nblat] - lat[1 : nblat - 1])) * pi / 180
-- 
GitLab


From 21744b259b253d061edce2053ed406c44b563fbb Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Mon, 2 Dec 2019 16:24:00 +0100
Subject: [PATCH 03/31] New function to interprete lons and regions

---
 R/Lon2Index.R | 134 ++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 134 insertions(+)
 create mode 100644 R/Lon2Index.R

diff --git a/R/Lon2Index.R b/R/Lon2Index.R
new file mode 100644
index 0000000..422ffce
--- /dev/null
+++ b/R/Lon2Index.R
@@ -0,0 +1,134 @@
+# Obtain the index of positions for a region in longitudes
+#
+#'@description This auxiliary function returns the index of position of a region of longitudes in a given vector of longitudes.
+#'
+#'@param lon vector of longitudes values.
+#'@param lonmin a numeric value indicating the minimum longitude of the region (understand as the left marging of the region).
+#'@param lonmax a numeric value indicating the maximum longitude of the region (understand as the right mariging of the region).
+#'
+#'@return the index of positions of all values inside the region in the vector lon.
+#'
+#'@examples
+#'
+#'lon <- 1 : 360
+#'pos <- Lon2Index(lon, lonmin = -20, lonmax = 20)
+#'lon[pos]
+#'pos <- Lon2Index(lon, lonmin = 340, lonmax = 20)
+#'lon[pos]
+#'pos <- Lon2Index(lon, lonmin = 20, lonmax = 340)
+#'lon[pos]
+#'pos <- Lon2Index(lon, lonmin = -40, lonmax = -20)
+#'lon[pos]
+#'pos <- Lon2Index(lon, lonmin = 320, lonmax = 340)
+#'lon[pos]
+#'pos <- Lon2Index(lon, lonmin = -220, lonmax = -170)
+#'lon[pos] # fails however it works
+#'pos <- Lon2Index(lon, lonmin = -350, lonmax = -300)
+#'lon[pos] # Fails no, it works
+#'pos <- Lon2Index(lon, lonmin = -400, lonmax = -370)
+#'lon[pos] # Fails
+#'pos <- Lon2Index(lon, lonmin = 340, lonmax = 380)
+#'lon[pos] #Fails however it works for now.
+#'lon <- -180 : 180
+#'pos <- Lon2Index(lon, lonmin = -20, lonmax = 20)
+#'lon[pos]
+#'pos <- Lon2Index(lon, lonmin = 340, lonmax = 20)
+#'lon[pos]
+#'pos <- Lon2Index(lon, lonmin = 20, lonmax = 340)
+#'lon[pos]
+#'pos <- Lon2Index(lon, lonmin = -40, lonmax = -20)
+#'lon[pos]
+#'pos <- Lon2Index(lon, lonmin = 320, lonmax = 340)
+#'lon[pos]
+#'pos <- Lon2Index(lon, lonmin = -220, lonmax = -170)
+#'lon[pos] # fails
+#'pos <- Lon2Index(lon, lonmin = -350, lonmax = -300)
+#'lon[pos] # Fails
+#'pos <- Lon2Index(lon, lonmin = -400, lonmax = -370)
+#'lon[pos] # Fails
+#'pos <- Lon2Index(lon, lonmin = 340, lonmax = 380)
+#'lon[pos] #Fails however it works for now.
+#'lon <- -360 : 360
+#'pos <- Lon2Index(lon, lonmin = -20, lonmax = 20)
+#'lon[pos]
+#'pos <- Lon2Index(lon, lonmin = 340, lonmax = 20)
+#'lon[pos] #Fails
+#'pos <- Lon2Index(lon, lonmin = 20, lonmax = 340)
+#'lon[pos]
+#'pos <- Lon2Index(lon, lonmin = -40, lonmax = -20)
+#'lon[pos]
+#'pos <- Lon2Index(lon, lonmin = 320, lonmax = 340)
+#'lon[pos]
+#'pos <- Lon2Index(lon, lonmin = -220, lonmax = -170)
+#'lon[pos] 
+#'pos <- Lon2Index(lon, lonmin = -350, lonmax = -300)
+#'lon[pos] 
+#'pos <- Lon2Index(lon, lonmin = -400, lonmax = -370)
+#'lon[pos] # Fails
+#'pos <- Lon2Index(lon, lonmin = 340, lonmax = 380)
+#'lon[pos] #Fails again it works
+#'
+#'@export
+Lon2Index <- function(lon, lonmin, lonmax) {
+    if (is.null(lon)) {
+        stop("Parameter 'lon' cannot be NULL.")
+    }
+    if (!is.numeric(lon)) {
+        stop("Parameter 'lon' must be numeric.")
+    }
+    if (!is.vector(lon)) {
+        stop("Parameter 'lon' must be a vector.")
+    }
+    if (!is.numeric(lonmin) | !is.numeric(lonmax)) {
+        stop("Parameter 'lonmin' and 'lonmax' must be numeric.")
+    }
+    if (!is.vector(lonmin) | !is.vector(lonmax)) {
+        stop("Parameter 'lonmin'and 'lonmax' must be a vector.")
+    }
+
+    vlonmax <- max(lon)
+    vlonmin <- min(lon)
+    if (vlonmin < 0 & !(vlonmax > 180)) { # -180 to 180
+        if (lonmin < -180) {
+            stop("Change parameter 'lonmin' to match longitudes ",
+                 "in the range -180 to 180.")
+        } else if (lonmin > 180) {
+            lonmin <- lonmin - 360
+        } 
+        if (lonmax < -180) {
+            stop("Change parameter 'lonmax' to match longitudes ",
+                 "in the range -180 to 180.")
+        } else if (lonmax > 180) {
+            lonmax <- lonmax - 360
+        }
+        if (lonmin > lonmax) {
+            index <- which(lon >= lonmin | lon <= lonmax)
+        } else {
+            index <- which(lon >= lonmin & lon <= lonmax)
+        }
+    } else if (vlonmin < 0 & vlonmax > 180) { # -360 to 360
+        if (lonmin > lonmax) {
+            index <- which(lon >= lonmin | lon <= lonmax)
+        } else {
+            index <- which(lon >= lonmin & lon <= lonmax)
+        }
+    } else { # 0 : 360
+        if (lonmin < 0) {
+            lonmin <- lonmin + 360
+        } else if (lonmin > 360) {
+            lonmin <- lonmin - 360
+        } 
+        if (lonmax < 0) {
+            lonmax <- lonmax + 360
+        } else if (lonmax > 360) {
+            lonmax <- lonmax - 360
+        }
+        if (lonmin > lonmax) {
+            index <- c(which(lon >= lonmin), which(lon <= lonmax))
+        } else {
+            index <- which(lon >= lonmin & lon <= lonmax)
+        }
+    }
+
+return(index)
+}
-- 
GitLab


From 0330e0cf42b3775e38733dae9c01d1f99b32ad62 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Wed, 11 Dec 2019 19:02:58 +0100
Subject: [PATCH 04/31] Reordering output

---
 NAMESPACE        |  1 +
 R/Lon2Index.R    |  4 +--
 man/Lon2Index.Rd | 82 ++++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 85 insertions(+), 2 deletions(-)
 create mode 100644 man/Lon2Index.Rd

diff --git a/NAMESPACE b/NAMESPACE
index 8c0c72b..8364b61 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -7,6 +7,7 @@ export(DTRIndicator)
 export(DTRRef)
 export(DailyAno)
 export(Extremes)
+export(Lon2Index)
 export(SeasonSelect)
 export(SelBox)
 export(Subset)
diff --git a/R/Lon2Index.R b/R/Lon2Index.R
index 422ffce..5b9ff88 100644
--- a/R/Lon2Index.R
+++ b/R/Lon2Index.R
@@ -102,13 +102,13 @@ Lon2Index <- function(lon, lonmin, lonmax) {
             lonmax <- lonmax - 360
         }
         if (lonmin > lonmax) {
-            index <- which(lon >= lonmin | lon <= lonmax)
+            index <- c(which(lon >= lonmin), which(lon <= lonmax))  
         } else {
             index <- which(lon >= lonmin & lon <= lonmax)
         }
     } else if (vlonmin < 0 & vlonmax > 180) { # -360 to 360
         if (lonmin > lonmax) {
-            index <- which(lon >= lonmin | lon <= lonmax)
+            index <- c(which(lon >= lonmin), which(lon <= lonmax))
         } else {
             index <- which(lon >= lonmin & lon <= lonmax)
         }
diff --git a/man/Lon2Index.Rd b/man/Lon2Index.Rd
new file mode 100644
index 0000000..4fecf2b
--- /dev/null
+++ b/man/Lon2Index.Rd
@@ -0,0 +1,82 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Lon2Index.R
+\name{Lon2Index}
+\alias{Lon2Index}
+\usage{
+Lon2Index(lon, lonmin, lonmax)
+}
+\arguments{
+\item{lon}{vector of longitudes values.}
+
+\item{lonmin}{a numeric value indicating the minimum longitude of the region (understand as the left marging of the region).}
+
+\item{lonmax}{a numeric value indicating the maximum longitude of the region (understand as the right mariging of the region).}
+}
+\value{
+the index of positions of all values inside the region in the vector lon.
+}
+\description{
+This auxiliary function returns the index of position of a region of longitudes in a given vector of longitudes.
+}
+\examples{
+
+lon <- 1 : 360
+pos <- Lon2Index(lon, lonmin = -20, lonmax = 20)
+lon[pos]
+pos <- Lon2Index(lon, lonmin = 340, lonmax = 20)
+lon[pos]
+pos <- Lon2Index(lon, lonmin = 20, lonmax = 340)
+lon[pos]
+pos <- Lon2Index(lon, lonmin = -40, lonmax = -20)
+lon[pos]
+pos <- Lon2Index(lon, lonmin = 320, lonmax = 340)
+lon[pos]
+pos <- Lon2Index(lon, lonmin = -220, lonmax = -170)
+lon[pos] # fails however it works
+pos <- Lon2Index(lon, lonmin = -350, lonmax = -300)
+lon[pos] # Fails no, it works
+pos <- Lon2Index(lon, lonmin = -400, lonmax = -370)
+lon[pos] # Fails
+pos <- Lon2Index(lon, lonmin = 340, lonmax = 380)
+lon[pos] #Fails however it works for now.
+lon <- -180 : 180
+pos <- Lon2Index(lon, lonmin = -20, lonmax = 20)
+lon[pos]
+pos <- Lon2Index(lon, lonmin = 340, lonmax = 20)
+lon[pos]
+pos <- Lon2Index(lon, lonmin = 20, lonmax = 340)
+lon[pos]
+pos <- Lon2Index(lon, lonmin = -40, lonmax = -20)
+lon[pos]
+pos <- Lon2Index(lon, lonmin = 320, lonmax = 340)
+lon[pos]
+pos <- Lon2Index(lon, lonmin = -220, lonmax = -170)
+lon[pos] # fails
+pos <- Lon2Index(lon, lonmin = -350, lonmax = -300)
+lon[pos] # Fails
+pos <- Lon2Index(lon, lonmin = -400, lonmax = -370)
+lon[pos] # Fails
+pos <- Lon2Index(lon, lonmin = 340, lonmax = 380)
+lon[pos] #Fails however it works for now.
+lon <- -360 : 360
+pos <- Lon2Index(lon, lonmin = -20, lonmax = 20)
+lon[pos]
+pos <- Lon2Index(lon, lonmin = 340, lonmax = 20)
+lon[pos] #Fails
+pos <- Lon2Index(lon, lonmin = 20, lonmax = 340)
+lon[pos]
+pos <- Lon2Index(lon, lonmin = -40, lonmax = -20)
+lon[pos]
+pos <- Lon2Index(lon, lonmin = 320, lonmax = 340)
+lon[pos]
+pos <- Lon2Index(lon, lonmin = -220, lonmax = -170)
+lon[pos] 
+pos <- Lon2Index(lon, lonmin = -350, lonmax = -300)
+lon[pos] 
+pos <- Lon2Index(lon, lonmin = -400, lonmax = -370)
+lon[pos] # Fails
+pos <- Lon2Index(lon, lonmin = 340, lonmax = 380)
+lon[pos] #Fails again it works
+
+}
+
-- 
GitLab


From 5a5973f613e90a060950f27dad07223c7feafe9a Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Thu, 12 Dec 2019 11:12:53 +0100
Subject: [PATCH 05/31] reorder and tests

---
 tests/testthat.R                |  5 ++
 tests/testthat/test-Lon2index.R | 98 +++++++++++++++++++++++++++++++++
 2 files changed, 103 insertions(+)
 create mode 100644 tests/testthat.R
 create mode 100644 tests/testthat/test-Lon2index.R

diff --git a/tests/testthat.R b/tests/testthat.R
new file mode 100644
index 0000000..6b1f95f
--- /dev/null
+++ b/tests/testthat.R
@@ -0,0 +1,5 @@
+library(testthat)
+library(ClimProjDiags)
+
+test_check("ClimProjDiags")
+
diff --git a/tests/testthat/test-Lon2index.R b/tests/testthat/test-Lon2index.R
new file mode 100644
index 0000000..8f53079
--- /dev/null
+++ b/tests/testthat/test-Lon2index.R
@@ -0,0 +1,98 @@
+context("Generic tests")
+test_that("Sanity checks", {
+  expect_error(Lon2Index(lon = NULL), "Parameter 'lon' cannot be NULL.")
+  expect_error(Lon2Index(lon = 'a'), "Parameter 'lon' must be numeric.")
+  expect_error(Lon2Index(lon = 1),
+               'argument "lonmin" is missing, with no default')
+  expect_error(Lon2Index(lon = 1, lonmin = 'a'),
+               'argument "lonmax" is missing, with no default')
+  expect_error(Lon2Index(lon = 1, lonmin = 'a', lonmax = 1),
+               "Parameter 'lonmin' and 'lonmax' must be numeric.")
+  expect_error(Lon2Index(lon = 1, lonmin = 1, lonmax = 'b'),
+               "Parameter 'lonmin' and 'lonmax' must be numeric.")
+  expect_equal(Lon2Index(lon = 1, lonmin = 1, lonmax = 1), 1)
+})
+
+test_that("Case 1 to 360", {
+lon <- 1 : 360
+  expect_equal(Lon2Index(lon, lonmin = -20, lonmax = 20),
+               c(340 : 360, 1 : 20))
+  expect_equal(Lon2Index(lon, lonmin = 340, lonmax = 20),
+               c(340 : 360, 1 : 20))
+  expect_equal(Lon2Index(lon, lonmin = 20, lonmax = 340),
+               20 : 340)
+  expect_equal(Lon2Index(lon, lonmin = -40, lonmax = -20),
+               320 : 340)
+  expect_equal(Lon2Index(lon, lonmin = 320, lonmax = 340),
+               320 : 340)
+  expect_equal(Lon2Index(lon, lonmin = -220, lonmax = -170),
+               140 : 190)
+  expect_equal(Lon2Index(lon, lonmin = -350, lonmax = -300),
+               10 : 60)
+  expect_equal(Lon2Index(lon, lonmin = -400, lonmax = -370), integer(0))
+  expect_equal(Lon2Index(lon, lonmin = 340, lonmax = 380),
+               c(340 : 360, 1 : 20))
+})
+
+test_that("Case -180 to 180", {
+lon <- -180 : 180
+  expect_equal(Lon2Index(lon, lonmin = -20, lonmax = 20),
+               161 : 201)
+  expect_equal(Lon2Index(lon, lonmin = 340, lonmax = 20),
+               161 : 201)
+  expect_equal(Lon2Index(lon, lonmin = 20, lonmax = 340),
+               c(201 : 361, 1 : 161))
+  expect_equal(Lon2Index(lon, lonmin = -40, lonmax = -20),
+               141 : 161)
+  expect_equal(Lon2Index(lon, lonmin = 320, lonmax = 340),
+               141 : 161)
+  expect_error(Lon2Index(lon, lonmin = -220, lonmax = -170),
+        "Change parameter 'lonmin' to match longitudes in the range -180 to 180.")
+  expect_error(Lon2Index(lon, lonmin = -350, lonmax = -300),
+        "Change parameter 'lonmin' to match longitudes in the range -180 to 180.")
+  expect_error(Lon2Index(lon, lonmin = -400, lonmax = -370),
+        "Change parameter 'lonmin' to match longitudes in the range -180 to 180.")
+  expect_equal(Lon2Index(lon, lonmin = 340, lonmax = 380),
+               161 : 201)
+})
+
+test_that("Case -360 to 360", {
+lon <- -360 : 360
+  expect_equal(Lon2Index(lon, lonmin = -20, lonmax = 20),
+               341 : 381)
+  expect_equal(Lon2Index(lon, lonmin = 340, lonmax = 20),
+               c(701 : 721, 1 : 381))
+  expect_equal(Lon2Index(lon, lonmin = 20, lonmax = 340),
+               381 : 701)
+  expect_equal(Lon2Index(lon, lonmin = -40, lonmax = -20),
+               321 : 341)
+  expect_equal(Lon2Index(lon, lonmin = 320, lonmax = 340),
+               681 : 701)
+  expect_equal(Lon2Index(lon, lonmin = -220, lonmax = -170),
+               141 : 191)
+  expect_equal(Lon2Index(lon, lonmin = -350, lonmax = -300),
+               11 : 61)
+  expect_equal(Lon2Index(lon, lonmin = -400, lonmax = -370), integer(0))
+  expect_equal(Lon2Index(lon, lonmin = 340, lonmax = 380),
+               701 : 721)
+})
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+  
-- 
GitLab


From 903f77d4f7f6e463117e20e9ba997969bb92fd98 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Thu, 12 Dec 2019 11:13:41 +0100
Subject: [PATCH 06/31] testthat in DESCRIPTION

---
 DESCRIPTION | 1 +
 1 file changed, 1 insertion(+)

diff --git a/DESCRIPTION b/DESCRIPTION
index d7691f7..e03c07b 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -30,5 +30,6 @@ LazyData: true
 RoxygenNote: 5.0.0
 Suggests:
     knitr,
+    testthat,
     rmarkdown
 VignetteBuilder: knitr
-- 
GitLab


From 992ddfc04d5e931aba2c1dcbe99cfead318eeed6 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Thu, 12 Dec 2019 11:14:52 +0100
Subject: [PATCH 07/31] Continuous integration yml

---
 .gitlab-ci.yml | 10 ++++++++++
 1 file changed, 10 insertions(+)
 create mode 100644 .gitlab-ci.yml

diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
new file mode 100644
index 0000000..d7d8605
--- /dev/null
+++ b/.gitlab-ci.yml
@@ -0,0 +1,10 @@
+stages:
+  - build
+build:
+  stage: build
+  script:
+    - module load R
+    - R CMD build --resave-data .
+    - R CMD check --as-cran --no-manual ClimProjDiags_*.tar.gz
+    - R -e 'covr::package_coverage()'
+
-- 
GitLab


From 3604abf9587adc4035d00b4821631e6f3cbb67a4 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Thu, 12 Dec 2019 11:26:31 +0100
Subject: [PATCH 08/31] Remove examples included in tests

---
 R/Lon2Index.R | 47 -----------------------------------------------
 1 file changed, 47 deletions(-)

diff --git a/R/Lon2Index.R b/R/Lon2Index.R
index 5b9ff88..78e2982 100644
--- a/R/Lon2Index.R
+++ b/R/Lon2Index.R
@@ -15,58 +15,11 @@
 #'lon[pos]
 #'pos <- Lon2Index(lon, lonmin = 340, lonmax = 20)
 #'lon[pos]
-#'pos <- Lon2Index(lon, lonmin = 20, lonmax = 340)
-#'lon[pos]
-#'pos <- Lon2Index(lon, lonmin = -40, lonmax = -20)
-#'lon[pos]
-#'pos <- Lon2Index(lon, lonmin = 320, lonmax = 340)
-#'lon[pos]
-#'pos <- Lon2Index(lon, lonmin = -220, lonmax = -170)
-#'lon[pos] # fails however it works
-#'pos <- Lon2Index(lon, lonmin = -350, lonmax = -300)
-#'lon[pos] # Fails no, it works
-#'pos <- Lon2Index(lon, lonmin = -400, lonmax = -370)
-#'lon[pos] # Fails
-#'pos <- Lon2Index(lon, lonmin = 340, lonmax = 380)
-#'lon[pos] #Fails however it works for now.
 #'lon <- -180 : 180
 #'pos <- Lon2Index(lon, lonmin = -20, lonmax = 20)
 #'lon[pos]
 #'pos <- Lon2Index(lon, lonmin = 340, lonmax = 20)
 #'lon[pos]
-#'pos <- Lon2Index(lon, lonmin = 20, lonmax = 340)
-#'lon[pos]
-#'pos <- Lon2Index(lon, lonmin = -40, lonmax = -20)
-#'lon[pos]
-#'pos <- Lon2Index(lon, lonmin = 320, lonmax = 340)
-#'lon[pos]
-#'pos <- Lon2Index(lon, lonmin = -220, lonmax = -170)
-#'lon[pos] # fails
-#'pos <- Lon2Index(lon, lonmin = -350, lonmax = -300)
-#'lon[pos] # Fails
-#'pos <- Lon2Index(lon, lonmin = -400, lonmax = -370)
-#'lon[pos] # Fails
-#'pos <- Lon2Index(lon, lonmin = 340, lonmax = 380)
-#'lon[pos] #Fails however it works for now.
-#'lon <- -360 : 360
-#'pos <- Lon2Index(lon, lonmin = -20, lonmax = 20)
-#'lon[pos]
-#'pos <- Lon2Index(lon, lonmin = 340, lonmax = 20)
-#'lon[pos] #Fails
-#'pos <- Lon2Index(lon, lonmin = 20, lonmax = 340)
-#'lon[pos]
-#'pos <- Lon2Index(lon, lonmin = -40, lonmax = -20)
-#'lon[pos]
-#'pos <- Lon2Index(lon, lonmin = 320, lonmax = 340)
-#'lon[pos]
-#'pos <- Lon2Index(lon, lonmin = -220, lonmax = -170)
-#'lon[pos] 
-#'pos <- Lon2Index(lon, lonmin = -350, lonmax = -300)
-#'lon[pos] 
-#'pos <- Lon2Index(lon, lonmin = -400, lonmax = -370)
-#'lon[pos] # Fails
-#'pos <- Lon2Index(lon, lonmin = 340, lonmax = 380)
-#'lon[pos] #Fails again it works
 #'
 #'@export
 Lon2Index <- function(lon, lonmin, lonmax) {
-- 
GitLab


From c410220fb42bae865c1b59c947e659db63b8158d Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Thu, 12 Dec 2019 11:28:49 +0100
Subject: [PATCH 09/31] documentation updated

---
 man/Lon2Index.Rd | 47 -----------------------------------------------
 1 file changed, 47 deletions(-)

diff --git a/man/Lon2Index.Rd b/man/Lon2Index.Rd
index 4fecf2b..605cba4 100644
--- a/man/Lon2Index.Rd
+++ b/man/Lon2Index.Rd
@@ -25,58 +25,11 @@ pos <- Lon2Index(lon, lonmin = -20, lonmax = 20)
 lon[pos]
 pos <- Lon2Index(lon, lonmin = 340, lonmax = 20)
 lon[pos]
-pos <- Lon2Index(lon, lonmin = 20, lonmax = 340)
-lon[pos]
-pos <- Lon2Index(lon, lonmin = -40, lonmax = -20)
-lon[pos]
-pos <- Lon2Index(lon, lonmin = 320, lonmax = 340)
-lon[pos]
-pos <- Lon2Index(lon, lonmin = -220, lonmax = -170)
-lon[pos] # fails however it works
-pos <- Lon2Index(lon, lonmin = -350, lonmax = -300)
-lon[pos] # Fails no, it works
-pos <- Lon2Index(lon, lonmin = -400, lonmax = -370)
-lon[pos] # Fails
-pos <- Lon2Index(lon, lonmin = 340, lonmax = 380)
-lon[pos] #Fails however it works for now.
 lon <- -180 : 180
 pos <- Lon2Index(lon, lonmin = -20, lonmax = 20)
 lon[pos]
 pos <- Lon2Index(lon, lonmin = 340, lonmax = 20)
 lon[pos]
-pos <- Lon2Index(lon, lonmin = 20, lonmax = 340)
-lon[pos]
-pos <- Lon2Index(lon, lonmin = -40, lonmax = -20)
-lon[pos]
-pos <- Lon2Index(lon, lonmin = 320, lonmax = 340)
-lon[pos]
-pos <- Lon2Index(lon, lonmin = -220, lonmax = -170)
-lon[pos] # fails
-pos <- Lon2Index(lon, lonmin = -350, lonmax = -300)
-lon[pos] # Fails
-pos <- Lon2Index(lon, lonmin = -400, lonmax = -370)
-lon[pos] # Fails
-pos <- Lon2Index(lon, lonmin = 340, lonmax = 380)
-lon[pos] #Fails however it works for now.
-lon <- -360 : 360
-pos <- Lon2Index(lon, lonmin = -20, lonmax = 20)
-lon[pos]
-pos <- Lon2Index(lon, lonmin = 340, lonmax = 20)
-lon[pos] #Fails
-pos <- Lon2Index(lon, lonmin = 20, lonmax = 340)
-lon[pos]
-pos <- Lon2Index(lon, lonmin = -40, lonmax = -20)
-lon[pos]
-pos <- Lon2Index(lon, lonmin = 320, lonmax = 340)
-lon[pos]
-pos <- Lon2Index(lon, lonmin = -220, lonmax = -170)
-lon[pos] 
-pos <- Lon2Index(lon, lonmin = -350, lonmax = -300)
-lon[pos] 
-pos <- Lon2Index(lon, lonmin = -400, lonmax = -370)
-lon[pos] # Fails
-pos <- Lon2Index(lon, lonmin = 340, lonmax = 380)
-lon[pos] #Fails again it works
 
 }
 
-- 
GitLab


From 4cb5bb135cc3509f64bb4c9c7fafbff1d02ac921 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Thu, 12 Dec 2019 11:36:01 +0100
Subject: [PATCH 10/31] Fix documentation

---
 R/Lon2Index.R    | 4 ++--
 man/Lon2Index.Rd | 1 +
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/R/Lon2Index.R b/R/Lon2Index.R
index 78e2982..1ab7197 100644
--- a/R/Lon2Index.R
+++ b/R/Lon2Index.R
@@ -1,5 +1,5 @@
-# Obtain the index of positions for a region in longitudes
-#
+#'Obtain the index of positions for a region in longitudes
+#'
 #'@description This auxiliary function returns the index of position of a region of longitudes in a given vector of longitudes.
 #'
 #'@param lon vector of longitudes values.
diff --git a/man/Lon2Index.Rd b/man/Lon2Index.Rd
index 605cba4..b2012a6 100644
--- a/man/Lon2Index.Rd
+++ b/man/Lon2Index.Rd
@@ -2,6 +2,7 @@
 % Please edit documentation in R/Lon2Index.R
 \name{Lon2Index}
 \alias{Lon2Index}
+\title{Obtain the index of positions for a region in longitudes}
 \usage{
 Lon2Index(lon, lonmin, lonmax)
 }
-- 
GitLab


From 585490b08876eae967e94a3a7dca7f858b7721f5 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Thu, 12 Dec 2019 11:44:57 +0100
Subject: [PATCH 11/31] Add format to as.PCICt to solve check to CRAN error in
 Extreme() example

---
 R/Extremes.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/Extremes.R b/R/Extremes.R
index 46ecbaf..2fcb7aa 100644
--- a/R/Extremes.R
+++ b/R/Extremes.R
@@ -114,7 +114,7 @@ Extremes <- function(data, threshold, op = ">", min.length = 6, spells.can.span.
   if (stop_error) {
     stop("Parameter 'dates' must be of the same length as the 'time' dimension of the parameter 'data'.")
   }
-  dates <- as.PCICt(dates, cal = calendar)
+  dates <- as.PCICt(dates, cal = calendar, format = "%Y%m%d")
   dates = as.character(dates)
   jdays <- as.numeric(strftime(dates, format = "%j"))
   if (calendar == "gregorian" | calendar == "standard" | calendar == "proleptic_gregorian") {
-- 
GitLab


From 4e43e5f23656c96e1c88514ef27a5366b03da1cf Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Thu, 12 Dec 2019 11:51:27 +0100
Subject: [PATCH 12/31] add format in as.PCICt

---
 R/Extremes.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/Extremes.R b/R/Extremes.R
index 2fcb7aa..870735c 100644
--- a/R/Extremes.R
+++ b/R/Extremes.R
@@ -94,7 +94,7 @@ Extremes <- function(data, threshold, op = ">", min.length = 6, spells.can.span.
       if (is.character(dates)) {
         as.POSIXct(dates, format = "%Y%m%d")
       } else {
-        as.POSIXct(dates)
+        as.POSIXct(dates, format = "%Y%m%d")
       }
     })
     if ('try-error' %in% class(dates) | sum(is.na(dates)) == length(dates)) {
-- 
GitLab


From 99991e8e161b6694a8689da6ccae811126aaa38f Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Thu, 12 Dec 2019 12:01:18 +0100
Subject: [PATCH 13/31] trying fixing the example

---
 R/Extremes.R    | 3 ++-
 man/Extremes.Rd | 3 ++-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/R/Extremes.R b/R/Extremes.R
index 870735c..ef81024 100644
--- a/R/Extremes.R
+++ b/R/Extremes.R
@@ -27,7 +27,8 @@
 #'##Example synthetic data:
 #'data <- 1:(2 * 3 * 372 * 1)
 #'dim(data) <- c(time = 372, lon = 2, lat = 3, model = 1)
-#'time <- as.POSIXct(paste(sort(rep(1900:1911, 31)), 1, 1:31, sep = "-"), tz = "CET")
+#'time <- as.POSIXct(paste(sort(rep(1900:1911, 31)), 1, 1:31, sep = "-"), tz = "CET",
+#'                   format = "%Y-%m-%d")
 #'metadata <- list(time = list(standard_name = 'time', long_name = 'time',  calendar = 'noleap', 
 #'                             units = 'days since 1970-01-01 00:00:00', prec = 'double', 
 #'                             dim = list(list(name = 'time', unlim = FALSE))))
diff --git a/man/Extremes.Rd b/man/Extremes.Rd
index 0939fa6..61b312a 100644
--- a/man/Extremes.Rd
+++ b/man/Extremes.Rd
@@ -46,7 +46,8 @@ This routine compares data to the thresholds using the given operator, generatin
 ##Example synthetic data:
 data <- 1:(2 * 3 * 372 * 1)
 dim(data) <- c(time = 372, lon = 2, lat = 3, model = 1)
-time <- as.POSIXct(paste(sort(rep(1900:1911, 31)), 1, 1:31, sep = "-"), tz = "CET")
+time <- as.POSIXct(paste(sort(rep(1900:1911, 31)), 1, 1:31, sep = "-"), tz = "CET",
+                  format = "\%Y-\%m-\%d")
 metadata <- list(time = list(standard_name = 'time', long_name = 'time',  calendar = 'noleap', 
                             units = 'days since 1970-01-01 00:00:00', prec = 'double', 
                             dim = list(list(name = 'time', unlim = FALSE))))
-- 
GitLab


From cd1de1fcceea051b5ee47ec5fabc145dd240a649 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Thu, 12 Dec 2019 12:15:34 +0100
Subject: [PATCH 14/31] trying tz parameter

---
 R/Extremes.R | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/R/Extremes.R b/R/Extremes.R
index ef81024..d80a4f3 100644
--- a/R/Extremes.R
+++ b/R/Extremes.R
@@ -27,8 +27,7 @@
 #'##Example synthetic data:
 #'data <- 1:(2 * 3 * 372 * 1)
 #'dim(data) <- c(time = 372, lon = 2, lat = 3, model = 1)
-#'time <- as.POSIXct(paste(sort(rep(1900:1911, 31)), 1, 1:31, sep = "-"), tz = "CET",
-#'                   format = "%Y-%m-%d")
+#'time <- as.POSIXct(paste(sort(rep(1900:1911, 31)), 1, 1:31, sep = "-"), tz = "CET")
 #'metadata <- list(time = list(standard_name = 'time', long_name = 'time',  calendar = 'noleap', 
 #'                             units = 'days since 1970-01-01 00:00:00', prec = 'double', 
 #'                             dim = list(list(name = 'time', unlim = FALSE))))
@@ -93,7 +92,7 @@ Extremes <- function(data, threshold, op = ">", min.length = 6, spells.can.span.
   if (!any(class(dates) %in% c('POSIXct'))) {
     dates <- try( {
       if (is.character(dates)) {
-        as.POSIXct(dates, format = "%Y%m%d")
+        as.POSIXct(dates, format = "%Y%m%d", tz = "CET")
       } else {
         as.POSIXct(dates, format = "%Y%m%d")
       }
-- 
GitLab


From 3ec076c08520cfbf2587b21956b0fbb5d8a33e3a Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Thu, 12 Dec 2019 12:17:51 +0100
Subject: [PATCH 15/31] trying tz parameter second place

---
 R/Extremes.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/Extremes.R b/R/Extremes.R
index d80a4f3..4b07c16 100644
--- a/R/Extremes.R
+++ b/R/Extremes.R
@@ -94,7 +94,7 @@ Extremes <- function(data, threshold, op = ">", min.length = 6, spells.can.span.
       if (is.character(dates)) {
         as.POSIXct(dates, format = "%Y%m%d", tz = "CET")
       } else {
-        as.POSIXct(dates, format = "%Y%m%d")
+        as.POSIXct(dates, format = "%Y%m%d", tz = "CET")
       }
     })
     if ('try-error' %in% class(dates) | sum(is.na(dates)) == length(dates)) {
-- 
GitLab


From 2dae68d6b581e7056ee8dc61e833e674dece9b74 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Thu, 12 Dec 2019 12:22:49 +0100
Subject: [PATCH 16/31] Try latest R version

---
 .gitlab-ci.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index d7d8605..3d9d3b6 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -3,7 +3,7 @@ stages:
 build:
   stage: build
   script:
-    - module load R
+    - module load R/3.6.1-foss-2015a-bare
     - R CMD build --resave-data .
     - R CMD check --as-cran --no-manual ClimProjDiags_*.tar.gz
     - R -e 'covr::package_coverage()'
-- 
GitLab


From 6f8cae7162f0a6d101703d4a9ddffd2bd7a43cc1 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Thu, 12 Dec 2019 12:33:58 +0100
Subject: [PATCH 17/31] remove as.character

---
 .gitlab-ci.yml | 2 +-
 R/Extremes.R   | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 3d9d3b6..d7d8605 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -3,7 +3,7 @@ stages:
 build:
   stage: build
   script:
-    - module load R/3.6.1-foss-2015a-bare
+    - module load R
     - R CMD build --resave-data .
     - R CMD check --as-cran --no-manual ClimProjDiags_*.tar.gz
     - R -e 'covr::package_coverage()'
diff --git a/R/Extremes.R b/R/Extremes.R
index 4b07c16..3b503a9 100644
--- a/R/Extremes.R
+++ b/R/Extremes.R
@@ -115,7 +115,7 @@ Extremes <- function(data, threshold, op = ">", min.length = 6, spells.can.span.
     stop("Parameter 'dates' must be of the same length as the 'time' dimension of the parameter 'data'.")
   }
   dates <- as.PCICt(dates, cal = calendar, format = "%Y%m%d")
-  dates = as.character(dates)
+#  dates = as.character(dates)
   jdays <- as.numeric(strftime(dates, format = "%j"))
   if (calendar == "gregorian" | calendar == "standard" | calendar == "proleptic_gregorian") {
     year <- as.numeric(strftime(dates, format = "%Y"))
-- 
GitLab


From 7315a243b5af979c6a321a7a2e7370590e83b5c8 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Thu, 12 Dec 2019 18:28:52 +0100
Subject: [PATCH 18/31] automatic test fails in gitlab but not in win-builder

---
 .gitlab-ci.yml | 10 ----------
 1 file changed, 10 deletions(-)
 delete mode 100644 .gitlab-ci.yml

diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
deleted file mode 100644
index d7d8605..0000000
--- a/.gitlab-ci.yml
+++ /dev/null
@@ -1,10 +0,0 @@
-stages:
-  - build
-build:
-  stage: build
-  script:
-    - module load R
-    - R CMD build --resave-data .
-    - R CMD check --as-cran --no-manual ClimProjDiags_*.tar.gz
-    - R -e 'covr::package_coverage()'
-
-- 
GitLab


From 6a52d012c4d4ee1b9811aaa6f46eb67d7bbd4160 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Mon, 20 Jan 2020 10:19:03 +0100
Subject: [PATCH 19/31] WaveDuration without climdex.pcict package

---
 R/WaveDuration.R | 156 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 156 insertions(+)

diff --git a/R/WaveDuration.R b/R/WaveDuration.R
index b46ef52..441fb3a 100644
--- a/R/WaveDuration.R
+++ b/R/WaveDuration.R
@@ -231,3 +231,159 @@ WaveDuration <- function(data, threshold, op = ">", spell.length = 6, by.seasons
                                                 spells.can.span.years = TRUE, 1)
   return(result)
 }
+#' @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.
+#' It then computes the lengths of the remaining sequences of TRUE values
+#' (spells) and sums their lengths.
+#' 
+#' The \code{spells.can.span.years} option controls whether spells must always
+#' terminate at the end of a period, or whether they may continue until the
+#' criteria ceases to be met or the end of the data is reached. The default for
+#' fclimdex is FALSE.
+#' 
+#' @param daily.temp Data to compute index on.
+#' @param date.factor Date factor to split by.
+#' @param jdays Timeseries of days of year.
+#' @param thresholds The thresholds to compare to.
+#' @param op The operator to use to compare data to threshold.
+#' @param min.length The minimum spell length to be considered.
+#' @param spells.can.span.years Whether spells can span years.
+#' @param max.missing.days Maximum number of NA values per time period.
+#' @return A timeseries of maximum spell lengths for each period.
+#' @seealso \code{\link{climdex.wsdi}}.
+#' @keywords ts climate
+#' @examples
+#' 
+#' prec.dat <- c(0.1, 3.0, 4.3, 1.9, 1.3, 6.0, 0, 0, 4.0, 1)
+#' phony.date.factor <- factor(rep(1:2, each=5))
+#' 
+#' ## With spells spanning years...
+#' alttedi <- threshold.exceedance.duration.index(prec.dat,
+#' phony.date.factor, rep(1:5, 2), rep(1, 5), ">=", 2, TRUE, 1)
+#' 
+#' ## Without spells spanning years...
+#' tedi <- threshold.exceedance.duration.index(prec.dat, phony.date.factor,
+#' rep(1:5, 2), rep(1, 5), ">=", 2, FALSE, 1)
+#' 
+#' @noRd
+threshold.exceedance.duration.index <- function(daily.temp, date.factor, jdays, thresholds, op=">", min.length=6, spells.can.span.years=TRUE, max.missing.days) {
+  stopifnot(is.numeric(c(daily.temp, thresholds, min.length)), is.factor(date.factor),
+            is.function(match.fun(op)),
+            min.length > 0)
+  f <- match.fun(op)
+  na.mask <- get.na.mask(is.na(daily.temp + thresholds[jdays]), date.factor, max.missing.days)
+
+  if(spells.can.span.years) {
+    periods <- select.blocks.gt.length(f(daily.temp, thresholds[jdays]), min.length - 1)
+    return(tapply.fast(periods, date.factor, sum) * na.mask)
+  } else {
+    ## fclimdex behaviour...
+    return(tapply.fast(1:length(daily.temp), date.factor, function(idx) { sum(select.blocks.gt.length(f(daily.temp[idx], thresholds[jdays[idx]]), min.length - 1)) } ) * na.mask)
+  }
+}
+## Get NA mask given threshold and split factor
+get.na.mask <- function(x, f, threshold) {
+  return(c(1, NA)[1 + as.numeric(tapply.fast(is.na(x), f, function(y) { return(sum(y) > threshold) } ))])
+}
+#' Select blocks of TRUE values of sufficient length.
+#' 
+#' Produces a sequence of booleans of the same length as input, with sequences
+#' of TRUE values shorter than n replaced with FALSE.
+#' 
+#' This function takes a series of booleans and returns a sequence of booleans
+#' of equal length, with all sequences of TRUE of length \code{n} or shorter
+#' replaced with sequences of FALSE. NA values are replaced with
+#' \code{na.value}.
+#' 
+#' @param d Sequence of booleans.
+#' @param n Longest sequence of TRUE to replace with FALSE.
+#' @param na.value Values to replace NAs with.
+#' @return A vector of booleans, with the length \code{n} or less sequences of
+#' TRUE replaced with FALSE.
+#' @keywords ts climate
+#' @examples
+#' 
+#' ## Return only the first sequence of TRUE... second sequence will be FALSE.
+#' foo <- select.blocks.gt.length(c(rep(TRUE, 4), FALSE, rep(TRUE, 3)), 3)
+#' 
+#' @noRd
+select.blocks.gt.length <- function(d, n, na.value=FALSE) {
+  stopifnot(is.logical(d), is.numeric(n))
+
+  if(n < 1)
+    return(d)
+
+  if(n >= length(d))
+    return(rep(FALSE, length(d)))
+
+  d[is.na(d)] <- na.value
+  
+  d2 <- Reduce(function(x, y) { return(c(rep(FALSE, y), d[1:(length(d) - y)]) & x) }, 1:n, d)
+  return(Reduce(function(x, y) { return(c(d2[(y + 1):length(d2)], rep(FALSE, y)) | x) }, 1:n, d2))
+}
+## Lower overhead version of tapply
+tapply.fast <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) {
+  FUN <- if (!is.null(FUN))
+    match.fun(FUN)
+  
+  if(!is.factor(INDEX))
+    stop("INDEX must be a factor.")
+  
+  if (length(INDEX) != length(X))
+    stop("arguments must have same length")
+  
+  if (is.null(FUN))
+    return(INDEX)
+  
+  namelist <- levels(INDEX)
+  ans <- lapply(split(X, INDEX), FUN, ...)
+  
+  ans <- unlist(ans, recursive = FALSE)
+  names(ans) <- levels(INDEX)
+  return(ans)
+}
+#' Select blocks of TRUE values of sufficient length.
+#' 
+#' Produces a sequence of booleans of the same length as input, with sequences
+#' of TRUE values shorter than n replaced with FALSE.
+#' 
+#' This function takes a series of booleans and returns a sequence of booleans
+#' of equal length, with all sequences of TRUE of length \code{n} or shorter
+#' replaced with sequences of FALSE. NA values are replaced with
+#' \code{na.value}.
+#' 
+#' @param d Sequence of booleans.
+#' @param n Longest sequence of TRUE to replace with FALSE.
+#' @param na.value Values to replace NAs with.
+#' @return A vector of booleans, with the length \code{n} or less sequences of
+#' TRUE replaced with FALSE.
+#' @keywords ts climate
+#' @examples
+#' 
+#' ## Return only the first sequence of TRUE... second sequence will be FALSE.
+#' foo <- select.blocks.gt.length(c(rep(TRUE, 4), FALSE, rep(TRUE, 3)), 3)
+#' 
+#' @noRd
+select.blocks.gt.length <- function(d, n, na.value=FALSE) {
+  stopifnot(is.logical(d), is.numeric(n))
+
+  if(n < 1)
+    return(d)
+
+  if(n >= length(d))
+    return(rep(FALSE, length(d)))
+
+  d[is.na(d)] <- na.value
+  
+  d2 <- Reduce(function(x, y) { return(c(rep(FALSE, y), d[1:(length(d) - y)]) & x) }, 1:n, d)
+  return(Reduce(function(x, y) { return(c(d2[(y + 1):length(d2)], rep(FALSE, y)) | x) }, 1:n, d2))
+}
+
-- 
GitLab


From fa151b6e73e223894ea6d6039f00d47c0995c3e9 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Mon, 20 Jan 2020 10:23:59 +0100
Subject: [PATCH 20/31] Threshold already independent of climdex.pcict

---
 R/Threshold.R | 1 -
 1 file changed, 1 deletion(-)

diff --git a/R/Threshold.R b/R/Threshold.R
index 311163a..181a5f0 100644
--- a/R/Threshold.R
+++ b/R/Threshold.R
@@ -13,7 +13,6 @@
 #'@return An array with similar dimensions as the \code{data} input, but without 'time' dimension, and a new 'jdays' dimension.
 #'
 #'@import multiApply
-#'@import climdex.pcic
 #'@import PCICt
 #'@importFrom stats quantile
 #'@examples
-- 
GitLab


From 7a79b0f0a859057ac9529ede54da21f23a231a8f Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Mon, 20 Jan 2020 10:29:42 +0100
Subject: [PATCH 21/31] remove import climdex.pcict

---
 R/Extremes.R     | 1 -
 R/WaveDuration.R | 1 -
 2 files changed, 2 deletions(-)

diff --git a/R/Extremes.R b/R/Extremes.R
index 46ecbaf..bedaa39 100644
--- a/R/Extremes.R
+++ b/R/Extremes.R
@@ -21,7 +21,6 @@
 #'@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. It then computes the lengths of the remaining sequences of TRUE values (spells) and sums their lengths. The \code{spells.can.spa .years} option controls whether spells must always terminate at the end of a period, or whether they may continue until the criteria ceases to be met or the end of the data is reached. The default for fclimdex is FALSE.
 #'
 #'@import multiApply
-#'@import climdex.pcic
 #'@import PCICt
 #'@examples
 #'##Example synthetic data:
diff --git a/R/WaveDuration.R b/R/WaveDuration.R
index 441fb3a..04f4f26 100644
--- a/R/WaveDuration.R
+++ b/R/WaveDuration.R
@@ -16,7 +16,6 @@
 #'  \item\code{$result}{An array with the same dimensions as the input \code{data}, but with the time dimension reduce from daily to monthly or seasonal resolution depending on the selected resolution in \code{by.season}.}
 #'  \item\code{$years}{A vector of the years and season/months corresponding to the resolution selected in \code{by.season} and temporal length of the input \code{data}}}
 #'
-#'@import climdex.pcic 
 #'@import multiApply   
 #'@import PCICt
 #'@examples
-- 
GitLab


From 50ed54fb981f7082f9e0c87d938385a99769fbdc Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Mon, 20 Jan 2020 10:46:54 +0100
Subject: [PATCH 22/31] Climdex independent of climdex.pcict

---
 R/Climdex.R | 271 +++++++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 270 insertions(+), 1 deletion(-)

diff --git a/R/Climdex.R b/R/Climdex.R
index e182979..a762f11 100644
--- a/R/Climdex.R
+++ b/R/Climdex.R
@@ -16,7 +16,6 @@
 #'  \item\code{$result} {An array with the same dimensions as the input array, except for the temporal dimension which is renamed to 'year', moved to the first dimension position and reduce to annual resolution.}
 #'  \item\code{$years}  {A vector of the corresponding years.}}
 #'
-#'@import climdex.pcic
 #'@import multiApply
 #'@import PCICt
 #'@examples 
@@ -243,3 +242,273 @@ Climdex <- function(data, metric, threshold = NULL, base.range = NULL, dates = N
   }
   return(result)
 }
+
+#' @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
+#' period (as defined by date.factor).
+#' 
+#' The \code{spells.can.span.years} option controls whether spells must always
+#' terminate at the end of a period, or whether they may continue until the
+#' criteria ceases to be met or the end of the data is reached. The default for
+#' fclimdex is TRUE.
+#' 
+#' @param daily.prec Data to compute index on.
+#' @param date.factor Date factor to split by.
+#' @param threshold The threshold to compare to.
+#' @param op The operator to use to compare data to threshold.
+#' @param spells.can.span.years Whether spells can span years.
+#' @return A timeseries of maximum spell lengths for each period.
+#' @seealso \code{\link{climdex.cdd}}.
+#' @keywords ts climate
+#' @examples
+#' 
+#' prec.dat <- c(0.1, 3.0, 4.3, 1.9, 1.3, 6.0, 0, 0, 4.0, 1)
+#' phony.date.factor <- factor(rep(1:2, each=5))
+#' 
+#' ## With spells spanning years...
+#' cwd <- spell.length.max(prec.dat, phony.date.factor, 1, ">=", TRUE)
+#' 
+#' ## Without spells spanning years...
+#' altcwd <- spell.length.max(prec.dat, phony.date.factor, 1, ">=", FALSE)
+#' 
+#' @noRd
+spell.length.max <- function(daily.prec, date.factor, threshold, op, spells.can.span.years) {
+  bools <- match.fun(op)(daily.prec, threshold)
+
+  if(spells.can.span.years) {
+    all.true <- tapply.fast(bools, date.factor, all)
+    max.spell <- tapply.fast(get.series.lengths.at.ends(bools), date.factor, max)
+    
+    ## Mask out values which are in the middle of a spell with NA
+    na.mask <- c(1, NA)[as.integer((max.spell == 0) & all.true) + 1]
+    return(max.spell * na.mask)
+  } else {
+    return(tapply.fast(bools, date.factor, function(x) { max(get.series.lengths.at.ends(x)) }))
+  }
+}
+#' 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.
+#' 
+#' It can often be useful to know how long a series of boolean values is. This
+#' function provides a method of knowing where and how long such sequences are.
+#' 
+#' @param x Sequence of booleans.
+#' @param na.value Value to replace NAs with.
+#' @return A vector consisting of the lengths of sequences of TRUE values at
+#' the location of the last TRUE value in the sequence, and zeroes elsewhere.
+#' @keywords ts climate
+#' @examples
+#' 
+#' ## Get lengths of sequences of TRUE values in a sequence
+#' series.lengths <- get.series.lengths.at.ends(c(TRUE, TRUE, TRUE, FALSE,
+#' TRUE, FALSE, TRUE, TRUE, TRUE, TRUE, FALSE))
+#' 
+#' 
+#' @noRd
+get.series.lengths.at.ends <- function(x, na.value=FALSE) {
+  stopifnot(is.logical(x) && is.logical(na.value))
+  n <- length(x)
+  if(n == 1)
+    return(as.numeric(x))
+
+  res <- rep(0, n)
+  x[is.na(x)] <- na.value
+
+  ## Compare series to lag-1 and lag+1 series; false added to trigger state transition from TRUE at ends of series
+  start <- which(x & !(c(FALSE, x[1:(n - 1)])))
+  end <- which(x & !(c(x[2:n], FALSE)))
+  res[end] <- end - start + 1
+  return(res)
+}
+#' @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
+#' precipitation and returns the maximum n-day total precipitation per unit
+#' time, as defined by \code{date.factor}.
+#' 
+#' @param daily.prec Daily timeseries of precipitation.
+#' @param date.factor Factor to aggregate by.
+#' @param ndays Number of days in the running window.
+#' @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 maximum n-day sum of precipitation 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(ec.1018935.tmax$MAX_TEMP,
+#' ec.1018935.tmin$MIN_TEMP, ec.1018935.prec$ONE_DAY_PRECIPITATION,
+#' tmax.dates, tmin.dates, prec.dates, base.range=c(1971, 2000))
+#' 
+#' ## Compute rx5day on a monthly basis.
+#' rx5day <- nday.consec.prec.max(ci@@data$prec, ci@@date.factors$monthly, 5)
+#' 
+#' @noRd
+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)))
+  }
+  ## 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 <- running.mean(daily.prec, ndays)
+  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 Running Mean of a Vector
+#'
+#' @description Calculates the running means of a vector with a shifting window
+#'
+#' @details Returns a new vector the same length as vec, where the ith
+#' element is the mean of the bin of elements centered at the ith element
+#' of the original vector. Means cannot be calculated for elements less
+#' than half the width of the bin from the beginning or end of the vector;
+#' the result vector has NA in those positions.
+#'
+#' @param vec A vector
+#' @param bin The number of entries to average over for each mean
+#'
+#' @return a vector containing the running mean of bin elements of vec
+#'
+#' @example
+#' \dontrun { 
+#' running.mean(c(1, 2, 3, 4, 5, 6), 2) 
+#' }
+#' \dontrun { 
+#' running.mean(c(5, 5, 5, 5, 5), 4) 
+#' }
+#'@noRd
+running.mean <- function(vec, bin){
+  vec = as.vector(vec)
+  len = length(vec)
+  if (bin<=1) {
+    return (vec)
+  }
+  if (bin > len) {
+    bin = len
+  }
+  left.bin = bin%/%2
+
+  means = double(len)
+
+  right.bin = bin - left.bin - 1
+  means = c( sum(vec[1:bin]), diff(vec,bin) ) # find the first sum and the differences from it
+  means = cumsum(means)/bin                  # apply precomputed differences
+  means = c(rep(NA,left.bin), means, rep(NA,right.bin))   # extend to original vector length
+  return(means)
+}
+#' Lengths of strings of TRUE values
+#' 
+#' Computes fraction of days above or below the baseline threshold for each
+#' day, and averages them using the date factor passed in.
+#' 
+#' This function computes fractions of days above or below baseline thresholds
+#' for each day, then aggregates them using \code{date.factor}. It is used to
+#' implement TN/TX 10/90p.
+#' 
+#' @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
+#' @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(ec.1018935.tmax$MAX_TEMP,
+#' ec.1018935.tmin$MIN_TEMP, ec.1018935.prec$ONE_DAY_PRECIPITATION,
+#' tmax.dates, tmin.dates, prec.dates, base.range=c(1971, 2000))
+#' 
+#' ## Compute monthly tx90p.
+#' tx90p <- percent.days.op.threshold(ci@@data$tmax, ci@@dates, ci@@jdays,
+#'                                    ci@@date.factors$monthly,
+#'                                    ci@@quantiles$tmax$outbase$q90,
+#'                                    ci@@quantiles$tmax$inbase$q90,
+#'                                    ci@@base.range, ">",
+#'                                    ci@@max.missing.days['monthly']) *
+#'          ci@@namasks$monthly$tmax
+#' 
+#' @noRd
+percent.days.op.threshold <- function(temp, dates, jdays, date.factor, threshold.outside.base, base.thresholds, base.range, op='<', max.missing.days) {
+  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 <- 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 - 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)
+  }
+  dat[is.nan(dat)] <- NA
+  if(missing(date.factor))
+    return(dat)
+  na.mask <- get.na.mask(dat, date.factor, max.missing.days)
+  ## FIXME: Need to monthly-ize the NA mask calculation, which will be ugly.
+  ret <- tapply.fast(dat, date.factor, mean, na.rm=TRUE) * 100 * na.mask
+  ret[is.nan(ret)] <- NA
+  return(ret)
+}
-- 
GitLab


From 91c4d7d8b2c120b1389ba0adff1c75fe4e5428fc Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Mon, 20 Jan 2020 10:48:59 +0100
Subject: [PATCH 23/31] documentation updated

---
 DESCRIPTION | 1 -
 NAMESPACE   | 1 -
 2 files changed, 2 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index d7691f7..0699598 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -18,7 +18,6 @@ Depends:
     R (>= 3.2.0)
 Imports:
     multiApply (>= 2.0.0),
-    climdex.pcic,
     PCICt,
     plyr,
     stats
diff --git a/NAMESPACE b/NAMESPACE
index 8c0c72b..e844fcb 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -14,7 +14,6 @@ export(Threshold)
 export(WaveDuration)
 export(WeightedMean)
 import(PCICt)
-import(climdex.pcic)
 import(multiApply)
 importFrom(plyr,aaply)
 importFrom(stats,quantile)
-- 
GitLab


From 5749a868c4e65901d816aa54fdd0c1578d20b6a9 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Mon, 20 Jan 2020 10:57:48 +0100
Subject: [PATCH 24/31] get.years function added

---
 R/Climdex.R | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/R/Climdex.R b/R/Climdex.R
index a762f11..8c6362e 100644
--- a/R/Climdex.R
+++ b/R/Climdex.R
@@ -512,3 +512,7 @@ percent.days.op.threshold <- function(temp, dates, jdays, date.factor, threshold
   ret[is.nan(ret)] <- NA
   return(ret)
 }
+## Get year
+get.years <- function(dates) {
+  return(as.POSIXlt(dates)$year + 1900)
+}
-- 
GitLab


From bb64be1a8fce56a82608104aa4f18c91e39be91a Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Mon, 20 Jan 2020 16:03:39 +0100
Subject: [PATCH 25/31] atomic functions starting with dot

---
 R/Climdex.R      | 41 ++++++++++++++++++-----------------
 R/Extremes.R     |  2 +-
 R/WaveDuration.R | 56 +++++++++---------------------------------------
 3 files changed, 32 insertions(+), 67 deletions(-)

diff --git a/R/Climdex.R b/R/Climdex.R
index 8c6362e..1944943 100644
--- a/R/Climdex.R
+++ b/R/Climdex.R
@@ -225,17 +225,17 @@ Climdex <- function(data, metric, threshold = NULL, base.range = NULL, dates = N
 .Climdex <- function(data, threshold, dates = dates, metric = metric,
                      jdays = jdays, date.factor = NULL, base.range = NULL) {
    if (metric == "cdd") {
-    result <- spell.length.max(daily.prec = as.vector(data), date.factor = dates, 1, "<", FALSE)
+    result <- .spell.length.max(daily.prec = as.vector(data), date.factor = dates, 1, "<", FALSE)
   } else if (metric == "rx5day") {
-    result <- nday.consec.prec.max(daily.prec = as.vector(data), date.factor = dates, 
+    result <- .nday.consec.prec.max(daily.prec = as.vector(data), date.factor = dates, 
                                    ndays = 5, center.mean.on.last.day = FALSE)
   } else if (metric == "t90p" || metric == "Wx") {
-    result <- percent.days.op.threshold(temp = data, dates = dates, jdays = jdays, 
+    result <- .percent.days.op.threshold(temp = data, dates = dates, jdays = jdays, 
                                         date.factor = date.factor, threshold.outside.base = threshold, 
                                         base.thresholds = threshold, base.range = base.range, 
                                         op = ">", 20)
   } else if (metric == "t10p") {
-    result <- percent.days.op.threshold(temp = data, dates = dates, jdays = jdays, 
+    result <- .percent.days.op.threshold(temp = data, dates = dates, jdays = jdays, 
                                         date.factor = date.factor, threshold.outside.base = threshold, 
                                         base.thresholds = threshold, base.range = base.range, 
                                         op = "<", 20)
@@ -280,18 +280,18 @@ Climdex <- function(data, metric, threshold = NULL, base.range = NULL, dates = N
 #' altcwd <- spell.length.max(prec.dat, phony.date.factor, 1, ">=", FALSE)
 #' 
 #' @noRd
-spell.length.max <- function(daily.prec, date.factor, threshold, op, spells.can.span.years) {
+.spell.length.max <- function(daily.prec, date.factor, threshold, op, spells.can.span.years) {
   bools <- match.fun(op)(daily.prec, threshold)
 
   if(spells.can.span.years) {
-    all.true <- tapply.fast(bools, date.factor, all)
-    max.spell <- tapply.fast(get.series.lengths.at.ends(bools), date.factor, max)
+    all.true <- .tapply.fast(bools, date.factor, all)
+    max.spell <- .tapply.fast(get.series.lengths.at.ends(bools), date.factor, max)
     
     ## Mask out values which are in the middle of a spell with NA
     na.mask <- c(1, NA)[as.integer((max.spell == 0) & all.true) + 1]
     return(max.spell * na.mask)
   } else {
-    return(tapply.fast(bools, date.factor, function(x) { max(get.series.lengths.at.ends(x)) }))
+    return(.tapply.fast(bools, date.factor, function(x) { max(.get.series.lengths.at.ends(x)) }))
   }
 }
 #' Get series length at ends
@@ -316,7 +316,7 @@ spell.length.max <- function(daily.prec, date.factor, threshold, op, spells.can.
 #' 
 #' 
 #' @noRd
-get.series.lengths.at.ends <- function(x, na.value=FALSE) {
+.get.series.lengths.at.ends <- function(x, na.value=FALSE) {
   stopifnot(is.logical(x) && is.logical(na.value))
   n <- length(x)
   if(n == 1)
@@ -371,19 +371,19 @@ get.series.lengths.at.ends <- function(x, na.value=FALSE) {
 #' rx5day <- nday.consec.prec.max(ci@@data$prec, ci@@date.factors$monthly, 5)
 #' 
 #' @noRd
-nday.consec.prec.max <- function(daily.prec, date.factor, ndays, center.mean.on.last.day=FALSE) {
+.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)))
   }
   ## 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 <- running.mean(daily.prec, ndays)
+  prec.runsum <- .running.mean(daily.prec, ndays)
   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)
+  return(.tapply.fast(prec.runsum, date.factor, max) * ndays)
 }
 #' @title Running Mean of a Vector
 #'
@@ -408,7 +408,7 @@ nday.consec.prec.max <- function(daily.prec, date.factor, ndays, center.mean.on.
 #' running.mean(c(5, 5, 5, 5, 5), 4) 
 #' }
 #'@noRd
-running.mean <- function(vec, bin){
+.running.mean <- function(vec, bin){
   vec = as.vector(vec)
   len = length(vec)
   if (bin<=1) {
@@ -478,7 +478,7 @@ running.mean <- function(vec, bin){
 #'          ci@@namasks$monthly$tmax
 #' 
 #' @noRd
-percent.days.op.threshold <- function(temp, dates, jdays, date.factor, threshold.outside.base, base.thresholds, base.range, op='<', max.missing.days) {
+.percent.days.op.threshold <- function(temp, dates, jdays, date.factor, threshold.outside.base, base.thresholds, base.range, op='<', max.missing.days) {
   f <- match.fun(op)
   dat <- f(temp, threshold.outside.base[jdays])
   
@@ -486,7 +486,7 @@ percent.days.op.threshold <- function(temp, dates, jdays, date.factor, threshold
   ## 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 <- get.years(dates[inset])
+    years.base <- .get.years(dates[inset])
 
     ## Get number of base years, subset temp data to base period only.
     temp.base <- temp[inset]
@@ -496,7 +496,7 @@ percent.days.op.threshold <- function(temp, dates, jdays, date.factor, threshold
     ## 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 - get.years(base.range)[1]) * bdim[1]
+    yday.byr.indices <- jdays.base + (years.base - .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])
 
@@ -506,13 +506,14 @@ percent.days.op.threshold <- function(temp, dates, jdays, date.factor, threshold
   dat[is.nan(dat)] <- NA
   if(missing(date.factor))
     return(dat)
-  na.mask <- get.na.mask(dat, date.factor, max.missing.days)
+  na.mask <- .get.na.mask(dat, date.factor, max.missing.days)
   ## FIXME: Need to monthly-ize the NA mask calculation, which will be ugly.
-  ret <- tapply.fast(dat, date.factor, mean, na.rm=TRUE) * 100 * na.mask
+  ret <- .tapply.fast(dat, date.factor, mean, na.rm=TRUE) * 100 * na.mask
   ret[is.nan(ret)] <- NA
   return(ret)
 }
 ## Get year
-get.years <- function(dates) {
+.get.years <- function(dates) {
   return(as.POSIXlt(dates)$year + 1900)
-}
+tr(a)
+
diff --git a/R/Extremes.R b/R/Extremes.R
index bedaa39..30dd569 100644
--- a/R/Extremes.R
+++ b/R/Extremes.R
@@ -211,7 +211,7 @@ Extremes <- function(data, threshold, op = ">", min.length = 6, spells.can.span.
 }
 .Extremes <- function(data, threshold, date.factor, jdays, op, min.length, spells.can.span.years,
                       max.missing.days) {
-  result <- threshold.exceedance.duration.index(data, date.factor, jdays, 
+  result <- .threshold.exceedance.duration.index(data, date.factor, jdays, 
                                                 threshold,op, min.length, 
                                                 spells.can.span.years, max.missing.days)
 }
diff --git a/R/WaveDuration.R b/R/WaveDuration.R
index 04f4f26..de0d573 100644
--- a/R/WaveDuration.R
+++ b/R/WaveDuration.R
@@ -223,7 +223,7 @@ WaveDuration <- function(data, threshold, op = ">", spell.length = 6, by.seasons
 }
 .WaveDuration <- function(data, threshold, date.factor = date.factor,
                           jdays = jdays,  op = op, spell.length = spell.length) {
-  result <- threshold.exceedance.duration.index(daily.temp = data, 
+  result <- .threshold.exceedance.duration.index(daily.temp = data, 
                                                 date.factor = date.factor, jdays = jdays,
                                                 thresholds = threshold,
                                                 op = op, min.length = spell.length, 
@@ -273,24 +273,24 @@ WaveDuration <- function(data, threshold, op = ">", spell.length = 6, by.seasons
 #' rep(1:5, 2), rep(1, 5), ">=", 2, FALSE, 1)
 #' 
 #' @noRd
-threshold.exceedance.duration.index <- function(daily.temp, date.factor, jdays, thresholds, op=">", min.length=6, spells.can.span.years=TRUE, max.missing.days) {
+.threshold.exceedance.duration.index <- function(daily.temp, date.factor, jdays, thresholds, op=">", min.length=6, spells.can.span.years=TRUE, max.missing.days) {
   stopifnot(is.numeric(c(daily.temp, thresholds, min.length)), is.factor(date.factor),
             is.function(match.fun(op)),
             min.length > 0)
   f <- match.fun(op)
-  na.mask <- get.na.mask(is.na(daily.temp + thresholds[jdays]), date.factor, max.missing.days)
+  na.mask <- .get.na.mask(is.na(daily.temp + thresholds[jdays]), date.factor, max.missing.days)
 
   if(spells.can.span.years) {
-    periods <- select.blocks.gt.length(f(daily.temp, thresholds[jdays]), min.length - 1)
-    return(tapply.fast(periods, date.factor, sum) * na.mask)
+    periods <- .select.blocks.gt.length(f(daily.temp, thresholds[jdays]), min.length - 1)
+    return(.tapply.fast(periods, date.factor, sum) * na.mask)
   } else {
     ## fclimdex behaviour...
-    return(tapply.fast(1:length(daily.temp), date.factor, function(idx) { sum(select.blocks.gt.length(f(daily.temp[idx], thresholds[jdays[idx]]), min.length - 1)) } ) * na.mask)
+    return(.tapply.fast(1:length(daily.temp), date.factor, function(idx) { sum(.select.blocks.gt.length(f(daily.temp[idx], thresholds[jdays[idx]]), min.length - 1)) } ) * na.mask)
   }
 }
 ## Get NA mask given threshold and split factor
-get.na.mask <- function(x, f, threshold) {
-  return(c(1, NA)[1 + as.numeric(tapply.fast(is.na(x), f, function(y) { return(sum(y) > threshold) } ))])
+.get.na.mask <- function(x, f, threshold) {
+  return(c(1, NA)[1 + as.numeric(.tapply.fast(is.na(x), f, function(y) { return(sum(y) > threshold) } ))])
 }
 #' Select blocks of TRUE values of sufficient length.
 #' 
@@ -314,7 +314,7 @@ get.na.mask <- function(x, f, threshold) {
 #' foo <- select.blocks.gt.length(c(rep(TRUE, 4), FALSE, rep(TRUE, 3)), 3)
 #' 
 #' @noRd
-select.blocks.gt.length <- function(d, n, na.value=FALSE) {
+.select.blocks.gt.length <- function(d, n, na.value=FALSE) {
   stopifnot(is.logical(d), is.numeric(n))
 
   if(n < 1)
@@ -329,7 +329,7 @@ select.blocks.gt.length <- function(d, n, na.value=FALSE) {
   return(Reduce(function(x, y) { return(c(d2[(y + 1):length(d2)], rep(FALSE, y)) | x) }, 1:n, d2))
 }
 ## Lower overhead version of tapply
-tapply.fast <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) {
+.tapply.fast <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) {
   FUN <- if (!is.null(FUN))
     match.fun(FUN)
   
@@ -349,40 +349,4 @@ tapply.fast <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) {
   names(ans) <- levels(INDEX)
   return(ans)
 }
-#' Select blocks of TRUE values of sufficient length.
-#' 
-#' Produces a sequence of booleans of the same length as input, with sequences
-#' of TRUE values shorter than n replaced with FALSE.
-#' 
-#' This function takes a series of booleans and returns a sequence of booleans
-#' of equal length, with all sequences of TRUE of length \code{n} or shorter
-#' replaced with sequences of FALSE. NA values are replaced with
-#' \code{na.value}.
-#' 
-#' @param d Sequence of booleans.
-#' @param n Longest sequence of TRUE to replace with FALSE.
-#' @param na.value Values to replace NAs with.
-#' @return A vector of booleans, with the length \code{n} or less sequences of
-#' TRUE replaced with FALSE.
-#' @keywords ts climate
-#' @examples
-#' 
-#' ## Return only the first sequence of TRUE... second sequence will be FALSE.
-#' foo <- select.blocks.gt.length(c(rep(TRUE, 4), FALSE, rep(TRUE, 3)), 3)
-#' 
-#' @noRd
-select.blocks.gt.length <- function(d, n, na.value=FALSE) {
-  stopifnot(is.logical(d), is.numeric(n))
-
-  if(n < 1)
-    return(d)
-
-  if(n >= length(d))
-    return(rep(FALSE, length(d)))
-
-  d[is.na(d)] <- na.value
-  
-  d2 <- Reduce(function(x, y) { return(c(rep(FALSE, y), d[1:(length(d) - y)]) & x) }, 1:n, d)
-  return(Reduce(function(x, y) { return(c(d2[(y + 1):length(d2)], rep(FALSE, y)) | x) }, 1:n, d2))
-}
 
-- 
GitLab


From 0eabf0021cfca368c15e14d07846e541fce1f6a7 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Mon, 20 Jan 2020 16:11:30 +0100
Subject: [PATCH 26/31] Updated automatic documentation

---
 R/Climdex.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/Climdex.R b/R/Climdex.R
index 1944943..261815c 100644
--- a/R/Climdex.R
+++ b/R/Climdex.R
@@ -515,5 +515,5 @@ Climdex <- function(data, metric, threshold = NULL, base.range = NULL, dates = N
 ## Get year
 .get.years <- function(dates) {
   return(as.POSIXlt(dates)$year + 1900)
-tr(a)
+}
 
-- 
GitLab


From b5e2168479ec39e8dd286e31e632c3d4aa19440e Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Mon, 20 Jan 2020 16:55:02 +0100
Subject: [PATCH 27/31] Climdex fixed and Selbox using Lon2Index

---
 R/Climdex.R |  2 +-
 R/SelBox.R  | 12 +++++++-----
 2 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/R/Climdex.R b/R/Climdex.R
index 261815c..20c3e43 100644
--- a/R/Climdex.R
+++ b/R/Climdex.R
@@ -373,7 +373,7 @@ Climdex <- function(data, metric, threshold = NULL, base.range = NULL, dates = N
 #' @noRd
 .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)))
+    return(suppressWarnings(.tapply.fast(daily.prec, date.factor, max, na.rm=TRUE)))
   }
   ## Ends of the data will be de-emphasized (padded with zero precip data); NAs replaced with 0
   daily.prec[is.na(daily.prec)] <- 0
diff --git a/R/SelBox.R b/R/SelBox.R
index 608e19e..5b7cd59 100644
--- a/R/SelBox.R
+++ b/R/SelBox.R
@@ -83,11 +83,13 @@ SelBox <- function(data, lon, lat, region, londim = NULL, latdim = NULL, mask =
   } else {
     LatIdx <- which(lat <= region[3] | lat >= region[4])
   }
-  if (region[1] <= region[2]) {
-    LonIdx <- which(lon >= region[1] & lon <= region[2])
-  } else {
-    LonIdx <- which(lon >= region[1] | lon <= region[2])
-  }
+  #if (region[1] <= region[2]) {
+  #  LonIdx <- which(lon >= region[1] & lon <= region[2])
+  #} else {
+  #  LonIdx <- which(lon >= region[1] | lon <= region[2])
+  #}
+  LonIdx <- Lon2Index(lon, lonmin = region[1], lonmax = region[2])
+
   data <- Subset(data, along = londim, indices = LonIdx, drop = "none")
   data <- Subset(data, along = latdim, indices = LatIdx, drop = "none")
   if (!is.null(mask)) {
-- 
GitLab


From aff21eb0a294746e2b8b21ea36541d81e0b534f4 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Tue, 21 Jan 2020 10:59:16 +0100
Subject: [PATCH 28/31] previous version Extremes function with working
 examples

---
 R/Extremes.R    | 9 +++++----
 man/Extremes.Rd | 3 +--
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/R/Extremes.R b/R/Extremes.R
index a609340..8fc5ca2 100644
--- a/R/Extremes.R
+++ b/R/Extremes.R
@@ -91,9 +91,9 @@ Extremes <- function(data, threshold, op = ">", min.length = 6, spells.can.span.
   if (!any(class(dates) %in% c('POSIXct'))) {
     dates <- try( {
       if (is.character(dates)) {
-        as.POSIXct(dates, format = "%Y%m%d", tz = "CET")
+        as.POSIXct(dates, format = "%Y%m%d")
       } else {
-        as.POSIXct(dates, format = "%Y%m%d", tz = "CET")
+        as.POSIXct(dates)
       }
     })
     if ('try-error' %in% class(dates) | sum(is.na(dates)) == length(dates)) {
@@ -113,8 +113,8 @@ Extremes <- function(data, threshold, op = ">", min.length = 6, spells.can.span.
   if (stop_error) {
     stop("Parameter 'dates' must be of the same length as the 'time' dimension of the parameter 'data'.")
   }
-  dates <- as.PCICt(dates, cal = calendar, format = "%Y%m%d")
-#  dates = as.character(dates)
+  dates <- as.PCICt(dates, cal = calendar)
+  dates = as.character(dates)
   jdays <- as.numeric(strftime(dates, format = "%j"))
   if (calendar == "gregorian" | calendar == "standard" | calendar == "proleptic_gregorian") {
     year <- as.numeric(strftime(dates, format = "%Y"))
@@ -215,3 +215,4 @@ Extremes <- function(data, threshold, op = ">", min.length = 6, spells.can.span.
                                                 threshold,op, min.length, 
                                                 spells.can.span.years, max.missing.days)
 }
+
diff --git a/man/Extremes.Rd b/man/Extremes.Rd
index 61b312a..0939fa6 100644
--- a/man/Extremes.Rd
+++ b/man/Extremes.Rd
@@ -46,8 +46,7 @@ This routine compares data to the thresholds using the given operator, generatin
 ##Example synthetic data:
 data <- 1:(2 * 3 * 372 * 1)
 dim(data) <- c(time = 372, lon = 2, lat = 3, model = 1)
-time <- as.POSIXct(paste(sort(rep(1900:1911, 31)), 1, 1:31, sep = "-"), tz = "CET",
-                  format = "\%Y-\%m-\%d")
+time <- as.POSIXct(paste(sort(rep(1900:1911, 31)), 1, 1:31, sep = "-"), tz = "CET")
 metadata <- list(time = list(standard_name = 'time', long_name = 'time',  calendar = 'noleap', 
                             units = 'days since 1970-01-01 00:00:00', prec = 'double', 
                             dim = list(list(name = 'time', unlim = FALSE))))
-- 
GitLab


From 4d85a6564b1133e3a3f3f388619935bddaacefcc Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Tue, 21 Jan 2020 13:16:16 +0100
Subject: [PATCH 29/31] reference to climdex.pcic and bumped version number

---
 DESCRIPTION    | 6 +++---
 R/Climdex.R    | 5 ++++-
 man/Climdex.Rd | 5 +++++
 3 files changed, 12 insertions(+), 4 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 8293274..3259032 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,11 +1,11 @@
 Package: ClimProjDiags
 Title: Set of Tools to Compute Various Climate Indices
-Version: 0.0.4
+Version: 0.1.0
 Authors@R: c(
     person("BSC-CNS", role = c("aut", "cph")),
-    person("Alasdair", "Hunter", , "alasdair.hunter@bsc.es", role = "aut"),   
-    person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = c("aut", "cre")),
+    person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-8568-3071")),
     person("Nicolau", "Manubens", , "nicolau.manubens@bsc.es", role = "ctb"),
+    person("Alasdair", "Hunter", , "alasdair.hunter@bsc.es", role = "aut"),
     person("Louis-Philippe", "Caron", , "louis-philippe.caron@bsc.es", role = "ctb"))
 Description: Set of tools to compute metrics and indices for climate analysis.
     The package provides functions to compute extreme indices, evaluate the
diff --git a/R/Climdex.R b/R/Climdex.R
index 20c3e43..e8203fb 100644
--- a/R/Climdex.R
+++ b/R/Climdex.R
@@ -47,6 +47,9 @@
 #'
 #'clim <- Climdex(data, metric = "t90p", threshold = thres)
 #'str(clim)
+#'@references David Bronaugh for the Pacific Climate Impacts Consortium (2015).
+#'  climdex.pcic: PCIC Implementation of Climdex Routines. R package
+#'  version 1.1-6. http://CRAN.R-project.org/package=climdex.pcic
 #'@export
 Climdex <- function(data, metric, threshold = NULL, base.range = NULL, dates = NULL, timedim = NULL, 
                     calendar = NULL, ncores = NULL) {
@@ -285,7 +288,7 @@ Climdex <- function(data, metric, threshold = NULL, base.range = NULL, dates = N
 
   if(spells.can.span.years) {
     all.true <- .tapply.fast(bools, date.factor, all)
-    max.spell <- .tapply.fast(get.series.lengths.at.ends(bools), date.factor, max)
+    max.spell <- .tapply.fast(.get.series.lengths.at.ends(bools), date.factor, max)
     
     ## Mask out values which are in the middle of a spell with NA
     na.mask <- c(1, NA)[as.integer((max.spell == 0) & all.true) + 1]
diff --git a/man/Climdex.Rd b/man/Climdex.Rd
index 3ff6bb7..669570a 100644
--- a/man/Climdex.Rd
+++ b/man/Climdex.Rd
@@ -63,4 +63,9 @@ str(thres)
 clim <- Climdex(data, metric = "t90p", threshold = thres)
 str(clim)
 }
+\references{
+David Bronaugh for the Pacific Climate Impacts Consortium (2015).
+ climdex.pcic: PCIC Implementation of Climdex Routines. R package
+ version 1.1-6. http://CRAN.R-project.org/package=climdex.pcic
+}
 
-- 
GitLab


From 7e3c1062192d1fa49a0bc645ccd9806a2c27e9b8 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Tue, 22 Dec 2020 12:10:59 +0100
Subject: [PATCH 30/31] Change version and add climdes.pcic dependency

---
 DESCRIPTION      |   5 +-
 NAMESPACE        |   1 +
 R/Climdex.R      | 282 +----------------------------------------------
 R/Extremes.R     |   3 +-
 R/WaveDuration.R | 122 +-------------------
 5 files changed, 13 insertions(+), 400 deletions(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 3259032..431d3bb 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: ClimProjDiags
 Title: Set of Tools to Compute Various Climate Indices
-Version: 0.1.0
+Version: 0.1.1
 Authors@R: c(
     person("BSC-CNS", role = c("aut", "cph")),
     person("Nuria", "Perez-Zanon", , "nuria.perez@bsc.es", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-8568-3071")),
@@ -20,8 +20,9 @@ Imports:
     multiApply (>= 2.0.0),
     PCICt,
     plyr,
+    climdex.pcic,
     stats
-License: LGPL-3
+License: Apache License 2.0
 URL: https://earth.bsc.es/gitlab/es/ClimProjDiags
 BugReports: https://earth.bsc.es/gitlab/es/ClimProjDiags/issues
 Encoding: UTF-8
diff --git a/NAMESPACE b/NAMESPACE
index 03fa0de..8364b61 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -15,6 +15,7 @@ export(Threshold)
 export(WaveDuration)
 export(WeightedMean)
 import(PCICt)
+import(climdex.pcic)
 import(multiApply)
 importFrom(plyr,aaply)
 importFrom(stats,quantile)
diff --git a/R/Climdex.R b/R/Climdex.R
index e8203fb..a13b15b 100644
--- a/R/Climdex.R
+++ b/R/Climdex.R
@@ -17,6 +17,7 @@
 #'  \item\code{$years}  {A vector of the corresponding years.}}
 #'
 #'@import multiApply
+#'@import climdex.pcic
 #'@import PCICt
 #'@examples 
 #'##Example synthetic data:
@@ -228,17 +229,17 @@ Climdex <- function(data, metric, threshold = NULL, base.range = NULL, dates = N
 .Climdex <- function(data, threshold, dates = dates, metric = metric,
                      jdays = jdays, date.factor = NULL, base.range = NULL) {
    if (metric == "cdd") {
-    result <- .spell.length.max(daily.prec = as.vector(data), date.factor = dates, 1, "<", FALSE)
+    result <- spell.length.max(daily.prec = as.vector(data), date.factor = dates, 1, "<", FALSE)
   } else if (metric == "rx5day") {
-    result <- .nday.consec.prec.max(daily.prec = as.vector(data), date.factor = dates, 
+    result <- nday.consec.prec.max(daily.prec = as.vector(data), date.factor = dates, 
                                    ndays = 5, center.mean.on.last.day = FALSE)
   } else if (metric == "t90p" || metric == "Wx") {
-    result <- .percent.days.op.threshold(temp = data, dates = dates, jdays = jdays, 
+    result <- percent.days.op.threshold(temp = data, dates = dates, jdays = jdays, 
                                         date.factor = date.factor, threshold.outside.base = threshold, 
                                         base.thresholds = threshold, base.range = base.range, 
                                         op = ">", 20)
   } else if (metric == "t10p") {
-    result <- .percent.days.op.threshold(temp = data, dates = dates, jdays = jdays, 
+    result <- percent.days.op.threshold(temp = data, dates = dates, jdays = jdays, 
                                         date.factor = date.factor, threshold.outside.base = threshold, 
                                         base.thresholds = threshold, base.range = base.range, 
                                         op = "<", 20)
@@ -246,277 +247,4 @@ Climdex <- function(data, metric, threshold = NULL, base.range = NULL, dates = N
   return(result)
 }
 
-#' @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
-#' period (as defined by date.factor).
-#' 
-#' The \code{spells.can.span.years} option controls whether spells must always
-#' terminate at the end of a period, or whether they may continue until the
-#' criteria ceases to be met or the end of the data is reached. The default for
-#' fclimdex is TRUE.
-#' 
-#' @param daily.prec Data to compute index on.
-#' @param date.factor Date factor to split by.
-#' @param threshold The threshold to compare to.
-#' @param op The operator to use to compare data to threshold.
-#' @param spells.can.span.years Whether spells can span years.
-#' @return A timeseries of maximum spell lengths for each period.
-#' @seealso \code{\link{climdex.cdd}}.
-#' @keywords ts climate
-#' @examples
-#' 
-#' prec.dat <- c(0.1, 3.0, 4.3, 1.9, 1.3, 6.0, 0, 0, 4.0, 1)
-#' phony.date.factor <- factor(rep(1:2, each=5))
-#' 
-#' ## With spells spanning years...
-#' cwd <- spell.length.max(prec.dat, phony.date.factor, 1, ">=", TRUE)
-#' 
-#' ## Without spells spanning years...
-#' altcwd <- spell.length.max(prec.dat, phony.date.factor, 1, ">=", FALSE)
-#' 
-#' @noRd
-.spell.length.max <- function(daily.prec, date.factor, threshold, op, spells.can.span.years) {
-  bools <- match.fun(op)(daily.prec, threshold)
-
-  if(spells.can.span.years) {
-    all.true <- .tapply.fast(bools, date.factor, all)
-    max.spell <- .tapply.fast(.get.series.lengths.at.ends(bools), date.factor, max)
-    
-    ## Mask out values which are in the middle of a spell with NA
-    na.mask <- c(1, NA)[as.integer((max.spell == 0) & all.true) + 1]
-    return(max.spell * na.mask)
-  } else {
-    return(.tapply.fast(bools, date.factor, function(x) { max(.get.series.lengths.at.ends(x)) }))
-  }
-}
-#' 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.
-#' 
-#' It can often be useful to know how long a series of boolean values is. This
-#' function provides a method of knowing where and how long such sequences are.
-#' 
-#' @param x Sequence of booleans.
-#' @param na.value Value to replace NAs with.
-#' @return A vector consisting of the lengths of sequences of TRUE values at
-#' the location of the last TRUE value in the sequence, and zeroes elsewhere.
-#' @keywords ts climate
-#' @examples
-#' 
-#' ## Get lengths of sequences of TRUE values in a sequence
-#' series.lengths <- get.series.lengths.at.ends(c(TRUE, TRUE, TRUE, FALSE,
-#' TRUE, FALSE, TRUE, TRUE, TRUE, TRUE, FALSE))
-#' 
-#' 
-#' @noRd
-.get.series.lengths.at.ends <- function(x, na.value=FALSE) {
-  stopifnot(is.logical(x) && is.logical(na.value))
-  n <- length(x)
-  if(n == 1)
-    return(as.numeric(x))
-
-  res <- rep(0, n)
-  x[is.na(x)] <- na.value
-
-  ## Compare series to lag-1 and lag+1 series; false added to trigger state transition from TRUE at ends of series
-  start <- which(x & !(c(FALSE, x[1:(n - 1)])))
-  end <- which(x & !(c(x[2:n], FALSE)))
-  res[end] <- end - start + 1
-  return(res)
-}
-#' @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
-#' precipitation and returns the maximum n-day total precipitation per unit
-#' time, as defined by \code{date.factor}.
-#' 
-#' @param daily.prec Daily timeseries of precipitation.
-#' @param date.factor Factor to aggregate by.
-#' @param ndays Number of days in the running window.
-#' @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 maximum n-day sum of precipitation 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(ec.1018935.tmax$MAX_TEMP,
-#' ec.1018935.tmin$MIN_TEMP, ec.1018935.prec$ONE_DAY_PRECIPITATION,
-#' tmax.dates, tmin.dates, prec.dates, base.range=c(1971, 2000))
-#' 
-#' ## Compute rx5day on a monthly basis.
-#' rx5day <- nday.consec.prec.max(ci@@data$prec, ci@@date.factors$monthly, 5)
-#' 
-#' @noRd
-.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)))
-  }
-  ## 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 <- .running.mean(daily.prec, ndays)
-  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 Running Mean of a Vector
-#'
-#' @description Calculates the running means of a vector with a shifting window
-#'
-#' @details Returns a new vector the same length as vec, where the ith
-#' element is the mean of the bin of elements centered at the ith element
-#' of the original vector. Means cannot be calculated for elements less
-#' than half the width of the bin from the beginning or end of the vector;
-#' the result vector has NA in those positions.
-#'
-#' @param vec A vector
-#' @param bin The number of entries to average over for each mean
-#'
-#' @return a vector containing the running mean of bin elements of vec
-#'
-#' @example
-#' \dontrun { 
-#' running.mean(c(1, 2, 3, 4, 5, 6), 2) 
-#' }
-#' \dontrun { 
-#' running.mean(c(5, 5, 5, 5, 5), 4) 
-#' }
-#'@noRd
-.running.mean <- function(vec, bin){
-  vec = as.vector(vec)
-  len = length(vec)
-  if (bin<=1) {
-    return (vec)
-  }
-  if (bin > len) {
-    bin = len
-  }
-  left.bin = bin%/%2
-
-  means = double(len)
-
-  right.bin = bin - left.bin - 1
-  means = c( sum(vec[1:bin]), diff(vec,bin) ) # find the first sum and the differences from it
-  means = cumsum(means)/bin                  # apply precomputed differences
-  means = c(rep(NA,left.bin), means, rep(NA,right.bin))   # extend to original vector length
-  return(means)
-}
-#' Lengths of strings of TRUE values
-#' 
-#' Computes fraction of days above or below the baseline threshold for each
-#' day, and averages them using the date factor passed in.
-#' 
-#' This function computes fractions of days above or below baseline thresholds
-#' for each day, then aggregates them using \code{date.factor}. It is used to
-#' implement TN/TX 10/90p.
-#' 
-#' @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
-#' @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(ec.1018935.tmax$MAX_TEMP,
-#' ec.1018935.tmin$MIN_TEMP, ec.1018935.prec$ONE_DAY_PRECIPITATION,
-#' tmax.dates, tmin.dates, prec.dates, base.range=c(1971, 2000))
-#' 
-#' ## Compute monthly tx90p.
-#' tx90p <- percent.days.op.threshold(ci@@data$tmax, ci@@dates, ci@@jdays,
-#'                                    ci@@date.factors$monthly,
-#'                                    ci@@quantiles$tmax$outbase$q90,
-#'                                    ci@@quantiles$tmax$inbase$q90,
-#'                                    ci@@base.range, ">",
-#'                                    ci@@max.missing.days['monthly']) *
-#'          ci@@namasks$monthly$tmax
-#' 
-#' @noRd
-.percent.days.op.threshold <- function(temp, dates, jdays, date.factor, threshold.outside.base, base.thresholds, base.range, op='<', max.missing.days) {
-  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 <- .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 - .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)
-  }
-  dat[is.nan(dat)] <- NA
-  if(missing(date.factor))
-    return(dat)
-  na.mask <- .get.na.mask(dat, date.factor, max.missing.days)
-  ## FIXME: Need to monthly-ize the NA mask calculation, which will be ugly.
-  ret <- .tapply.fast(dat, date.factor, mean, na.rm=TRUE) * 100 * na.mask
-  ret[is.nan(ret)] <- NA
-  return(ret)
-}
-## Get year
-.get.years <- function(dates) {
-  return(as.POSIXlt(dates)$year + 1900)
-}
 
diff --git a/R/Extremes.R b/R/Extremes.R
index 8fc5ca2..f3e7b16 100644
--- a/R/Extremes.R
+++ b/R/Extremes.R
@@ -22,6 +22,7 @@
 #'
 #'@import multiApply
 #'@import PCICt
+#'@import climdex.pcic
 #'@examples
 #'##Example synthetic data:
 #'data <- 1:(2 * 3 * 372 * 1)
@@ -211,7 +212,7 @@ Extremes <- function(data, threshold, op = ">", min.length = 6, spells.can.span.
 }
 .Extremes <- function(data, threshold, date.factor, jdays, op, min.length, spells.can.span.years,
                       max.missing.days) {
-  result <- .threshold.exceedance.duration.index(data, date.factor, jdays, 
+  result <- threshold.exceedance.duration.index(data, date.factor, jdays, 
                                                 threshold,op, min.length, 
                                                 spells.can.span.years, max.missing.days)
 }
diff --git a/R/WaveDuration.R b/R/WaveDuration.R
index de0d573..23c2d97 100644
--- a/R/WaveDuration.R
+++ b/R/WaveDuration.R
@@ -18,6 +18,7 @@
 #'
 #'@import multiApply   
 #'@import PCICt
+#'@import climdex.pcic
 #'@examples
 #'##Example synthetic data:
 #'data <- 1:(2 * 3 * 31 * 5)
@@ -223,130 +224,11 @@ WaveDuration <- function(data, threshold, op = ">", spell.length = 6, by.seasons
 }
 .WaveDuration <- function(data, threshold, date.factor = date.factor,
                           jdays = jdays,  op = op, spell.length = spell.length) {
-  result <- .threshold.exceedance.duration.index(daily.temp = data, 
+  result <- threshold.exceedance.duration.index(daily.temp = data, 
                                                 date.factor = date.factor, jdays = jdays,
                                                 thresholds = threshold,
                                                 op = op, min.length = spell.length, 
                                                 spells.can.span.years = TRUE, 1)
   return(result)
 }
-#' @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.
-#' It then computes the lengths of the remaining sequences of TRUE values
-#' (spells) and sums their lengths.
-#' 
-#' The \code{spells.can.span.years} option controls whether spells must always
-#' terminate at the end of a period, or whether they may continue until the
-#' criteria ceases to be met or the end of the data is reached. The default for
-#' fclimdex is FALSE.
-#' 
-#' @param daily.temp Data to compute index on.
-#' @param date.factor Date factor to split by.
-#' @param jdays Timeseries of days of year.
-#' @param thresholds The thresholds to compare to.
-#' @param op The operator to use to compare data to threshold.
-#' @param min.length The minimum spell length to be considered.
-#' @param spells.can.span.years Whether spells can span years.
-#' @param max.missing.days Maximum number of NA values per time period.
-#' @return A timeseries of maximum spell lengths for each period.
-#' @seealso \code{\link{climdex.wsdi}}.
-#' @keywords ts climate
-#' @examples
-#' 
-#' prec.dat <- c(0.1, 3.0, 4.3, 1.9, 1.3, 6.0, 0, 0, 4.0, 1)
-#' phony.date.factor <- factor(rep(1:2, each=5))
-#' 
-#' ## With spells spanning years...
-#' alttedi <- threshold.exceedance.duration.index(prec.dat,
-#' phony.date.factor, rep(1:5, 2), rep(1, 5), ">=", 2, TRUE, 1)
-#' 
-#' ## Without spells spanning years...
-#' tedi <- threshold.exceedance.duration.index(prec.dat, phony.date.factor,
-#' rep(1:5, 2), rep(1, 5), ">=", 2, FALSE, 1)
-#' 
-#' @noRd
-.threshold.exceedance.duration.index <- function(daily.temp, date.factor, jdays, thresholds, op=">", min.length=6, spells.can.span.years=TRUE, max.missing.days) {
-  stopifnot(is.numeric(c(daily.temp, thresholds, min.length)), is.factor(date.factor),
-            is.function(match.fun(op)),
-            min.length > 0)
-  f <- match.fun(op)
-  na.mask <- .get.na.mask(is.na(daily.temp + thresholds[jdays]), date.factor, max.missing.days)
-
-  if(spells.can.span.years) {
-    periods <- .select.blocks.gt.length(f(daily.temp, thresholds[jdays]), min.length - 1)
-    return(.tapply.fast(periods, date.factor, sum) * na.mask)
-  } else {
-    ## fclimdex behaviour...
-    return(.tapply.fast(1:length(daily.temp), date.factor, function(idx) { sum(.select.blocks.gt.length(f(daily.temp[idx], thresholds[jdays[idx]]), min.length - 1)) } ) * na.mask)
-  }
-}
-## Get NA mask given threshold and split factor
-.get.na.mask <- function(x, f, threshold) {
-  return(c(1, NA)[1 + as.numeric(.tapply.fast(is.na(x), f, function(y) { return(sum(y) > threshold) } ))])
-}
-#' Select blocks of TRUE values of sufficient length.
-#' 
-#' Produces a sequence of booleans of the same length as input, with sequences
-#' of TRUE values shorter than n replaced with FALSE.
-#' 
-#' This function takes a series of booleans and returns a sequence of booleans
-#' of equal length, with all sequences of TRUE of length \code{n} or shorter
-#' replaced with sequences of FALSE. NA values are replaced with
-#' \code{na.value}.
-#' 
-#' @param d Sequence of booleans.
-#' @param n Longest sequence of TRUE to replace with FALSE.
-#' @param na.value Values to replace NAs with.
-#' @return A vector of booleans, with the length \code{n} or less sequences of
-#' TRUE replaced with FALSE.
-#' @keywords ts climate
-#' @examples
-#' 
-#' ## Return only the first sequence of TRUE... second sequence will be FALSE.
-#' foo <- select.blocks.gt.length(c(rep(TRUE, 4), FALSE, rep(TRUE, 3)), 3)
-#' 
-#' @noRd
-.select.blocks.gt.length <- function(d, n, na.value=FALSE) {
-  stopifnot(is.logical(d), is.numeric(n))
-
-  if(n < 1)
-    return(d)
-
-  if(n >= length(d))
-    return(rep(FALSE, length(d)))
-
-  d[is.na(d)] <- na.value
-  
-  d2 <- Reduce(function(x, y) { return(c(rep(FALSE, y), d[1:(length(d) - y)]) & x) }, 1:n, d)
-  return(Reduce(function(x, y) { return(c(d2[(y + 1):length(d2)], rep(FALSE, y)) | x) }, 1:n, d2))
-}
-## Lower overhead version of tapply
-.tapply.fast <- function (X, INDEX, FUN = NULL, ..., simplify = TRUE) {
-  FUN <- if (!is.null(FUN))
-    match.fun(FUN)
-  
-  if(!is.factor(INDEX))
-    stop("INDEX must be a factor.")
-  
-  if (length(INDEX) != length(X))
-    stop("arguments must have same length")
-  
-  if (is.null(FUN))
-    return(INDEX)
-  
-  namelist <- levels(INDEX)
-  ans <- lapply(split(X, INDEX), FUN, ...)
-  
-  ans <- unlist(ans, recursive = FALSE)
-  names(ans) <- levels(INDEX)
-  return(ans)
-}
 
-- 
GitLab


From 835129f9583a9470eadf62231edec263ddc197e4 Mon Sep 17 00:00:00 2001
From: nperez <nperez@bscearth372.int.bsc.es>
Date: Thu, 7 Jan 2021 09:59:39 +0100
Subject: [PATCH 31/31] url

---
 DESCRIPTION | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/DESCRIPTION b/DESCRIPTION
index 431d3bb..2956cce 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -24,12 +24,13 @@ Imports:
     stats
 License: Apache License 2.0
 URL: https://earth.bsc.es/gitlab/es/ClimProjDiags
-BugReports: https://earth.bsc.es/gitlab/es/ClimProjDiags/issues
+BugReports: https://earth.bsc.es/gitlab/es/ClimProjDiags/-/issues
 Encoding: UTF-8
 LazyData: true
 RoxygenNote: 5.0.0
 Suggests:
     knitr,
     testthat,
+    markdown,
     rmarkdown
 VignetteBuilder: knitr
-- 
GitLab