Merge pull request #69 from geco-bern/fix-CRU-multisite-downscaling

Fix CRU multisite downscaling
geco-bern · Dec 6, 2024 · 1b6b19e · 1b6b19e
2 parents 6f40e07 + b8b79b0
commit 1b6b19e
Show file tree

Hide file tree

Showing 2 changed files with 73 additions and 11 deletions.
diff --git a/R/ingest_globalfields.R b/R/ingest_globalfields.R
@@ -346,8 +346,9 @@ ingest_globalfields <- function(
       # expand monthly to daily data
 
       if (length(cruvars)>0){
-        df_out <- expand_clim_cru_monthly( mdf, cruvars ) %>%
-          right_join( df_out, by = "date" )
+        df_out <- left_join(df_out,
+                            expand_clim_cru_monthly( mdf, cruvars ),
+                            by = c("date", "sitename") )
       }
 
       if ("vpd" %in% getvars){
@@ -924,15 +925,12 @@ ingest_globalfields_cru_byvar <- function( siteinfo, dir, varnam ){
 
 expand_clim_cru_monthly <- function( mdf, cruvars ){
 
-  # ensure this function is always called with a single site only
-  stopifnot(length(unique(mdf$sitename)) == 1)
-  # for multiple sites the code would need to be adapted, e.g.:
-  # ddf2 <- mdf |>
-  #   group_split(sitename, year) |>
-  #   purrr::map(\(df) expand_clim_cru_monthly_byyr(first(df$year), df, cruvars))
-
-  ddf <- purrr::map(as.list(unique(mdf$year)),
-      ~expand_clim_cru_monthly_byyr( ., mdf, cruvars ) ) %>%
+  ddf <- mdf |>
+    # apply it separately for each site and each year
+    group_split(sitename, year) |>
+    purrr::map(\(df) expand_clim_cru_monthly_byyr(first(df$year), df, cruvars) |>
+                 mutate(sitename = first(df$sitename)) #ensure to keep sitename
+                 ) |>
     bind_rows()
 
   return( ddf )

diff --git a/tests/testthat/test_CRU_WFDEI_NDEP.R b/tests/testthat/test_CRU_WFDEI_NDEP.R
@@ -122,6 +122,70 @@ test_that("test CRU data (monthly and downscaled daily)", {
 
 })
 
+
+test_that("test CRU data multisite downscaling (monthly and downscaled daily)", {
+  skip_on_cran()
+  library(dplyr)
+  library(tidyr)
+  library(ingestr)
+  library(testthat)
+
+  siteinfo_test <- tibble(
+    sitename   = c("Reichetal_Colorado", "Reichetal_New_Mexico", "Reichetal_Venezuela", "Reichetal_Wisconsin", "Lulea"),
+    lon        = c(-105.60, -107.00, -67.05, -90.00, 22.15),
+    lat        = c(  40.05,   34.00,   1.93,  42.50, 65.59),
+    elv        = c(  3360L,   1620L,   120L,   275L,    0L),
+    year_start = c(2010, 2010, 2010, 2010, 2010),
+    year_end   = c(2015, 2015, 2015, 2015, 2015))
+
+  # Daily:
+  # site-separate downscaling
+  df_cru_daily_separate <- siteinfo_test[c(2,3,5),] |> 
+    rowwise() |> group_split() |> # do it separately for each row
+    lapply(function(curr_site_inf){
+      ingestr::ingest(
+        siteinfo = curr_site_inf,
+        source = "cru",
+        getvars = c("temp", "ccov"), # , "prec" NOTE: we can't test prec, unless we specify a seed somewhere
+        dir = "/data/archive/cru_harris_2024/data", 
+        timescale = "d")}) |> bind_rows() |> ungroup() |> unnest(data)
+
+  # site-combined downscaling
+  df_cru_daily_combined <- siteinfo_test[c(2,3,5),] |> 
+    ungroup() |> group_split() |> # do it together for all rows
+    lapply(function(curr_site_inf){
+      ingestr::ingest(
+        siteinfo = curr_site_inf,
+        source = "cru",
+        getvars = c("temp", "ccov"), # , "prec" NOTE: we can't test prec, unless we specify a seed somewhere
+        dir = "/data/archive/cru_harris_2024/data", 
+        timescale = "d")}) |> bind_rows() |> ungroup() |> unnest(data)
+
+  testthat::expect_equal(df_cru_daily_separate, 
+                         df_cru_daily_combined)
+
+  # # Monthly:
+  # # site-combined monthly
+  # df_cru_monthly_combined <- siteinfo_test[c(2,3,5),] |> 
+  #   ungroup() |> group_split() |> # do it together for all rows
+  #   lapply(function(curr_site_inf){
+  #     ingestr::ingest(
+  #       siteinfo = curr_site_inf,
+  #       source = "cru",
+  #       getvars = c("temp", "ccov"), # , "prec" NOTE: we can't test prec, unless we specify a seed somewhere
+  #       dir = "/data/archive/cru_harris_2024/data", 
+  #       timescale = "m")}) |> bind_rows() |> ungroup() |> unnest(data)
+  # 
+  # # Illustration of failing test
+  # library(ggplot2)
+  # p1 <- ggplot(df_cru_monthly_combined, aes(y=temp, x=date, linetype=sitename, color=sitename)) + ggtitle("df_cru_monthly_combined") + geom_point() # CORRECT
+  # p2 <- ggplot(df_cru_daily_separate, aes(y=temp, x=date, linetype=sitename, color=sitename))   + ggtitle("df_cru_daily_separate")   + geom_line()   # CORRECT DOWNSCALING
+  # p3 <- ggplot(df_cru_daily_combined, aes(y=temp, x=date, linetype=sitename, color=sitename))   + ggtitle("df_cru_daily_combined")   + geom_line()   # WRONG DOWNSCALING: Gives same time series for each site
+  # gridExtra::grid.arrange(p1, p2, p3) 
+
+})
+
+
 test_that("test WATCH_WFDEI data (daily)", {
   skip_on_cran()