Merge branch 'implement_user_requests' into 'main'

implement user requests

See merge request lpjml/lpjmlstats!16

.buildlibrary

@@ -1,4 +1,4 @@
-ValidationKey: '617303'
+ValidationKey: '796880'
 AutocreateReadme: yes
 AcceptedWarnings:
 - 'Warning: package ''.*'' was built under R version'

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: 'lpjmlstats: Statistical tools for LPJmL data analysis'
-version: 0.3.1
-date-released: '2024-07-09'
+version: 0.4.0
+date-released: '2024-07-18'
 abstract: This package provides statistical tools for LPJmL data analysis to be used
   for benchmarking LPJmL outputs.
 authors:

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: lpjmlstats
 Title: Statistical tools for LPJmL data analysis
-Version: 0.3.1
+Version: 0.4.0
 Authors@R: c(
     person("David","Hötten", , "[email protected]", role = c("aut", "cre")),
     person("Jannes","Breier", , "[email protected]", role = c("aut"), comment = c(ORCID = "0000-0002-9055-6904"))
@@ -9,7 +9,7 @@ Description: This package provides statistical tools for LPJmL data analysis
     to be used for benchmarking LPJmL outputs.
 License: AGPL-3
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.2
 Encoding: UTF-8
 URL: https://github.com/PIK-LPJmL/lpjmlstats
 BugReports: https://github.com/PIK-LPJmL/lpjmlstats/issues
@@ -18,7 +18,6 @@ Imports:
     units,
     lpjmlkit,
     Matrix,
-    abind,
     methods,
     memoise,
     rlang,
@@ -34,7 +33,8 @@ Imports:
     cli,
     raster,
     tidyselect,
-    R6
+    R6,
+    rnaturalearth
 Suggests: 
     testthat (>= 3.0.0),
     maps,
@@ -44,5 +44,4 @@ Config/testthat/edition: 3
 VignetteBuilder: knitr
 Depends:
     R (>= 3.5.0)
-Date: 2024-07-09
-
+Date: 2024-07-18

NAMESPACE

@@ -37,7 +37,6 @@ importFrom(Matrix,Matrix)
 importFrom(Matrix,colSums)
 importFrom(Matrix,sparseMatrix)
 importFrom(R6,R6Class)
-importFrom(abind,abind)
 importFrom(dplyr,"%>%")
 importFrom(jsonlite,fromJSON)
 importFrom(memoise,forget)
@@ -46,6 +45,7 @@ importFrom(methods,as)
 importFrom(rlang,":=")
 importFrom(rlang,.data)
 importFrom(rlang,hash)
+importFrom(rnaturalearth,ne_countries)
 importFrom(tidyselect,matches)
 importFrom(units,as_units)
 importFrom(units,deparse_unit)

R/LPJmLDataCalc.R

@@ -2,7 +2,6 @@
 #'
 #' @importFrom units set_units as_units drop_units deparse_unit
 #' @import lpjmlkit
-#' @importFrom abind abind
 #' @importFrom Matrix sparseMatrix Matrix colSums
 #' @importFrom R6 R6Class
 #' @description
@@ -34,6 +33,14 @@ LPJmLDataCalc <- R6::R6Class( # nolint:object_linter_name
       private$.__aggregate__(ref_area, ...)
     },
 
+    #' @description
+    #' Add a band to the object by applying a function to the band vector for each spacial and temporal unit
+    #' @param band_name Name of band
+    #' @param fun function
+    add_band = function(band_name, fun) {
+      private$.__add_band__(band_name, fun)
+    },
+
     #' @description
     #' Get the reference area of the LPJmLDataCalc object.
     #' For an area density variable the reference area should
@@ -195,10 +202,6 @@ LPJmLDataCalc$set( # nolint: object_name_linter.
     #check internal consistency of grid or region data object
     ncells_meta <- private$.meta$ncell
     if (inherits(private$.grid, "LPJmLGridData")) {
-      if (!(private$.grid$meta$ncell == ncells_meta)) {
-        stop("Number of cells in grid meta data
-             is inconsistent withLPJmLDataCalc meta data")
-      }
       if (!(private$.grid$meta$ncell == dim(private$.grid$data)[1])) {
         stop("Number of cells in grid meta data
              is inconsistent grid data array")
@@ -325,6 +328,14 @@ LPJmLDataCalc$set(
   "private",
   ".__apply_operator__",
   function(sec_operand, operator) {
+    # check for matching bandnames
+    dimnames_to_match <- which(dim(sec_operand) > 1)
+    if (length(dimnames_to_match) > 0)
+      if (!identical(dimnames(sec_operand)[dimnames_to_match],
+                     dimnames(self$data)[dimnames_to_match]))
+        stop("Dimnames of second operand do not match first operator.")
+
+    # the dimensions of "self" should stay
     tar_dim <- dim(private$.data)
 
     # this function is used to expand the second operand
@@ -513,7 +524,7 @@ LPJmLDataCalc$set("private", ".initialize",  function(lpjml_data) {
   if (!inherits(lpjml_data, "LPJmLData")) {
     stop("Expected an LPJmLData object")
   }
-  if (!methods::is(lpjml_data$meta, "LPJmLMetaData")) {
+  if (!inherits(lpjml_data$meta, "LPJmLMetaData")) {
     stop("Meta data is missing")
   }
   if (!lpjml_data$meta$._space_format_ == "cell") {
@@ -585,13 +596,14 @@ LPJmLDataCalc$set("private", ".__plot_aggregated__", # nolint: object_name_linte
                       lapply(1:private$.meta$nbands, function(band) {
                         as.array(region_matrix %*% self$data[, , band])
                       })
-                    disaggr_data <- abind(list_of_disaggr_bands, along = 3)
+                    first_band <- list_of_disaggr_bands[[1]]
+                    disaggr_data <- unlist(list_of_disaggr_bands)
 
                     # recover dims and dimnames
                     dim(disaggr_data) <- c(
-                      cell = unname(dim(disaggr_data)[1]),
-                      time = unname(dim(disaggr_data)[2]),
-                      band = unname(dim(disaggr_data)[3])
+                      cell = unname(dim(first_band)[1]),
+                      time = unname(dim(first_band)[2]),
+                      band = private$.meta$nbands
                     )
                     dimnames(disaggr_data) <-
                       list(
@@ -691,6 +703,32 @@ subset.LPJmLDataCalc <- function(x, ...) {
   return(.as_LPJmLDataCalc(lpjml_dat))
 }
 
+LPJmLDataCalc$set("private", ".__add_band__",
+                  function(band_name, fun) {
+                    # create new larger array and copy content
+                    old_dimnames <- dimnames(private$.data)
+                    old_dims <- dim(private$.data)
+                    new_array <- array(NA, dim = c(old_dims[1], old_dims[2], old_dims[3] + 1))
+                    new_array[, , -(old_dims[3] + 1)] <- private$.data
+
+                    # insert new band
+                    new_array[, , old_dims[3] + 1] <- apply(private$.data, MARGIN = c(1, 2), FUN = fun)
+
+                    # update unit
+                    resulting_unit <- units(fun(private$.data[1, 1, ])) # test case to get unit
+                    new_array <- units::set_units(new_array, resulting_unit)
+
+                    # update dimnames
+                    old_dimnames[["band"]] <- c(old_dimnames[["band"]], band_name)
+                    dimnames(new_array) <- old_dimnames
+
+                    private$.data <- new_array
+
+                    # meta data update
+                    private$.meta$.__set_attribute__("nbands", private$.meta$nbands + 1)
+                    private$copy_unit_array2meta()
+                    private$.meta$.__set_attribute__("band_names", c(self$.meta$band_names, band_name))
+                  })
 
 
 # ----- add_grid -----

R/LPJmLDataCalc_aggregate.R

@@ -215,6 +215,7 @@ LPJmLDataCalc$set(
     # the grid is also needed to calculate the supporting cell areas
     # and construct dynamic regions
     self$add_grid()
+    private$.grid$subset(cell = dimnames(self$data)[["cell"]]) # subset the grid
 
     # select the correct LPJmLRegionData object as aggregation unit
     if (inherits(spatial_agg_units, "LPJmLRegionData")) {
@@ -375,13 +376,14 @@ LPJmLDataCalc$set(
       lapply(1:private$.meta$nbands, function(band) {
         as.array(region_matrix %*% self$data[, , band]) # core aggregation step
       })
-    aggr_data <- abind(list_of_aggr_bands, along = 3)
+    aggr_data <- unlist(list_of_aggr_bands)
+    first_band <- list_of_aggr_bands[[1]]
 
     # recover names of dimension vector
-    region_names <- dimnames(aggr_data)[[1]]
-    dim(aggr_data) <- c(region = dim(aggr_data)[1],
-                        time = dim(aggr_data)[2],
-                        band = dim(aggr_data)[3])
+    region_names <- dimnames(first_band)[[1]]
+    dim(aggr_data) <- c(region = dim(first_band)[1],
+                        time = dim(first_band)[2],
+                        band = private$.meta$nbands)
 
     # recover dimnames
     dimnames(aggr_data) <- list(
@@ -422,11 +424,12 @@ LPJmLDataCalc$set(
       stop("The ref_area must be either 'terr_area' or
            'cell_area'")
     }
+    cell_areas <- subset(cell_areas, cell = dimnames(self$data)[["cell"]])
     return(cell_areas)
   }
 )
 
-read_file <- function(searchdir, name, add_grid = TRUE) {
+read_file <- function(searchdir, name, add_grid = TRUE, ...) {
   # find path of file
   filename <- find_varfile(searchdir, name)
 
@@ -436,7 +439,7 @@ read_file <- function(searchdir, name, add_grid = TRUE) {
     " read from ",
     sQuote(basename(filename))
   ))
-  lpjml_calc <- read_io(filename)
+  lpjml_calc <- read_io(filename, ...)
 
   if (add_grid)
     lpjml_calc$add_grid()
@@ -461,7 +464,7 @@ calc_cellarea_wrapper <- function(lpjml_grid) {
   dim(cell_areas) <- c(cell = ncell,
                        time = 1,
                        band = 1)
-  dimnames(cell_areas) <- list(cell = 1:ncell,
+  dimnames(cell_areas) <- list(cell = dimnames(lpjml_grid$data)[["cell"]],
                                time = 1,
                                band = 1)
 

R/LPJmLRegionData.R

@@ -2,7 +2,6 @@
 #'
 #' @import lpjmlkit
 #' @importFrom Matrix sparseMatrix Matrix colSums
-#' @importFrom abind abind
 #' @importFrom utils read.csv
 #' @importFrom methods as
 #' @importFrom jsonlite fromJSON
@@ -160,7 +159,7 @@ read_countrymap <- function(file_path) {
 build_global_region <- function(grid) {
   region_matrix <- Matrix::Matrix(1, nrow = 1, ncol = grid$meta$ncell,
                                   sparse = TRUE)
-  region_matrix <- as(as(region_matrix, "generalMatrix"), "CsparseMatrix")
+  region_matrix <- methods::as(methods::as(region_matrix, "generalMatrix"), "CsparseMatrix")
   dimnames(region_matrix) <- list(c("global"), NULL)
   return(LPJmLRegionData$new(grid, region_matrix))
 }

R/MetricVarGrp.R

@@ -101,7 +101,7 @@ Metric <- R6::R6Class( # nolint: cyclocomp_linter object_linter_name
     #' @param var Variable name
     #' @param type Type of data ("baseline" or "under_test")
     capture_summary = function(lpjml_calc, var, type) {
-      if (!is.null(self$m_options$year_range)) {
+      if (!is.null(self$m_options$year_subset)) {
         subset_years <- function(lpjml_calc, years) {
           lpjml_calc %>%
             transform("year_month_day") %>%
@@ -110,7 +110,7 @@ Metric <- R6::R6Class( # nolint: cyclocomp_linter object_linter_name
         }
         lpjml_calc <-
           keep_units_lpjml_calc(lpjml_calc,
-                                function(x) subset_years(x, self$m_options$year_range))
+                                function(x) subset_years(x, self$m_options$year_subset))
       }
       summary <- self$summarize(lpjml_calc)
       self$store_summary(summary, var, type)
@@ -185,6 +185,7 @@ Metric <- R6::R6Class( # nolint: cyclocomp_linter object_linter_name
     generate_report_content = function() {
       self$print_metric_header()
       self$print_metric_description()
+      self$print_year_subset()
       plotlist <- self$plot()
       self$arrange_plots(plotlist)
     },
@@ -201,6 +202,15 @@ Metric <- R6::R6Class( # nolint: cyclocomp_linter object_linter_name
     #' Function to print the metric description.
     print_metric_description = function() {
       cat(self$description, "\n")
+    },
+
+    #' @description
+    #' !Package internal method!
+    #' Function to print the year_subset metric option.
+    print_year_subset = function() {
+      if (!is.null(self$m_options$year_subset)) {
+        pretty_print_year_subset(self$m_options$year_subset)
+      }
     }
   )
 )
@@ -214,28 +224,34 @@ VarGrp <- # nolint:object_linter_name
       # Function to retrieve the minimum and maximum of contained data
       # for the different types of data contained
       get_limits = function(type = "all", quantiles = c(0, 1)) {
-        get_quantiles_of_lpjml_calc_list <- function(list, quantiles) {
-          data <- sapply(c(unlist(list)), function(x) x$data) # nolint
-          lower_lim <- quantile(data, quantiles[1], na.rm = TRUE)
-          upper_lim <- quantile(data, quantiles[2], na.rm = TRUE)
+        get_quantiles_of_lpjml_calc <- function(lpjml_calc, quantiles) {
+          lower_lim <- Inf
+          upper_lim <- -Inf
+          for (band in dimnames(lpjml_calc$data)[["band"]]) {
+            data <- subset(lpjml_calc, band = band)$data
+            data[data == 0] <- NA # we assume that a zero means "no exisitng data" for the quantiles
+            if (all(is.na(data))) {
+              low <- 0
+              up <- 0
+            } else {
+              low <- unlist(quantile(data, quantiles[1], na.rm = TRUE))
+              up <- unlist(quantile(data, quantiles[2], na.rm = TRUE))
+            }
+            if (low < lower_lim) lower_lim <- low
+            if (up > upper_lim) upper_lim <- up
+          }
           return(c(lower_lim, upper_lim))
         }
-        if (type == "all") {
-          data <- list(self$under_test, self$compare, self$baseline)
-        } else {
-          data <- self[[type]]
-        }
-        limits <- get_quantiles_of_lpjml_calc_list(data, quantiles)
-        return(limits)
+        var_grp_lims <- unlist(self$apply_to_lpjml_calcs(get_quantiles_of_lpjml_calc, quantiles))
+        return(c(min(var_grp_lims), max(var_grp_lims)))
       },
 
       get_band_names = function() {
-        if (!is.null(self$baseline))
-          return(dimnames(self$baseline)[["band"]])
-        else if (!is.null(self$under_test[[1]]))
-          return(dimnames(self$under_test[[1]])[["band"]])
-        else if (!is.null(self$compare[[1]][[1]]))
-          return(dimnames(self$compare[[1]][[1]])[["band"]])
+        self$apply_to_any_lpjml_calc(function(x) dimnames(x$data)[["band"]])
+      },
+
+      get_var_name = function() {
+        self$apply_to_any_lpjml_calc(function(x) x$meta$variable)
       },
 
       # Function applies the function `fun`
@@ -342,3 +358,16 @@ VarGrp <- # nolint:object_linter_name
       compare = NULL
     )
   )
+
+# Utility functions
+pretty_print_year_subset <- function(years) {
+  numeric_years <- as.numeric(years)
+  vec_range <- range(numeric_years)
+  cat("Data subset: ")
+  if (identical(as.integer(numeric_years), as.integer(vec_range[1]:vec_range[2])))
+    cat(vec_range[1], "-", vec_range[2])
+  else
+    cat(numeric_years)
+  cat(" (years)")
+  cat("\n\n")
+}