Potential create_mosaic_raster for spatialplanr

[ric325]

#' Creates a single terra raster layer with a mosaic of patches with the specified number of categories and approximate percentage area of each based on relative cover values provided
#'
#' The area of each category in the mosaic will sum approximately to its specified relative cover provided
#' The patch generation process works by using a random initial patch centre
#' Thus, higher cover numbers lead to more initial patch centres and
#' The number of patch centres for each category is specified by the relative cover of each category
#' Patches then grow by assigning each pixel to the nearest random patch centre, resulting in patches for each category
#' Note that the number of patches depends on how dispersed the initial random centres are and this is more likely with higher covers
#'
#' @param nrows Integer specifying the number of columns in the raster
#' @param ncols Integer specifying the number of columns in the raster
#' @param setseed Boolean specifying whether the mosaic raster should be reproducible (1) or random (0)
#' @param cover A named vector of category cover (does not need to sum to 100)
#'
#' @return A terra raster stack of a single terra raster layer with a mosaic of patches with the specified number of categories and approximate % area of each
#'
#' @examples
#' \dontrun{
#' mosaic_cover <- create_mosaic_raster(
#' nrows = 100, ncols = 100, setseed = 1,
#' cover = c("cocoa plantation" = 40,
#' "disturbed forest" = 20,
#' "undisturbed forest" = 10,
#' "other areas" = 20,
#' "protected areas" = 10)
#' )
#' plot(mosaic_cover)
#'
#'#' \dontrun{
#' mosaic_cover <- create_mosaic_raster(
#' nrows = 10, ncols = 10, setseed = 0,
#' cover = c("Zone1" = 4,
#' "Zone2" = 2)
#' )
#' plot(mosaic_cover)

#'#' \dontrun{
#' mosaic_cover <- create_mosaic_raster(
#' nrows = 10, ncols = 10, setseed = 0,
#' cover = c("Zone1" = 40,
#' "Zone2" = 20)
#' )
#' plot(mosaic_cover)

create_mosaic_raster <- function(nrows = 100,
ncols = 100,
setseed = 1,
cover = c("Category1" = 50,
"Category2" = 30,
"Category3" = 20,
)){

Assign integer codes to each category and store their names

categories <- seq_along(cover)
names_vec <- names(cover)
total_cells <- nrows * ncols

Generate random patch centres for each category

The number of patch centres for each category is set by its cover value

if(setseed == 1) set.seed(123) # for reproducibility

patch_centres <- tibble::tibble(
x = sample(1:ncols, sum(cover), replace = TRUE), # x-coordinates of patch centers
y = sample(1:nrows, sum(cover), replace = TRUE), # y-coordinates of patch centers
cat = rep(categories, cover) # category code for each patch center
)

Create a grid of coordinates for all pixels in the raster

pixel_coords <- expand.grid(x = 1:ncols, y = 1:nrows)

For each pixel, find the nearest patch center and assign its category

This ensures pixels close to a patch center are grouped together (clumped)

assign_category <- function(px, py) {
dists <- sqrt((patch_centres$x - px)^2 + (patch_centres$y - py)^2) # Euclidean distance to all patch centers
nearest <- which.min(dists) # Find the nearest patch center
patch_centres$cat[nearest] # Assign the category of that patch center
}

Use purrr::map2_int to efficiently assign categories to all pixels

assigned <- purrr::map2_int(pixel_coords$x, pixel_coords$y, assign_category)

Convert the assigned categories into a matrix and then to a raster

mat <- matrix(assigned, nrow = nrows, ncol = ncols, byrow = FALSE)
r <- terra::rast(mat)

Add category names as levels for easier plotting and interpretation

levels(r) <- data.frame(id = categories, cover = names_vec)

Return the patchy raster

return(r)
}