Nothing
#' randomize_raster
#'
#' @description Randomized-habitats procedure
#'
#' @param raster SpatRaster with discrete habitat classes.
#' @param n_random Integer with number of randomizations.
#' @param directions Interger with cells neighbourhood rule: 4 (rook's case), 8 (queen's case).
#' @param return_input Logical if the original input data is returned.
#' @param simplify Logical if only the raster will be returned if \code{n_random = 1}
#' and \code{return_input = FALSE}.
#' @param verbose Logical if progress report is printed.
#'
#' @details
#' The function randomizes a habitat map with discrete classes (as SpatRaster) as proposed
#' by Harms et al. (2001) as “randomized-habitats procedure”. The algorithm starts with an
#' empty habitat map and starts to assign random neighbouring cells to each habitat
#' (in increasing order of abundance in observed map). We modified the procedure
#' slightly by increasing a probability to jump to a non-neighbouring cell as the
#' current patch becomes larger.
#'
#' In case the SpatRaster contains NA cells, this needs to be reflected in the observation
#' window of the point pattern as well (i.e., no point locations possible in these areas).
#'
#' @seealso
#' \code{\link{translate_raster}}
#'
#' @return rd_ras
#'
#' @examples
#' \dontrun{
#' landscape_classified <- classify_habitats(terra::rast(landscape), n = 5, style = "fisher")
#' landscape_random <- randomize_raster(landscape_classified, n_random = 19)
#' }
#'
#' @references
#' Harms, K.E., Condit, R., Hubbell, S.P., Foster, R.B., 2001. Habitat associations of
#' trees and shrubs in a 50-ha neotropical forest plot. Journal of Ecology 89, 947–959.
#' <https://doi.org/10.1111/j.1365-2745.2001.00615.x>
#'
#'@export
randomize_raster <- function(raster,
n_random = 1,
directions = 4,
return_input = TRUE,
simplify = FALSE,
verbose = TRUE){
# warning if NA are present
if (anyNA(terra::values(raster, mat = FALSE))) {
warning("NA values present. Please make sure the observation window of the point pattern reflects this.", call. = FALSE)
}
# check if n_random is >= 1
if (!n_random >= 1) {
stop("n_random must be >= 1.", call. = FALSE)
}
# set names of randomization randomized_1 ... randomized_n
names_randomization <- paste0("randomized_", seq_len(n_random))
habitats <- sort(table(terra::values(raster, mat = FALSE))) # get table of habitats
# print warning if more than 10 classes are present
if (length(habitats) > 10) {
warning("The raster has more than 10 classes. Please make sure discrete classes are provided.",
call. = FALSE)
}
n_cells <- sum(habitats) # number of cells
# create n_random rasters
result_list <- lapply(seq_len(n_random), function(current_raster) {
random_matrix <- terra::as.matrix(raster, wide = TRUE, na.rm = FALSE) # new raster without values
random_matrix[!is.na(random_matrix)] <- -999 # set all non-NAs to unique number
# loop through habitats but last one (all remaining cells will be last)
for (current_habitat in seq_len(length(habitats) - 1)) {
k <- 0 # counter since last jump
habitat_id <- as.numeric(names(habitats[current_habitat])) # get value of current habitat
random_cell <- sample(x = which(random_matrix == -999), size = 1) # random cell which is still -999
random_matrix[random_cell] <- habitat_id # assign habitat to cell
# loop until same number as number in original raster are assigned (break criterion within loop)
for (i in seq_len(habitats[current_habitat] - 1)) {
# increases as loop continious and increases prob to jump to non-neighbouring cell
ratio <- k / n_cells
r <- stats::runif(n = 1, min = 0, max = 1) # random number
# assign value to neighbouring patch
if (r >= ratio) {
# cells already assigned to habitat
cells_habitat <- which(random_matrix == habitat_id,
arr.ind = TRUE, useNames = FALSE)
# get neighbour cells
neighbours <- create_neighbourhood(cells = cells_habitat, matrix = random_matrix,
directions = directions)
# all neighbouring cells that are -999
empty_neighbours <- which(random_matrix[neighbours] == -999,
arr.ind = TRUE, useNames = FALSE)
# neighbours without habitat and inside plot present
if (length(empty_neighbours) > 0) {
# sample random neighbour
random_neighbour <- sample(x = empty_neighbours, size = 1)
# get matrix index of sampled neighbour
random_neighbour <- matrix(neighbours[random_neighbour, ], ncol = 2)
# assign cell to habitat
random_matrix[random_neighbour] <- habitat_id
k <- k + 1 # count since laste time jumped
# no neighbour with habitat and inside plot present
} else{
# random cell which is still -999
random_cell <- sample(x = which(random_matrix == -999), size = 1)
# assign habitat to cell
random_matrix[random_cell] <- habitat_id
k <- 0 # set counter since last jump zero
}
# jump to random starting cell
} else {
# random cell which is still -999
random_cell <- sample(x = which(random_matrix == -999), size = 1)
# assign habitat to cell
random_matrix[random_cell] <- habitat_id
k <- 0 # set counter since last jump zero
}
# print progess
if (verbose) {
message("\r> Progress: n_random: ", current_raster, "/", n_random,
" || habitats: " , current_habitat, "/", length(habitats), "\t\t",
appendLF = FALSE) # add habitat and number empty cells
}
}
}
# cells not assigned to any habitat yet
empty_cells <- which(random_matrix == -999, arr.ind = TRUE, useNames = FALSE)
# assign all still empty cells to last habitat
random_matrix[empty_cells] <- as.numeric(names(habitats[length(habitats)]))
# convert back to raster
random_raster <- terra::setValues(x = raster, values = random_matrix)
names(random_raster) <- "layer"
return(random_raster)
})
names(result_list) <- names_randomization
# combine to one list
randomization <- list(randomized = result_list, observed = raster,
method = "randomize_raster()")
# set class of result
class(randomization) <- "rd_ras"
# remove input if return_input = FALSE
if (!return_input) {
# set observed to NA
randomization$observed <- NA
# check if output should be simplified
if (simplify) {
# not possible if more than one raster is present
if (n_random > 1) {
warning("'simplify = TRUE' not possible for 'n_random > 1'.",
call. = FALSE)
# only one random raster is present that should be returend
} else if (n_random == 1) {
randomization <- randomization$randomized[[1]]
}
}
# return input if return_input = TRUE
} else {
# return warning if simply = TRUE because not possible if return_input = TRUE (only verbose = TRUE)
if (simplify) {
warning("'simplify = TRUE' not possible for 'return_input = TRUE'.", call. = FALSE)
}
}
# write result in new line if progress was printed
if (verbose) {
message("\r")
}
return(randomization)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.