--- title: "Landscape Metric Covariates" author: "Bill Peterman" output: rmarkdown::html_vignette: number_sections: true toc: true vignette: > %\VignetteIndexEntry{Landscape Metric Covariates} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} library(multiScaleR) library(terra) knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 6.5, fig.height = 4.25, eval = FALSE, warning = FALSE, message = FALSE ) pkg_extdata <- function(...) { src_path <- normalizePath(file.path("..", "inst", "extdata", ...), winslash = "/", mustWork = FALSE) if (file.exists(src_path)) { return(src_path) } inst_path <- system.file("extdata", ..., package = "multiScaleR") if (nzchar(inst_path) && file.exists(inst_path)) { return(inst_path) } stop("Could not locate required extdata file: ", paste(..., collapse = "/")) } vignette_cache <- readRDS(pkg_extdata("vignette_cache.rds")) ``` # Why explicit covariate specifications? Traditional `multiScaleR` workflows assume that each raster layer becomes one optimized, kernel-weighted covariate. That remains the default behavior. For landscape metrics, however, the relationship between raster layers and model covariates is no longer one-to-one. For example, the same binary forest raster might be used to calculate: - a distance-weighted proportion of forest cover, - edge density within a fixed buffer, and - edge density at an optimized scale. Those are different covariates, even though they come from the same source raster. The `msr_vars()` interface lets you name those derived covariates explicitly. --- # Landscape metric reference `landscape_var()` supports 25 metrics grouped into five families. All metrics are computed within a circular buffer centered on each sample point (or on each raster cell when projecting with `kernel_scale.raster()`). Metrics that operate on categorical rasters expect integer class values; continuous rasters should be classified before use. Most metrics describe the whole landscape within the buffer. Three are *class-level*: they describe a single category and require a `class` argument naming the focal class (covered in [Class-level metrics](#class-level-metrics)). All of the categorical metrics reproduce the values returned by the `landscapemetrics` package, which itself follows the FRAGSTATS conventions. ## Composition metrics These metrics describe the diversity and abundance of patch types within the buffer. They require a categorical (integer-valued) raster. | Code | Name | Description | |------|------|-------------| | `shdi` | Shannon Diversity Index | $-\sum p_i \ln p_i$, where $p_i$ is the proportion of class $i$. Increases with richness and evenness. The `base` argument controls the logarithm base (default: natural log). | | `shei` | Shannon Evenness Index | SHDI / ln(richness). Ranges from 0 (one dominant class) to 1 (perfectly even). | | `sidi` | Simpson Diversity Index | $1 - \sum p_i^2$. Probability that two randomly selected cells belong to different classes. Ranges from 0 to 1. | | `siei` | Simpson Evenness Index | SIDI / (1 - 1/richness). Standardizes SIDI against maximum possible diversity. Requires at least 2 classes and can be uninformative for optimized scales when local evenness is uniform across sample points. | | `msidi` | Modified Simpson Diversity | $-\ln \sum p_i^2$. Logarithmic transformation of Simpson concentration; unbounded above. | | `msiei` | Modified Simpson Evenness | MSIDI / ln(richness). Evenness analogue of MSIDI. Can be uninformative for optimized scales when local evenness is uniform across sample points. | | `pr` | Patch Richness | Count of distinct classes present in the buffer. | | `prd` | Patch Richness Density | PR per 100 ha of buffer area. | | `rpr` | Relative Patch Richness | (PR / `classes_max`) × 100. Requires the `classes_max` argument. | | `ta` | Total Area | Total area of the buffer in hectares. | ## Edge metrics These metrics quantify the amount and complexity of boundaries between classes. They require a categorical raster and use cell adjacency information. | Code | Name | Description | |------|------|-------------| | `ed` | Edge Density | Total internal edge length (m) per hectare of buffer area. Higher values indicate a more fragmented landscape. | | `te` | Total Edge | Total internal edge length in metres. | | `lsi` | Landscape Shape Index | Total perimeter / minimum possible perimeter for the same number of cells. LSI = 1 for a perfectly compact (circular) patch; increases with irregularity. | ## Adjacency metrics These metrics describe how often similar or dissimilar classes occur next to one another. They require a categorical raster and use cell-pair counts (the class-by-class co-occurrence matrix). | Code | Name | Description | |------|------|-------------| | `ai` | Aggregation Index | Area-weighted aggregation of like adjacencies within classes. Higher values indicate that classes are more aggregated. | | `pladj` | Proportion of Like Adjacencies | Percentage of cell-pair boundaries shared between cells of the *same* class. High values indicate a clumped or aggregated landscape. | | `contag` | Contagion | Measures the degree of clumping across all class pairs. Ranges from 0 (maximum interspersion) to 100 (single class). Requires at least 2 classes. | | `iji` | Interspersion & Juxtaposition | Evenness of adjacencies *among different classes* (like adjacencies are ignored). Ranges from near 0 (one class pair dominates the between-class boundaries) to 100 (every class is equally adjacent to every other). Requires at least 3 classes and enough between-class adjacencies in each buffer; otherwise it returns `NA`. | ## Information theory metrics These metrics come from treating the class-pair co-occurrence matrix as a bivariate probability distribution, following Nowosad and Stepinski (2019). Their value is that they separate *composition* (how many classes, how evenly mixed) from *configuration* (how the classes are arranged), which the diversity and contagion metrics confound. They require a categorical raster. The `base` argument sets the logarithm base; use `base = 2` for bits, matching `landscapemetrics`. | Code | Name | Description | |------|------|-------------| | `ent` | Marginal entropy, H(x) | Diversity of the cell values entering adjacencies. Increases with the number and evenness of classes; this is the composition signal. | | `joinent` | Joint entropy, H(x, y) | Diversity of the class-pair adjacencies themselves: the total complexity of the pattern. | | `condent` | Conditional entropy, H(y \| x) | Configurational complexity remaining once composition is accounted for: the uncertainty about a neighbor's class given the focal cell's class. | | `mutinf` | Mutual information, I(y, x) | How much a cell's class tells you about its neighbor's. Rises as like classes cluster; this is configuration separated from composition. | | `relmutinf` | Relative mutual information | `mutinf / ent`, rescaled to 0–1 so configuration can be compared across landscapes of differing composition. Equal to 1 when `mutinf` is 0. Unaffected by `base`. | ## Class-level metrics {#class-level-metrics} The metrics above summarize the whole landscape within the buffer. These three instead describe a single focal class, supplied through the `class` argument of `landscape_var()` (an integer matching one of the raster's class codes). Asking for one of these metrics without a `class`, or supplying `class` to any other metric, is an error. Where the focal class is absent from a buffer, `pland` and `ca` return 0 and `clumpy` returns `NA`. | Code | Name | Description | |------|------|-------------| | `clumpy` | Clumpiness Index | Aggregation of the focal class, corrected for its own abundance. Ranges from -1 (maximally disaggregated) through 0 (random) to 1 (maximally clumped). | | `pland` | Percentage of Landscape | The focal class's share of the valid buffer cells, as a percentage. | | `ca` | Class Area | Area of the focal class within the buffer, in hectares. | ```{r class-level-example} ## Clumpiness of land-cover class 2, and that class's percentage cover, ## both within a fixed 500 m radius class_vars <- msr_vars( cover2_clumpy = landscape_var("landcover", metric = "clumpy", radius = 500, class = 2), cover2_pland = landscape_var("landcover", metric = "pland", radius = 500, class = 2) ) class_vars ``` --- # Load example data This example uses the simulated count data and survey locations included with `multiScaleR`. The data contain 100 spatial sample points, which is large enough to demonstrate the workflow without making the vignette expensive to build. ```{r load-data} data("landscape_counts") dat <- landscape_counts data("surv_pts") pts <- vect(surv_pts) land_rast <- rast(pkg_extdata("landscape.tif")) ``` The raster layer `land1` is a binary habitat layer. To demonstrate a categorical landscape-diversity metric, we also create a simple three-class land-cover layer from the continuous `land2` raster. In a real analysis, this would usually be a mapped categorical land-cover product. ```{r build-landcover} landcover <- classify( land_rast$land2, rcl = matrix( c(-Inf, -0.5, 1, -0.5, 0.5, 2, 0.5, Inf, 3), ncol = 3, byrow = TRUE ) ) names(landcover) <- "landcover" metric_rasters <- c(land_rast$land1, landcover) ``` ```{r plot-inputs, fig.cap = "***Source rasters used to derive model covariates.***"} plot(metric_rasters) plot(metric_rasters$land1, main = "land1 with survey points") points(pts, pch = 19, cex = 0.6, col = "red") ``` --- # Define derived covariates Use `msr_vars()` to define model covariates as transformations of source raster layers. Here, `land1_prop` is optimized as a kernel-weighted covariate, while `land1_ed_500` and `landcover_shdi_500` are fixed-radius landscape metrics calculated within 500 m of each point. ```{r define-vars} scale_vars <- msr_vars( land1_prop = kernel_var("land1"), land1_ed_500 = landscape_var("land1", metric = "ed", radius = 500), landcover_shdi_500 = landscape_var("landcover", metric = "shdi", radius = 500) ) scale_vars ``` This is the key distinction: - `land1` and `landcover` are source raster layers. - `land1_prop`, `land1_ed_500`, and `landcover_shdi_500` are model covariates. - Only `land1_prop` receives an optimized sigma parameter in this example. --- # Prepare model data Pass the covariate specification to `kernel_prep()`. The returned `kernel_dat` data frame contains the derived covariates, not merely the original raster layer names. ```{r kernel-prep} kernel_inputs <- kernel_prep( pts = pts, raster_stack = metric_rasters, max_D = 1200, kernel = "gaussian", scale_vars = scale_vars, verbose = FALSE ) head(kernel_inputs$kernel_dat) ``` Combine the derived covariates with the response data and fit the starting model. ```{r fit-model} df <- data.frame(dat, kernel_inputs$kernel_dat) mod <- glm( counts ~ site + land1_prop + land1_ed_500 + landcover_shdi_500, family = poisson(), data = df ) ``` --- # Optimize the kernel covariate Because the two landscape metrics were defined with fixed radii, only `land1_prop` is optimized. This keeps the example fast while still showing how kernel-scaled and landscape-metric covariates can be included in the same model. ```{r optimize} opt <- multiScale_optim( fitted_mod = mod, kernel_inputs = kernel_inputs, verbose = FALSE ) ``` ```{r summary, eval = TRUE} opt <- vignette_cache$landscape$opt summary(opt) ``` --- # Project derived covariates to rasters `kernel_scale.raster()` also understands `scale_vars`. When a fitted `multiScaleR` object is supplied, optimized covariates use their fitted scale and fixed-radius landscape metrics use their specified radius. ```{r project, fig.cap = "***Projected derived covariates.***", eval = TRUE} projected <- terra::unwrap(vignette_cache$landscape$projected) plot(projected) ``` The resulting raster names match the model covariates: ```{r projected-names, eval = TRUE} names(projected) ``` --- # Producing landscape metric rasters without a fitted model `kernel_scale.raster()` can generate landscape metric raster surfaces directly, without having fitted a model at all. This is useful when you want to explore how a metric varies across the landscape at a given scale, compare metrics side-by-side before committing to a model structure, or simply produce a covariate map from a classified raster. Pass a `scale_vars` object in which all covariates have a fixed `radius` (i.e., `radius` is not `NULL` in `landscape_var()`). No `sigma` argument is needed because no scale is being optimized. ```{r standalone-metric-rasters, fig.cap = "***Landscape metric surfaces computed directly from the source rasters at a fixed 500 m radius.***", fig.height = 7, eval = FALSE} ## Define several metrics at a fixed 500 m radius standalone_vars <- msr_vars( land1_ai = landscape_var("land1", metric = "ai", radius = 500), land1_ed = landscape_var("land1", metric = "ed", radius = 500), land1_lsi = landscape_var("land1", metric = "lsi", radius = 500), land1_pladj = landscape_var("land1", metric = "pladj", radius = 500), cover_shdi = landscape_var("landcover", metric = "shdi", radius = 500), cover_pr = landscape_var("landcover", metric = "pr", radius = 500), cover_sidi = landscape_var("landcover", metric = "sidi", radius = 500) ) ## Apply directly: no multiScaleR object or sigma values required metric_surfaces <- kernel_scale.raster( raster_stack = metric_rasters, scale_vars = standalone_vars, verbose = FALSE ) plot(metric_surfaces) ``` Each output layer is named after the covariate defined in `msr_vars()`: ```{r standalone-names, eval = FALSE} names(metric_surfaces) ``` You can also compute a single metric in one line without defining a `msr_vars()` object. Simply wrap a single `landscape_var()` call: ```{r standalone-single, fig.cap = "***Edge density within a 300 m radius.***", eval = FALSE} ed_300 <- kernel_scale.raster( raster_stack = metric_rasters["land1"], scale_vars = msr_vars(ed_300m = landscape_var("land1", metric = "ed", radius = 300)), verbose = FALSE ) plot(ed_300, main = "Edge density (m/ha) within 300 m") ``` Different radii produce notably different surfaces for the same metric. The example below computes edge density at three scales and plots them together for comparison: ```{r standalone-multiscale, fig.cap = "***Edge density at three spatial scales.***", fig.height = 3.5, eval = FALSE} ed_multi <- kernel_scale.raster( raster_stack = metric_rasters["land1"], scale_vars = msr_vars( ed_200m = landscape_var("land1", metric = "ed", radius = 200), ed_500m = landscape_var("land1", metric = "ed", radius = 500), ed_900m = landscape_var("land1", metric = "ed", radius = 900) ), verbose = FALSE ) plot(ed_multi, main = c("ED: 200 m", "ED: 500 m", "ED: 900 m")) ``` ## Separating composition from configuration The information-theory metrics are most useful when you want to ask whether an effect is driven by *what* is in the neighborhood (its composition) or by *how* it is arranged (its configuration). The surfaces below project three metrics for the three-class `landcover` layer at a fixed 500 m radius. ```{r info-theory-surfaces, fig.cap = "***Composition and configuration surfaces for the same landscape.***", fig.height = 3.5, eval = FALSE} info_surfaces <- kernel_scale.raster( raster_stack = metric_rasters["landcover"], scale_vars = msr_vars( H_x = landscape_var("landcover", metric = "ent", radius = 500, base = 2), MI = landscape_var("landcover", metric = "mutinf", radius = 500, base = 2), relMI = landscape_var("landcover", metric = "relmutinf", radius = 500) ), verbose = FALSE ) plot(info_surfaces, main = c("Marginal entropy (bits)", "Mutual information (bits)", "Relative mutual information")) ``` **How to read these surfaces:** marginal entropy (`ent`, here in bits because `base = 2`) is a composition measure: it is high where the buffer holds several classes in roughly equal proportion, and near 0 where one class dominates. Mutual information (`mutinf`) is a configuration measure: it is high where like classes cluster into contiguous patches and low where the classes are finely interspersed, *independent* of how many classes are present. Relative mutual information (`relmutinf`) divides the two, putting configuration on a 0 to 1 scale so it is comparable between neighborhoods that differ in composition. Two buffers can share the same `ent` yet have very different `mutinf`: same ingredients, different spatial arrangement. That separation is the reason to prefer these metrics over contagion when arrangement is the hypothesis of interest. For a single class, `clumpy` answers the same configuration question on the -1 to 1 scale familiar from FRAGSTATS, while `pland` reports that class's cover: ```{r clumpy-surface, fig.cap = "***Clumpiness and cover of land-cover class 2.***", fig.height = 3.5, eval = FALSE} class2_surfaces <- kernel_scale.raster( raster_stack = metric_rasters["landcover"], scale_vars = msr_vars( clumpy2 = landscape_var("landcover", metric = "clumpy", radius = 500, class = 2), pland2 = landscape_var("landcover", metric = "pland", radius = 500, class = 2) ), verbose = FALSE ) plot(class2_surfaces, main = c("CLUMPY (class 2)", "PLAND (class 2, %)")) ``` --- # Optimizing landscape metric radii To optimize the scale of a landscape metric, omit `radius` from `landscape_var()`. For example: ```{r optimized-landscape-vars, eval = FALSE} scale_vars_opt <- msr_vars( land1_prop = kernel_var("land1"), land1_ed = landscape_var("land1", metric = "ed"), landcover_shdi = landscape_var("landcover", metric = "shdi") ) ``` This is statistically meaningful, but it is more computationally expensive than optimizing a kernel-weighted mean. At each candidate radius, the landscape metric must be recalculated for every sample point. For moderate or large data sets, consider first testing the workflow with a subset of points or using precomputed results in reproducible documents. Before optimizing `iji`, `siei`, or `msiei`, inspect a pilot `kernel_prep()` result or compare a few fixed radii. These metrics can be all `NA` or nearly constant when buffers contain too few classes, too few between-class adjacencies, or nearly identical evenness at all sample points. In that case, use a larger radius, increase `max_D`, or choose a more stable diversity, entropy, or edge metric. --- # Practical guidance Use fixed-radius landscape metrics when you have a biologically motivated buffer distance or when the metric is expensive to recompute. Use optimized landscape metric radii when the scale itself is a central question and computation time is acceptable. At present, optimized landscape metrics are not supported with `kernel = "expow"` because the exponential-power shape parameter does not apply cleanly to fixed-buffer landscape metrics. --- # See also - `vignette("surface_metric_covariates", package = "multiScaleR")` for the continuous-surface analogue of these metrics (roughness, gradient, and texture of a continuous raster). - `vignette("multiScaleR_Guide", package = "multiScaleR")` for the core kernel-weighting workflow. - The help page `?msr_vars` lists every supported metric with its definition. # References Nowosad, J., & Stepinski, T. F. (2019). Information theory as a consistent framework for quantification and classification of landscape patterns. *Landscape Ecology*, 34(9), 2091–2101.