R/spatial.R

In metaweave: Spatial Ecological Network Inference

#' Create a standard raster grid
#' @param extent Numeric vector `c(xmin, xmax, ymin, ymax)`.
#' @param resolution Cell resolution.
#' @param crs Coordinate reference system.
#' @return A `terra::SpatRaster` with the requested extent, resolution, and coordinate reference system.
#' @examples
#' create_standard_grid(c(0, 2, 0, 2), resolution = 1)
#' @export
create_standard_grid <- function(extent = c(-180, 180, -90, 90), resolution = 1, crs = "EPSG:4326") {
  terra::rast(xmin = extent[1], xmax = extent[2], ymin = extent[3], ymax = extent[4], resolution = resolution, crs = crs)
}

#' Create a bounding-box polygon
#' @param extent Numeric vector `c(xmin, xmax, ymin, ymax)`.
#' @param crs Coordinate reference system.
#' @return A polygon `terra::SpatVector` representing the bounding box.
#' @export
make_bbox <- function(extent, crs = "EPSG:4326") {
  terra::as.polygons(terra::ext(extent), crs = crs)
}

#' List ESRI shapefiles
#' @param path Directory.
#' @param recursive Search recursively.
#' @return A character vector of shapefile paths, empty if no matches are found.
#' @export
list_shapefiles <- function(path, recursive = FALSE) {
  list.files(path, pattern = "\\.shp$", full.names = TRUE, recursive = recursive, ignore.case = TRUE)
}

#' Rasterize one range file
#' @param shp_file Shapefile path.
#' @param template Raster template.
#' @param touches Count all touched cells.
#' @return A one-layer `terra::SpatRaster`, named from the file stem, with 1 for presence and `NA` for background.
#' @export
rasterize_range_file <- function(shp_file, template, touches = TRUE) {
  v <- terra::project(terra::vect(shp_file), terra::crs(template))
  r <- terra::rasterize(v, template, field = 1, background = NA, touches = touches)
  names(r) <- tools::file_path_sans_ext(basename(shp_file))
  r
}

#' Rasterize rows grouped by species
#' @param x An sf or SpatVector object.
#' @param template Raster template.
#' @param species_col Species-name column.
#' @param species_keep Optional species subset.
#' @param touches Count all touched cells.
#' @return A `terra::SpatRaster` with one named layer per retained species, using 1 for presence and `NA` for background.
#' @export
rasterize_range_rows <- function(x, template, species_col, species_keep = NULL, touches = TRUE) {
  v <- terra::vect(x)
  vals <- terra::values(v)
  if (!species_col %in% names(vals)) stop("Species column not found.")
  spp <- as.character(vals[[species_col]])
  if (!is.null(species_keep)) spp <- intersect(unique(spp), species_keep) else spp <- unique(spp)
  out <- lapply(spp, function(sp) {
    one <- v[as.character(terra::values(v)[[species_col]]) == sp, ]
    one <- terra::project(one, terra::crs(template))
    r <- terra::rasterize(one, template, field = 1, background = NA, touches = touches)
    names(r) <- sp
    r
  })
  if (!length(out)) stop("No species remained after filtering.")
  terra::rast(out)
}

#' Build a species range stack
#' @param source Directory, vector file, sf, or SpatVector.
#' @param template Raster template.
#' @param species_col Species column for multi-species data.
#' @param species_keep Optional species subset.
#' @param region Optional crop region.
#' @param touches Count all touched cells.
#' @return A `terra::SpatRaster` with named species layers, optionally cropped and filtered to the region.
#' @export
build_range_stack <- function(source, template, species_col = NULL, species_keep = NULL, region = NULL, touches = TRUE) {
  if (is.character(source) && length(source) == 1L && dir.exists(source)) {
    files <- list_shapefiles(source)
    if (!length(files)) stop("No shapefiles found.")
    out <- terra::rast(lapply(files, rasterize_range_file, template = template, touches = touches))
  } else {
    if (is.character(source)) source <- terra::vect(source)
    if (is.null(species_col)) stop("`species_col` is required for multi-species data.")
    out <- rasterize_range_rows(source, template, species_col, species_keep, touches)
  }
  if (!is.null(region)) out <- crop_filter_stack(out, region)
  out
}

#' Crop a stack and remove species absent from the region
#' @param stack Species raster stack.
#' @param region Crop region.
#' @param return_cropped Return cropped rather than original layers.
#' @return A `terra::SpatRaster` containing layers with nonmissing cells in the crop region. Zero values count as nonmissing; this helper expects presence/NA range rasters.
#' @export
crop_filter_stack <- function(stack, region, return_cropped = TRUE) {
  cropped <- terra::crop(stack, region)
  present <- as.numeric(terra::global(!is.na(cropped), "sum", na.rm = TRUE)[, 1]) > 0
  if (return_cropped) cropped[[present]] else stack[[present]]
}

#' Prepare palm and mammal range stacks
#' @param palm_source,mammal_source Range sources accepted by `build_range_stack()`.
#' @param mammal_species_col Species-name column in the mammal source.
#' @param mammal_species_keep Optional mammal species subset.
#' @param output_dir Optional directory in which to save the range stacks.
#' @inheritParams build_range_stack
#' @return A list with `palms` and `mammals`, each a `terra::SpatRaster`. If `output_dir` is supplied, RDS files are also written there.
#' @export
prepare_range_stacks <- function(palm_source, mammal_source, mammal_species_col, mammal_species_keep = NULL,
                                 template = create_standard_grid(), region = NULL, output_dir = NULL, touches = TRUE) {
  palms <- build_range_stack(palm_source, template, region = region, touches = touches)
  mammals <- build_range_stack(mammal_source, template, mammal_species_col, mammal_species_keep, region, touches)
  if (!is.null(output_dir)) {
    dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
    saveRDS(palms, file.path(output_dir, "palm_ranges_grid_stack.RDS"))
    saveRDS(mammals, file.path(output_dir, "mammal_ranges_grid_stack.RDS"))
  }
  list(palms = palms, mammals = mammals)
}

#' Assemble local communities from aligned distribution rasters
#' @param distributions Nonempty named list of `SpatRaster` or
#'   `distribution_collection` objects with unique group and species names,
#'   aligned geometry, and nonnegative values. Zero and `NA` indicate absence.
#' @param min_species Positive integer minimum species per group.
#' @return A list of `assemblage` objects, one per cell meeting `min_species` in every group. Nonzero, nonmissing raster values indicate presence; the list is empty if no cell qualifies.
#' @export
assemble_communities <- function(distributions, min_species = 1L) {
  if (!is.list(distributions) || !length(distributions) ||
      is.null(names(distributions)) || anyNA(names(distributions)) ||
      any(!nzchar(names(distributions))) || anyDuplicated(names(distributions))) {
    stop("`distributions` must be a non-empty list with unique group names.")
  }
  if (!is.numeric(min_species) || length(min_species) != 1L ||
      !is.finite(min_species) || min_species < 1 || min_species != floor(min_species)) {
    stop("`min_species` must be one positive integer.")
  }
  rasters <- lapply(distributions, function(x) if (inherits(x, "distribution_collection")) x$data else x)
  if (any(vapply(rasters, function(x) !inherits(x, "SpatRaster") ||
                 !length(x) || is.null(names(x)) || anyNA(names(x)) ||
                 any(!nzchar(names(x))) || anyDuplicated(names(x)) > 0L,
                 logical(1)))) {
    stop("Each group must have a non-empty raster with unique species layer names.")
  }
  ref <- rasters[[1]]
  ok <- vapply(rasters[-1], function(x) terra::compareGeom(ref, x, stopOnError = FALSE), logical(1))
  if (length(ok) && !all(ok)) stop("All distribution rasters must share geometry.")
  values <- lapply(rasters, terra::values, mat = TRUE)
  if (any(vapply(values, function(v) any(is.infinite(v)) ||
                 any(v < 0, na.rm = TRUE), logical(1)))) {
    stop("Distribution raster values must be non-negative and not infinite.")
  }
  valid <- Reduce(`&`, lapply(values, function(v) rowSums(!is.na(v) & v != 0) >= min_species))
  cells <- which(valid)
  xy <- terra::xyFromCell(ref, cells)
  lapply(seq_along(cells), function(k) {
    i <- cells[k]
    spp <- Map(function(v, r) names(r)[which(!is.na(v[i, ]) & v[i, ] != 0)], values, rasters)
    new_assemblage(spp, cells[k], c(x = xy[k, 1], y = xy[k, 2]))
  })
}