Merge branch 'dev_getcellindex' into 'master'

Dev getcellindex See merge request lpjml/lpjmlkit!86
PIK-LPJmL · Mar 27, 2024 · 8341481 · 8341481
2 parents d1cdf36 + 459b26c
commit 8341481
Show file tree

Hide file tree

Showing 24 changed files with 888 additions and 320 deletions.
diff --git a/.buildlibrary b/.buildlibrary
@@ -1,4 +1,4 @@
-ValidationKey: '2634331'
+ValidationKey: '2773260'
 AutocreateReadme: yes
 AcceptedWarnings:
 - 'Warning: package ''.*'' was built under R version'

diff --git a/CITATION.cff b/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: 'lpjmlkit: Toolkit for Basic LPJmL Handling'
-version: 1.3.3
-date-released: '2024-03-25'
+version: 1.4.0
+date-released: '2024-03-27'
 abstract: A collection of basic functions to facilitate the work with the Dynamic
   Global Vegetation Model (DGVM) Lund-Potsdam-Jena managed Land (LPJmL) hosted at
   the Potsdam Institute for Climate Impact Research (PIK). It provides functions for

diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: lpjmlkit
 Type: Package
 Title: Toolkit for Basic LPJmL Handling
-Version: 1.3.3
+Version: 1.4.0
 Authors@R: c(
     person("Jannes", "Breier", , "[email protected]", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-9055-6904")),
     person("Sebastian","Ostberg", , "[email protected]", role = "aut", comment = c(ORCID = "0000-0002-2368-7015")),
@@ -54,4 +54,4 @@ Suggests:
     sf
 Config/testthat/edition: 3
 VignetteBuilder: knitr
-Date: 2024-03-25
+Date: 2024-03-27
diff --git a/NAMESPACE b/NAMESPACE
@@ -21,12 +21,14 @@ export(calc_cellarea)
 export(check_config)
 export(create_header)
 export(detect_io_type)
+export(get_cellindex)
 export(get_datatype)
 export(get_header_item)
 export(get_headersize)
 export(make_lpjml)
 export(plot.LPJmLData)
 export(read_config)
+export(read_grid)
 export(read_header)
 export(read_io)
 export(read_meta)

diff --git a/R/LPJmLData.R b/R/LPJmLData.R
@@ -254,7 +254,7 @@ LPJmLData <- R6::R6Class( # nolint:object_name_linter
       # Summary
       cat(col_var("$summary()\n"))
       print(self$summary(cutoff = TRUE))
-
+      # default handling of LPJmLData objects else inherited like LPJmLGridData
       if (class(self)[1] == "LPJmLData") {
         cat(
           col_note("Note: summary is not weighted by grid area.\n")

diff --git a/R/LPJmLData_as.R b/R/LPJmLData_as.R
@@ -69,11 +69,14 @@ as_array <- function(x,
 }
 
 # as_array method roxygen documentation in LPJmlData.R
-LPJmLData$set("private",
-              ".as_array",
-              function(subset = NULL,
-                       aggregate = NULL,
-                       ...) {
+LPJmLData$set(
+  "private",
+  ".as_array",
+  function(
+    subset = NULL,
+    aggregate = NULL,
+    ...
+  ) {
 
     # Initiate clone to be returned on which following methods are executed
     data_subset <- self$clone(deep = TRUE)
@@ -150,12 +153,15 @@ as_tibble.LPJmLData <- function(x,
 }
 
 # as_tibble method roxygen documentation in LPJmlData.R
-LPJmLData$set("private",
-              ".as_tibble",
-              function(subset = NULL,
-                       aggregate = NULL,
-                       value_name = "value",
-                       ...) {
+LPJmLData$set(
+  "private",
+  ".as_tibble",
+  function(
+    subset = NULL,
+    aggregate = NULL,
+    value_name = "value",
+    ...
+  ) {
 
     data <- self$as_array(subset, aggregate, ...)
 
@@ -236,11 +242,14 @@ as_raster <- function(x,
 }
 
 # as_raster method roxygen documentation in LPJmlData.R
-LPJmLData$set("private",
-              ".as_raster",
-              function(subset = NULL,
-                       aggregate = NULL,
-                       ...) {
+LPJmLData$set(
+  "private",
+  ".as_raster",
+  function(
+    subset = NULL,
+    aggregate = NULL,
+    ...
+  ) {
 
     data_subset <- private$.subset_raster_data(self, subset, aggregate, ...)
 
@@ -258,9 +267,20 @@ LPJmLData$set("private",
     # dimensions.
     if (data_subset$meta$._space_format_ == "lon_lat") {
 
-      data_subset$transform(to = "cell")
+      # default handling of LPJmLData objects else inherited like LPJmLGridData
+      if (class(data_subset)[1] == "LPJmLData") {
 
-      multi_dims <- get_multidims(data_subset)
+        data_subset$transform(to = "cell")
+
+        multi_dims <- get_multidims(data_subset)
+
+      } else {
+
+        tmp_raster <- raster::setValues(
+          tmp_raster,
+          t(data_subset$data)
+        )
+      }
     }
 
     # For space_format "cell" allow one additional dimension
@@ -297,13 +317,24 @@ LPJmLData$set("private",
       } else {
         tmp_data <- data_subset$data
       }
-      tmp_raster[
-        raster::cellFromXY(
-          tmp_raster,
-          cbind(subset_array(data_subset$grid$data, list(band = "lon")),
-                subset_array(data_subset$grid$data, list(band = "lat")))
-        )
-      ] <- tmp_data
+      # default handling of LPJmLData objects else inherited like LPJmLGridData
+      if (class(data_subset)[1] == "LPJmLData") {
+        tmp_raster[
+          raster::cellFromXY(
+            tmp_raster,
+            cbind(subset_array(data_subset$grid$data, list(band = "lon")),
+                  subset_array(data_subset$grid$data, list(band = "lat")))
+          )
+        ] <- tmp_data
+      } else {
+        tmp_raster[
+          raster::cellFromXY(
+            tmp_raster,
+            cbind(subset_array(data_subset$data, list(band = "lon")),
+                  subset_array(data_subset$data, list(band = "lat")))
+          )
+        ] <- tmp_data
+      }
     }
 
     return(tmp_raster)
@@ -356,22 +387,21 @@ LPJmLData$set("private",
 #' @md
 #' @aliases as_rast as_SpatRaster
 #' @export
-as_terra <- function(x,
-                     subset = NULL,
-                     aggregate = NULL,
-                      ...) {
-  y <- x$as_terra(subset,
-                  aggregate,
-                  ...)
+as_terra <- function(x, subset = NULL, aggregate = NULL, ...) {
+
+  y <- x$as_terra(subset, aggregate, ...)
   y
 }
 
 # as_terra method roxygen documentation in LPJmlData.R
-LPJmLData$set("private",
-              ".as_terra",
-              function(subset = NULL,
-                       aggregate = NULL,
-                       ...) {
+LPJmLData$set(
+  "private",
+  ".as_terra",
+  function(
+    subset = NULL,
+    aggregate = NULL,
+    ...
+  ) {
 
     data_subset <- private$.subset_raster_data(self, subset, aggregate, ...)
 
@@ -388,9 +418,20 @@ LPJmLData$set("private",
     # dimensions.
     if (data_subset$meta$._space_format_ == "lon_lat") {
 
-      data_subset$transform(to = "cell")
+      # default handling of LPJmLData objects else inherited like LPJmLGridData
+      if (class(data_subset)[1] == "LPJmLData") {
+
+        data_subset$transform(to = "cell")
+
+        multi_dims <- get_multidims(data_subset)
+
+      } else {
 
-      multi_dims <- get_multidims(data_subset)
+        tmp_rast <- terra::setValues(
+          tmp_rast,
+          t(data_subset$data)
+        )
+      }
     }
 
     # For space_format "cell" allow one additional dimension
@@ -410,7 +451,8 @@ LPJmLData$set("private",
         multi_layer <- multi_dims[which(multi_dims != "cell")]
 
         tmp_rast <- terra::rast(tmp_rast,
-                                nl = dim(data_subset$data)[multi_layer])
+                                nl = dim(data_subset$data)[multi_layer],
+                                vals = NA)
 
         names(tmp_rast) <- dimnames(data_subset$data)[[multi_layer]]
 
@@ -444,13 +486,24 @@ LPJmLData$set("private",
       } else {
         tmp_data <- data_subset$data
       }
-      tmp_rast[
-        terra::cellFromXY(
-          tmp_rast,
-          cbind(subset_array(data_subset$grid$data, list(band = "lon")),
-                subset_array(data_subset$grid$data, list(band = "lat")))
-        )
-      ] <- tmp_data
+      # default handling of LPJmLData objects else inherited like LPJmLGridData
+      if (class(data_subset)[1] == "LPJmLData") {
+        tmp_rast[
+          terra::cellFromXY(
+            tmp_rast,
+            cbind(subset_array(data_subset$grid$data, list(band = "lon")),
+                  subset_array(data_subset$grid$data, list(band = "lat")))
+          )
+        ] <- tmp_data
+      } else {
+        tmp_rast[
+          terra::cellFromXY(
+            tmp_rast,
+            cbind(subset_array(data_subset$data, list(band = "lon")),
+                  subset_array(data_subset$data, list(band = "lat")))
+          )
+        ] <- tmp_data
+      }
     }
 
     # Assign units (meta data)
@@ -461,16 +514,22 @@ LPJmLData$set("private",
 )
 
 
-LPJmLData$set("private",
-              ".subset_raster_data",
-              function(self,
-                       subset = NULL,
-                       aggregate = NULL,
-                       ...) {
+LPJmLData$set(
+  "private",
+  ".subset_raster_data",
+  function(
+    self,
+    subset = NULL,
+    aggregate = NULL,
+    ...
+  ) {
 
     # Support lazy loading of grid for meta files. This throws an error if no
     # suitable grid file is detected.
-    self$add_grid()
+    # default handling of LPJmLData objects else inherited like LPJmLGridData
+    if (class(self)[1] == "LPJmLData") {
+      self$add_grid()
+    }
 
     # Workflow adjusted for subsetted grid (via cell)
     data_subset <- self$clone(deep = TRUE)
@@ -506,9 +565,14 @@ create_tmp_raster <- function(data_subset, is_terra = FALSE) {
 
   # Calculate grid extent from range to span raster
   if (data_subset$meta$._space_format_ == "cell") {
-    data_extent <- rbind(min = apply(data_subset$grid$data, "band", min),
-                         max = apply(data_subset$grid$data, "band", max))
-
+    # default handling of LPJmLData objects else inherited like LPJmLGridData
+    if (class(data_subset)[1] == "LPJmLData") {
+      data_extent <- rbind(min = apply(data_subset$grid$data, "band", min),
+                           max = apply(data_subset$grid$data, "band", max))
+    } else {
+      data_extent <- rbind(min = apply(data_subset$data, "band", min),
+                           max = apply(data_subset$data, "band", max))
+    }
   } else {
     data_extent <- cbind(
       lon = range(as.numeric(dimnames(data_subset$data)[["lon"]])),
@@ -532,7 +596,8 @@ create_tmp_raster <- function(data_subset, is_terra = FALSE) {
       xmax = grid_extent[2, 1],
       ymin = grid_extent[1, 2],
       ymax = grid_extent[2, 2],
-      crs = "EPSG:4326"
+      crs = "EPSG:4326",
+      vals = NA
     )
 
   } else {

diff --git a/R/LPJmLData_subset.R b/R/LPJmLData_subset.R
@@ -11,13 +11,13 @@
 #'   `cell = c(27411:27416)`, `band = -c(14:16, 19:32)`, or character vectors if
 #'   the dimension has a dimnames attribute, e.g.
 #'   `band = c("rainfed rice", "rainfed maize")`.\
-#'   Coordinate pairs of individual cells can be selected by providing a tibble
-#'   in the form of `coords = tibble(lon = ..., lat =...)`. Coordinate values
-#'   in the tibble need to be supplied as character vectors. The argument can
-#'   also be called `coordinates`. When coordinates are supplied as character
-#'   vectors to subset either along the `lon` or `lat` dimension or to subset
-#'   by coordinate pair, the function matches the grid cells closest to the
-#'   supplied coordinate value.
+#'   Coordinate pairs of individual cells can be selected by providing a list or
+#'   tibble in the form of `coords = list(lon = ..., lat =...)`. Coordinate
+#'   values need to be supplied as character vectors. The argument
+#'   can also be called `coordinates`. When coordinates are supplied as
+#'   character vectors to subset either along the `lon` or `lat` dimension or to
+#'   subset by coordinate pair, the function matches the grid cells closest to
+#'   the supplied coordinate value.
 #'
 #' @return An [`LPJmLData`] object with dimensions resulting from the selection
 #'   in `subset`. Meta data are updated as well.
@@ -95,9 +95,10 @@ LPJmLData$set(
         subset_array_pair(x = self$data,
                           pair = subset_list[[coords]])
       )
-
-      # Subset grid with coordinates and update corresponding grid meta data
-      private$.grid$.__subset_space__(subset_list[coords])
+      if (!is.null(private$.grid)) {
+        # Subset grid with coordinates and update corresponding grid meta data
+        do.call(private$.grid$subset, subset_list[coords])
+      }
 
     } else {
       # Avoid errors when subsetting list with coords
@@ -122,10 +123,9 @@ LPJmLData$set(
                    subset_list[names(subset_list) != coords],
                    drop = FALSE)
     )
-
     # Subset grid with space dimensions and update corresponding grid meta data
     if (!is.null(private$.grid) && !is.null(subset_space_dim)) {
-      private$.grid$.__subset_space__(subset_list[subset_space_dim])
+      do.call(private$.grid$subset, subset_list[subset_space_dim])
     }
 
     if ("time" %in% names(subset_list)) {
@@ -153,11 +153,17 @@ LPJmLData$set(
 
       } else {
 
-        if (is.null(private$.grid)) {
+        if (is.null(private$.grid) && class(self)[1] == "LPJmLData") {
           stop("Missing $grid attribute. Add via $add_grid()")
         }
-        cell_dimnames <- sort(private$.grid$data) %>%
-          format(trim = TRUE, scientific = FALSE, justify = "none")
+        # default handling of LPJmLData objects else inherited like LPJmLGridData
+        if (class(self)[1] == "LPJmLData") {
+          cell_dimnames <- sort(private$.grid$data) %>%
+            format(trim = TRUE, scientific = FALSE, justify = "none")
+        } else {
+          cell_dimnames <- sort(self$data) %>%
+            format(trim = TRUE, scientific = FALSE, justify = "none")
+        }
       }
 
     } else {