Merge branch 'update_get_cellindex' into 'master'

Added shape file support and changed return value of get_cellindex

Closes #63

See merge request lpjml/lpjmlkit!100

.buildlibrary

@@ -1,4 +1,4 @@
-ValidationKey: '3437592'
+ValidationKey: '3458270'
 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: 'lpjmlkit: Toolkit for Basic LPJmL Handling'
-version: 1.7.2
-date-released: '2024-09-20'
+version: 1.7.3
+date-released: '2024-09-24'
 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

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: lpjmlkit
 Type: Package
 Title: Toolkit for Basic LPJmL Handling
-Version: 1.7.2
+Version: 1.7.3
 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")),
@@ -55,4 +55,4 @@ Suggests:
     sf
 Config/testthat/edition: 3
 VignetteBuilder: knitr
-Date: 2024-09-20
+Date: 2024-09-24

R/get_cellindex.R

@@ -12,10 +12,13 @@
 #' @param extent A numeric vector (lonmin, lonmax, latmin, latmax) containing the
 #'  longitude and latitude boundaries between which values included in the subset.
 #' @param coordinates A list of two named (lon, lat) numeric vectors representing the coordinates.
+#' @param shape A terra SpatVector object in the WGS 84 coordinate reference system
+#' representing the shape to subset the grid.
+#' @param simplify A logical indicating whether to simplify the output to a vector.
 #'
 #'
-#' @return The cell index from the grid file based on the provided extent or
-#' coordinates.
+#' @return Either an LPJmLData object containing the grid or a vector subsetted
+#' to the provided extent, coordinates or shape.
 #' @export
 #'
 #' @examples
@@ -43,28 +46,36 @@
 #' should be a list of two character vectors representing the longitude and
 #' latitude values as for [`subset()`].
 #'
-#' If both `extent` and `coordinates` are provided, the function will stop and
-#' ask for only one of them. If neither `extent` nor `coordinates` are provided,
-#' the function will return the cell numbers for all cells in the grid.
+#' If a shape is provided as a SpatVector object, the function will return the
+#' cell index for the cells that intersect with the shape.
+#'
+#' If more than on of `extent`, `coordinates` `shape` are provided, the function
+#' will stop and ask for only one of them. If neither `extent` nor `coordinates`
+#' nor `shape` are provided, the function will return the cell numbers for all
+#' cells in the grid.
 #'
 #' The function also includes checks for input types and values, and gives
 #' specific error messages for different error conditions. For example, it
 #' checks if the `grid_filename` exists, if the `extent` vector has the correct
 #' length, and if the `coordinates` list contains two vectors of equal length.
-get_cellindex <- function(grid_filename, extent = NULL, coordinates = NULL) {
+get_cellindex <- function(grid_filename, extent = NULL, coordinates = NULL, shape = NULL, simplify = TRUE) {
   # check if filepath is valid
   check_filepath(grid_filename)
   # check if (only) one of extent or coordinates is provided
-  check_extent_and_coordinates(extent, coordinates)
+  check_multiple(extent, coordinates, shape)
 
-  grid_lonlat <- read_io(filename = grid_filename) %>%
+  grid_lonlat <- read_io(filename = grid_filename) |>
     LPJmLGridData$new()
 
   if (!is.null(extent)) {
-    extent <- check_extent(extent) %>%
+    extent <- check_extent(extent) |>
       correct_extent()
   } else if (!is.null(coordinates)) {
     coordinates <- check_coordinates(coordinates)
+  } else if (!is.null(shape)) {
+    if (!("SpatVector" %in% class(shape))) {
+      stop("shape must be a SpatVector object.")
+    }
   }
 
   # Read the grid file and create a data frame
@@ -76,7 +87,7 @@ get_cellindex <- function(grid_filename, extent = NULL, coordinates = NULL) {
 
   # Check if extent values are within the longitude and latitude range in the cells
   if (!is.null(extent)) {
-    cells$cellindex <- as.numeric(row.names(cells)) + 1
+    cells$cellindex <- as.numeric(row.names(cells))
 
     out_of_bounds_lon <- extent[c(1, 2)][extent[c(1, 2)] < lon_range[1] |
                                            extent[c(1, 2)] > lon_range[2]]
@@ -94,7 +105,14 @@ get_cellindex <- function(grid_filename, extent = NULL, coordinates = NULL) {
 
     cells <- cells[cells$lon >= extent[1] &
                      cells$lon <= extent[2] &
-                     cells$lat >= extent[3] & cells$lat <= extent[4], ]$cellindex
+                     cells$lat >= extent[3] & cells$lat <= extent[4], ]
+
+    grid_cell <- transform(grid_lonlat, "lon_lat")
+
+    cells <- grid_cell$subset(coordinates = lapply(X = list(lon = cells$lon,
+                                                            lat = cells$lat),
+                                                   FUN = as.character))
+
   }
 
   # Check if coordinates are within the longitude and latitude range in the cells
@@ -120,9 +138,8 @@ get_cellindex <- function(grid_filename, extent = NULL, coordinates = NULL) {
 
     grid_cell <- transform(grid_lonlat, "lon_lat")
 
-    grid_cell$subset(coordinates = lapply(X = coordinates, FUN = as.character))
-
-    cells <- c(stats::na.omit(c(grid_cell$data + 1)))
+    cells <- grid_cell$subset(coordinates = lapply(X = coordinates,
+                                                   FUN = as.character))
 
     message(
       col_note(
@@ -131,6 +148,22 @@ get_cellindex <- function(grid_filename, extent = NULL, coordinates = NULL) {
     )
   }
 
+  if (!is.null(shape)) {
+    cell_coords <- grid_lonlat |>
+      as_terra() |>
+      terra::mask(shape) |>
+      terra::as.data.frame(xy = TRUE) |>
+      dplyr::select("x", "y")
+
+    cells <- grid_lonlat$transform("lon_lat") |>
+      subset(coordinates = lapply(list(lon = cell_coords$x, lat = cell_coords$y),
+                                  FUN = as.character))
+  }
+
+  if (simplify) {
+    cells <- c(stats::na.omit(c(cells$data + 1)))
+  }
+
   cells
 }
 
@@ -167,14 +200,16 @@ check_extent <- function(extent) {
 
 
 # Check if both extent and coordinates are provided
-check_extent_and_coordinates <- function(extent, coordinates) {
-  if (is.null(extent) && is.null(coordinates)) {
-    warning("Neither extent or coordinates provided. Full grid will be returned.")
+check_multiple <- function(extent, coordinates, shape) {
+  if (is.null(extent) && is.null(coordinates) && is.null(shape)) {
+    warning("Neither extent, coordinates or shape provided. Full grid will be returned.")
   }
-  if (!is.null(extent) && !is.null(coordinates)) {
+  if ((!is.null(extent) && !is.null(coordinates)) ||
+        (!is.null(extent) && !is.null(shape)) ||
+        (!is.null(coordinates) && !is.null(shape))) {
     stop(
-      "Both extent and coordinates are provided.",
-      " Please provide only one of them."
+      "Multiple subset options provided.",
+      " Please provide only one of coordinates, extent and shape."
     )
   }
 }

README.md

@@ -1,6 +1,6 @@
 # <a href=''><img src='inst/img/logo.png' align='right' alt='logo' height=139 /></a> Toolkit for Basic LPJmL Handling
 
-R package **lpjmlkit**, version **1.7.2**
+R package **lpjmlkit**, version **1.7.3**
 
 [![CRAN status](https://www.r-pkg.org/badges/version/lpjmlkit)](https://cran.r-project.org/package=lpjmlkit) [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.7773134.svg)](https://doi.org/10.5281/zenodo.7773134) [![R build status](https://github.com/PIK-LPJmL/lpjmlkit/workflows/check/badge.svg)](https://github.com/PIK-LPJmL/lpjmlkit/actions) [![codecov](https://codecov.io/gh/PIK-LPJmL/lpjmlkit/branch/master/graph/badge.svg)](https://app.codecov.io/gh/PIK-LPJmL/lpjmlkit) [![r-universe](https://pik-piam.r-universe.dev/badges/lpjmlkit)](https://pik-piam.r-universe.dev/builds)
 
@@ -76,7 +76,7 @@ In case of questions / problems please contact Jannes Breier <jannesbr@pik-potsd
 
 To cite package **lpjmlkit** in publications use:
 
-Breier J, Ostberg S, Wirth S, Minoli S, Stenzel F, Hötten D, Müller C (2024). _lpjmlkit: Toolkit for Basic LPJmL Handling_. doi:10.5281/zenodo.7773134 <https://doi.org/10.5281/zenodo.7773134>, R package version 1.7.2, <https://github.com/PIK-LPJmL/lpjmlkit>.
+Breier J, Ostberg S, Wirth S, Minoli S, Stenzel F, Hötten D, Müller C (2024). _lpjmlkit: Toolkit for Basic LPJmL Handling_. doi:10.5281/zenodo.7773134 <https://doi.org/10.5281/zenodo.7773134>, R package version 1.7.3, <https://github.com/PIK-LPJmL/lpjmlkit>.
 
 A BibTeX entry for LaTeX users is
 
@@ -85,7 +85,7 @@ A BibTeX entry for LaTeX users is
   title = {lpjmlkit: Toolkit for Basic LPJmL Handling},
   author = {Jannes Breier and Sebastian Ostberg and Stephen Björn Wirth and Sara Minoli and Fabian Stenzel and David Hötten and Christoph Müller},
   year = {2024},
-  note = {R package version 1.7.2},
+  note = {R package version 1.7.3},
   url = {https://github.com/PIK-LPJmL/lpjmlkit},
   doi = {10.5281/zenodo.7773134},
 }

tests/testthat/test-get_cellindex.R

@@ -26,7 +26,7 @@ test_that("get_cellindex handles both extent and coordinates provided", {
       extent = c(1.25, 2.75, 3.25, 4.75),
       coordinates = list(c(1.25, 2.75), c(1.25, 2.75))
     ),
-    "Both extent and coordinates are provided. Please provide only one of them."
+    "Multiple subset options provided. Please provide only"
   )
 })
 
@@ -62,3 +62,11 @@ test_that("get_cellindex returns correct cell index for given coordinates", {
   )
   expect_true(length(result) == 2 && result[1] == 10001 && result[2] == 10002)
 })
+
+test_that("get_cellindex returns correct cell index for a given shape", {
+  result <- get_cellindex("../testdata/output/grid.bin.json",
+                          shape = terra::vect(terra::ext(c(-87.25, -87.25,
+                                                           55.25, 55.75)),
+                                              crs = "+proj=longlat +datum=WGS84"))
+  expect_true(length(result) == 2 && result[1] == 10002 && result[2] == 10001)
+})