National forest raster summaries by ecozone

What this is

LandR provides three functions for summarizing a national forest raster over a set of polygons:

• calc_raster_counts() — the value-frequency table of raster values within each polygon;
• calc_raster_stats() — min, mean, max, quartiles and the proportion of zeroes, per polygon;
• plot_raster_stats() — a per-polygon figure combining a histogram, a map, and the statistics.

This vignette applies them to a national stand age map, summarized by ecozone.

How this vignette is built

The rasters these functions are designed for are large: SCANFI at 30 m is 2.1 × 10¹⁰ pixels, and a single pass reads tens of GB. That is not something a CI runner can do — and the ecostratification polygons are served from a host that throttles CI IP ranges — so this vignette is not computed from the full raster when the package is built.

Instead the work is cut at its natural seam. calc_raster_counts() is the only step that touches every pixel, and its output is small: a few thousand rows. That table is computed once on a workstation by data-raw/precompute-vignettes.R and committed, together with a coarsened raster used only for drawing maps and a simplified copy of the ecozone polygons.

Everything below — every statistic, histogram and map — is then recomputed from those artifacts when this vignette is built, in seconds and offline. So the numbers you are reading cannot drift away from the code in R/: they are produced by it, right now.

What can drift is the committed counts table, if calc_raster_counts() itself changes. That is guarded: the manifest records a hash of the normalized source of the functions that produced the table, and tests/testthat/test-vignette-artifacts.R fails when the two diverge, naming the script to re-run.

Note that the summary statistics are deliberately not part of that hash. They are an exact function of the counts table, so a change to, say, the quantile definition simply shows up in the next build — it must not force a multi-hour regeneration of a counts table that is still perfectly valid.

Loading the artifacts

library(LandR)

extdata <- function(f) system.file("extdata", f, package = "LandR")

## exact value-frequency table, computed from the full-resolution raster
counts <- utils::read.csv(extdata(manifest$CountsFile))

## coarsened raster -- used for maps only, never for statistics
display <- terra::rast(extdata(manifest$DisplayFile))

ecozones <- sf::st_read(extdata(manifest$PolygonFile), quiet = TRUE)

str(counts)
#> 'data.frame':    4493 obs. of  3 variables:
#>  $ ID   : chr  "Arctic Cordillera" "Arctic Cordillera" "Arctic Cordillera" "Arctic Cordillera" ...
#>  $ value: int  0 1 2 3 4 5 6 7 8 9 ...
#>  $ count: int  971 211 818 375 451 258 12 24 17 23 ...

Summary statistics

Every statistic is an exact function of the value-frequency table, so these are the true national numbers at full resolution — not an approximation from the coarsened raster used for the maps below.

stats <- calc_raster_stats(raster = NULL, counts_df = counts)

knitr::kable(stats, digits = 1, row.names = FALSE)

Note the quartiles use the inverse-ECDF definition (the smallest value whose cumulative proportion reaches p). For an integer raster with a large point mass at zero, an interpolated definition would report a non-zero median even where the majority of cells are zero.

Age distributions

library(ggplot2)

binned <- counts |>
  subset(!is.na(value)) |>
  transform(bin = floor(value / 10) * 10) |>
  aggregate(count ~ ID + bin, data = _, FUN = sum)

ggplot(binned, aes(x = bin, y = count)) +
  geom_col(fill = "steelblue") +
  facet_wrap(~ID, scales = "free_y", ncol = 4) +
  theme_minimal(base_size = 9) +
  labs(x = "stand age (years)", y = "pixel count")

Per-ecozone figures

plot_raster_stats() combines the histogram, a map and the statistics into one figure per polygon. Here it is driven from the precomputed tables — passing counts_df and stats_df skips the expensive recomputation — and from the coarsened display raster, so the whole set renders in seconds.

## Write the figures where knitr puts its own, rather than under `tempdir()`:
## pkgdown knits each article in a separate process and then only rewrites
## absolute `<img src>` paths relative to the site -- it never copies the files
## -- so a temp path is both unreachable and already deleted by then, and the
## site build dies with "Failed to realize ... no such file or directory".
## A path under `fig.path` is carried to the output by every knitr front end:
## base64-embedded into the self-contained vignette by `R CMD build`, and
## deployed with the article by pkgdown.
fig_dir <- paste0(knitr::opts_chunk$get("fig.path"), "ecozones")

plot_raster_stats(
  raster = display,
  polygons = ecozones,
  polygon_id = "ZONE_NAME",
  counts_df = counts,
  stats_df = stats,
  aggregate_factor = 1, ## already coarsened
  raster_label = "stand age (years)",
  inset_canada = TRUE,
  ## dissolving the ecozones gives the same national outline the inset would
  ## otherwise download from GADM, so this build stays offline
  inset_polygons = sf::st_union(ecozones),
  bin_width = 10,
  output_dir = fig_dir,
  ## web-sized rather than the print-sized 12 x 7.5 in @ 300 dpi default
  fig_width = 10,
  fig_height = 6.25,
  fig_dpi = 96
)

Reproducing this at full resolution

To run the same summary against the raster at native resolution — or against ecoprovinces or ecoregions instead of ecozones — see data-raw/precompute-vignettes.R. It is parameterized by environment variable:

LANDR_VIGNETTE_SOURCE=SCANFI LANDR_VIGNETTE_YEAR=2020 LANDR_VIGNETTE_VERSION=V2 \
  Rscript data-raw/precompute-vignettes.R

SCANFI is distributed through a Google Drive folder that prepInputsStandAgeMap() can only reach with credentials for it. If you already have a local copy, skip the download by pointing the script straight at it:

LANDR_VIGNETTE_RASTER=~/data/SCANFI/v2/2020/SCANFI_age_median_2020_v2_20260119.tif \
  Rscript data-raw/precompute-vignettes.R

Process large rasters on disk with bounded RAM, and point tempdir at fast local storage:

terra::terraOptions(memmax = 16, todisk = TRUE, tempdir = "/fast/local/scratch")

Prefer memmax over memfrac: memfrac is a fraction of free RAM rather than a hard cap, and memfrac = 0.0 crashes (terra’s block size collapses to zero). Never point tempdir at NFS, and never at a tmpfs such as /tmp on systems where that is RAM-backed — a raster “copied to disk” there really sits in RAM.