Merge branch 'water_consumption' into 'master'

Implement new approach for bluewater PB based on consumptive water use

See merge request tess/boundaries!18

.buildlibrary

@@ -1,4 +1,4 @@
-ValidationKey: '2082496'
+ValidationKey: '2205390'
 AutocreateReadme: yes
 AcceptedWarnings:
 - 'Warning: package ''.*'' was built under R version'
@@ -9,3 +9,4 @@ allowLinterWarnings: no
 AddInReadme: inst/README.md
 AddLogoReadme: inst/img/logo.png
 enforceVersionUpdate: no
+skipCoverage: no

.github/workflows/check.yaml

@@ -63,6 +63,6 @@ jobs:
         shell: Rscript {0}
         run: |
           nonDummyTests <- setdiff(list.files("./tests/testthat/"), c("test-dummy.R", "_snaps"))
-          if(length(nonDummyTests) > 0) covr::codecov(quiet = FALSE)
+          if(length(nonDummyTests) > 0 && !lucode2:::loadBuildLibraryConfig()[["skipCoverage"]]) covr::codecov(quiet = FALSE)
         env:
           NOT_CRAN: "true"

CITATION.cff

@@ -2,8 +2,8 @@ cff-version: 1.2.0
 message: If you use this software, please cite it using the metadata from this file.
 type: software
 title: 'boundaries: Planetary Boundary Status based on LPJmL simulations'
-version: 1.0.4
-date-released: '2024-10-28'
+version: 1.1.0
+date-released: '2024-11-22'
 abstract: A systematic approach to quantify the status of the terrestrial planetary
   boundaries based on the Dynamic Global Vegetation Model (DGVM) Lund-Potsdam-Jena
   managed Land (LPJmL) hosted at the Potsdam Institute for Climate Impact Research
@@ -14,18 +14,22 @@ authors:
   given-names: Johanna
   email: braun@pik-potsdam.de
   orcid: https://orcid.org/0000-0002-8809-1044
+  orcid: https://orcid.org/0000-0002-8809-1044
 - family-names: Breier
   given-names: Jannes
   email: jannesbr@pik-potsdam.de
   orcid: https://orcid.org/0000-0002-9055-6904
+  orcid: https://orcid.org/0000-0002-9055-6904
 - family-names: Stenzel
   given-names: Fabian
   email: stenzel@pik-potsdam.de
   orcid: https://orcid.org/0000-0002-5109-0048
+  orcid: https://orcid.org/0000-0002-5109-0048
 - family-names: Vanelli
   given-names: Caterina
   email: caterina@pik-potsdam.de
   orcid: https://orcid.org/0000-0001-6552-0320
+  orcid: https://orcid.org/0000-0001-6552-0320
 license: AGPL-3.0
 repository-code: https://github.com/PIK-tess/boundaries
 doi: 10.5281/zenodo.11550559

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: boundaries
 Title: Planetary Boundary Status based on LPJmL simulations
-Version: 1.0.4
+Version: 1.1.0
 Authors@R: c(
     person("Johanna", "Braun", , "braun@pik-potsdam.de", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-8809-1044")),
     person("Jannes", "Breier", , "jannesbr@pik-potsdam.de", role = c("aut"), comment = c(ORCID = "0000-0002-9055-6904")),
@@ -40,9 +40,7 @@ Imports:
     methods,
     tidyselect,
     matrixStats,
-    biospheremetrics
-Suggests: 
-    testthat (>= 3.0.0),
+    biospheremetrics,
     cowplot,
     ggh4x,
     ggnewscale,
@@ -57,8 +55,10 @@ Suggests:
     ggrepel,
     scales,
     tidyterra
+Suggests: 
+    testthat (>= 3.0.0)
 Remotes:
     github::stenzelf/biospheremetrics,
     github::yjunechoe/ggtrace
 Config/testthat/edition: 3
-Date: 2024-10-28
+Date: 2024-11-22

R/as_risk_level.R

@@ -54,22 +54,30 @@ as_risk_level <- function(
       if (methods::is(thresholds, "list")) {
         risk_level <- (control_variable - thresholds[["holocene"]]) /
           (thresholds[["pb"]] - thresholds[["holocene"]])
-        attr(risk_level, "thresholds") <-
-          list(
-            holocene = 0, pb = 1,
-            highrisk = (thresholds[["highrisk"]] - thresholds[["holocene"]]) /
-              (thresholds[["pb"]] - thresholds[["holocene"]])
-          )
+        risk_level <- set_attributes(
+          risk_level,
+          thresholds =
+            list(
+              holocene = 0, pb = 1,
+              highrisk = (thresholds[["highrisk"]] - thresholds[["holocene"]]) /
+                (thresholds[["pb"]] - thresholds[["holocene"]])
+            ),
+          overwrite = TRUE
+        )
 
       } else if (methods::is(thresholds, "array")) {
         risk_level <- (control_variable - thresholds[, , "holocene"]) /
           (thresholds[, , "pb"] - thresholds[, , "holocene"])
-        attr(risk_level, "thresholds") <-
-          list(
-            holocene = 0, pb = 1,
-            highrisk = (thresholds[, , "highrisk"] - thresholds[, , "holocene"]) / # nolint
-              (thresholds[, , "pb"] - thresholds[, , "holocene"])
-          )
+        risk_level <- set_attributes(
+          risk_level,
+          thresholds =
+            list(
+              holocene = 0, pb = 1,
+              highrisk = (thresholds[, , "highrisk"] - thresholds[, , "holocene"]) / # nolint
+                (thresholds[, , "pb"] - thresholds[, , "holocene"])
+            ),
+          overwrite = TRUE
+        )
       }
 
     } else if (normalize == "increasing risk") {
@@ -93,11 +101,14 @@ as_risk_level <- function(
           1 - (control_variable[control_variable < thresholds[["pb"]] &
                      is.na(control_variable) == FALSE] - thresholds[["pb"]]) / # nolint
           (thresholds[["holocene"]] - thresholds[["pb"]])
-        attr(risk_level, "thresholds") <-
-          list(
-            holocene = 0,
-            pb = 1, highrisk = 2
-          )
+        risk_level <- set_attributes(
+          risk_level,
+          thresholds =
+            list(
+              holocene = 0, pb = 1, highrisk = 2
+            ),
+          overwrite = TRUE
+        )
         # alternative, if no additional normalization from holocene to pb:
         # holocene = (thresholds[["holocene"]] -
         #                   thresholds[["pb"]]) /
@@ -122,11 +133,14 @@ as_risk_level <- function(
           1 - (control_variable[subset_safe] -
                thresh_pb[subset_safe]) /
           (thresh_holocene[subset_safe] - thresh_pb[subset_safe])
-        attr(risk_level, "thresholds") <-
-          list(
-            holocene = 0,
-            pb = 1, highrisk = 2
-          )
+        risk_level <- set_attributes(
+          risk_level,
+          thresholds =
+            list(
+              holocene = 0, pb = 1, highrisk = 2
+            ),
+          overwrite = TRUE
+        )
 
       } else {
         stop("thresholds must be a list or array")

R/biosphere_status.R

@@ -272,14 +272,15 @@ biosphere_status <- function(
     x = control_variable_raw * 100,
     time_series_avg = time_series_avg
   )
-  attr(control_variable, "spatial_scale") <- spatial_scale
-  attr(control_variable, "thresholds") <- thresholds
-  attr(control_variable, "unit") <- list_unit("biosphere", approach,
-                                              spatial_scale)
-  attr(control_variable, "control_variable") <- "BioCol"
-  attr(control_variable, "long_name") <- list_long_name("biosphere")
 
-  class(control_variable) <- c("control_variable")
+  control_variable <- set_attributes(
+    control_variable,
+    approach,
+    "biosphere",
+    spatial_scale,
+    thresholds
+  )
+
   return(control_variable)
 
 } # end of biosphere_status

R/bluewater_status.R

@@ -82,7 +82,7 @@ bluewater_status <- function(
   # verify available methods and resolution
   approach <- match.arg(
     approach,
-    c("gerten2020", "wang-erlandsson2022", "porkka2024")
+    c("gerten2020", "wang-erlandsson2022", "porkka2024", "rockstroem2009")
   )
   spatial_scale <- match.arg(spatial_scale, c("global", "subglobal", "grid"))
 
@@ -108,7 +108,40 @@ bluewater_status <- function(
     )
 
   } else if (spatial_scale %in% c("subglobal", "global")) {
-    if (!approach %in% c("wang-erlandsson2022", "porkka2024")) {
+    if (approach %in% c("wang-erlandsson2022", "porkka2024")) {
+      control_variable <- calc_water_deviations(
+        files_scenario = files_scenario,
+        files_reference = files_reference,
+        spatial_scale = spatial_scale,
+        time_span_scenario = time_span_scenario,
+        time_span_reference =  time_span_reference,
+        approach = approach,
+        time_series_avg = time_series_avg,
+        config_args = config_args,
+        thresholds = thresholds,
+        variable = "discharge"
+      )
+    } else if (approach == "rockstroem2009" && spatial_scale == "global") {
+      # calculater bluewater consumption
+      bw_consumption <- NULL
+      bw_consumption <- calc_bw_consumption(
+        files_scenario = files_scenario,
+        files_reference = files_reference,
+        time_span_scenario = time_span_scenario,
+        time_series_avg = time_series_avg,
+        config_args = config_args
+      )
+
+      # aggregate to global value (conversion to km3/yr)
+      dim_remain <- names(dim(bw_consumption))[
+        names(dim(bw_consumption)) != "cell"
+      ]
+      # conversion from l/yr to km3/yr
+      control_variable <- apply(bw_consumption,
+                                dim_remain, sum, na.rm = TRUE) * 10^-12
+
+
+    } else {
       stop(
         "Approach \"",
         approach,
@@ -117,21 +150,15 @@ bluewater_status <- function(
         "."
       )
     }
-
-    control_variable <- calc_water_deviations(
-      files_scenario = files_scenario,
-      files_reference = files_reference,
-      spatial_scale = spatial_scale,
-      time_span_scenario = time_span_scenario,
-      time_span_reference =  time_span_reference,
-      approach = approach,
-      time_series_avg = time_series_avg,
-      config_args = config_args,
-      thresholds = thresholds,
-      variable = "discharge"
-    )
   }
-  attr(control_variable, "long_name") <- list_long_name("bluewater")
+
+  control_variable <- set_attributes(
+    control_variable,
+    approach,
+    "bluewater",
+    spatial_scale,
+    thresholds
+  )
 
   return(control_variable)
 }
@@ -254,14 +281,85 @@ calc_bluewater_efrs <- function(
   # if ratio is above >75% transgression (red)
   # define PB thresholds as attributes
 
-  attr(control_variable, "control_variable") <-
-    "EFR transgression to uncertainty ratio"
-  attr(control_variable, "thresholds") <- thresholds
-  attr(control_variable, "spatial_scale") <- spatial_scale
-  attr(control_variable, "unit") <- list_unit("bluewater", approach,
-                                              spatial_scale)
-  attr(control_variable, "long_name") <- list_long_name("bluewater")
-
-  class(control_variable) <- c("control_variable")
   return(control_variable)
 }
+
+calc_irrig_consumption <- function(
+  files_scenario,
+  files_reference,
+  time_span_scenario = time_span_scenario,
+  config_args = list()
+) {
+
+  # irrigation
+  irrig <- NULL
+  irrig %<-% read_io_format(
+    files_scenario$irrig,
+    time_span_scenario,
+    aggregate = list(band = sum, month = sum),
+    spatial_subset = config_args$spatial_subset
+  )
+  # evaporative conveyance losses
+  conv_loss_evap <- NULL
+  conv_loss_evap %<-% read_io_format(
+    files_scenario$conv_loss_evap,
+    time_span_scenario,
+    aggregate = list(band = sum, month = sum),
+    spatial_subset = config_args$spatial_subset
+  )
+  # bluewater return flow (from irrigation)
+  return_flow_b <- NULL
+  return_flow_b %<-% read_io_format(
+    files_scenario$return_flow_b,
+    time_span_scenario,
+    aggregate = list(band = sum, month = sum),
+    spatial_subset = config_args$spatial_subset
+  )
+  # calculate terrestrial area
+  terr_area <- lpjmlkit::read_io(
+    files_scenario$terr_area
+  ) %>%
+    conditional_subset(config_args$spatial_subset) %>%
+    lpjmlkit::as_array()
+
+  terr_area <- terr_area[, , 1]
+
+  # calculate bluewater consumption for irrigation, in l/yr
+  consumption_irrig <- (irrig + conv_loss_evap - return_flow_b) * terr_area
+
+  return(consumption_irrig)
+}
+
+calc_bw_consumption <- function(
+  files_scenario,
+  files_reference,
+  time_span_scenario = time_span_scenario,
+  time_series_avg = NULL,
+  config_args = list()
+) {
+  irrig_consumption <- NULL
+  irrig_consumption <- calc_irrig_consumption(
+    files_scenario = files_scenario,
+    files_reference = files_reference,
+    time_span_scenario = time_span_scenario,
+    config_args = config_args
+  )
+  # read in water consumption for HIL (houshoulds, industry, livestock)
+  consumption_hil <- NULL
+  consumption_hil %<-% read_io_format(
+    files_scenario$wateruse_hil,
+    time_span_scenario,
+    aggregate = list(band = sum, month = sum),
+    spatial_subset = config_args$spatial_subset
+  )
+
+  # calculate total bluewater consumption in l/yr
+  total_consumption <- irrig_consumption + consumption_hil
+
+  # average over time
+  avg_total_consumption <- aggregate_time(
+    x = total_consumption,
+    time_series_avg = time_series_avg
+  )
+  return(avg_total_consumption)
+}