diff --git a/DESCRIPTION b/DESCRIPTION index 2b14576..d5a6797 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -20,6 +20,7 @@ Imports: terra, tibble Suggests: - testthat (>= 3.0.0) + testthat (>= 3.0.0), + withr Config/testthat/edition: 3 LazyData: true diff --git a/R/check_type_and_length.R b/R/check_type_and_length.R new file mode 100644 index 0000000..2bec7df --- /dev/null +++ b/R/check_type_and_length.R @@ -0,0 +1,51 @@ +check_type_and_length <- function(..., + call = rlang::caller_env(), + env = rlang::caller_env()) { + dots <- list(...) + if (length(dots) == 0) { + return(invisible(TRUE)) + } + + if (is.null(names(dots)) || any(names(dots) == "")) { + rlang::abort("All arguments to `check_type_and_length()` must be named.") + } + + problem_args <- character() + for (dot in names(dots)) { + arg <- get(dot, envir = env) + if (is.null(arg)) next + + arg_class <- class(arg) + dot_class <- class(dots[[dot]]) + if (arg_class != dot_class) { + if ("integer" %in% arg_class && "numeric" %in% dot_class) { + next # Purposefully doing nothing -- rely on implicit conversion + } else if ("integer" %in% dot_class && rlang::is_integerish(arg)) { + next # Purposefully doing nothing -- rely on implicit conversion + } else { + problem_args <- c( + problem_args, + glue::glue("{dot} should be a {class(dots[[dot]])}, but is a {class(arg)}.") + ) + } + } + if (length(dots[[dot]]) && length(arg) != length(dots[[dot]])) { + problem_args <- c( + problem_args, + glue::glue("{dot} should be of length {length(dots[[dot]])}, but is length {length(arg)}.") + ) + } + } + + if (length(problem_args)) { + names(problem_args) <- rep("*", length(problem_args)) + rlang::abort( + c( + "Some input arguments weren't the right class or length:", + problem_args + ), + call = call + ) + } + return(invisible(TRUE)) +} diff --git a/R/get_stac_data.R b/R/get_stac_data.R index e4d8dcc..b8da1ca 100644 --- a/R/get_stac_data.R +++ b/R/get_stac_data.R @@ -1,11 +1,15 @@ -#' Retrieve images (or composite images) from STAC endpoints +#' Retrieve and composite images from STAC endpoints #' -#' This function retrieves composites of satellite images from STAC endpoints. +#' These function retrieves composites of satellite images from STAC endpoints. +#' `get_stac_data()` is a generic function which should be able to download and +#' composite images from a variety of data sources, while the other helper +#' functions have useful defaults for downloading common data sets from standard +#' STAC sources. #' #' @section Usage Tips: #' It's often useful to buffer your `aoi` object slightly, on the order of 1-2 #' cell widths, in order to ensure that data is downloaded for your entire AOI -#' even after accounting for any reprojection necessary to compare your AOI to +#' even after accounting for any reprojection needed to compare your AOI to #' the data on the STAC server. #' #' Setting the following GDAL configuration variables via `Sys.setenv()` or in @@ -75,17 +79,11 @@ #' @param mask_function A function that takes a raster and returns a boolean #' raster, where `TRUE` pixels will be preserved and `FALSE` or `NA` pixels will #' be masked out. See [sentinel2_mask_function()]. -#' @param output_filename The filename to write the output raster to. If -#' `composite_function` is `NULL`, something will be appended to the end of the -#' filename (before the extension) to create a unique file for each downloaded -#' image. Currently this means sequential numbers are added to each file, but -#' this may change in the future -- do not depend on the resulting filenames -#' staying the same across versions! -#' @param composite_function The (quoted) name of a function from -#' `terra` (for instance, [terra::median]) used to combine downloaded images -#' into a single composite (i.e., to aggregate pixel values from multiple images -#' into a single value). Set to `NULL` to not composite images, but instead -#' save each image as a separate file. +#' @param output_filename The filename to write the output raster to. +#' @param composite_function Character of length 1: The (quoted) name of a +#' function from `terra` (for instance, [terra::median]) used to combine +#' downloaded images into a single composite (i.e., to aggregate pixel values +#' from multiple images into a single value). #' @inheritParams rstac::stac_search #' @param gdalwarp_options Options passed to `gdalwarp` through the `options` #' argument of [sf::gdal_utils()]. The same set of options are used for all @@ -139,6 +137,29 @@ get_stac_data <- function(aoi, "-co", "PREDICTOR=2", "-co", "NUM_THREADS=ALL_CPUS" )) { + + if (!(inherits(aoi, "sf") || inherits(aoi, "sfc"))) { + rlang::abort( + "`aoi` must be an sf or sfc object.", + class = "rsi_aoi_not_sf" + ) + } + + check_type_and_length( + start_date = character(1), + end_date = character(1), + pixel_x_size = numeric(1), + pixel_y_size = numeric(1), + stac_source = character(1), + collection = character(1), + rescale_bands = logical(1), + mask_band = character(1), + output_filename = character(1), + composite_function = character(1), + limit = numeric(1), + gdalwarp_options = character() + ) + gdalwarp_options <- process_gdalwarp_options( gdalwarp_options = gdalwarp_options, aoi = aoi, @@ -248,54 +269,52 @@ get_stac_data <- function(aoi, ) } - if (!is.null(composite_function)) { - download_dir <- file.path(tempdir(), "composite_dir") - if (!dir.exists(download_dir)) dir.create(download_dir) - - composited_bands <- vapply( - names(downloaded_bands), - function(band_name) { - out_file <- file.path(download_dir, paste0(toupper(band_name), ".tif")) - - do.call( - getFromNamespace(composite_function, "terra"), - list( - x = terra::rast(downloaded_bands[[band_name]]), - na.rm = TRUE, - filename = out_file, - overwrite = TRUE - ) + download_dir <- file.path(tempdir(), "composite_dir") + if (!dir.exists(download_dir)) dir.create(download_dir) + + composited_bands <- vapply( + names(downloaded_bands), + function(band_name) { + out_file <- file.path(download_dir, paste0(toupper(band_name), ".tif")) + + do.call( + getFromNamespace(composite_function, "terra"), + list( + x = terra::rast(downloaded_bands[[band_name]]), + na.rm = TRUE, + filename = out_file, + overwrite = TRUE ) + ) - if (rescale_bands) { - out_file <- rescale_band(downloaded_bands, band_name, out_file) - } + if (rescale_bands) { + out_file <- rescale_band(downloaded_bands, band_name, out_file) + } - out_file - }, - character(1) - ) - on.exit(file.remove(composited_bands), add = TRUE) - - out_vrt <- tempfile(fileext = ".vrt") - invisible( - stack_rasters( - composited_bands, - out_vrt, - band_names = remap_band_names(names(items_urls), asset_names) - ) - ) - on.exit(file.remove(out_vrt), add = TRUE) + out_file + }, + character(1) + ) + on.exit(file.remove(composited_bands), add = TRUE) - sf::gdal_utils( - "warp", + out_vrt <- tempfile(fileext = ".vrt") + invisible( + stack_rasters( + composited_bands, out_vrt, - output_filename, - options = gdalwarp_options + band_names = remap_band_names(names(items_urls), asset_names) ) + ) + on.exit(file.remove(out_vrt), add = TRUE) - output_filename - } + sf::gdal_utils( + "warp", + out_vrt, + output_filename, + options = gdalwarp_options + ) + + output_filename } #' @rdname get_stac_data diff --git a/R/spectral_indices.R b/R/spectral_indices.R index f8bb026..054206f 100644 --- a/R/spectral_indices.R +++ b/R/spectral_indices.R @@ -50,7 +50,7 @@ spectral_indices <- function(..., url = spectral_indices_url(), update_cache = N if (update_cache) { tryCatch( - download_indices(url), + update_cached_indices(url), error = function(e) { rlang::warn(c( "Failed to update the cache of indices.", @@ -99,7 +99,7 @@ download_indices <- function(url = spectral_indices_url()) { update_cached_indices <- function(url = spectral_indices_url()) { if (!dir.exists(tools::R_user_dir("rsi"))) { - dir.create(tools::R_user_dir("rsi")) + dir.create(tools::R_user_dir("rsi"), recursive = TRUE) } indices_path <- file.path( tools::R_user_dir("rsi"), diff --git a/R/stack_rasters.R b/R/stack_rasters.R index 7a3137f..e9f174c 100644 --- a/R/stack_rasters.R +++ b/R/stack_rasters.R @@ -102,31 +102,31 @@ stack_rasters <- function(rasters, ref_crs <- terra::crs(ref_rast) - if (missing(band_names) || is.function(band_names)) { - var_names <- unlist( - lapply( - rasters, - function(r) { - r <- terra::rast(r) - # this is the only place we instantiate these rasters, so may as well - # check CRS alignment while we're here... - if (terra::crs(r) != ref_crs) { - rlang::abort(c( - "Rasters do not all share the reference raster's CRS.", - i = "Reproject rasters to all share the same CRS." - ), - class = "rsi_multiple_crs" - ) - } - names(r) + var_names <- unlist( + lapply( + rasters, + function(r) { + r <- terra::rast(r) + # this is the only place we instantiate these rasters, so may as well + # check CRS alignment while we're here... + if (terra::crs(r) != ref_crs) { + rlang::abort(c( + "Rasters do not all share the reference raster's CRS.", + i = "Reproject rasters to all share the same CRS." + ), + class = "rsi_multiple_crs", + call = rlang::caller_env() + ) } - ) + names(r) + } ) - } + ) + if (!missing(band_names) && is.function(band_names)) { var_names <- band_names(var_names) } - if (is.character(band_names)) { + if (!missing(band_names) && is.character(band_names)) { var_names <- band_names } diff --git a/data-raw/band_mappings.R b/data-raw/band_mappings.R index ff0d8d8..83bd520 100644 --- a/data-raw/band_mappings.R +++ b/data-raw/band_mappings.R @@ -1,3 +1,4 @@ +# devtools::install() library(rsi) sentinel2_band_mapping <- list( aws_v0 = structure( diff --git a/man/get_stac_data.Rd b/man/get_stac_data.Rd index 1e49451..2d796e3 100644 --- a/man/get_stac_data.Rd +++ b/man/get_stac_data.Rd @@ -6,7 +6,7 @@ \alias{get_sentinel2_imagery} \alias{get_landsat_imagery} \alias{get_dem} -\title{Retrieve images (or composite images) from STAC endpoints} +\title{Retrieve and composite images from STAC endpoints} \usage{ get_stac_data( aoi, @@ -169,18 +169,12 @@ STAC source to use to mask the data. Set to \code{NULL} to not mask. See the raster, where \code{TRUE} pixels will be preserved and \code{FALSE} or \code{NA} pixels will be masked out. See \code{\link[=sentinel2_mask_function]{sentinel2_mask_function()}}.} -\item{output_filename}{The filename to write the output raster to. If -\code{composite_function} is \code{NULL}, something will be appended to the end of the -filename (before the extension) to create a unique file for each downloaded -image. Currently this means sequential numbers are added to each file, but -this may change in the future -- do not depend on the resulting filenames -staying the same across versions!} +\item{output_filename}{The filename to write the output raster to.} \item{composite_function}{The (quoted) name of a function from \code{terra} (for instance, \link[terra:summarize-generics]{terra::median}) used to combine downloaded images into a single composite (i.e., to aggregate pixel values from multiple images -into a single value). Set to \code{NULL} to not composite images, but instead -save each image as a separate file.} +into a single value).} \item{limit}{an \code{integer} defining the maximum number of results to return. If not informed, it defaults to the service implementation.} @@ -195,13 +189,17 @@ varying data types.} \code{output_filename}, unchanged. } \description{ -This function retrieves composites of satellite images from STAC endpoints. +These function retrieves composites of satellite images from STAC endpoints. +\code{get_stac_data()} is a generic function which should be able to download and +composite images from a variety of data sources, while the other helper +functions have useful defaults for downloading common data sets from standard +STAC sources. } \section{Usage Tips}{ It's often useful to buffer your \code{aoi} object slightly, on the order of 1-2 cell widths, in order to ensure that data is downloaded for your entire AOI -even after accounting for any reprojection necessary to compare your AOI to +even after accounting for any reprojection needed to compare your AOI to the data on the STAC server. Setting the following GDAL configuration variables via \code{Sys.setenv()} or in diff --git a/tests/testthat/test-calculate_indices.R b/tests/testthat/test-calculate_indices.R index 18e3d3e..3029192 100644 --- a/tests/testthat/test-calculate_indices.R +++ b/tests/testthat/test-calculate_indices.R @@ -17,3 +17,25 @@ test_that("Index calculation is stable", { terra::values(terra::rast(system.file("rasters/dpdd.tif", package = "rsi"))) ) }) + +test_that("Index calculations fail when missing a column", { + expect_error( + calculate_indices( + system.file("rasters/example_sentinel1.tif", package = "rsi"), + filter_platforms(platforms = "Sentinel-1 (Dual Polarisation VV-VH)")["formula"], + index_out + ), + class = "rsi_missing_column" + ) +}) + +test_that("Index calculations fail when missing bands", { + expect_error( + calculate_indices( + system.file("rasters/example_sentinel1.tif", package = "rsi"), + filter_platforms(platforms = "Landsat-OLI"), + index_out + ), + class = "rsi_missing_indices" + ) +}) diff --git a/tests/testthat/test-get_stac_data.R b/tests/testthat/test-get_stac_data.R new file mode 100644 index 0000000..e136e12 --- /dev/null +++ b/tests/testthat/test-get_stac_data.R @@ -0,0 +1,97 @@ +test_that("get_landsat_imagery() is stable", { + skip_on_cran() + aoi <- sf::st_point(c(-74.912131, 44.080410)) |> + sf::st_sfc() |> + sf::st_set_crs(4326) |> + sf::st_transform(3857) |> + sf::st_buffer(1000) + + expect_no_error( + out <- get_landsat_imagery( + aoi, + "2022-06-01", + "2022-06-30", + output_filename = tempfile(fileext = ".tif") + ) + ) + expect_no_error(terra::rast(out)) + expect_true( + all( + names(terra::rast(out)) %in% + as.vector(landsat_band_mapping$planetary_computer_v1) + ) + ) +}) + +test_that("get_sentinel1_imagery() is stable", { + skip_on_cran() + aoi <- sf::st_point(c(-74.912131, 44.080410)) |> + sf::st_sfc() |> + sf::st_set_crs(4326) |> + sf::st_transform(3857) |> + sf::st_buffer(1000) + + expect_no_error( + out <- get_sentinel1_imagery( + aoi, + "2022-06-01", + "2022-06-30", + output_filename = tempfile(fileext = ".tif") + ) + ) + expect_no_error(terra::rast(out)) + expect_true( + all( + names(terra::rast(out)) %in% + as.vector(sentinel1_band_mapping$planetary_computer_v1) + ) + ) +}) + +test_that("get_sentinel2_imagery() is stable", { + skip_on_cran() + aoi <- sf::st_point(c(-74.912131, 44.080410)) |> + sf::st_sfc() |> + sf::st_set_crs(4326) |> + sf::st_transform(3857) |> + sf::st_buffer(1000) + + expect_no_error( + out <- get_sentinel2_imagery( + aoi, + "2022-06-01", + "2022-06-30", + output_filename = tempfile(fileext = ".tif") + ) + ) + expect_no_error(terra::rast(out)) + expect_true( + all( + names(terra::rast(out)) %in% + as.vector(sentinel2_band_mapping$planetary_computer_v1) + ) + ) +}) + +test_that("get_dem() is stable", { + skip_on_cran() + aoi <- sf::st_point(c(-74.912131, 44.080410)) |> + sf::st_sfc() |> + sf::st_set_crs(4326) |> + sf::st_transform(3857) |> + sf::st_buffer(1000) + + expect_no_error( + out <- get_dem( + aoi, + output_filename = tempfile(fileext = ".tif") + ) + ) + expect_no_error(terra::rast(out)) + expect_true( + all( + names(terra::rast(out)) %in% + as.vector(dem_band_mapping$planetary_computer_v1) + ) + ) +}) diff --git a/tests/testthat/test-spectral_indices.R b/tests/testthat/test-spectral_indices.R index 7545717..fd4ae27 100644 --- a/tests/testthat/test-spectral_indices.R +++ b/tests/testthat/test-spectral_indices.R @@ -18,6 +18,27 @@ test_that("spectral_indices() works", { # we should be able to update our cache: skip_on_cran() expect_no_warning(spectral_indices()) + expect_no_warning(spectral_indices(update_cache = TRUE)) }) +test_that("spectral_indices_url() respects options", { + skip_if_not_installed("withr") + expect_identical( + withr::with_options( + list("rsi_url" = "example"), + spectral_indices_url() + ), + "example" + ) +}) +test_that("spectral_indices_url() respects environment variables", { + skip_if_not_installed("withr") + expect_identical( + withr::with_envvar( + list("rsi_url" = "example"), + spectral_indices_url() + ), + "example" + ) +}) diff --git a/tests/testthat/test-stack_rasters.R b/tests/testthat/test-stack_rasters.R index d3f66d0..25836ea 100644 --- a/tests/testthat/test-stack_rasters.R +++ b/tests/testthat/test-stack_rasters.R @@ -12,3 +12,78 @@ test_that("stack_rasters works", { terra::values() ) }) + +test_that("stack_rasters fails when reference_raster isn't in the vector", { + expect_error( + stack_rasters( + list( + system.file("rasters/example_sentinel1.tif", package = "rsi") + ), + tempfile(fileext = ".vrt"), + reference_raster = 2 + ), + class = "rsi_not_in_vec" + ) + + expect_error( + stack_rasters( + list( + system.file("rasters/example_sentinel1.tif", package = "rsi") + ), + tempfile(fileext = ".vrt"), + reference_raster = "some_raster" + ), + class = "rsi_not_in_vec" + ) +}) + +test_that("stack_rasters fails when extent isn't four numbers", { + expect_error( + stack_rasters( + list( + system.file("rasters/example_sentinel1.tif", package = "rsi") + ), + tempfile(fileext = ".vrt"), + extent = 20 + ), + class = "rsi_bad_extent" + ) +}) + +test_that("stack_rasters fails when resolution isn't 1-2 numbers", { + expect_error( + stack_rasters( + list( + system.file("rasters/example_sentinel1.tif", package = "rsi") + ), + tempfile(fileext = ".vrt"), + resolution = c(20, 30, 40) + ), + class = "rsi_bad_resolution" + ) +}) + +test_that("stack_rasters fails when rasters don't share a CRS", { + s1 <- tempfile(fileext = ".tif") + terra::writeRaster( + terra::project( + terra::rast( + system.file("rasters/example_sentinel1.tif", package = "rsi") + ), + "EPSG:4326" + ), + s1 + ) + + expect_error( + stack_rasters( + list( + system.file("rasters/example_sentinel1.tif", package = "rsi"), + s1 + ), + tempfile(fileext = ".vrt") + ), + class = "rsi_multiple_crs" + ) +}) +