R/spatialize_lsm.R

In landscapemetrics: Landscape Metrics for Categorical Map Patterns

#' spatialize_lsm
#'
#' @description Spatialize landscape metric values
#'
#' @param landscape Raster* Layer, Stack, Brick, SpatRaster (terra), stars, or a list of rasterLayers.
#' @param level Level of metrics. Either 'patch', 'class' or 'landscape' (or vector with combination).
#' @param metric Abbreviation of metrics (e.g. 'area').
#' @param name Full name of metrics (e.g. 'core area')
#' @param type Type according to FRAGSTATS grouping (e.g. 'aggregation metrics').
#' @param what Selected level of metrics: either "patch", "class" or "landscape".
#' It is also possible to specify functions as a vector of strings, e.g. `what = c("lsm_c_ca", "lsm_l_ta")`.
#' @param directions The number of directions in which patches should be
#' connected: 4 (rook's case) or 8 (queen's case).
#' @param progress Print progress report.
#' @param to_disk If TRUE raster will be saved to disk.
#' @param ... Arguments passed on to \code{calculate_lsm()}.
#'
#' @details
#' The functions returns a nested list with \code{RasterLayer}s. The first level
#' contains each input layer (only one element if \code{RasterLayer} was provided).
#' The second level contains a \code{RasterLayer} for each selected metric
#' (see \code{list_lsm} for details) where each cell has the landscape metric
#' value of the patch it belongs to. Only patch level metrics are allowed.
#'
#' @seealso
#' \code{\link{list_lsm}} \cr
#' \code{\link{show_lsm}}
#'
#' @return list
#'
#' @examples
#' spatialize_lsm(landscape, what = "lsm_p_area")
#'
#' @aliases spatialize_lsm
#'
#' @rdname spatialize_lsm
#'
#' @export
spatialize_lsm <- function(landscape,
                                level = "patch",
                                metric = NULL,
                                name = NULL,
                                type = NULL,
                                what = NULL,
                                directions = 8,
                                progress = FALSE,
                                to_disk = getOption("to_disk", default = FALSE),
                                ...) {

    landscape <- landscape_as_list(landscape)

    result <- lapply(X = seq_along(landscape), FUN = function(x) {

        if (progress) {

            cat("\r> Progress nlayers: ", x , "/", length(landscape))
        }

        spatialize_lsm_internal(landscape = landscape[[x]],
                                level = level,
                                metric = metric,
                                name = name,
                                type = type,
                                what = what,
                                directions = directions,
                                progress = FALSE,
                                to_disk = to_disk,
                                ...)
    })

    if (progress) {cat("\n")}

    names(result) <- paste0("layer_", 1:length(result))

    return(result)
}

spatialize_lsm_internal <- function(landscape,
                                    level, metric, name, type, what,
                                    directions,
                                    progress,
                                    to_disk,
                                    ...) {

    # get name of metrics
    metrics <- list_lsm(level = level,
                        metric = metric,
                        name = name,
                        type = type,
                        what = what,
                        simplify = TRUE,
                        verbose = FALSE)

    # how many metrics need to be calculated?
    number_metrics <- length(metrics)

    # error if no patch level metrics are provided
    if (!all(metrics %in% list_lsm(level = "patch", simplify = TRUE))) {
        stop("'spatialize_lsm()' only takes patch level metrics.",
             call. = FALSE)
    }

    # get CRS of input
    crs_input <- raster::crs(landscape)

    # get patches
    landscape_labeled <- get_patches(landscape,
                                     class = "all",
                                     directions = directions,
                                     to_disk = to_disk,
                                     return_raster = TRUE)[[1]]

    # continuous, unique patch id
    for (i in seq_len(length(landscape_labeled) - 1)) {

        max_id <- max(raster::values(landscape_labeled[[i]]), na.rm = TRUE)

        landscape_labeled[[i + 1]] <- landscape_labeled[[i + 1]] + max_id
    }

    # get dataframe with patch ID and coordinates to merge with result of metric
    # MH: Do we really want to remove NA?
    patches_tibble <- raster::as.data.frame(sum(raster::stack(landscape_labeled),
                                                na.rm = TRUE),
                                            xy = TRUE)

    # modify names
    names(patches_tibble) <- c("x", "y", "id")

    # replace all 0 values for NA
    patches_tibble$id <- replace(patches_tibble$id,
                                 patches_tibble$id == 0,
                                 NA)

    # create object for warning messages
    warning_messages <- character(0)

    # loop through metrics and return raster with value for each patch
    result <- withCallingHandlers(expr = {lapply(seq_along(metrics), function(x) {

        # print progress using the non-internal name
        if (progress) {

            cat("\r> Progress metrics: ", x, "/", number_metrics)
        }

        # get metric value
        fill_value <- calculate_lsm(landscape,
                                    what = metrics[[x]],
                                    progress = FALSE,
                                    ...)

        # merge with coords data frame
        fill_value <- merge(x = patches_tibble,
                            y = fill_value,
                            by = "id",
                            all.x = TRUE)

        if (to_disk) {

            # order fill_value by x
            index <- order(fill_value$x)
            fill_value <- fill_value[index, ]

            # split by y
            fill_value <- rev(split(x = fill_value, f = fill_value$y))

            # create empty raster
            result <- raster::raster(landscape)

            # get block size
            block_size <- raster::blockSize(result)

            # starting to write values in raster
            result <- raster::writeStart(x = result,
                                         filename = raster::rasterTmpFile(),
                                         overwrite = TRUE)

            # loop through all block sizes
            for (i in 1:block_size$n) {

                # start and end row of current block
                start_row <- block_size$row[i]
                end_row <- block_size$row[i] + (block_size$nrows[i] - 1)

                # get values of current rows and combine to df
                values_temp <- do.call("rbind", fill_value[start_row:end_row])

                # write current block
                result <- raster::writeValues(x = result,
                                              v = values_temp$value,
                                              start = block_size$row[i])
            }

            # close writing connections
            result <- raster::writeStop(result)

            return(result)
        } else {

            # convert to raster (wrap)
            result <- raster::rasterFromXYZ(fill_value[, c(2, 3, 8)],
                                            crs = crs_input)

            return(result)
        }})}, warning = function(cond) {

        warning_messages <<- c(warning_messages, conditionMessage(cond))

        invokeRestart("muffleWarning")}
    )

    # using metrics to name list
    names(result) <- metrics

    if (progress) {cat("\n")}

    # warnings present
    if (length(warning_messages) > 0) {

        # only unique warnings
        warning_messages <- unique(warning_messages)

        # print warnings
        lapply(warning_messages, function(x){ warning(x, call. = FALSE)})
    }

    return(result)
}