From 459b26c73c87d0e1c00e717037148f58477e66f1 Mon Sep 17 00:00:00 2001 From: Stephen Wirth Date: Wed, 27 Mar 2024 12:07:36 +0100 Subject: [PATCH] Dev getcellindex --- .buildlibrary | 2 +- CITATION.cff | 4 +- DESCRIPTION | 4 +- NAMESPACE | 2 + R/LPJmLData.R | 2 +- R/LPJmLData_as.R | 181 +++++++++++----- R/LPJmLData_subset.R | 36 ++-- R/LPJmLData_transform.R | 53 +++-- R/LPJmLGridData.R | 67 +----- R/LPJmLMetaData.R | 18 +- R/get_cellindex.R | 204 ++++++++++++++++++ R/read_grid.R | 20 ++ R/read_io.R | 4 +- R/subset_array.R | 3 +- README.md | 6 +- man/LPJmLGridData.Rd | 114 ++-------- man/get_cellindex.Rd | 64 ++++++ man/read_grid.Rd | 28 +++ man/subset.LPJmLData.Rd | 14 +- .../test-LPJmGridLData_subset_transform.R | 148 +++++++++++++ .../test-LPJmLData_subset_transform.R | 35 +-- tests/testthat/test-LPJmLGridData.R | 19 +- tests/testthat/test-LPJmLGridData_as.R | 116 ++++++++++ tests/testthat/test-get_cellindex.R | 64 ++++++ 24 files changed, 888 insertions(+), 320 deletions(-) create mode 100644 R/get_cellindex.R create mode 100644 R/read_grid.R create mode 100644 man/get_cellindex.Rd create mode 100644 man/read_grid.Rd create mode 100644 tests/testthat/test-LPJmGridLData_subset_transform.R create mode 100644 tests/testthat/test-LPJmLGridData_as.R create mode 100644 tests/testthat/test-get_cellindex.R diff --git a/.buildlibrary b/.buildlibrary index b845ea4..346238f 100644 --- 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 index b6f7112..1f45f0c 100644 --- 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 index 4559e0e..2eef9ae 100644 --- 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", , "jannesbr@pik-potsdam.de", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-9055-6904")), person("Sebastian","Ostberg", , "ostberg@pik-potsdam.de", 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 index 625c8c4..f99800e 100644 --- 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 index 695dc81..4e7c197 100644 --- 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 index 6045fd2..a5966a7 100644 --- 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 index 02ca7de..71744b9 100644 --- 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 { diff --git a/R/LPJmLData_transform.R b/R/LPJmLData_transform.R index 89a3a6f..ba88d9c 100644 --- a/R/LPJmLData_transform.R +++ b/R/LPJmLData_transform.R @@ -56,9 +56,10 @@ transform <- function(x, } # transform method roxygen documentation in LPJmlData.R -LPJmLData$set("private", - ".transform", - function(to) { +LPJmLData$set( + "private", + ".transform", + function(to) { # Transform time format if (any(to %in% private$.meta$._dimension_map_$time)) { @@ -104,9 +105,10 @@ LPJmLData$set("private", # Method to transform space dimension of data array from "cell" to "lon_lat" or # the other way around. If required add_grid to LPJmLData along the way. -LPJmLData$set("private", - ".transform_space", - function(to) { +LPJmLData$set( + "private", + ".transform_space", + function(to) { # If to equals current format return directly if (private$.meta$._space_format_ == to) { @@ -125,9 +127,9 @@ LPJmLData$set("private", # Case 1: Transformation from cell dimension to lon, lat dimensions if (private$.meta$._space_format_ == "cell" && - to == "lon_lat") { + to == "lon_lat") { - private$.grid$.__transform_space__(to = to) + private$.grid$transform(to = to) # Matrix with ilon and ilat indices of cells in new array ilonilat <- arrayInd( @@ -186,9 +188,8 @@ LPJmLData$set("private", # Insert data from source array into target array target_array[index_target] <- self$data - # Case 2: Transformation between lon, lat dimensions and cell dimension - } else if (private$.meta$._space_format_ == "lon_lat" && - to == "cell") { + # Case 2: Transformation between lon, lat dimensions and cell dimension + } else if (private$.meta$._space_format_ == "lon_lat" && to == "cell") { # Matrix with ilon and ilat indices of cells in new array ilonilat <- arrayInd(match(sort(private$.grid$data), private$.grid$data), @@ -199,7 +200,7 @@ LPJmLData$set("private", band = c("lon", "lat")) # Transform grid to target space format - private$.grid$.__transform_space__(to = to) + private$.grid$transform(to = to) new_dimnames <- append( other_dimnames, @@ -230,13 +231,13 @@ LPJmLData$set("private", MoreArgs = list(ilonilat = ilonilat) ) %>% unlist(recursive = TRUE, use.names = FALSE) %>% - matrix(nrow = nrow(index_source)) + matrix(nrow = nrow(index_source)) rm(index_source, ilonilat) gc(full = TRUE) target_array <- array(self$data[index_target], dim = new_dims, - dimnames = new_dimnames) + dimnames = new_dimnames) } else { return(invisible(self)) @@ -255,20 +256,20 @@ LPJmLData$set("private", # Method to transform the time dimension of the data array from "time" to # "year_month_day" or the other way around -LPJmLData$set("private", - ".transform_time", - function(to) { +LPJmLData$set( + "private", + ".transform_time", + function(to) { - # If to equals current format return directly + # If to equals current format return directly if (private$.meta$._time_format_ == to) { return(invisible(self)) } # Case 1: Transformation from "time" dimension to "year", "month", "day" # dimensions (if available) - if (private$.meta$._time_format_ == "time" && - to == "year_month_day") { + if (private$.meta$._time_format_ == "time" && to == "year_month_day") { # Possible ndays of months ndays_in_month <- c(31, 30, 28) @@ -283,17 +284,15 @@ LPJmLData$set("private", } # Remove month dimension if there is only one month in data -> annual - if (length(time_dimnames$month) == 1 && - is.null(time_dimnames[["day"]])) { + if (length(time_dimnames$month) == 1 && is.null(time_dimnames[["day"]])) { time_dimnames[["month"]] <- NULL } new_time_dim <- "year_month_day" - # Case 2: Transformation from dimensions "year", "month", "day" - # (if available) to "time" dimension - } else if (private$.meta$._time_format_ == "year_month_day" && - to == "time") { + # Case 2: Transformation from dimensions "year", "month", "day" + # (if available) to "time" dimension + } else if (private$.meta$._time_format_ == "year_month_day" && to == "time") { # nolint:line_length_linter # Convert time dimnames back to time pre_dimnames <- self$dimnames() %>% @@ -309,7 +308,7 @@ LPJmLData$set("private", new_time_dim <- "time" - # Return directly if no transformation is necessary + # Return directly if no transformation is necessary } else { return(invisible(self)) } diff --git a/R/LPJmLGridData.R b/R/LPJmLGridData.R index 30c0eeb..fa8060f 100644 --- a/R/LPJmLGridData.R +++ b/R/LPJmLGridData.R @@ -32,70 +32,20 @@ LPJmLGridData <- R6::R6Class( # nolint:object_name_linter stop("Not allowed for an object of class LPJmLGridData") }, - - #' @description - #' ! Not allowed to use dimension names of `LPJmLGridData$data` - #' array directly to subset each dimension to match the supplied vectors. - #' - #' @param ... See [`subset.LPJmLData()`] - subset = function(...) { - - stop("Not allowed for an object of class LPJmLGridData") - }, - - - #' @description - #' ! Not allowed to transform inner `LPJmLGridData$data` array - #' into another space or time format. - #' - #' @param ... See [`transform()`]. - transform = function(...) { - - stop("Not allowed for an object of class LPJmLGridData") - }, - - # export methods --------------------------------------------------------- # - - #' @description - #' ! Not allowed to coerce (convert) an `LPJmLGridData` object into a - #' \link[raster]{raster} or \link[raster]{brick} object. + #' ! No plot function available for `LPJmLGridData` object. Use + #' [`as_raster()`] or [`as_terra()`] (and `plot()`) to visualize the grid. #' - #' @param ... See [`as_raster()`]. - as_raster = function(...) { + #' @param ... See [`plot()`]. + plot = function(...) { - stop("Not allowed for an object of class LPJmLGridData") - }, - - - #' @description - #' ! Not allowed to coerce (convert) an `LPJmLGridData` object into a - #' \link[terra]{rast} object. - #' - #' @param ... See [`as_terra()`]. - as_terra = function(...) { - - stop("Not allowed for an object of class LPJmLGridData") - }, - - - #' @description - #' ! Internal method only to be used for package development. - #' - #' @param ... See [`subset()`]. - .__subset_space__ = function(...) { - private$.subset_space(...) + stop( + "No plot function available for LPJmLGridData object. Use ", + "as_raster() or as_terra() (and plot()) to visualize grid data." + ) }, - #' @description - #' ! Internal method only to be used for package development. - #' - #' @param ... See [`transform()`]. - .__transform_space__ = function(...) { - private$.transform_space(...) - }, - # Create a new LPJmLGridData object only to be used internally or explicitly #' @description @@ -113,7 +63,6 @@ LPJmLGridData <- R6::R6Class( # nolint:object_name_linter # Assign LPJmLData data private$.data <- x$data - if (!is.null(private$.meta$variable)) { if (private$.meta$variable %in% c("grid", "cellid")) { diff --git a/R/LPJmLMetaData.R b/R/LPJmLMetaData.R index 6ad277b..5c03d10 100644 --- a/R/LPJmLMetaData.R +++ b/R/LPJmLMetaData.R @@ -144,12 +144,18 @@ LPJmLMetaData <- R6::R6Class( # nolint # Update fields_set private$.fields_set <- private$.fields_set[ - -na.omit(match(c("nyear", - "firstyear", - "lastyear", - "nstep", - "timestep"), - private$.fields_set)) + -stats::na.omit( + match( + c( + "nyear", + "firstyear", + "lastyear", + "nstep", + "timestep" + ), + private$.fields_set + ) + ) ] }, diff --git a/R/get_cellindex.R b/R/get_cellindex.R new file mode 100644 index 0000000..3d88639 --- /dev/null +++ b/R/get_cellindex.R @@ -0,0 +1,204 @@ +#' @title Get Cell Index +#' +#' @description This function returns the cell index from a grid file based on +#' the provided extent or coordinates. If neither extent nor coordinates are +#' provided, the full grid will be returned. If both extent and coordinates are +#' provided, the function will stop and ask for only one of them. The extent +#' should be a vector of length 4 in the form c(lonmin, lonmax, latmin, latmax). +#' If the extent is not in the correct form, the function will swap the values +#' to correct it. +#' +#' @param grid_filename A string representing the grid file name. +#' @param extent A numeric vector of length 4 representing the extent +#' (lonmin, lonmax, latmin, latmax). +#' @param coordinates A list of two numeric vectors representing the coordinates. +#' +#' +#' @return The cell index from the grid file based on the provided extent or +#' coordinates. +#' @export +#' +#' @examples +#' \dontrun{ +#' get_cellindex( +#' grid_filename = "my_grid.bin.json", +#' extent = c(-123.25, -122.75, 49.25, 49.75) +#' ) +#' get_cellindex( +#' grid_filename = "my_grid.bin.json", +#' coordinates = list(c(-123.25, -122.75), c(49.25, 49.75)) +#' ) +#' } +#' @details +#' The function reads a grid file specified by `grid_filename` and creates a +#' data frame with columns for longitude, latitude, and cell number. The cell +#' number is a sequence from 1 to the number of rows in the data frame. +#' +#' If an `extent` is provided, the function filters the cells to include only +#' those within the specified longitude and latitude range. The `extent` should +#' be a numeric vector of length 4 in the form c(lonmin, lonmax, latmin, latmax). +#' +#' If a list of `coordinates` is provided, the function filters the cells to +#' include only those that match the specified coordinates. The `coordinates` +#' 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. +#' +#' 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) { + # 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) + + grid_lonlat <- read_io(filename = grid_filename) %>% + LPJmLGridData$new() + + if (!is.null(extent)) { + extent <- check_extent(extent) %>% + correct_extent() + } else if (!is.null(coordinates)) { + check_coordinates_length(coordinates) + } + + # Read the grid file and create a data frame + cells <- as.data.frame(grid_lonlat$data) + + # Get the range of longitude and latitude in the cells + lon_range <- range(cells$lon) + lat_range <- range(cells$lat) + + # 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 + + out_of_bounds_lon <- extent[c(1, 2)][extent[c(1, 2)] < lon_range[1] | + extent[c(1, 2)] > lon_range[2]] + out_of_bounds_lat <- extent[c(3, 4)][extent[c(3, 4)] < lat_range[1] | + extent[c(3, 4)] > lat_range[2]] + + if (length(out_of_bounds_lon) > 0 || length(out_of_bounds_lat) > 0) { + warning(paste( + "Extent values out of bounds: longitude", + paste(out_of_bounds_lon, collapse = ", "), + "latitude", + paste(out_of_bounds_lat, collapse = ", ") + )) + } + + cells <- cells[cells$lon >= extent[1] & + cells$lon <= extent[2] & + cells$lat >= extent[3] & cells$lat <= extent[4], ]$cellindex + } + + # Check if coordinates are within the longitude and latitude range in the cells + if (!is.null(coordinates)) { + out_of_bounds <- coordinates[[1]] < lon_range[1] | + coordinates[[1]] > lon_range[2] | + coordinates[[2]] < lat_range[1] | + coordinates[[2]] > lat_range[2] + + out_of_bounds_coords <- paste0( + "(", + coordinates[[1]][out_of_bounds], + ", ", + coordinates[[2]][out_of_bounds], + ")" + ) + if (any(out_of_bounds)) { + warning(paste( + "Coordinates out of bounds:", + paste(out_of_bounds_coords, collapse = ", ") + )) + } + + grid_cell <- transform(grid_lonlat, "lon_lat") + + grid_cell$subset(coordinates = coordinates) + + cells <- c(stats::na.omit(c(grid_cell$data + 1))) + + message( + col_note( + "Note: Possible duplicates of cell ids are removed." + ) + ) + } + + cells +} + +# Check if the input is a valid file path +check_filepath <- function(grid_filename) { + if (!is.character(grid_filename) || grid_filename == "") { + stop("grid_filename must be a string representing a valid file path.") + } + if (!file.exists(grid_filename)) { + stop("grid_filename does not exist.") + } +} + +# Check if the extent is a numeric vector of length 4 +check_extent <- function(extent) { + + if (!is.null(extent) && (length(extent) != 4)) { + stop("extent must be a numeric vector of length 4.") + } + + # check if character vector present and convert to numeric (to be consistent with coordinates) # nolint + if (is.character(extent)) { + extent <- as.numeric(extent) + } + + # Check if any NA's are present in extent + if (any(is.na(extent))) { + stop("NA's (induced by conversion) are not allowed in extent.") + } + + extent +} + + +# Check if the coordinates are a list of two numeric vectors of equal length + +# 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.") + } + if (!is.null(extent) && !is.null(coordinates)) { + stop( + "Both extent and coordinates are provided.", + " Please provide only one of them." + ) + } +} + +# Check and correct the extent if necessary +correct_extent <- function(extent) { + if (extent[1] > extent[2] || extent[3] > extent[4]) { + warning("Extent should be in the form c(lonmin, lonmax, latmin, latmax).") + if (extent[1] > extent[2]) { + extent[1:2] <- sort(extent[1:2]) + warning("Swapped values of lonmin and lonmax.") + } + if (extent[3] > extent[4]) { + extent[3:4] <- sort(extent[3:4]) + warning("Swapped values of latmin and latmax.") + } + } + return(extent) +} + +# Check the length of coordinates +check_coordinates_length <- function(coordinates) { + if (length(coordinates[[1]]) != length(coordinates[[2]])) { + stop("The two vectors in coordinates must have the same length.") + } +} diff --git a/R/read_grid.R b/R/read_grid.R new file mode 100644 index 0000000..0de5685 --- /dev/null +++ b/R/read_grid.R @@ -0,0 +1,20 @@ +#' @title Read LPJmL input and output grid files +#' +#' @description Generic function to read LPJmL input & output files in different +#' formats. Depending on the format, arguments can be automatically detected +#' or have to be passed as individual arguments. +#' +#' @param ... See [read_io] for further arguments. +#' @return An [LPJmLGridData] object. +#' @examples +#' \dontrun{ +#' my_grid <- read_io("grid.bin.json") +#' +#' } +#' @details +#' See [read_io] for more details. +#' @export +read_grid <- function(...) { + read_io(...) %>% + LPJmLGridData$new() +} diff --git a/R/read_io.R b/R/read_io.R index 039de49..57293d9 100644 --- a/R/read_io.R +++ b/R/read_io.R @@ -204,7 +204,7 @@ read_io <- function( # nolint:cyclocomp_linter. valid_dim_names <- c("cell", "time", "band") if (!all(dim_order %in% valid_dim_names)) { stop( - "Invalid dim_order prodided: c(", + "Invalid dim_order provided: c(", toString(sQuote(dim_order)), ")\n", "dim_order can be in any order but must include all of the following ", "band names: c(", toString(sQuote(sort(valid_dim_names))), ")" @@ -212,7 +212,7 @@ read_io <- function( # nolint:cyclocomp_linter. } if (!all(valid_dim_names %in% dim_order)) { stop( - "Invalid dim_order prodided: c(", + "Invalid dim_order provided: c(", toString(sQuote(dim_order)), ")\n", "dim_order can be in any order but must include all of the following ", "band names: c(", toString(sQuote(sort(valid_dim_names))), ")" diff --git a/R/subset_array.R b/R/subset_array.R index 4721faf..21b2245 100644 --- a/R/subset_array.R +++ b/R/subset_array.R @@ -208,7 +208,8 @@ subset_array_pair <- function(x, pair_names <- names(pair) # Ensure that coordinates are character vectors. - if (!all(sapply(pair, is.character) | sapply(pair, is.factor), na.rm = TRUE)) { # nolint + if (length(pair) == 0 || + !all(sapply(pair, is.character) | sapply(pair, is.factor), na.rm = TRUE)) { # nolint stop("Values for coordinate pairs must be supplied as strings") } diff --git a/README.md b/README.md index b16a223..397259b 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ # Toolkit for Basic LPJmL Handling -R package **lpjmlkit**, version **1.3.3** +R package **lpjmlkit**, version **1.4.0** [![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 . +Breier J, Ostberg S, Wirth S, Minoli S, Stenzel F, Müller C (2024). _lpjmlkit: Toolkit for Basic LPJmL Handling_. doi: 10.5281/zenodo.7773134 (URL: https://doi.org/10.5281/zenodo.7773134), R package version 1.4.0, . 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 Christoph Müller}, year = {2024}, - note = {R package version 1.3.3}, + note = {R package version 1.4.0}, doi = {10.5281/zenodo.7773134}, url = {https://github.com/PIK-LPJmL/lpjmlkit}, } diff --git a/man/LPJmLGridData.Rd b/man/LPJmLGridData.Rd index dd3bb88..98b354b 100644 --- a/man/LPJmLGridData.Rd +++ b/man/LPJmLGridData.Rd @@ -22,12 +22,7 @@ the meta data via \verb{$meta}. \subsection{Public methods}{ \itemize{ \item \href{#method-LPJmLGridData-add_grid}{\code{LPJmLGridData$add_grid()}} -\item \href{#method-LPJmLGridData-subset}{\code{LPJmLGridData$subset()}} -\item \href{#method-LPJmLGridData-transform}{\code{LPJmLGridData$transform()}} -\item \href{#method-LPJmLGridData-as_raster}{\code{LPJmLGridData$as_raster()}} -\item \href{#method-LPJmLGridData-as_terra}{\code{LPJmLGridData$as_terra()}} -\item \href{#method-LPJmLGridData-.__subset_space__}{\code{LPJmLGridData$.__subset_space__()}} -\item \href{#method-LPJmLGridData-.__transform_space__}{\code{LPJmLGridData$.__transform_space__()}} +\item \href{#method-LPJmLGridData-plot}{\code{LPJmLGridData$plot()}} \item \href{#method-LPJmLGridData-new}{\code{LPJmLGridData$new()}} \item \href{#method-LPJmLGridData-print}{\code{LPJmLGridData$print()}} \item \href{#method-LPJmLGridData-clone}{\code{LPJmLGridData$clone()}} @@ -39,12 +34,15 @@ the meta data via \verb{$meta}.
  • lpjmlkit::LPJmLData$.__set_data__()
  • lpjmlkit::LPJmLData$.__set_grid__()
  • lpjmlkit::LPJmLData$as_array()
    • +
  • lpjmlkit::LPJmLData$as_raster()
    • +
  • lpjmlkit::LPJmLData$as_terra()
  • lpjmlkit::LPJmLData$as_tibble()
  • lpjmlkit::LPJmLData$dim()
  • lpjmlkit::LPJmLData$dimnames()
  • lpjmlkit::LPJmLData$length()
    • -
  • lpjmlkit::LPJmLData$plot()
    • +
  • lpjmlkit::LPJmLData$subset()
  • lpjmlkit::LPJmLData$summary()
    • +
  • lpjmlkit::LPJmLData$transform()
    • }} @@ -66,107 +64,19 @@ the meta data via \verb{$meta}. } } \if{html}{\out{
    }} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-LPJmLGridData-subset}{}}} -\subsection{Method \code{subset()}}{ -! Not allowed to use dimension names of \code{LPJmLGridData$data} -array directly to subset each dimension to match the supplied vectors. +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-LPJmLGridData-plot}{}}} +\subsection{Method \code{plot()}}{ +! No plot function available for \code{LPJmLGridData} object. Use +\code{\link[=as_raster]{as_raster()}} or \code{\link[=as_terra]{as_terra()}} (and \code{plot()}) to visualize the grid. \subsection{Usage}{ -\if{html}{\out{
    }}\preformatted{LPJmLGridData$subset(...)}\if{html}{\out{
    }} +\if{html}{\out{
    }}\preformatted{LPJmLGridData$plot(...)}\if{html}{\out{
    }} } \subsection{Arguments}{ \if{html}{\out{
    }} \describe{ -\item{\code{...}}{See \code{\link[=subset.LPJmLData]{subset.LPJmLData()}}} -} -\if{html}{\out{
    }} -} -} -\if{html}{\out{
    }} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-LPJmLGridData-transform}{}}} -\subsection{Method \code{transform()}}{ -! Not allowed to transform inner \code{LPJmLGridData$data} array -into another space or time format. -\subsection{Usage}{ -\if{html}{\out{
    }}\preformatted{LPJmLGridData$transform(...)}\if{html}{\out{
    }} -} - -\subsection{Arguments}{ -\if{html}{\out{
    }} -\describe{ -\item{\code{...}}{See \code{\link[=transform]{transform()}}.} -} -\if{html}{\out{
    }} -} -} -\if{html}{\out{
    }} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-LPJmLGridData-as_raster}{}}} -\subsection{Method \code{as_raster()}}{ -! Not allowed to coerce (convert) an \code{LPJmLGridData} object into a -\link[raster]{raster} or \link[raster]{brick} object. -\subsection{Usage}{ -\if{html}{\out{
    }}\preformatted{LPJmLGridData$as_raster(...)}\if{html}{\out{
    }} -} - -\subsection{Arguments}{ -\if{html}{\out{
    }} -\describe{ -\item{\code{...}}{See \code{\link[=as_raster]{as_raster()}}.} -} -\if{html}{\out{
    }} -} -} -\if{html}{\out{
    }} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-LPJmLGridData-as_terra}{}}} -\subsection{Method \code{as_terra()}}{ -! Not allowed to coerce (convert) an \code{LPJmLGridData} object into a -\link[terra]{rast} object. -\subsection{Usage}{ -\if{html}{\out{
    }}\preformatted{LPJmLGridData$as_terra(...)}\if{html}{\out{
    }} -} - -\subsection{Arguments}{ -\if{html}{\out{
    }} -\describe{ -\item{\code{...}}{See \code{\link[=as_terra]{as_terra()}}.} -} -\if{html}{\out{
    }} -} -} -\if{html}{\out{
    }} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-LPJmLGridData-.__subset_space__}{}}} -\subsection{Method \code{.__subset_space__()}}{ -! Internal method only to be used for package development. -\subsection{Usage}{ -\if{html}{\out{
    }}\preformatted{LPJmLGridData$.__subset_space__(...)}\if{html}{\out{
    }} -} - -\subsection{Arguments}{ -\if{html}{\out{
    }} -\describe{ -\item{\code{...}}{See \code{\link[=subset]{subset()}}.} -} -\if{html}{\out{
    }} -} -} -\if{html}{\out{
    }} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-LPJmLGridData-.__transform_space__}{}}} -\subsection{Method \code{.__transform_space__()}}{ -! Internal method only to be used for package development. -\subsection{Usage}{ -\if{html}{\out{
    }}\preformatted{LPJmLGridData$.__transform_space__(...)}\if{html}{\out{
    }} -} - -\subsection{Arguments}{ -\if{html}{\out{
    }} -\describe{ -\item{\code{...}}{See \code{\link[=transform]{transform()}}.} +\item{\code{...}}{See \code{\link[=plot]{plot()}}.} } \if{html}{\out{
    }} } diff --git a/man/get_cellindex.Rd b/man/get_cellindex.Rd new file mode 100644 index 0000000..a2dec19 --- /dev/null +++ b/man/get_cellindex.Rd @@ -0,0 +1,64 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/get_cellindex.R +\name{get_cellindex} +\alias{get_cellindex} +\title{Get Cell Index} +\usage{ +get_cellindex(grid_filename, extent = NULL, coordinates = NULL) +} +\arguments{ +\item{grid_filename}{A string representing the grid file name.} + +\item{extent}{A numeric vector of length 4 representing the extent +(lonmin, lonmax, latmin, latmax).} + +\item{coordinates}{A list of two numeric vectors representing the coordinates.} +} +\value{ +The cell index from the grid file based on the provided extent or +coordinates. +} +\description{ +This function returns the cell index from a grid file based on +the provided extent or coordinates. If neither extent nor coordinates are +provided, the full grid will be returned. If both extent and coordinates are +provided, the function will stop and ask for only one of them. The extent +should be a vector of length 4 in the form c(lonmin, lonmax, latmin, latmax). +If the extent is not in the correct form, the function will swap the values +to correct it. +} +\details{ +The function reads a grid file specified by \code{grid_filename} and creates a +data frame with columns for longitude, latitude, and cell number. The cell +number is a sequence from 1 to the number of rows in the data frame. + +If an \code{extent} is provided, the function filters the cells to include only +those within the specified longitude and latitude range. The \code{extent} should +be a numeric vector of length 4 in the form c(lonmin, lonmax, latmin, latmax). + +If a list of \code{coordinates} is provided, the function filters the cells to +include only those that match the specified coordinates. The \code{coordinates} +should be a list of two character vectors representing the longitude and +latitude values as for \code{\link[=subset]{subset()}}. + +If both \code{extent} and \code{coordinates} are provided, the function will stop and +ask for only one of them. If neither \code{extent} nor \code{coordinates} 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 \code{grid_filename} exists, if the \code{extent} vector has the correct +length, and if the \code{coordinates} list contains two vectors of equal length. +} +\examples{ +\dontrun{ +get_cellindex( + grid_filename = "my_grid.bin.json", + extent = c(-123.25, -122.75, 49.25, 49.75) +) +get_cellindex( + grid_filename = "my_grid.bin.json", + coordinates = list(c(-123.25, -122.75), c(49.25, 49.75)) +) +} +} diff --git a/man/read_grid.Rd b/man/read_grid.Rd new file mode 100644 index 0000000..1467a29 --- /dev/null +++ b/man/read_grid.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/read_grid.R +\name{read_grid} +\alias{read_grid} +\title{Read LPJmL input and output grid files} +\usage{ +read_grid(...) +} +\arguments{ +\item{...}{See \link{read_io} for further arguments.} +} +\value{ +An \link{LPJmLGridData} object. +} +\description{ +Generic function to read LPJmL input & output files in different +formats. Depending on the format, arguments can be automatically detected +or have to be passed as individual arguments. +} +\details{ +See \link{read_io} for more details. +} +\examples{ +\dontrun{ +my_grid <- read_io("grid.bin.json") + +} +} diff --git a/man/subset.LPJmLData.Rd b/man/subset.LPJmLData.Rd index a98d2fc..7ac70a6 100644 --- a/man/subset.LPJmLData.Rd +++ b/man/subset.LPJmLData.Rd @@ -15,13 +15,13 @@ dimensions. Subsets may either specify integer indices, e.g. \code{cell = c(27411:27416)}, \code{band = -c(14:16, 19:32)}, or character vectors if the dimension has a dimnames attribute, e.g. \code{band = c("rainfed rice", "rainfed maize")}.\ -Coordinate pairs of individual cells can be selected by providing a tibble -in the form of \code{coords = tibble(lon = ..., lat =...)}. Coordinate values -in the tibble need to be supplied as character vectors. The argument can -also be called \code{coordinates}. When coordinates are supplied as character -vectors to subset either along the \code{lon} or \code{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 \code{coords = list(lon = ..., lat =...)}. Coordinate +values need to be supplied as character vectors. The argument +can also be called \code{coordinates}. When coordinates are supplied as +character vectors to subset either along the \code{lon} or \code{lat} dimension or to +subset by coordinate pair, the function matches the grid cells closest to +the supplied coordinate value.} } \value{ An \code{\link{LPJmLData}} object with dimensions resulting from the selection diff --git a/tests/testthat/test-LPJmGridLData_subset_transform.R b/tests/testthat/test-LPJmGridLData_subset_transform.R new file mode 100644 index 0000000..3fe2244 --- /dev/null +++ b/tests/testthat/test-LPJmGridLData_subset_transform.R @@ -0,0 +1,148 @@ +# Utility function to test data integrity +# tests designed for data to still have sequential order (c(1,2,3) NOT c(1,3)) +# latter would not work with following simplified tests +test_integrity <- function(output) { + + # Do call dim and dimnames only once + dim_data <- dim(output$data) + dimnames_data <- dimnames(output$data) + + # Check for two cases "cell" or "lon_lat" + if ("cell" %in% names(dim_data)) { + # Test for equal length of cell in data and meta data (ncell) + testthat::expect_equal(dim_data[["cell"]], output$meta$ncell) + # Test for equal dimnames of cell in data and those constructed by meta data + testthat::expect_equal( + dimnames_data$cell, + format( + seq(output$meta$firstcell, length.out = output$meta$ncell), + trim = TRUE, scientific = FALSE + ) + ) + } + + if (output$meta$._space_format_ == "cell") { + # Test for equal length of bands in data and meta data (nbands) + testthat::expect_equal(dim_data[["band"]], output$meta$nbands) + } +} + + +# test subset method +test_that("test subset method", { + file_name <- "../testdata/output/grid.bin.json" + output <- read_grid(filename = file_name) + + output_sub <- subset( + output, + cell = 1:2 + ) + test_integrity(output_sub) +}) + + +# test transform_time method +test_that("test transform (space) method", { + file_name <- "../testdata/output/grid.bin.json" + output <- read_grid(filename = file_name) + + test_integrity(output) + + output2 <- read_grid( + filename = file_name, + dim_order = c("band", "cell", "time") + ) + + output2$transform(to = "lon_lat") + + test_integrity(output2) +}) + + +# test transform and subset method +test_that("test transform (space) method", { + file_name <- "../testdata/output/grid.bin.json" + output <- read_grid(filename = file_name) + + # Explicitly load grid + output$transform(to = c("lon_lat")) + output$subset(lat = c("55.25", "55.75", "56.25", "56.75")) + output$transform(to = "cell") + test_integrity(output) + + output2 <- read_grid( + filename = file_name, + dim_order = c("band", "cell", "time") + ) + + output2$transform(to = c("lon_lat")) + output2$subset(lat = c("55.25", "55.75", "56.25", "56.75")) + output2$transform(to = "cell") + test_integrity(output2) + expect_identical( + output2$data, + aperm(output$data, names(dim(output2))) + ) +}) + + +# coordinates located within the last two cells +coordinates <- list( + lat = c("55.9", "63.7"), + lon = c("-87.3", "-87.1") +) +coordinates_numeric <- list( + lon = as.numeric(coordinates$lon), + lat = as.numeric(coordinates$lat) +) + +# test subset method for coordinates (pair) +test_that("test transform (space) method", { + + file_name <- "../testdata/output/grid.bin.json" + output <- read_grid(filename = file_name) + + output_trans <- transform(output, to = c("lon_lat")) + output_trans$subset(coordinates = coordinates) + output_back <- transform(output_trans, to = "cell") + test_integrity(output_trans) + test_integrity(output_back) + if (!anyNA(output$data)) { + expect_false(anyNA(output_trans$data)) + } + + output2 <- read_grid( + filename = file_name, + dim_order = c("band", "cell", "time") + ) + + output_trans2 <- transform(output2, to = "lon_lat") + output_trans2$subset(coordinates = coordinates) + output_back2 <- transform(output_trans2, to = "cell") + test_integrity(output_trans2) + test_integrity(output_back2) + if (!anyNA(output2$data)) { + expect_false(anyNA(output_trans2$data)) + } + expect_identical( + output_trans2$data, + aperm(output_trans$data, names(dim(output_trans2))) + ) + expect_identical( + output_back2$data, + aperm(output_back$data, names(dim(output_back2))) + ) + + # For coordinate subsetting output has to be transformed first + expect_error( + output$subset(coordinates = coordinates), + "convert into suitable format" + ) + + # Coordinate values must be strings, not numerical values + output_trans <- transform(output, to = c("lon_lat")) + expect_error( + output_trans$subset(coordinates = coordinates_numeric), + "Values for coordinate pairs must be supplied as strings" + ) +}) diff --git a/tests/testthat/test-LPJmLData_subset_transform.R b/tests/testthat/test-LPJmLData_subset_transform.R index 8c93ecb..fb62bb0 100644 --- a/tests/testthat/test-LPJmLData_subset_transform.R +++ b/tests/testthat/test-LPJmLData_subset_transform.R @@ -11,9 +11,9 @@ test_integrity <- function(output) { # Check for two cases "cell" or "lon_lat" if ("cell" %in% names(dim_data)) { # Test for equal length of cell in data and meta data (ncell) - expect_equal(dim_data[["cell"]], output$meta$ncell) + testthat::expect_equal(dim_data[["cell"]], output$meta$ncell) # Test for equal dimnames of cell in data and those constructed by meta data - expect_equal( + testthat::expect_equal( dimnames_data$cell, format( seq(output$meta$firstcell, length.out = output$meta$ncell), @@ -22,20 +22,20 @@ test_integrity <- function(output) { ) } else { # Test for equal dimnames of lat, lon in data and those of underlying grid - expect_equal(dimnames_data$lat, - dimnames(output$grid)$lat) - expect_equal(dimnames_data$lon, - dimnames(output$grid)$lon) + testthat::expect_equal(dimnames_data$lat, + dimnames(output$grid)$lat) + testthat::expect_equal(dimnames_data$lon, + dimnames(output$grid)$lon) } # Check for two cases "time" or "year_month_day" if ("time" %in% names(dim_data)) { # Test for equal length of time steps in data and meta data (nyear * nstep) - expect_equal(dim_data[["time"]], - output$meta$nyear * output$meta$nstep) + testthat::expect_equal(dim_data[["time"]], + output$meta$nyear * output$meta$nstep) # Test for equal dimnames of time steps in data and those constructed by # meta data with create_time_names function (nstep, firstyear, nyear) - expect_equal( + testthat::expect_equal( dimnames_data$time, create_time_names( output$meta$nstep, @@ -45,10 +45,10 @@ test_integrity <- function(output) { ) } else { # Test for equal length of years in data and meta data (nyear) - expect_equal(dim_data[["year"]], output$meta$nyear) + testthat::expect_equal(dim_data[["year"]], output$meta$nyear) # Test for equal dimnames of years in data and those constructed by # meta data (firstyear, nyear) - expect_equal( + testthat::expect_equal( dimnames_data$year, format( seq(output$meta$firstyear, by = output$meta$timestep, @@ -60,19 +60,19 @@ test_integrity <- function(output) { # For month there is no meta data available (like nmonth, firstmonth) # following tests only via hardcoded pre defined month to be tested if ("month" %in% names(dim_data) && !output$meta$subset) { - expect_equal(dimnames_data$month, as.character(1:12)) + testthat::expect_equal(dimnames_data$month, as.character(1:12)) } else if ("month" %in% names(dim_data) && output$meta$subset) { - expect_equal(dimnames_data$month, as.character(6:9)) + testthat::expect_equal(dimnames_data$month, as.character(6:9)) } } # Test for equal length of bands in data and meta data (nbands) - expect_equal(dim_data[["band"]], output$meta$nbands) + testthat::expect_equal(dim_data[["band"]], output$meta$nbands) # Check if band dimension > 1 -> then has band_names if (output$meta$nbands > 1) { # Test for equal dimnames of band in data and those constructed by meta data # (band_names) - expect_equal(dimnames_data$band, output$meta$band_names) + testthat::expect_equal(dimnames_data$band, output$meta$band_names) } # check for grid @@ -84,7 +84,7 @@ test_integrity <- function(output) { if ("cell" %in% names(dimnames_grid)) { # Test for equal dimnames of cell in grid data and those constructed by # output meta data - expect_equal( + testthat::expect_equal( dimnames_grid$cell, format( seq(output$meta$firstcell, length.out = output$meta$ncell), @@ -94,7 +94,7 @@ test_integrity <- function(output) { } else { # Test to match data of grid (cell numbers) and cell numbers constructed # by meta data of output - expect_true( + testthat::expect_true( all(as.vector(stats::na.omit(output$grid$data)) %in% seq(output$meta$firstcell, length.out = output$meta$ncell)) ) @@ -112,6 +112,7 @@ test_that("test subset method", { output$add_grid(), "grid.bin.json" ) + output_sub <- subset( output, cell = 1:2, diff --git a/tests/testthat/test-LPJmLGridData.R b/tests/testthat/test-LPJmLGridData.R index a91d7e4..66c5d97 100644 --- a/tests/testthat/test-LPJmLGridData.R +++ b/tests/testthat/test-LPJmLGridData.R @@ -97,23 +97,8 @@ test_that("test LPJmLGridData methods", { ) expect_error( - grid$subset(), - "Not allowed for an object of class LPJmLGridData" - ) - - expect_error( - grid$transform(), - "Not allowed for an object of class LPJmLGridData" - ) - - expect_error( - grid$as_raster(), - "Not allowed for an object of class LPJmLGridData" - ) - - expect_error( - grid$as_terra(), - "Not allowed for an object of class LPJmLGridData" + grid$plot(), + "No plot function available for LPJmLGridData object." ) # Only variable "grid" (and "LPJGRID") are valid to init LPJmLGridData diff --git a/tests/testthat/test-LPJmLGridData_as.R b/tests/testthat/test-LPJmLGridData_as.R new file mode 100644 index 0000000..f892ce7 --- /dev/null +++ b/tests/testthat/test-LPJmLGridData_as.R @@ -0,0 +1,116 @@ +# Test as_array method (basic way) +test_that("basic array export", { + file_name <- "../testdata/output/grid.bin.json" + output <- read_grid(filename = file_name) + output_array <- output$as_array() + + output$transform(to = "lon_lat") + + expect_true( + all( + output_array %in% as.numeric(unlist(dimnames(output$data))) + ) + ) +}) + + +# test as_tibble method (basic way) +test_that("basic tibble export", { + file_name <- "../testdata/output/grid.bin.json" + output <- read_grid(filename = file_name) + output_tibble <- tibble::as_tibble(output) + + # test for equality with test_tibble + expect_equal(output_tibble$value, as.vector(output$data)) +}) + + +# test as_raster method with lon lat dimensions +test_that("raster export lon_lat", { + file_name <- "../testdata/output/grid.bin.json" + output <- read_grid(filename = file_name) + + output_raster <- as_raster(output) + + output2 <- transform(output, to = "lon_lat") + replace_array <- output2$data + # create tmp_raster with expected dimensions and coordinate ref system + test_raster <- raster::raster(res = 0.5, + crs = raster_crs_fallback("EPSG:4326"), + xmn = -87.5, + xmx = -87, + ymn = 55, + ymx = 64) %>% + raster::brick(nl = 2) + + # replace values in tmp_raster with expected values of longitude + replace_array[which(!is.na(replace_array))] <- rev( + subset(output, band = "lon")$data + ) + test_raster <- raster::setValues( + test_raster, + values = replace_array, + layer = 1 + ) + + # replace values in tmp_raster with expected values of latitude + replace_array[which(!is.na(replace_array))] <- rev( + subset(output, band = "lat")$data + ) + test_raster <- raster::setValues( + test_raster, + values = replace_array, + layer = 2 + ) + + raster::crs(test_raster) <- raster_crs_fallback("EPSG:4326") + # add variable as layer name + names(test_raster) <- names(dim(output2)) + + # test for equality with test_raster + expect_equal(output_raster, test_raster) + + + test_raster2 <- output2$as_raster() + + # test raster matrix values with those of transformed grid + expect_equal(as.vector(test_raster2[]), as.vector(output2$data)) +}) + + +# test as_terra method with lon lat dimensions +test_that("terra export lon_lat", { + file_name <- "../testdata/output/grid.bin.json" + output <- read_grid(filename = file_name) + + output_rast <- as_terra(output) + + output2 <- transform(output, to = "lon_lat") + replace_array <- array(NA, dim = c(length(output2$data), 2)) + + + test_rast <- terra::rast( + crs = "EPSG:4326", + extent = terra::ext(-87.5, -87, 55, 64), + resolution = 0.5, + nlyrs = 2 + ) + + # replace values in tmp_rast with expected values of longitude and latitude + replace_array[which(!is.na(output2$data)), 1] <- rev(subset(output, band = "lon")$data) + replace_array[which(!is.na(output2$data)), 2] <- rev(subset(output, band = "lat")$data) + test_rast <- terra::setValues(test_rast, replace_array) + + # add variable as layer name + names(test_rast) <- names(dim(output2)) + terra::units(test_rast) <- output$meta$unit + + # test for equality with test_rast + expect_equal(output_rast[], test_rast[]) + + test_raster2 <- output2$as_raster() + + # test raster matrix values with those of transformed grid + expect_equal(as.vector(test_raster2[]), as.vector(output2$data)) + +}) diff --git a/tests/testthat/test-get_cellindex.R b/tests/testthat/test-get_cellindex.R new file mode 100644 index 0000000..1795651 --- /dev/null +++ b/tests/testthat/test-get_cellindex.R @@ -0,0 +1,64 @@ +test_that("get_cellindex handles invalid file path", { + + expect_error( + get_cellindex(""), + "grid_filename must be a string representing a valid file path." + ) + + expect_error( + get_cellindex("../testdata/output/invalid_file_path"), + "grid_filename does not exist." + ) +}) + +test_that("get_cellindex handles invalid coordinates", { + expect_error( + get_cellindex("../testdata/output/grid.bin.json", + coordinates = list(c(1.25, 2.75), c(1.25, 2.75, 3.25)) + ), + "The two vectors in coordinates must have the same length" + ) +}) + +test_that("get_cellindex handles both extent and coordinates provided", { + expect_error( + get_cellindex("../testdata/output/grid.bin.json", + 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." + ) +}) + +test_that("get_cellindex handles extent values out of order", { + expect_warning( + get_cellindex("../testdata/output/grid.bin.json", + extent = c(-88.25, -87.25, 25.75, 55.25) + ), + "Extent values out of bounds:" + ) +}) + +test_that("get_cellindex handles valid coordinates", { + expect_error( + get_cellindex( + "../testdata/output/grid.bin.json", + coordinates = list(c(-87.25, -87.25), c(55.25, 55.75)) + ), + "Values for coordinate pairs must be supplied as strings" + ) +}) + +test_that("get_cellindex returns correct cell index for given extent", { + result <- get_cellindex("../testdata/output/grid.bin.json", + extent = c(-87.25, -87.25, 55.25, 55.75) + ) + expect_true(length(result) == 2 && result[1] == 10001 && result[2] == 10002) +}) + +test_that("get_cellindex returns correct cell index for given coordinates", { + result <- get_cellindex("../testdata/output/grid.bin.json", + coordinates = list(lon = c("-87.25", "-87.25"), lat = c("55.25", "55.75")) + ) + expect_true(length(result) == 2 && result[1] == 10001 && result[2] == 10002) +})