---
title: "National forest raster summaries by ecozone"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{National forest raster summaries by ecozone}
  %\VignetteDepends{dplyr, ggplot2, knitr, LandR, sf, terra}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.align = "center",
  out.width = "100%",
  dpi = 96,
  fig.retina = 1
)

## Suggested packages this vignette needs, plus the `--as-cran` "no suggests"
## check, which hides them from the build even when they are installed. Same
## guard the other PredictiveEcology packages use.
SuggestedPkgsNeeded <- c("dplyr", "knitr")
hasSuggests <- all(vapply(SuggestedPkgsNeeded, requireNamespace, logical(1), quietly = TRUE))
useSuggests <- !identical(tolower(Sys.getenv("_R_CHECK_DEPENDS_ONLY_")), "true")

manifest <- LandR:::.read_artifact_manifest()

## The manifest existing is not proof its files do -- check each one, so a
## partial checkout degrades to the notice below rather than erroring mid-build.
artifact_files <- if (is.null(manifest)) {
  character(0)
} else {
  vapply(
    c(manifest$CountsFile, manifest$DisplayFile, manifest$PolygonFile),
    function(f) system.file("extdata", f, package = "LandR"),
    character(1)
  )
}

have_artifacts <- length(artifact_files) == 3L && all(nzchar(artifact_files))

knitr::opts_chunk$set(eval = hasSuggests && useSuggests && have_artifacts)
```

```{r notice, echo = FALSE, results = "asis", eval = TRUE}
if (!hasSuggests || !useSuggests) {
  cat(
    "> **This vignette was built without its suggested packages**, so the results",
    "below are not shown. It needs", paste(SuggestedPkgsNeeded, collapse = ", "),
    "-- install them and rebuild to see the summaries.\n"
  )
} else if (!have_artifacts) {
  cat(
    "> **The precomputed artifacts are not present in this build**, so the",
    "results below are not shown. They are generated by",
    "`data-raw/precompute-vignettes.R`, which needs the national forest",
    "rasters. See *How this vignette is built*, below.\n"
  )
}
```

## 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^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`](https://github.com/PredictiveEcology/LandR/blob/development/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.

### Provenance of the committed artifacts

```{r provenance, echo = FALSE}
knitr::kable(
  data.frame(
    field = names(manifest),
    value = unlist(manifest, use.names = FALSE)
  ),
  col.names = c("", "")
)
```

## Loading the artifacts

```{r load}
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)
```

## 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.

```{r stats}
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

```{r histograms, fig.width = 9, fig.height = 7}
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.

```{r maps, results = "hide", message = FALSE, warning = FALSE}
## 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
)
```

```{r show-maps, echo = FALSE, out.width = "100%"}
## relative paths -- `include_graphics()` resolves them against the document
knitr::include_graphics(sort(list.files(fig_dir, pattern = "[.]png$", full.names = TRUE)))
```

## 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:

```sh
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:

```sh
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:

```{r terra-options, eval = FALSE}
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.