Merge pull request #141 from CeresBarros/135-calc_-functions-should-o…

…utput-nas-is

135 `calc_*` functions should output NAs is there is no available data

.github/workflows/R-CMD-check.yaml

@@ -54,6 +54,7 @@ jobs:
           extra-packages: |
             any::rcmdcheck
             any::Rcpp
+            any::withr
 
       - uses: r-lib/actions/check-r-package@v2
         with:

DESCRIPTION

@@ -35,7 +35,8 @@ Imports:
 Suggests: 
     parallel,
     rmarkdown,
-    testthat (>= 3.0.0)
+    testthat (>= 3.0.0),
+    withr
 Depends: 
     R (>= 4.0)
 Config/testthat/edition: 3

R/append_clim_vars.R

@@ -667,7 +667,9 @@ append_clim_vars <- function(dt, vars) {
   )
 
   # Append vars except default one
-  for (var in vars[!vars %in% sprintf(c("PPT%02d", "Tmax%02d", "Tmin%02d"), sort(rep(1:12, 3)))]) {
+  vars2 <- vars[!vars %in% sprintf(c("PPT%02d", "Tmax%02d", "Tmin%02d"), sort(rep(1:12, 3)))]
+  vars2 <- vars2[order(match(vars2, names(appenders)))]   ## run functions in the order of appenders
+  for (var in vars2) {
     f(var)
   }
 

R/calc_DD.R

@@ -10,6 +10,8 @@ calc_DD_m_above <- function(tm, k, a, b, t0, beta, c) {
   i <- which(tm >= k)
   DD_m[i] <- c + beta * tm[i]
 
+  DD_m[is.na(tm)] <- tm[is.na(tm)] ## use tm[is.na(tm)] to respect NA type
+
   return(DD_m)
 }
 
@@ -25,6 +27,8 @@ calc_DD_m_below <- function(tm, k, a, b, t0, beta, c) {
   i <- which(tm <= k)
   DD_m[i] <- c + beta * tm[i]
 
+  DD_m[is.na(tm)] <- tm[is.na(tm)] ## use tm[is.na(tm)] to respect NA type
+
   return(DD_m)
 }
 

R/calc_EMT_EXT.R

@@ -11,7 +11,6 @@ calc_EMT <- function(t_min_list, td) {
 
   year_min <- do.call(pmin, t_min_list)
 
-
   -23.02164 + 0.77908 * tmin1 + 0.67048 * tmin12 + 0.01075 * year_min^2 + 0.11565 * td
 }
 
@@ -30,6 +29,5 @@ calc_EXT <- function(t_max_list, td) {
 
   year_max <- do.call(pmax, t_max_list)
 
-
   10.64245 + -1.92005 * tmax7 + 0.04816 * tmax7^2 + 2.51176 * tmax8 - 0.03088 * tmax8^2 - 0.01311 * year_max^2 + 0.33167 * td - 0.001 * td^2
 }

R/calc_Eref_CMD.R

@@ -17,6 +17,9 @@ calc_Eref <- function(m, tmin, tmax, latitude) {
     calc_S0_I(day_julian[m], tmean[i], latitude[i]) *
     (tmean[i] + 17.8) * sqrt(tmax[i] - tmin[i]) *
     (1.18 - 0.0065 * latitude[i])
+
+  Eref[is.na(tmax)] <- tmax[is.na(tmax)] ## use tmax[is.na(tmax)] to respect NA type
+
   return(Eref)
 }
 
@@ -28,6 +31,10 @@ calc_CMD <- function(Eref, PPT) {
   CMD <- numeric(length(Eref))
   i <- which(Eref > PPT)
   CMD[i] <- Eref[i] - PPT[i]
+
+  ## return 0s to NaNs if missing values 
+  CMD[is.na(Eref)] <- Eref[is.na(Eref)]   ## use Eref[is.na(Eref)] to respect NA type
+
   return(CMD)
 }
 

R/helper-testInit.R

@@ -0,0 +1,22 @@
+# loads and libraries indicated
+
+testInit <- function(libraries = character()) {
+
+  pf <- parent.frame()
+
+  if (length(libraries)) {
+    libraries <- unique(libraries)
+    loadedAlready <- vapply(libraries, function(pkg)
+      any(grepl(paste0("package:", pkg), search())), FUN.VALUE = logical(1))
+    libraries <- libraries[!loadedAlready]
+
+    if (length(libraries)) {
+      pkgsLoaded <- unlist(lapply(libraries, requireNamespace, quietly = TRUE))
+      if (!all(pkgsLoaded)) {
+        lapply(libraries[!pkgsLoaded], skip_if_not_installed)
+      }
+      suppressWarnings(lapply(libraries, withr::local_package, .local_envir = pf))
+    }
+  }
+}
+

README.md

@@ -30,7 +30,7 @@ To do so, it uses a remote PostGIS database, and optionally caches data locally.
 
 -   Historical observational time series (1902-2022), currently limited to the ClimateNA time series (Wang et al., 2016)
 
--   Multiple historical (1851-2010) and future (2001-2100) climate model simulations for each of 13 CMIP6 global climate models, in monthly time series and 20-year normals
+-   Multiple historical (1851-2014) and future (2015-2100) climate model simulations for each of 13 CMIP6 global climate models, in monthly time series and 20-year normals
 
 -   User selection of single or multiple climate variables, with derived variables following the ClimateNA methodology of Wang et al. (2016).
 

tests/testthat/setup.R

@@ -0,0 +1,3 @@
+# Preload all Suggests so that examples don't take tons of time
+requireNamespace("parallel", quietly = TRUE)
+

tests/testthat/test-calcfuns.R

@@ -1,31 +1,87 @@
 test_that("calc_* functions work", {
-  climr:::calc_DD_below_0(2, -14)
+  expect_identical(round(climr:::calc_DD_below_0(2, -14), 4), 393.9186)
+  expect_identical(climr:::calc_DD_below_0(2, NA_real_), NA_real_)
+
+  expect_identical(round(climr:::calc_DD_above_5(2, -14, "All"), 4), 0.2144)
+  expect_identical(climr:::calc_DD_above_5(2, NA, "All"), NA_real_)
 
-  ## TODO: what are the correct outputs?
+  expect_identical(round(climr:::calc_DD_below_18(2, -14), 4), 892.0826)
+  expect_identical(climr:::calc_DD_below_18(2, NA_real_), NA_real_)
 
-  # expect_???
-
-  climr:::calc_DD_above_5(2, -14, "All")
-
-  climr:::calc_DD_below_18(2, -14)
-
-  climr:::calc_DD_above_18(2, -14, "All")
+  expect_identical(round(climr:::calc_DD_above_18(2, -14, "All"), 4), 0.0001)
+  expect_identical(climr:::calc_DD_above_18(2, NA, "All"), NA_real_)
 
   t_min_list <- list(
     "1" = -35, "2" = -32, "3" = -25, "4" = -10,
     "5" = -5, "6" = 3, "7" = 15, "8" = 17, "9" = 10, "10" = -5,
     "11" = -20, "12" = -30
   )
 
-  climr:::calc_bFFP(td = 30, NFFD = 10, t_min_list = t_min_list)
+  expect_identical(round(climr:::calc_bFFP(td = 30, NFFD = 10, t_min_list = t_min_list), 4),
+                   214.5964)
+  expect_identical(climr:::calc_bFFP(td = 30, NFFD = NA, t_min_list = t_min_list), NA_real_)
+  expect_identical(climr:::calc_bFFP(td = NA, NFFD = 10, t_min_list = t_min_list), NA_real_)
 
-  climr:::calc_eFFP(NFFD = 10, t_min_list = t_min_list)
+  expect_identical(round(climr:::calc_eFFP(NFFD = 10, t_min_list = t_min_list), 4),
+                   265.4581)
+  expect_identical(climr:::calc_eFFP(NFFD = NA, t_min_list = t_min_list), NA_real_)
 
-  climr:::calc_FFP(bFFP = 214.5964, eFFP = 265.4581)
+  expect_identical(round(climr:::calc_FFP(bFFP = 214.5964, eFFP = 265.4581), 4),
+                   50.8617)
+  expect_identical(climr:::calc_FFP(bFFP = NA, eFFP = 265.4581), NA_real_)
+  expect_identical(climr:::calc_FFP(bFFP = 214.5964, eFFP = NA), NA_real_)
 
-  climr:::calc_NFFD(3, 2.05)
+  expect_identical(round(climr:::calc_NFFD(3, 2.05), 4), 21.1018)
+  expect_identical(climr:::calc_NFFD(3, NA_real_), NA_real_)
+
+  expect_identical(round(climr:::calc_PAS(4, 2, 600), 4), 308.4204)
+  expect_identical(climr:::calc_PAS(4, NA, 600), NA_real_)
+  expect_identical(climr:::calc_PAS(4, 2, NA_real_), NA_real_)
+
+  expect_identical(round(climr:::calc_RH(tmin = 10, tmax = 40), 4), 28.5378)
+  expect_identical(climr:::calc_RH(tmin = NA, tmax = 40), NA_real_)
+  expect_identical(climr:::calc_RH(tmin = 10, tmax = NA), NA_real_)
+})
 
-  climr:::calc_PAS(4, 2, 600)
 
-  climr:::calc_RH(tmin = 10, tmax = 40)
+test_that("calc_* give sensible outputs", {
+  library(pool)
+  library(data.table)
+  library(terra)
+
+  dbCon <- data_connect()
+
+  ## the following includes NAs for the test
+  xyz <- data.frame(lon = c(-128, -125, -128, -125), lat = c(50, 50, 48, 48), elev = runif(4))
+  thebb <- get_bb(xyz)
+
+  normalbc <- normal_input(dbCon = dbCon, normal = "normal_bc", bbox = thebb, cache = TRUE) 
+
+  gcm <- gcm_input(dbCon, bbox = thebb, 
+                   gcm = c("BCC-CSM2-MR"), 
+                   ssp = c("ssp126"), 
+                   period = "2041_2060",
+                   max_run = 0,
+                   cache = TRUE)
+
+  set.seed(678)  ## a situation with known NAs (ocean where elev is 0)
+  n <- 20
+  sample_xyz <- data.frame(lon = runif(n, xmin(normalbc$dem2_WNA), xmax(normalbc$dem2_WNA)), 
+                           lat = runif(n, ymin(normalbc$dem2_WNA), ymax(normalbc$dem2_WNA)), 
+                           elev = NA)
+  sample_xyz[, 3] <- terra::extract(normalbc$dem2_WNA, sample_xyz[, 1:2], method = "bilinear")[, -1L]
+
+  ds_res_bc <- downscale(sample_xyz, normal = normalbc, gcm = gcm, var = list_variables())
+
+  sample_xyz$ID <- 1:nrow(sample_xyz)
+  ds_res_bc <- as.data.table(sample_xyz)[ds_res_bc, on = .(ID)]
+  ds_res_bc[, .(lat, lon, elev, Tmax, PPT01, CMD, Eref)]
+
+  ## if elevation of inout climate data are NA, downscaled variables should be as well.
+  expect_true(all(is.na(ds_res_bc[is.na(elev), .SD, .SDcols = list_variables()])))
+  expect_true(all(is.na(ds_res_bc[is.na(PPT01), .SD, .SDcols = list_variables()])))
+  ## conversely, if there is input data, there should be no NAs in downscaled vars
+  expect_true(all(!is.na(ds_res_bc[!is.na(PPT01), .SD, .SDcols = list_variables()])))
+
+  ## TODO: more sanity checks?
 })