get_layer(): bbox

[rellimylime]

#' Retrieve a WRI layer for a bounding box
#'
#' Streams one or more WRI layers from cloud storage and returns a
#' \code{SpatRaster} cropped to a bounding box.
#'
#' @param bbox numeric vector. Bounding box in the format
#' \code{c(xmin, ymin, xmax, ymax)} in WGS84 (EPSG:4326).
#' @param layer_id character. One or more layer IDs to retrieve. Use
#' \code{wri_overview()} to find valid IDs. Defaults to \code{"WRI_score"}.
#'
#' @return A \code{terra::SpatRaster} cropped to the bounding box.
#' Multi-layer if multiple IDs provided.
#'
#' @examples
#' # Bounding box around Butte County, CA
#' bbox <- c(-122.0, 39.5, -121.2, 40.2)
#' get_layer(bbox)
#'
#' # Specific layer
#' get_layer(bbox, layer_id = "water_domain_score")
#'
#' # Multiple layers
#' get_layer(bbox, layer_id = c("water_domain_score", "air_domain_score"))
#'
#' @export
get_layer <- function(bbox, layer_id = "WRI_score") {}

[FlukeAndFeather]

I suggest referring to a layer rather than a tif, since the former is more general and the latter is a particular file format, which the user doesn't necessarily need to know about.

[rellimylime]

How to implement #17

check <- check_extent_overlap(bbox, id)
status <- attr(check, "status")
msg <- attr(check, "message")

if (status == "none") {
  stop(msg, call. = FALSE)
}

if (status == "partial") {
  warning(msg, call. = FALSE)
}

[kaimorquecho]

Function

get_layer <- function(id, bbox = NULL) {

  # Validate inputs
  if (missing(id) || !is.character(id) || length(id) != 1 || !nzchar(id)) {
    stop("`id` must be a single non-empty character string.", call. = FALSE)
  }

  if (!is.null(bbox)) {
    if (!is.numeric(bbox) || length(bbox) != 4) {
      stop("`bbox` must be numeric: c(xmin, ymin, xmax, ymax)", call. = FALSE)
    }
  }

  # Get layer metadata (asset-level)
  layer <- layer_info(id)

  # Check layer exists
  layer <- .wri_pick_raster_asset(layer)

  if (nrow(layer) == 0) {
    stop(
      "No raster asset found for layer '", id, "'.",
      call. = FALSE
    )
  }

  # Ensure it's hosted
  is_hosted <- tolower(layer$is_hosted[1]) == "true"

  if (!is_hosted) {
    stop(
      "Layer '", id, "' is not hosted and cannot be streamed.",
      call. = FALSE
    )
  }

  href <- layer$asset_href[1]

  # Build GDAL VSI path (COG streaming)
  vsi_path <- paste0("/vsicurl/", href)

  tryCatch({

    rast <- terra::rast(vsi_path)

    if (!is.null(bbox)) {

      bbox_ext <- terra::ext(bbox[1], bbox[3], bbox[2], bbox[4])
      bbox_vect <- terra::as.polygons(bbox_ext, crs = "EPSG:4326")

      bbox_transformed <- terra::project(bbox_vect, terra::crs(rast))

      rast <- terra::crop(rast, bbox_transformed)
    }

    rast

  }, error = function(e) {
    stop(
      "Failed to retrieve layer '", id, "'.\n",
      "URL: ", href, "\n",
      "Error: ", conditionMessage(e),
      call. = FALSE
    )
  })
}

Example output:

get_layer("WRI_score")
class       : SpatRaster 
dimensions  : 57865, 52355, 1  (nrow, ncol, nlyr)
resolution  : 90, 90  (x, y)
extent      : -5216640, -504689.7, 991231.7, 6199082  (xmin, xmax, ymin, ymax)
coord. ref. : NAD83 / Conus Albers (EPSG:5070) 
source      : WRI_score.tif 
name        :     mean 
min value   : 48.31954 
max value   : 90.27865