From f5faf715c5ce8d4a6b0e6bcf7bcf0cbc52bcb088 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Mon, 28 Sep 2026 18:10:57 -0700 Subject: [PATCH 01/10] Initialize PWF baseline for #81 Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- planning/active/.gitkeep | 0 planning/active/findings.md | 45 +++++++++++++++++++++ planning/active/progress.md | 8 ++++ planning/active/task_plan.md | 77 ++++++++++++++++++++++++++++++++++++ 4 files changed, 130 insertions(+) delete mode 100644 planning/active/.gitkeep create mode 100644 planning/active/findings.md create mode 100644 planning/active/progress.md create mode 100644 planning/active/task_plan.md diff --git a/planning/active/.gitkeep b/planning/active/.gitkeep deleted file mode 100644 index e69de29..0000000 diff --git a/planning/active/findings.md b/planning/active/findings.md new file mode 100644 index 0000000..ee08f8d --- /dev/null +++ b/planning/active/findings.md @@ -0,0 +1,45 @@ +# Findings — Accuracy assessment for change maps (#81) + +## Issue context + +## Problem + +drift produces land-cover change maps (`dft_rast_transition()`, `dft_rast_break_class()`, `dft_transition_vectors()`), but it has no way to say **how right they are**. Every area it reports is the mapped area, and mapped area is a biased estimator whenever the map has errors. For a change map it is usually badly biased, because a single wrong label on either date manufactures a transition. + +The standard remedy is a stratified sample of reference labels and estimators that correct area for the measured error. That method is generic to any classified map, so it belongs in drift, not in each driver. floodplains#93 is the first consumer: NECR floodplain change, fire/harvest attribution, and wetland change. + +## Scope + +Follow Olofsson et al. 2014, *Good practices for estimating area and assessing accuracy of land change* (Remote Sensing of Environment 148:42–57), and Stehman's stratified estimators it builds on. + +1. **Stratified sample design.** Given a strata raster or polygons (the map classes, or caller-defined strata such as "change attributed / unattributed / stable"), draw random **points** per stratum. Sample points rather than patches: patch sampling over-weights large patches, and area is the quantity being estimated. + - Allocation can be equal per stratum or a caller-supplied `n`. + - The output is a reproducible sample: seeded, with a stable point id, the stratum, and each stratum's mapped area (weights) recorded alongside. + - A helper to size the sample from a pilot's variance, or from target user's accuracies. +2. **Accuracy and error-adjusted area from labels.** Given the sample with its reference label and the stratum weights, return: + - an area-weighted error matrix + - overall, user's and producer's accuracy, with standard errors + - **error-adjusted area per class, with confidence intervals** +3. **A tidy label contract.** Document the columns a reference-label table must carry (point id, stratum, map class, reference class, and optionally confidence and reviewer), so any review tool can feed (2). drift does not store labels; the caller does. +4. **Train and test separation.** If labels are ever reused to train a classifier, the accuracy sample must be held out. Provide a documented split, or at minimum a flag that makes the estimator refuse training points. + +Out of scope: the review UI (drift#79 covers the imagery and the map layers), and choosing strata or storing labels, which are the driver's job. + +## Acceptance + +- Estimators reproduce the worked example in Olofsson et al. 2014 (the published error matrix and its error-adjusted areas with CIs) to the published precision. That is the external reference; a self-generated fixture would only test the code against itself. +- A must-fail test: dropping the stratum weights (unweighted accuracy) gives a different, wrong answer on the paper's example. +- Sample draws are reproducible from a seed and stable across terra versions for the same inputs. + +Relates: floodplains#93, drift#79 + +## Plan-gate decisions (2026-09-28) + +- Reference values come from the PDFs. The user is adding Olofsson 2014 and Stehman 2014 to Zotero; neither was in the library (checked zotero.sqlite read-only). The Zotero MCP fails on missing `ZOTERO_LIBRARY_ID` / `ZOTERO_API_KEY`. +- The estimator takes the general Stehman 2014 form (strata may differ from map classes), because floodplains#93 strata are not map classes. +- No existing implementation in the org: swept the exports of 19 NGE packages plus `gh search code org:NewGraphEnvironment olofsson`. `mapaccuracy`, `survey` and `sampling` are not installed and not needed. + +## Errors Encountered + +| Error | Resolution | +|-------|------------| diff --git a/planning/active/progress.md b/planning/active/progress.md new file mode 100644 index 0000000..873181b --- /dev/null +++ b/planning/active/progress.md @@ -0,0 +1,8 @@ +# Progress — Accuracy assessment for change maps (#81) + +## Session 2026-09-28 + +- Plan-mode exploration; phases approved by user +- Created branch `81-accuracy-assessment-for-change-maps-stra` off main +- Scaffolded PWF baseline from issue #81 with approved phases +- Next: Phase 2 estimator code and tests; published-value pins wait on the PDFs (Phase 1) diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md new file mode 100644 index 0000000..4c92e32 --- /dev/null +++ b/planning/active/task_plan.md @@ -0,0 +1,77 @@ +# Task: Accuracy assessment for change maps: stratified sampling and error-adjusted area with CIs (Olofsson et al. 2014) (#81) + + +drift produces land-cover change maps (`dft_rast_transition()`, `dft_rast_break_class()`, `dft_transition_vectors()`), but it has no way to say **how right they are**. Every area it reports is the mapped area, and mapped area is a biased estimator whenever the map has errors. For a change map it is usually badly biased, because a single wrong label on either date manufactures a transition. + +The standard remedy is a stratified sample of reference labels and estimators that correct area for the measured error. That method is generic to any classified map, so it belongs in drift, not in each driver. floodplains#93 is the first consumer: NECR floodplain change, fire/harvest attribution, and wetland change. + +## Design + +A new `dft_accuracy_*` family, one function per file (noun-first, so the family autocompletes together): + +| function | input → output | +|---|---| +| `dft_accuracy_sample(strata, n, seed, n_min)` | strata SpatRaster (integer codes) → `list(points = sf, strata = tibble)` | +| `dft_accuracy_size(weights, ua, se_target, …)` | stratum weights + target / pilot user's accuracies → total n, plus allocation | +| `dft_accuracy_estimate(labels, strata, …)` | label table + `$strata` → `list(matrix, accuracy, area)` with SEs and CIs | +| `dft_accuracy_labels(labels)` | validates the label contract; returns it invisibly or errors naming the fault | + +**Sampler, stable across terra versions.** `terra::spatSample()` is not stable across versions, so it is not used. Instead: +1. Count cells per stratum with `terra::freq()` (C++, streaming). +2. Draw within-stratum ranks with base `sample.int(N_h, n_h)` under a local seed. The caller's RNG state is saved and restored, and `sample.kind = "Rejection"` is set explicitly. +3. Resolve each rank to a cell id in one block-wise pass (`terra::readStart` / `readValues` by row block, cumulative counts per stratum). + +The result depends only on cell order, which is row-major and fixed, and on base R's RNG. Memory stays bounded at BULK's 169M cells. Points are cell centres (`terra::xyFromCell`) and carry `point_id` (stable, zero-padded, stratum-then-rank order), `stratum`, `cell`, and `x`/`y`. `$strata` carries `stratum`, `n_cells`, `area` (ha), `weight` (`n_cells / Σ n_cells`) and `n`. Other details: +- Allocation is either equal per stratum (`n` as a scalar) or caller-supplied (a named vector). `n_h > N_h` is refused. +- `dft_check_crs()` guards the area. +- `NA` cells are outside the population. +- Raster strata only. Polygon strata are rasterised first with `terra::rasterize()` onto the map grid, documented in the roxygen. It is one line, and it keeps the grid (and so the weights) the caller's explicit choice. + +**Label contract** (columns the caller's table must carry; drift stores nothing): +- required: `point_id`, `stratum`, `map_class`, `ref_class` +- optional: `confidence`, `reviewer`, `use` (`"accuracy"` / `"training"`) + +**Train/test.** `dft_accuracy_estimate()` refuses any row with `use == "training"` and names the count. It has no silent drop: dropping those rows would bias the design, because the training rows were not drawn as part of the probability sample. + +**Estimators** (Stehman 2014 general form: per-stratum sample means of indicator variables, weighted by `N_h`): +- An area-weighted error matrix of estimated proportions p̂_ij. +- Overall, user's and producer's accuracy. UA and PA are ratio estimators, with their variances from Stehman eq. 26–28. +- Error-adjusted area per reference class = `A_total · p̂_·j`, with its SE and a 95% CI (`z = 1.96`, as the argument `level`). +- When strata = map classes, these are algebraically Olofsson eq. 2–11. A test asserts that equality on the paper's table. + +**Sizing.** Olofsson eq. 13: `n = (Σ W_i S_i / S(Ô))²` with `S_i = sqrt(U_i(1−U_i))`. Two allocation helpers: `"equal"`, and `"proportional_min"` (proportional with `n_min` per rare stratum, Olofsson §5.1.1). Pilot variance comes in as `ua` taken from a pilot's `dft_accuracy_estimate()$accuracy`, documented with an example. + +### Phase 1: Reference values in hand +- [ ] User adds the Olofsson 2014 and Stehman 2014 PDFs (decided at the plan gate). Until then, Phases 2–4 code and non-published tests proceed, and the published-value pins wait. +- [ ] Transcribe the worked examples into `tests/testthat/helper-accuracy.R`, each value cited to its page and table number +- [ ] `findings.md`: the estimator equations with numbers, and a check of the Olofsson example by hand arithmetic (deforestation 21,158 ha is reproducible from the row counts; confirm against the PDF) + +### Phase 2: Estimator (tests first) +- [ ] `test-dft_accuracy_estimate.R`: Olofsson Table 8/9 (the error matrix, UA/PA/OA with SEs, and the adjusted areas with 95% CIs, to the published precision); the Stehman 2014 example (strata ≠ map classes); a must-fail test showing unweighted accuracy differs from the published values on the Olofsson example; `use == "training"` refused; a stratum with n_h = 1 (variance undefined) refused; a class absent from the reference handled +- [ ] `R/dft_accuracy_estimate.R` + `R/dft_accuracy_labels.R` (the contract validator, which the estimator calls) +- [ ] Restore-the-bug check: remove the weights from the estimator and confirm the published-value tests go red + +### Phase 3: Sampler (tests first) +- [ ] `test-dft_accuracy_sample.R`: pinned golden cell ids for a seeded draw on `example_2017.tif`; the same seed gives an identical draw; the caller's `.Random.seed` is untouched; within-stratum uniformity (chi-square over a large draw); weights sum to 1 and areas match `dft_rast_summarize()`; `n_h > N_h`, a lonlat raster and a non-integer raster are refused; NA cells are never drawn; the block-wise resolver agrees with a brute-force `which()` on the small tile +- [ ] `R/dft_accuracy_sample.R` +- [ ] Scale test on BULK `classified_2017.tif` (169M cells) with the RSS sampler; record the time and peak RSS in the PR body + +### Phase 4: Sizing +- [ ] `test-dft_accuracy_size.R`: reproduce Olofsson §5.1.1's sample-size example (n and allocation) from the paper; edge cases (UA = 1, a zero weight) +- [ ] `R/dft_accuracy_size.R` + +### Phase 5: Docs and release +- [ ] A runnable `@examples` block on every function (the estimator example uses Olofsson's published counts; the sampler uses the bundled tile) +- [ ] `devtools::document()`, `lintr::lint_package()`, `pkgdown::check_pkgdown()` (the new exports go in the reference index) +- [ ] NEWS.md, then version 0.19.0 as the final commit +- [ ] CLAUDE.md Core Pipeline: add the accuracy block +- [ ] Update the floodplains#93 body's "What lives where" table if the delivered API differs from what it assumes + +No vignette. One made with fabricated labels would illustrate a number nobody measured. The worked example belongs in floodplains#93 once real labels exist. + +## Validation + +- [ ] Tests pass +- [ ] `/code-check` clean on each commit +- [ ] PWF checkboxes match landed work +- [ ] `/planning-archive` on completion From 44d1564a07e7ce7b56d30c7b4c41afe19312c9de Mon Sep 17 00:00:00 2001 From: almac2022 Date: Mon, 28 Sep 2026 18:17:58 -0700 Subject: [PATCH 02/10] Fold plan review into #81 PWF Two blockers. freq() and readValues() disagree on factor strata (confirmed by probe). Sizing assumed strata equal map classes. Also: per-stratum seeds so a pilot extends, a census for tiny strata under an fpc argument, map extraction at sampled cells, and a census oracle that needs no PDFs. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- planning/active/review-plan.md | 27 +++++++++ planning/active/task_plan.md | 104 +++++++++++++++++++++------------ 2 files changed, 95 insertions(+), 36 deletions(-) create mode 100644 planning/active/review-plan.md diff --git a/planning/active/review-plan.md b/planning/active/review-plan.md new file mode 100644 index 0000000..7767fbc --- /dev/null +++ b/planning/active/review-plan.md @@ -0,0 +1,27 @@ +# Plan review — #81 (Plan agent, 2026-09-28) + +Read-only reviewer, so its findings came back as reply text; I transcribed them here. Every finding is folded into `task_plan.md` except where the disposition says otherwise. + +| id | finding | disposition | +|---|---|---| +| B1 | `terra::freq()` returns labels on factor strata while `readValues()` returns codes; `freq()` also rounds floats. Count and resolve steps could not join | **Confirmed by probe** (`Water -> Water` vs `1001`). Two passes on one chunked reader; `stratum_label` from `levels()` | +| B2 | Sizing through `ua` assumes strata = map classes; #93 needs per-stratum SD for an area target | Primary form takes `s_h` from the estimator's `$stratum`; `ua` kept as a convenience | +| B3 | FPC undecided; `n_h > N_h` refusal breaks on tiny strata (the tile has 2- and 3-cell strata) | `fpc = TRUE` by default, and the Olofsson pin runs with `fpc = FALSE`; a stratum with `n_h ≥ N_h` becomes a census, with a message | +| G1 | The sampler does not emit `map_class` | `map =` argument, extracted at each cell | +| G2 | One `blocks()` chunk on the tile, so the carry logic is untested | Explicit row chunks, capped at about 1e7 cells; chunk-invariance test | +| G3 | RNG save/restore under-specified (L'Ecuyer, absent seed) | Pin all three kinds; restore, or remove; three tests | +| G4 | A single seed across strata cannot extend a pilot | Per-stratum sub-seeds; `point_id` in draw order, independent of n. Prefix stability **confirmed by probe** | +| G5 | Nonresponse and coverage rules | Refusals listed in the contract | +| G6 | Change targets and union SEs | Composition recipe and "recode, then estimate" in the roxygen, plus a test | +| G7 | Per-stratum reporting | `$stratum` output | +| G8 | Estimator edge cases | Listed in the Phase 2 tests | +| A1 | Equation and table numbers were from memory | Verified against the PDFs in Phase 1 | +| A2 | `level` vs `z`, and relative tolerance | `qnorm`, and absolute tolerance on the pins | +| A3 | The train/test rationale was backwards | Rewritten | +| A4 | `x`/`y` duplicate the geometry; `cell` is grid-bound | `x`/`y` dropped; the grid identity is recorded in `$design` | +| A5 | Polygon rasterise traps (background written as 0, mask to footprint) | In the roxygen recipe | +| O1 | A census oracle validates the strata ≠ map-classes path without the PDFs | Added to Phase 2 | +| O2 | Freeze the output shapes in Phase 2 | Added | +| O3 | `_pkgdown.yml` has no reference index | Noted; out of scope | +| AC1–4 | The must-fail test must run package code; `rep()` expansion; BULK on the factor path | Added | +| S2 | A provenance record | `$design` | diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index 4c92e32..b783131 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -7,65 +7,97 @@ The standard remedy is a stratified sample of reference labels and estimators th ## Design -A new `dft_accuracy_*` family, one function per file (noun-first, so the family autocompletes together): +Revised 2026-09-28 after the plan review (`planning/active/review-plan.md`). The review found two blockers, both confirmed or folded in here: B1 (the `freq()` step reports labels while `readValues()` reports codes on factor strata; measured) and B2 (sizing assumed strata = map classes). + +A new `dft_accuracy_*` family, one function per file: | function | input → output | |---|---| -| `dft_accuracy_sample(strata, n, seed, n_min)` | strata SpatRaster (integer codes) → `list(points = sf, strata = tibble)` | -| `dft_accuracy_size(weights, ua, se_target, …)` | stratum weights + target / pilot user's accuracies → total n, plus allocation | -| `dft_accuracy_estimate(labels, strata, …)` | label table + `$strata` → `list(matrix, accuracy, area)` with SEs and CIs | -| `dft_accuracy_labels(labels)` | validates the label contract; returns it invisibly or errors naming the fault | - -**Sampler, stable across terra versions.** `terra::spatSample()` is not stable across versions, so it is not used. Instead: -1. Count cells per stratum with `terra::freq()` (C++, streaming). -2. Draw within-stratum ranks with base `sample.int(N_h, n_h)` under a local seed. The caller's RNG state is saved and restored, and `sample.kind = "Rejection"` is set explicitly. -3. Resolve each rank to a cell id in one block-wise pass (`terra::readStart` / `readValues` by row block, cumulative counts per stratum). - -The result depends only on cell order, which is row-major and fixed, and on base R's RNG. Memory stays bounded at BULK's 169M cells. Points are cell centres (`terra::xyFromCell`) and carry `point_id` (stable, zero-padded, stratum-then-rank order), `stratum`, `cell`, and `x`/`y`. `$strata` carries `stratum`, `n_cells`, `area` (ha), `weight` (`n_cells / Σ n_cells`) and `n`. Other details: -- Allocation is either equal per stratum (`n` as a scalar) or caller-supplied (a named vector). `n_h > N_h` is refused. -- `dft_check_crs()` guards the area. -- `NA` cells are outside the population. -- Raster strata only. Polygon strata are rasterised first with `terra::rasterize()` onto the map grid, documented in the roxygen. It is one line, and it keeps the grid (and so the weights) the caller's explicit choice. - -**Label contract** (columns the caller's table must carry; drift stores nothing): +| `dft_accuracy_sample(strata, n, seed, map = NULL)` | strata SpatRaster → `list(points = sf, strata = tibble, design = list)` | +| `dft_accuracy_estimate(labels, strata, level = 0.95, fpc = TRUE)` | label table + `$strata` → `list(matrix, accuracy, area, stratum)` | +| `dft_accuracy_size(weights, s_h, se_target, allocation, n_min)` | stratum weights + per-stratum SD + target SE → total n and allocation | +| `dft_accuracy_labels(labels, strata)` | validates the label contract; returns it invisibly or errors naming the fault | + +**Sampler.** `terra::spatSample()` and `terra::freq()` are not used: `spatSample()` is not stable across versions, and `freq()` returns labels on factor rasters and rounds floats. +1. **Pass 1** reads the layer in explicit row chunks with `readStart`/`readValues`, capped at about 1e7 cells, with `on.exit(readStop())`. It counts `N_h` per integer code (`tabulate`/`table`), and on the same pass refuses non-integer values and `nlyr != 1`. +2. **Draw:** each stratum gets its own sub-stream: `set.seed(, kind = "Mersenne-Twister", normal.kind = "Inversion", sample.kind = "Rejection")`, then `sample.int(N_h, n_h)`. `sample.int` keeps its prefix when n grows, and per-stratum seeds mean that raising one stratum's n, or adding a stratum, leaves every other draw unchanged. So a pilot extends into a superset of itself. The caller's `.Random.seed` is restored, or removed if it was absent. `seed` is required. +3. **Pass 2** runs the same chunked reader, turning ranks into cell ids. It ends by asserting the cumulative counts equal `N_h`. +4. **Allocation** is equal per stratum (a scalar) or a named vector. A named vector that omits a non-empty stratum is refused. Where `n_h ≥ N_h` the stratum is a **census** (`n_h = N_h`) and a message names it; tiny transition strata are normal (the bundled tile has a 3-cell one). +5. **Output:** + - Factor strata keep the integer code as `stratum` and carry `levels()` as `stratum_label`. + - Points are cell centres with `point_id = "_"` in draw order, which is independent of total n. They also carry `stratum`, `stratum_label` and `cell`, with no `x`/`y` columns: the geometry is authoritative. + - Optional `map =` takes a SpatRaster or a named list of years on the same grid (`terra::compareGeom`) and extracts `map_class` (or `map_`) at each cell. + - `$strata` carries `stratum`, `stratum_label`, `n_cells`, `area` (ha), `weight` and `n`. + - `$design` records the seed, RNG kinds, R and terra versions, allocation, and the grid dims, extent and CRS. + - `dft_check_crs()` guards area. NA cells are outside the population. + - Raster strata only. The roxygen gives the polygon recipe: rasterise in memory (not `filename =` with an integer datatype, which writes the background as 0, not NA), then `terra::mask()` to the map's NA footprint. + +**Label contract**: - required: `point_id`, `stratum`, `map_class`, `ref_class` -- optional: `confidence`, `reviewer`, `use` (`"accuracy"` / `"training"`) +- optional: `confidence`, `reviewer`, `use` (`"accuracy"` / `"training"` / NA) + +`dft_accuracy_labels()` refuses: +- `ref_class` NA (nonresponse is the caller's decision, and a silent drop changes `n_h`) +- duplicate `point_id` +- a stratum absent from `$strata` +- a `$strata` row with `N_h > 0` and no labels +- any other `use` value + +Its roxygen states that filtering by `confidence` changes the design. For change accuracy, `map_class` is the transition code and `ref_class = ref_from * 1000 + ref_to`, the `dft_rast_transition()` id scheme. Union targets (tree loss, the unattributed residual) are "recode, then estimate": the SE of a union is not the sum of the SEs. -**Train/test.** `dft_accuracy_estimate()` refuses any row with `use == "training"` and names the count. It has no silent drop: dropping those rows would bias the design, because the training rows were not drawn as part of the probability sample. +**Train/test.** `use == "training"` rows are refused with their count. Including them is the harm: points that trained a classifier cannot measure it. Reusing a non-random subset of accuracy points for training biases the estimate too. The documented-split option in the issue is not provided. That meets "at minimum a flag"; an `n_train` split at draw time is a possible follow-up. -**Estimators** (Stehman 2014 general form: per-stratum sample means of indicator variables, weighted by `N_h`): -- An area-weighted error matrix of estimated proportions p̂_ij. -- Overall, user's and producer's accuracy. UA and PA are ratio estimators, with their variances from Stehman eq. 26–28. -- Error-adjusted area per reference class = `A_total · p̂_·j`, with its SE and a 95% CI (`z = 1.96`, as the argument `level`). -- When strata = map classes, these are algebraically Olofsson eq. 2–11. A test asserts that equality on the paper's table. +**Estimators** (Stehman 2014 general form: stratified means of indicator variables, weighted by `N_h`): +- The error matrix is long format over the union of map and reference classes, with estimated proportions. When strata ≠ map classes, the row totals are **estimated**, not the known `W_i`; the roxygen says so. +- OA, UA and PA with SEs. UA and PA are ratio estimators; PA is NA where `p̂_·j = 0`. +- Adjusted area per reference class = `A_total · p̂_·j`, with SE and a Wald CI, `z = qnorm(1 − (1 − level)/2)`. The interval is not truncated at 0, and the roxygen says so. +- `fpc = TRUE` applies `(1 − n_h/N_h)`, so a census stratum contributes zero variance. With `fpc = FALSE` the estimator is algebraically Olofsson eq. 2–11; the Olofsson pin runs that way. +- `$stratum` reports per-stratum `n_h`, `N_h`, agreement mean and SE, and the SD of the OA indicator and of each reference-class indicator. The sizer consumes it. +- A stratum with `n_h = 1` is refused (its variance is undefined), unless it is a census and `fpc = TRUE`. -**Sizing.** Olofsson eq. 13: `n = (Σ W_i S_i / S(Ô))²` with `S_i = sqrt(U_i(1−U_i))`. Two allocation helpers: `"equal"`, and `"proportional_min"` (proportional with `n_min` per rare stratum, Olofsson §5.1.1). Pilot variance comes in as `ua` taken from a pilot's `dft_accuracy_estimate()$accuracy`, documented with an example. +**Sizing.** The primary form is `n = (Σ W_h S_h)² / SE_target²`, with `S_h` per **stratum** taken from a pilot's `$stratum` for a named quantity: OA or a class's area proportion. `ua =` is a convenience that holds only when strata = map classes, where `S_i = sqrt(U_i(1−U_i))` (Olofsson eq. 13). There are two allocations: `"equal"`, and `"proportional_min"` (Olofsson §5.1.1). ### Phase 1: Reference values in hand - [ ] User adds the Olofsson 2014 and Stehman 2014 PDFs (decided at the plan gate). Until then, Phases 2–4 code and non-published tests proceed, and the published-value pins wait. -- [ ] Transcribe the worked examples into `tests/testthat/helper-accuracy.R`, each value cited to its page and table number +- [ ] Transcribe the worked examples into `tests/testthat/helper-accuracy.R`, each value cited to its page and table number, with the equation numbers verified against the PDF and both papers checked for errata - [ ] `findings.md`: the estimator equations with numbers, and a check of the Olofsson example by hand arithmetic (deforestation 21,158 ha is reproducible from the row counts; confirm against the PDF) ### Phase 2: Estimator (tests first) -- [ ] `test-dft_accuracy_estimate.R`: Olofsson Table 8/9 (the error matrix, UA/PA/OA with SEs, and the adjusted areas with 95% CIs, to the published precision); the Stehman 2014 example (strata ≠ map classes); a must-fail test showing unweighted accuracy differs from the published values on the Olofsson example; `use == "training"` refused; a stratum with n_h = 1 (variance undefined) refused; a class absent from the reference handled -- [ ] `R/dft_accuracy_estimate.R` + `R/dft_accuracy_labels.R` (the contract validator, which the estimator calls) +- [ ] `test-dft_accuracy_estimate.R`, published: the Olofsson example (error matrix, UA/PA/OA with SEs, adjusted areas with CIs), run with `fpc = FALSE` and **absolute** tolerance at the published precision; the Stehman 2014 example (strata ≠ map classes). The fixture expands counts to per-point rows with base `rep()` +- [ ] Must-fail: run the **estimator** with equal `n_cells` on the Olofsson counts and assert it differs from the published values +- [ ] Census oracle (independent truth, no PDF needed): map = 2017 and "reference" = 2023 on the bundled tile, with the true error matrix and areas from `terra::crosstab`. Run about 500 stratified draws under map-class strata and under a changed/stable split; check bias ≈ 0, empirical SD ≈ mean SE, and CI coverage ≈ level. Skippable if slow +- [ ] Perfect labels (`ref = map`) give UA = PA = OA = 1, SE = 0, and adjusted area equal to mapped area. Recoding to a 2-class union gives an SE that is not the sum of the SEs +- [ ] Contract refusals: training rows, NA `ref_class`, duplicate ids, an unknown stratum, a stratum with no labels, `n_h = 1`. A reference-only class appears in the matrix, and PA is NA at `p̂_·j = 0` +- [ ] `R/dft_accuracy_estimate.R` + `R/dft_accuracy_labels.R`. Freeze the `$strata` and `$stratum` shapes here - [ ] Restore-the-bug check: remove the weights from the estimator and confirm the published-value tests go red ### Phase 3: Sampler (tests first) -- [ ] `test-dft_accuracy_sample.R`: pinned golden cell ids for a seeded draw on `example_2017.tif`; the same seed gives an identical draw; the caller's `.Random.seed` is untouched; within-stratum uniformity (chi-square over a large draw); weights sum to 1 and areas match `dft_rast_summarize()`; `n_h > N_h`, a lonlat raster and a non-integer raster are refused; NA cells are never drawn; the block-wise resolver agrees with a brute-force `which()` on the small tile +- [ ] `test-dft_accuracy_sample.R`: + - golden `point_id`s and cells for a seeded draw on `example_2017.tif`, with an allocation the tile can satisfy (classes 4 and 9 have 2 cells) + - chunk invariance: identical results at row chunks 1, 7, 50 and full, matching brute-force `which()` + - pilot extension: the first 30 per stratum at n = 30 are identical at n = 50, and adding a stratum leaves the others unchanged + - RNG: an existing seed is unchanged, an absent one is still absent, and a caller's L'Ecuyer kind is restored + - the factor transition raster: codes plus labels, joined correctly + - census with message + - within-stratum uniformity (chi-square) + - weights sum to 1 and areas match `dft_rast_summarize()` + - refusals: lonlat, non-integer, multi-layer, and a named allocation that omits a stratum + - NA cells are never drawn + - `map =` extraction, including a grid-mismatch refusal - [ ] `R/dft_accuracy_sample.R` -- [ ] Scale test on BULK `classified_2017.tif` (169M cells) with the RSS sampler; record the time and peak RSS in the PR body +- [ ] Sampler → estimator integration: draw, fake labels from a reference raster, estimate (covered by the census oracle once both exist) +- [ ] Scale test on BULK: `classified_2017.tif` and its `dft_rast_transition()` factor output (the #93 input), with an RSS sampler. Record pass-1 and pass-2 time and peak RSS in the PR body ### Phase 4: Sizing -- [ ] `test-dft_accuracy_size.R`: reproduce Olofsson §5.1.1's sample-size example (n and allocation) from the paper; edge cases (UA = 1, a zero weight) +- [ ] `test-dft_accuracy_size.R`: reproduce Olofsson §5.1.1's sample-size example (n and allocation); the `s_h` form from a pilot `$stratum` agrees with the `ua` form when strata = map classes; edge cases (UA = 1, a zero weight) - [ ] `R/dft_accuracy_size.R` ### Phase 5: Docs and release -- [ ] A runnable `@examples` block on every function (the estimator example uses Olofsson's published counts; the sampler uses the bundled tile) -- [ ] `devtools::document()`, `lintr::lint_package()`, `pkgdown::check_pkgdown()` (the new exports go in the reference index) +- [ ] A runnable `@examples` block on every function (the estimator example uses Olofsson's published counts; the sampler uses the bundled tile with a satisfiable allocation) +- [ ] `devtools::document()`, `lintr::lint_package()`, `pkgdown::check_pkgdown()`. `_pkgdown.yml` has no `reference:` index, so the check cannot catch an omission; it stays that way (out of scope) - [ ] NEWS.md, then version 0.19.0 as the final commit -- [ ] CLAUDE.md Core Pipeline: add the accuracy block -- [ ] Update the floodplains#93 body's "What lives where" table if the delivered API differs from what it assumes +- [ ] CLAUDE.md Core Pipeline: add the accuracy block, and correct the bundled tile's size (314 x 326, not 600 x 600) +- [ ] Update the floodplains#93 body's "What lives where" table: the `map =` argument, the `ref_class` composition recipe, the pilot-extension rule, and per-stratum sizing No vignette. One made with fabricated labels would illustrate a number nobody measured. The worked example belongs in floodplains#93 once real labels exist. From 9d3885d3f492fb66ad1b8f50813ffe151fba6690 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Mon, 28 Sep 2026 22:57:47 -0700 Subject: [PATCH 03/10] Record mapaccuracy/mapac and Olofsson 2014 PA typos in #81 findings Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- planning/active/findings.md | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/planning/active/findings.md b/planning/active/findings.md index ee08f8d..68ee179 100644 --- a/planning/active/findings.md +++ b/planning/active/findings.md @@ -39,6 +39,22 @@ Relates: floodplains#93, drift#79 - The estimator takes the general Stehman 2014 form (strata may differ from map classes), because floodplains#93 strata are not map classes. - No existing implementation in the org: swept the exports of 19 NGE packages plus `gh search code org:NewGraphEnvironment olofsson`. `mapaccuracy`, `survey` and `sampling` are not installed and not needed. +## Existing implementations (reported by a soul session, 2026-09-28) + +My org-only sweep missed these. Both are outside NewGraphEnvironment. + +- **`mapaccuracy`** (CRAN 0.1.2, 2024-04-03, Hugo Costa, MIT). Imports only `stats`. It implements Olofsson 2014 and Stehman 2014 (`olofsson()`, `stehman2014()`), and its docs check it against published examples from Olofsson 2013 (two), Olofsson 2014 and Stehman 2014. Its docs record a confirmed typo in Olofsson 2013 (a CI lower bound). +- **`mapac`** (Dirk Pflugmacher, HU Berlin GitLab, v0.31, 91 commits 2020–2026). Not on CRAN. It covers stratified and Stehman-2014 estimators, allocation, and report tables. Its tests are thin: one file, which checks only Stehman 2014. + +**The Olofsson 2014 PDF is in the NGE Zotero group** as `olofsson_etal2014Goodpractices`; its md5 was verified against the published PDF. Run through `mapac`, Tables 8–9 match on: +- all four areas with their CIs (deforestation 21,158 ± 6,158 ha) +- all four user's accuracies +- OA 0.947 ± 0.018, which the paper prints rounded as 0.95 ± 0.02 + +**Two producer's-accuracy CIs in the paper look like typos.** Forest gain is printed ±0.23 but Eq. 7 gives ±0.254; stable non-forest is printed ±0.01 but Eq. 7 gives ±0.018. The soul session recomputed both independently of `mapac`, and no erratum is registered. So those two pins take the Eq. 7 value, cite the discrepancy, and do not assert the printed figure. + +Stehman 2014 and Olofsson 2013 are paywalled and not yet saved. + ## Errors Encountered | Error | Resolution | From 75d7b1eea34cfc9c4f5692c1bc950b0d9d7b7b99 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Mon, 28 Sep 2026 22:59:35 -0700 Subject: [PATCH 04/10] Plan #81 estimator as a wrapper over mapaccuracy::stehman2014() The user chose to import mapaccuracy (CRAN, stats-only) over a second implementation. drift keeps the contract, CIs in ha, the tidy matrix and the per-stratum table the sizer needs. The FPC argument goes: the package always applies it. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- planning/active/progress.md | 1 + planning/active/review-plan.md | 2 +- planning/active/task_plan.md | 24 ++++++++++++++---------- 3 files changed, 16 insertions(+), 11 deletions(-) diff --git a/planning/active/progress.md b/planning/active/progress.md index 873181b..d2b51fb 100644 --- a/planning/active/progress.md +++ b/planning/active/progress.md @@ -6,3 +6,4 @@ - Created branch `81-accuracy-assessment-for-change-maps-stra` off main - Scaffolded PWF baseline from issue #81 with approved phases - Next: Phase 2 estimator code and tests; published-value pins wait on the PDFs (Phase 1) +- Decision: the estimator wraps `mapaccuracy::stehman2014()` (user, 2026-09-28) instead of reimplementing it; plan revised diff --git a/planning/active/review-plan.md b/planning/active/review-plan.md index 7767fbc..18afe39 100644 --- a/planning/active/review-plan.md +++ b/planning/active/review-plan.md @@ -6,7 +6,7 @@ Read-only reviewer, so its findings came back as reply text; I transcribed them |---|---|---| | B1 | `terra::freq()` returns labels on factor strata while `readValues()` returns codes; `freq()` also rounds floats. Count and resolve steps could not join | **Confirmed by probe** (`Water -> Water` vs `1001`). Two passes on one chunked reader; `stratum_label` from `levels()` | | B2 | Sizing through `ua` assumes strata = map classes; #93 needs per-stratum SD for an area target | Primary form takes `s_h` from the estimator's `$stratum`; `ua` kept as a convenience | -| B3 | FPC undecided; `n_h > N_h` refusal breaks on tiny strata (the tile has 2- and 3-cell strata) | `fpc = TRUE` by default, and the Olofsson pin runs with `fpc = FALSE`; a stratum with `n_h ≥ N_h` becomes a census, with a message | +| B3 | FPC undecided; `n_h > N_h` refusal breaks on tiny strata (the tile has 2- and 3-cell strata) | A stratum with `n_h ≥ N_h` becomes a census, with a message. FPC is now fixed on by `mapaccuracy`, so there is no `fpc` argument (superseded 2026-09-28) | | G1 | The sampler does not emit `map_class` | `map =` argument, extracted at each cell | | G2 | One `blocks()` chunk on the tile, so the carry logic is untested | Explicit row chunks, capped at about 1e7 cells; chunk-invariance test | | G3 | RNG save/restore under-specified (L'Ecuyer, absent seed) | Pin all three kinds; restore, or remove; three tests | diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index b783131..faa2078 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -14,7 +14,7 @@ A new `dft_accuracy_*` family, one function per file: | function | input → output | |---|---| | `dft_accuracy_sample(strata, n, seed, map = NULL)` | strata SpatRaster → `list(points = sf, strata = tibble, design = list)` | -| `dft_accuracy_estimate(labels, strata, level = 0.95, fpc = TRUE)` | label table + `$strata` → `list(matrix, accuracy, area, stratum)` | +| `dft_accuracy_estimate(labels, strata, level = 0.95)` | label table + `$strata` → `list(matrix, accuracy, area, stratum)` | | `dft_accuracy_size(weights, s_h, se_target, allocation, n_min)` | stratum weights + per-stratum SD + target SE → total n and allocation | | `dft_accuracy_labels(labels, strata)` | validates the label contract; returns it invisibly or errors naming the fault | @@ -47,13 +47,15 @@ Its roxygen states that filtering by `confidence` changes the design. For change **Train/test.** `use == "training"` rows are refused with their count. Including them is the harm: points that trained a classifier cannot measure it. Reusing a non-random subset of accuracy points for training biases the estimate too. The documented-split option in the issue is not provided. That meets "at minimum a flag"; an `n_train` split at draw time is a possible follow-up. -**Estimators** (Stehman 2014 general form: stratified means of indicator variables, weighted by `N_h`): -- The error matrix is long format over the union of map and reference classes, with estimated proportions. When strata ≠ map classes, the row totals are **estimated**, not the known `W_i`; the roxygen says so. -- OA, UA and PA with SEs. UA and PA are ratio estimators; PA is NA where `p̂_·j = 0`. -- Adjusted area per reference class = `A_total · p̂_·j`, with SE and a Wald CI, `z = qnorm(1 − (1 − level)/2)`. The interval is not truncated at 0, and the roxygen says so. -- `fpc = TRUE` applies `(1 − n_h/N_h)`, so a census stratum contributes zero variance. With `fpc = FALSE` the estimator is algebraically Olofsson eq. 2–11; the Olofsson pin runs that way. -- `$stratum` reports per-stratum `n_h`, `N_h`, agreement mean and SE, and the SD of the OA indicator and of each reference-class indicator. The sizer consumes it. -- A stratum with `n_h = 1` is refused (its variance is undefined), unless it is a census and `fpc = TRUE`. +**Estimators: wrap `mapaccuracy::stehman2014()`** (decided 2026-09-28, after a soul session surfaced the package; see findings). `mapaccuracy` is on CRAN, MIT-licensed, imports only `stats`, and is checked against the published Olofsson 2013/2014 and Stehman 2014 examples. It becomes an Import. `dft_accuracy_estimate()` owns only what the package does not: +- the label contract, and refusing training rows (via `dft_accuracy_labels()`) +- passing `stratum` as character codes, because `stehman2014()` matches stratum names by **regex** (`grep(paste0("^", nm, "$"))`), so a label such as `"Trees -> Water"` or one with `.`/`+`/`(` is unsafe +- area in hectares = `A_total · area`, with SE, and a Wald CI with `z = qnorm(1 − (1 − level)/2)`, not truncated at 0 +- a tidy long-format error matrix. The package returns zero cells as NA; they become 0, and the conversion is documented +- a `$stratum` table (`n_h`, `N_h`, agreement mean and SE, and the SD of the OA indicator and of each reference-class indicator), which `mapaccuracy` does not return and the sizer needs +- a refusal for a stratum with `n_h = 1` unless it is a census. The package only warns + +The package always applies the FPC `(1 − n_h/N_h)`, so a census stratum contributes zero variance and there is no `fpc` argument. Olofsson 2014 omits the FPC; at its pixel-scale `N_h` the difference falls below the published precision, and the test asserts that. When strata ≠ map classes, the matrix row totals are estimated, not the known `W_i`; the roxygen says so. Any upstream defect found is filed on its tracker and shown to the user, not worked around (karpathy §8). **Sizing.** The primary form is `n = (Σ W_h S_h)² / SE_target²`, with `S_h` per **stratum** taken from a pilot's `$stratum` for a named quantity: OA or a class's area proportion. `ua =` is a convenience that holds only when strata = map classes, where `S_i = sqrt(U_i(1−U_i))` (Olofsson eq. 13). There are two allocations: `"equal"`, and `"proportional_min"` (Olofsson §5.1.1). @@ -63,13 +65,15 @@ Its roxygen states that filtering by `confidence` changes the design. For change - [ ] `findings.md`: the estimator equations with numbers, and a check of the Olofsson example by hand arithmetic (deforestation 21,158 ha is reproducible from the row counts; confirm against the PDF) ### Phase 2: Estimator (tests first) -- [ ] `test-dft_accuracy_estimate.R`, published: the Olofsson example (error matrix, UA/PA/OA with SEs, adjusted areas with CIs), run with `fpc = FALSE` and **absolute** tolerance at the published precision; the Stehman 2014 example (strata ≠ map classes). The fixture expands counts to per-point rows with base `rep()` +- [ ] `mapaccuracy` in Imports (DESCRIPTION) +- [ ] `test-dft_accuracy_estimate.R`, published: Olofsson 2014 Tables 8–9 through `dft_accuracy_estimate()`, with **absolute** tolerance at the published precision. The forest-gain and stable-non-forest PA CIs pin the Eq. 7 values (±0.254, ±0.018), not the printed ±0.23 / ±0.01, and cite the discrepancy (findings). Add the Stehman 2014 example once its PDF is in hand. The fixture expands counts to per-point rows with base `rep()` - [ ] Must-fail: run the **estimator** with equal `n_cells` on the Olofsson counts and assert it differs from the published values - [ ] Census oracle (independent truth, no PDF needed): map = 2017 and "reference" = 2023 on the bundled tile, with the true error matrix and areas from `terra::crosstab`. Run about 500 stratified draws under map-class strata and under a changed/stable split; check bias ≈ 0, empirical SD ≈ mean SE, and CI coverage ≈ level. Skippable if slow - [ ] Perfect labels (`ref = map`) give UA = PA = OA = 1, SE = 0, and adjusted area equal to mapped area. Recoding to a 2-class union gives an SE that is not the sum of the SEs - [ ] Contract refusals: training rows, NA `ref_class`, duplicate ids, an unknown stratum, a stratum with no labels, `n_h = 1`. A reference-only class appears in the matrix, and PA is NA at `p̂_·j = 0` +- [ ] Measure `stehman2014()` runtime at #93 scale (about 1,000 points × about 80 transition classes; it builds `classes²` indicator columns). File upstream and report if it is impractical - [ ] `R/dft_accuracy_estimate.R` + `R/dft_accuracy_labels.R`. Freeze the `$strata` and `$stratum` shapes here -- [ ] Restore-the-bug check: remove the weights from the estimator and confirm the published-value tests go red +- [ ] Restore-the-bug check: pass unweighted `N_h` inside the wrapper and confirm the published-value tests go red ### Phase 3: Sampler (tests first) - [ ] `test-dft_accuracy_sample.R`: From b32d063b3a430fd09c41077759272babfbcab1ca Mon Sep 17 00:00:00 2001 From: almac2022 Date: Tue, 29 Sep 2026 06:46:11 -0700 Subject: [PATCH 05/10] Transcribe Olofsson 2014 reference values for #81 Tables 8-9, section 5.1.1 and 5.2, with page and table cited. Two PA half-widths and one area SE are printed inconsistently with the paper's own equations; the equation values are kept and the printed ones recorded. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- planning/active/findings.md | 22 ++++++++ planning/active/progress.md | 1 + planning/active/task_plan.md | 6 +- tests/testthat/helper-accuracy.R | 94 ++++++++++++++++++++++++++++++++ 4 files changed, 120 insertions(+), 3 deletions(-) create mode 100644 tests/testthat/helper-accuracy.R diff --git a/planning/active/findings.md b/planning/active/findings.md index 68ee179..ec54691 100644 --- a/planning/active/findings.md +++ b/planning/active/findings.md @@ -55,6 +55,28 @@ My org-only sweep missed these. Both are outside NewGraphEnvironment. Stehman 2014 and Olofsson 2013 are paywalled and not yet saved. +## Olofsson 2014 reference values (read from the PDF, 2026-09-29) + +The values are transcribed into `tests/testthat/helper-accuracy.R` with page and table. From Tables 8–9 (p. 55) and §5.2 (p. 54), through `mapaccuracy::stehman2014()`: +- Areas match to the hectare (21,157.8 / 11,686.2 / 285,769.9 / 581,386.2 against the printed 21,158 / 11,686 / 285,770 / 581,386). +- The error matrix matches Table 9 to 4 dp. The package returns zero cells as NA. +- UA, the two PA values that agree with the text, and OA (0.9465 ± 0.0185) all match at 2 dp. + +**The half-widths depend on two choices.** Stehman's FPC and `qnorm(.975)` versus the paper's no-FPC and 1.96 move the half-widths by up to 1.1 ha: stable non-forest is 16,280.9 against the printed 16,282, while no-FPC with 1.96 gives 16,281.7. So the area half-width pin uses an absolute tolerance of 1.5 ha, and the proportion pins use 0.005 (the half-step at 2 dp). + +**Three printed values contradict the paper's own equations:** +- PA half-width, forest gain: printed 0.23, Eq. 7 gives 0.254. +- PA half-width, stable non-forest: printed 0.01, Eq. 7 gives 0.018. +- "S(Â₁) … = 34,097 pixels": 1.96 × 34,097 = 66,830, not the 68,418 printed next to it. 68,418 / 1.96 = 34,907. + +**Sample size:** Eq. 13 with Table 5's W and U and a target SE(O) of 0.01 gives (0.25312 / 0.01)² = 640.7, so n = 641. The Equal (160) and Prop (13 / 10 / 205 / 413, by rounding n·W) columns reproduce. **Alloc1–3 do not:** the stated rule (100 per change stratum, remainder proportional to the stable classes) gives 146 / 295, where the table has 149 / 292. They are not asserted. + +**`mapaccuracy` internals worth knowing:** +- `stehman2014()` matches stratum names by regex (`grep(paste0("^", nm, "$"), ...)`). drift passes internal ids `s1..sH` to it. +- Its `order =` default is `sort(union(r, m))` on character, so "10" sorts before "2". drift passes the order explicitly. +- It only *warns* on a stratum with one observation. +- It applies the FPC in eq. 25 and eq. 28. + ## Errors Encountered | Error | Resolution | diff --git a/planning/active/progress.md b/planning/active/progress.md index d2b51fb..a0e5725 100644 --- a/planning/active/progress.md +++ b/planning/active/progress.md @@ -7,3 +7,4 @@ - Scaffolded PWF baseline from issue #81 with approved phases - Next: Phase 2 estimator code and tests; published-value pins wait on the PDFs (Phase 1) - Decision: the estimator wraps `mapaccuracy::stehman2014()` (user, 2026-09-28) instead of reimplementing it; plan revised +- Phase 1 done: Olofsson 2014 read from the local Zotero PDF; values in `helper-accuracy.R`; three printed values contradict the paper's own equations (findings). Issue body revised. Stehman 2014 is not needed (mapaccuracy tests it) diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index faa2078..3da7cf9 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -60,9 +60,9 @@ The package always applies the FPC `(1 − n_h/N_h)`, so a census stratum contri **Sizing.** The primary form is `n = (Σ W_h S_h)² / SE_target²`, with `S_h` per **stratum** taken from a pilot's `$stratum` for a named quantity: OA or a class's area proportion. `ua =` is a convenience that holds only when strata = map classes, where `S_i = sqrt(U_i(1−U_i))` (Olofsson eq. 13). There are two allocations: `"equal"`, and `"proportional_min"` (Olofsson §5.1.1). ### Phase 1: Reference values in hand -- [ ] User adds the Olofsson 2014 and Stehman 2014 PDFs (decided at the plan gate). Until then, Phases 2–4 code and non-published tests proceed, and the published-value pins wait. -- [ ] Transcribe the worked examples into `tests/testthat/helper-accuracy.R`, each value cited to its page and table number, with the equation numbers verified against the PDF and both papers checked for errata -- [ ] `findings.md`: the estimator equations with numbers, and a check of the Olofsson example by hand arithmetic (deforestation 21,158 ha is reproducible from the row counts; confirm against the PDF) +- [x] User adds the Olofsson 2014 and Stehman 2014 PDFs (decided at the plan gate). Until then, Phases 2–4 code and non-published tests proceed, and the published-value pins wait. +- [x] Transcribe the worked examples into `tests/testthat/helper-accuracy.R`, each value cited to its page and table number, with the equation numbers verified against the PDF and both papers checked for errata +- [x] `findings.md`: the estimator equations with numbers, and a check of the Olofsson example by hand arithmetic (deforestation 21,158 ha is reproducible from the row counts; confirm against the PDF) ### Phase 2: Estimator (tests first) - [ ] `mapaccuracy` in Imports (DESCRIPTION) diff --git a/tests/testthat/helper-accuracy.R b/tests/testthat/helper-accuracy.R new file mode 100644 index 0000000..b44b672 --- /dev/null +++ b/tests/testthat/helper-accuracy.R @@ -0,0 +1,94 @@ +# Published reference values for the accuracy-assessment tests (#81). +# +# Olofsson, P., Foody, G.M., Herold, M., Stehman, S.V., Woodcock, C.E., +# Wulder, M.A. (2014). Good practices for estimating area and assessing +# accuracy of land change. Remote Sensing of Environment 148: 42-57. +# doi:10.1016/j.rse.2014.02.015. Zotero: olofsson_etal2014Goodpractices. +# +# Every value below is transcribed from the PDF, with its page and table. Where +# the printed value contradicts the paper's own equation, the equation value is +# kept and the printed one is recorded beside it -- see olofsson_printed_typos. + +olofsson_classes <- c("deforestation", "forest_gain", "stable_forest", + "stable_nonforest") + +# Table 8 (p. 55): sample counts n_ij, rows = map (= stratum), cols = reference +olofsson_counts <- matrix( + c(66, 0, 5, 4, + 0, 55, 8, 12, + 1, 0, 153, 11, + 2, 1, 9, 313), + nrow = 4, byrow = TRUE, dimnames = list(olofsson_classes, olofsson_classes) +) + +# Table 8 (p. 55): mapped area in 30 m pixels. 1 pixel = 0.09 ha. +olofsson_pixels <- c(200000, 150000, 3200000, 6450000) +names(olofsson_pixels) <- olofsson_classes + +# Table 9 (p. 55): estimated area proportions p_ij, to 4 dp +olofsson_p <- matrix( + c(0.0176, 0, 0.0013, 0.0011, + 0, 0.0110, 0.0016, 0.0024, + 0.0019, 0, 0.2967, 0.0213, + 0.0040, 0.0020, 0.0179, 0.6212), + nrow = 4, byrow = TRUE, dimnames = list(olofsson_classes, olofsson_classes) +) + +# Section 5.2.1 (p. 54): estimate +/- 95% half-width, to 2 dp +olofsson_user <- c(0.88, 0.73, 0.93, 0.96) +olofsson_user_hw <- c(0.07, 0.10, 0.04, 0.02) +olofsson_prod <- c(0.75, 0.85, 0.93, 0.96) +# Printed as 0.21, 0.23, 0.03, 0.01. Eq. (7) gives 0.254 for forest gain and +# 0.018 for stable non-forest (recomputed independently of any package, and +# matching mapac and mapaccuracy), so those two are kept at the equation value. +olofsson_prod_hw <- c(0.21, 0.25, 0.03, 0.02) +olofsson_overall <- 0.95 +olofsson_overall_hw <- 0.02 + +# Section 5.2.2 (p. 54): error-adjusted area and 95% half-width, ha +olofsson_area <- c(21158, 11686, 285770, 581386) +olofsson_area_hw <- c(6158, 3756, 15510, 16282) + +# Printed values that disagree with the paper's own arithmetic, kept for the +# record and asserted nowhere +olofsson_printed_typos <- list( + producer_hw_forest_gain = c(printed = 0.23, eq7 = 0.254), + producer_hw_stable_nonforest = c(printed = 0.01, eq7 = 0.018), + # "S(A_1) = 0.0035 x 10,000,000 = 34,097 pixels" -- 1.96 x 34,097 is 66,830, + # not the 68,418 printed next; 68,418 / 1.96 = 34,907 + se_area_deforestation_px = c(printed = 34097, implied = 34907) +) + +# Section 5.1.1 and Table 5 (p. 53): sample size planning +olofsson_plan_weights <- c(0.020, 0.015, 0.320, 0.645) +olofsson_plan_ua <- c(0.70, 0.60, 0.90, 0.95) +olofsson_plan_n <- 641 +olofsson_plan_equal <- c(160, 160, 160, 160) +olofsson_plan_prop <- c(13, 10, 205, 413) +# Table 5's Alloc1-3 columns (e.g. 100/100/149/292) are not reproducible from +# the rule the paper states ("allocate the remainder proportionally to the +# stable classes" gives 146/295), so they are not asserted. + +# One row per sample point, expanded from Table 8, in the label contract shape +olofsson_labels <- function() { + idx <- which(olofsson_counts > 0, arr.ind = TRUE) + map <- rep(olofsson_classes[idx[, "row"]], olofsson_counts[idx]) + ref <- rep(olofsson_classes[idx[, "col"]], olofsson_counts[idx]) + data.frame( + point_id = sprintf("p%04d", seq_along(map)), + stratum = map, + map_class = map, + ref_class = ref, + stringsAsFactors = FALSE + ) +} + +# The strata table as dft_accuracy_sample() emits it +olofsson_strata <- function(n_cells = olofsson_pixels) { + tibble::tibble( + stratum = olofsson_classes, + n_cells = unname(n_cells), + area = unname(n_cells) * 0.09, + weight = unname(n_cells) / sum(n_cells) + ) +} From ccdafc369dd3df80aa0e8f33b1fd0b1d48473b64 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Tue, 29 Sep 2026 07:18:10 -0700 Subject: [PATCH 06/10] Add stratified accuracy sampling and error-adjusted area (#81) dft_accuracy_sample() draws stratified random points from any strata raster in two chunked passes, with per-stratum seeded streams, so a larger n extends a pilot and draws do not depend on terra's sampler. dft_accuracy_estimate() wraps mapaccuracy::stehman2014() (strata may differ from the map classes) and adds CIs in ha, a long matrix and a per-stratum table. dft_accuracy_labels() is the label contract and dft_accuracy_size() is Olofsson Eq. 13 plus allocation. Tests reproduce Olofsson et al. 2014 Tables 8-9, check against a census of the bundled tiles, and pin a golden draw. Three code-check rounds, the last ended by enumerating every identity site. Relates to NewGraphEnvironment/sred-2025-2026#16 Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- DESCRIPTION | 1 + NAMESPACE | 4 + R/dft_accuracy_estimate.R | 276 +++++++++++++++ R/dft_accuracy_labels.R | 191 ++++++++++ R/dft_accuracy_sample.R | 368 ++++++++++++++++++++ R/dft_accuracy_size.R | 179 ++++++++++ man/dft_accuracy_estimate.Rd | 137 ++++++++ man/dft_accuracy_labels.Rd | 89 +++++ man/dft_accuracy_sample.Rd | 123 +++++++ man/dft_accuracy_size.Rd | 107 ++++++ planning/active/findings.md | 28 ++ planning/active/progress.md | 2 + planning/active/review-round1.md | 50 +++ planning/active/review-round2.md | 41 +++ planning/active/review-round3.md | 151 ++++++++ planning/active/task_plan.md | 32 +- tests/testthat/helper-accuracy.R | 23 +- tests/testthat/test-dft_accuracy_estimate.R | 278 +++++++++++++++ tests/testthat/test-dft_accuracy_sample.R | 294 ++++++++++++++++ tests/testthat/test-dft_accuracy_size.R | 85 +++++ 20 files changed, 2439 insertions(+), 20 deletions(-) create mode 100644 R/dft_accuracy_estimate.R create mode 100644 R/dft_accuracy_labels.R create mode 100644 R/dft_accuracy_sample.R create mode 100644 R/dft_accuracy_size.R create mode 100644 man/dft_accuracy_estimate.Rd create mode 100644 man/dft_accuracy_labels.Rd create mode 100644 man/dft_accuracy_sample.Rd create mode 100644 man/dft_accuracy_size.Rd create mode 100644 planning/active/review-round1.md create mode 100644 planning/active/review-round2.md create mode 100644 planning/active/review-round3.md create mode 100644 tests/testthat/test-dft_accuracy_estimate.R create mode 100644 tests/testthat/test-dft_accuracy_sample.R create mode 100644 tests/testthat/test-dft_accuracy_size.R diff --git a/DESCRIPTION b/DESCRIPTION index 03ecb95..dd3cfcb 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -25,6 +25,7 @@ Imports: cli, digest, dplyr, + mapaccuracy, parallel, rappdirs, rlang, diff --git a/NAMESPACE b/NAMESPACE index c78362a..e6a766b 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,5 +1,9 @@ # Generated by roxygen2: do not edit by hand +export(dft_accuracy_estimate) +export(dft_accuracy_labels) +export(dft_accuracy_sample) +export(dft_accuracy_size) export(dft_break_category) export(dft_break_strength) export(dft_cache_clear) diff --git a/R/dft_accuracy_estimate.R b/R/dft_accuracy_estimate.R new file mode 100644 index 0000000..c1acb04 --- /dev/null +++ b/R/dft_accuracy_estimate.R @@ -0,0 +1,276 @@ +#' Accuracy and error-adjusted area from a stratified reference sample +#' +#' Turn reference labels at stratified random sample points into the numbers +#' a change map should be reported with: an error matrix in estimated +#' proportions of area, overall / user's / producer's accuracy, and +#' **error-adjusted area per class with a confidence interval** -- following +#' Olofsson et al. (2014) and Stehman (2014). +#' +#' Mapped area is a biased estimator whenever the map has errors, and for a +#' change map the bias is usually large: a single wrong label on either date +#' manufactures a transition. The error-adjusted area is the area the +#' *reference* labels imply, weighted back to the population by each +#' stratum's size. +#' +#' @param labels A data frame meeting the [dft_accuracy_labels()] contract: +#' `point_id`, `stratum`, `map_class`, `ref_class`. +#' @param strata The `$strata` table from [dft_accuracy_sample()], or any data +#' frame with `stratum`, `n_cells` (the stratum's size in cells, `n_cells_h`) and +#' `area` (its mapped area in hectares). +#' @param level Confidence level for the intervals. Default `0.95`. +#' +#' @return A list: +#' - `matrix` -- tibble, long: `map_class`, `ref_class`, `proportion` (the +#' estimated share of total area), one row per pair of classes, zeros +#' included. +#' - `accuracy` -- tibble: `measure` (`"overall"`, `"user"`, `"producer"`), +#' `class` (`NA` for overall), `estimate`, `se`, `lower`, `upper`. +#' - `area` -- tibble, one row per class: `proportion`, `proportion_se`, +#' `area`, `area_se`, `lower`, `upper` (hectares). +#' - `stratum` -- tibble, long, per stratum and `target`: `n`, `n_cells`, +#' `weight`, `mean` and `sd` of an indicator within the stratum. `target` is +#' `"agreement"` (map equals reference) or a reference class, written in +#' full as a string (`"100000"`, never `"1e+05"`). `sd` is what +#' [dft_accuracy_size()] takes to size a full sample from a pilot. +#' - `level`, `area_total` (ha). +#' +#' @section Estimators: +#' The estimates come from [mapaccuracy::stehman2014()], which implements +#' Stehman (2014)'s estimators for stratified random sampling. Those hold +#' **whether or not the strata are the map classes** -- strata such as "change +#' attributed to fire", "change unattributed" and "stable" are fine. When the +#' strata are the map classes they give Olofsson et al. (2014)'s estimates, +#' and the package test suite reproduces that paper's worked example. +#' +#' - The **finite population correction** `(1 - n_h / n_cells_h)` is always +#' applied, so a stratum sampled in full (a census) contributes no variance. +#' Olofsson et al. omit it; at pixel-scale `n_cells_h` the difference is below +#' anything reported. +#' - Intervals are Wald intervals, `estimate +/- z * se` with +#' `z = qnorm(1 - (1 - level) / 2)`, and are **not** truncated at 0 or 1: a +#' lower bound below 0 says the class is too rare for the sample to bound +#' away from zero, and truncating it would hide that. +#' - **A rare class hiding in a large stratum makes the interval too narrow +#' at small samples.** If 1\% of a big stratum is really class *j* and +#' the stratum gets 25 points, most samples see none of it, estimate that +#' stratum's variance for *j* as 0, and report an interval that misses. On +#' the bundled tiles (98 of 7,127 "Trees" cells reference Water) Water's +#' 95\% interval covered 61\% of the time at 25 points per stratum, 83\% +#' at 75 and 89\% at 150; the estimate itself stayed unbiased. Omission +#' hides in the large stable strata, so do not starve them. +#' - When the strata are not the map classes, the error matrix's row totals +#' (the map-class shares) are **estimated** from the sample, not the known +#' stratum weights. That surprises readers used to Olofsson's tables. +#' - Producer's accuracy is `NA` for a class no reference label falls in. +#' - Run time grows with roughly the 2.5th power of the number of classes: +#' about 13 s for 1,000 points over 80 transition classes. Collapse classes +#' you will not report before estimating. +#' +#' @section Targets that are unions of classes: +#' "Tree loss" is every Trees -> non-Trees transition; an "unattributed +#' residual" is loss outside any fire or harvest. Recode `map_class` and +#' `ref_class` to the target (say `"loss"` / `"other"`) and estimate again. The +#' standard error of a union is **not** the sum of its members' standard +#' errors, because the members' estimates are correlated. +#' +#' @section Training points: +#' A row with `use == "training"` is refused. Points that trained a +#' classifier cannot also measure it -- the estimate would be optimistic by +#' construction. Keep the accuracy sample separate from the start; a subset of +#' accuracy points chosen for training by judgement also stops being a random +#' sample of its stratum. +#' +#' @references +#' Olofsson, P., Foody, G.M., Herold, M., Stehman, S.V., Woodcock, C.E. and +#' Wulder, M.A. (2014). Good practices for estimating area and assessing +#' accuracy of land change. *Remote Sensing of Environment* 148, 42-57. +#' \doi{10.1016/j.rse.2014.02.015} +#' +#' Stehman, S.V. (2014). Estimating area and map accuracy for stratified +#' random sampling when the strata are different from the map classes. +#' *International Journal of Remote Sensing* 35(13), 4923-4939. +#' \doi{10.1080/01431161.2014.930207} +#' +#' @seealso [dft_accuracy_sample()] to draw the points, +#' [dft_accuracy_labels()] for the label contract, +#' [dft_accuracy_size()] to size a sample from a pilot. +#' +#' @export +#' @examples +#' # Olofsson et al. (2014), Table 8: 640 labelled points in four strata that +#' # are the map classes, over 10 million 30 m pixels +#' cls <- c("deforestation", "forest_gain", "stable_forest", "stable_nonforest") +#' counts <- matrix(c(66, 0, 5, 4, 0, 55, 8, 12, 1, 0, 153, 11, 2, 1, 9, 313), +#' nrow = 4, byrow = TRUE) +#' idx <- which(counts > 0, arr.ind = TRUE) +#' map <- rep(cls[idx[, 1]], counts[idx]) +#' ref <- rep(cls[idx[, 2]], counts[idx]) +#' labels <- data.frame(point_id = seq_along(map), stratum = map, +#' map_class = map, ref_class = ref) +#' pixels <- c(200000, 150000, 3200000, 6450000) +#' strata <- data.frame(stratum = cls, n_cells = pixels, area = pixels * 0.09) +#' +#' res <- dft_accuracy_estimate(labels, strata) +#' res$area # deforestation: 21,158 ha +/- 6,157 -- mapped was 18,000 +#' res$accuracy +dft_accuracy_estimate <- function(labels, strata, level = 0.95) { + dft_accuracy_labels(labels, strata) + if (!"area" %in% names(strata)) { + stop("`strata` needs an `area` column (hectares) for error-adjusted area.", + call. = FALSE) + } + level_ok <- is.numeric(level) && length(level) == 1L && !is.na(level) && + level > 0 && level < 1 + if (!level_ok) { + stop("`level` must be a single number between 0 and 1.", call. = FALSE) + } + if ("use" %in% names(labels)) { + n_train <- sum(as.character(labels$use) %in% "training") + if (n_train > 0L) { + stop(n_train, " row(s) have `use == \"training\"`. Points that trained ", + "a classifier cannot measure its accuracy; estimate from the ", + "held-out accuracy sample only.", call. = FALSE) + } + } + + strata <- strata[strata$n_cells > 0, , drop = FALSE] + cell_area <- strata$area / strata$n_cells + grid_ok <- all(is.finite(cell_area)) && + max(abs(cell_area / cell_area[1] - 1)) <= 1e-6 + if (!grid_ok) { + stop("`strata$area` / `strata$n_cells` differs between strata: the ", + "weights (cells) and the areas describe different grids.", + call. = FALSE) + } + area_total <- sum(strata$area) + + key <- accuracy_key(strata$stratum) + lab_key <- accuracy_key(labels$stratum) + n_h <- as.vector(table(factor(lab_key, levels = key))) + n_cells_h <- strata$n_cells + if (any(n_h > n_cells_h)) { + bad <- key[n_h > n_cells_h] + stop("More labels than cells in stratum/strata: ", + paste(bad, collapse = ", "), ".", call. = FALSE) + } + census <- n_h == n_cells_h + thin <- n_h < 2L & !census + if (any(thin)) { + stop("Stratum/strata with a single labelled point: ", + paste(key[thin], collapse = ", "), ". Its variance is undefined; ", + "label at least two points (or all of a stratum's cells).", + call. = FALSE) + } + + # one normalisation, accuracy_key(), for the class set, the labels and every + # comparison: c() on a factor and a non-factor falls back to the factor's + # integer codes, and as.character() writes 100000 as "1e+05" for a double + # but not an integer -- either would split one class in two + map_chr <- accuracy_key(labels$map_class) + ref_chr <- accuracy_key(labels$ref_class) + classes <- accuracy_class_order(c(map_chr, ref_chr)) + classes_numeric <- is.numeric(labels$map_class) && + is.numeric(labels$ref_class) + + # stehman2014() matches stratum names by regex, so pass ids that cannot + # carry a metacharacter; its default class order sorts "10" before "2", so + # pass the order too + sid <- paste0("s", seq_along(key)) + nh_strata <- stats::setNames(n_cells_h, sid) + s_lab <- sid[match(lab_key, key)] + est <- withCallingHandlers( + mapaccuracy::stehman2014(s = s_lab, r = ref_chr, m = map_chr, + Nh_strata = nh_strata, margins = FALSE, + order = classes), + # a single-point stratum is refused above unless it is a census, where + # the finite population correction makes its variance exactly 0 + warning = function(w) { + if (grepl("include only one observation", conditionMessage(w))) { + invokeRestart("muffleWarning") + } + } + ) + + z <- stats::qnorm(1 - (1 - level) / 2) + cls_out <- if (classes_numeric) as.numeric(classes) else classes + + m <- est$matrix + m[is.na(m)] <- 0 + matrix_long <- tibble::tibble( + map_class = rep(cls_out, times = length(classes)), + ref_class = rep(cls_out, each = length(classes)), + proportion = as.vector(m) + ) + + acc_row <- function(measure, class, e, se) { + tibble::tibble(measure = measure, class = class, estimate = unname(e), + se = unname(se), lower = unname(e - z * se), + upper = unname(e + z * se)) + } + accuracy <- rbind( + acc_row("overall", cls_out[NA_integer_], est$OA, est$SEoa), + acc_row("user", cls_out, est$UA[classes], est$SEua[classes]), + acc_row("producer", cls_out, est$PA[classes], est$SEpa[classes]) + ) + + area <- tibble::tibble( + class = cls_out, + proportion = unname(est$area[classes]), + proportion_se = unname(est$SEa[classes]), + area = unname(est$area[classes]) * area_total, + area_se = unname(est$SEa[classes]) * area_total + ) + area$lower <- area$area - z * area$area_se + area$upper <- area$area + z * area$area_se + + list( + matrix = matrix_long, + accuracy = accuracy, + area = area, + stratum = accuracy_stratum_table(labels, strata, classes), + level = level, + area_total = area_total + ) +} + +# Keys (from accuracy_key()) in a stable, locale-free order: numeric order when +# every class reads as a number (transition ids 1001 < 11011), radix (C-locale) +# order otherwise +accuracy_class_order <- function(x) { + u <- unique(x) + num <- suppressWarnings(as.numeric(u)) + if (!anyNA(num)) u[order(num)] else sort(u, method = "radix") +} + +accuracy_sd <- function(v) { + if (length(v) > 1L) stats::sd(v) else NA_real_ +} + +# Per-stratum mean and SD of the agreement indicator and of each reference +# class indicator -- the S_h a pilot hands to dft_accuracy_size() +accuracy_stratum_table <- function(labels, strata, classes) { + key <- accuracy_key(strata$stratum) + lab_key <- accuracy_key(labels$stratum) + ref_chr <- accuracy_key(labels$ref_class) + agree <- accuracy_key(labels$map_class) == ref_chr + weight <- strata$n_cells / sum(strata$n_cells) + targets <- c("agreement", classes) + rows <- lapply(seq_along(key), function(h) { + in_h <- lab_key == key[h] + ind <- c(list(agree[in_h]), + lapply(classes, function(cl) ref_chr[in_h] == cl)) + tibble::tibble( + stratum = strata$stratum[h], + n = sum(in_h), + n_cells = strata$n_cells[h], + weight = weight[h], + target = targets, + mean = vapply(ind, mean, numeric(1)), + # a census stratum has no sampling variance, whatever its sd -- and a + # one-cell census has no sd at all, which the sizer would refuse + sd = if (sum(in_h) == strata$n_cells[h]) 0 else + vapply(ind, accuracy_sd, numeric(1)) + ) + }) + do.call(rbind, rows) +} diff --git a/R/dft_accuracy_labels.R b/R/dft_accuracy_labels.R new file mode 100644 index 0000000..f16175e --- /dev/null +++ b/R/dft_accuracy_labels.R @@ -0,0 +1,191 @@ +#' Check a reference-label table against the accuracy-assessment contract +#' +#' drift does not store reference labels; the caller does, in whatever review +#' tool they use. This is the contract such a table must meet before +#' [dft_accuracy_estimate()] will read it, so a review tool can be checked +#' against it directly. +#' +#' @param labels A data frame with one row per labelled sample point. +#' @param strata Optional. The `$strata` table from [dft_accuracy_sample()], or +#' any data frame with columns `stratum` and `n_cells`. When given, the +#' labels are also checked for coverage of it. +#' +#' @return `labels`, invisibly. Every fault is an error that names it. +#' +#' @section Columns: +#' \tabular{lll}{ +#' `point_id` \tab required \tab unique, non-missing sample point id \cr +#' `stratum` \tab required \tab the stratum the point was drawn from \cr +#' `map_class` \tab required \tab the map's class at the point \cr +#' `ref_class` \tab required \tab the reference class the reviewer assigned \cr +#' `confidence` \tab optional \tab reviewer confidence, carried and not read \cr +#' `reviewer` \tab optional \tab who labelled it, carried and not read \cr +#' `use` \tab optional \tab `"accuracy"`, `"training"` or `NA` +#' } +#' +#' `stratum` and `map_class` come from the design, not from the reviewer: +#' [dft_accuracy_sample()] writes both (`map_class` through its `map` +#' argument). Only `ref_class`, and the optional columns, are the reviewer's. +#' +#' @section What is refused, and why: +#' - **A missing or blank `ref_class`** (`read.csv()` reads an empty cell as +#' `""`, not `NA`; both are refused). A point the reviewer could not label (cloud, no +#' imagery) is nonresponse, and dropping it changes the stratum's sample +#' size. That is a design decision, so the caller makes it explicitly -- +#' remove the rows and say so -- rather than having it happen silently here. +#' - **A duplicated `point_id`**, which would count one point twice. +#' - **A stratum absent from `strata`**, a stratum with no cells that +#' nevertheless has labels, or a stratum with cells but no labels: either +#' way the labels do not cover the population the weights describe, and no +#' estimator can repair that. +#' - **A `use` value other than `"accuracy"`, `"training"` or `NA`.** +#' +#' Filtering by `confidence` before estimating is possible and is also a +#' design change: the kept points are no longer a random sample of their +#' stratum when low confidence is not random (it rarely is -- edges, mixed +#' cells). Report it if you do it. +#' +#' @section Change maps: +#' For accuracy of a transition map, `map_class` is the transition id and +#' `ref_class` is composed from the reviewer's two endpoint labels with the +#' same scheme [dft_rast_transition()] uses: `ref_from * 1000 + ref_to`. The +#' contract needs no endpoint columns; carry them alongside if useful. +#' +#' @seealso [dft_accuracy_estimate()], which calls this; +#' [dft_accuracy_sample()], which writes `point_id`, `stratum` and +#' `map_class`. +#' +#' @export +#' @examples +#' labels <- data.frame( +#' point_id = c("1_00001", "1_00002", "2_00001", "2_00002"), +#' stratum = c(1, 1, 2, 2), +#' map_class = c(1, 1, 2, 2), +#' ref_class = c(1, 2, 2, 2) +#' ) +#' dft_accuracy_labels(labels) +#' +#' # a point nobody could label is refused, not dropped +#' labels$ref_class[2] <- NA +#' try(dft_accuracy_labels(labels)) +dft_accuracy_labels <- function(labels, strata = NULL) { + if (!is.data.frame(labels)) { + stop("`labels` must be a data frame, not ", class(labels)[1], ".", + call. = FALSE) + } + cols_required <- c("point_id", "stratum", "map_class", "ref_class") + missing_cols <- setdiff(cols_required, names(labels)) + if (length(missing_cols) > 0L) { + stop("`labels` is missing required column(s): ", + paste(missing_cols, collapse = ", "), ".", call. = FALSE) + } + if (nrow(labels) == 0L) { + stop("`labels` has no rows.", call. = FALSE) + } + for (col in cols_required) { + n_na <- sum(accuracy_is_missing(labels[[col]])) + if (n_na > 0L) { + hint <- if (col == "ref_class") { + paste0(" A point that could not be labelled is nonresponse; remove ", + "it deliberately (it changes the stratum's sample size) ", + "rather than passing NA.") + } else { + "" + } + stop("`", col, "` has ", n_na, " missing or blank value(s).", hint, + call. = FALSE) + } + } + # ids are compared as given: numeric normalisation would make "01" and "1" + # the same point + ids <- if (is.factor(labels$point_id)) as.character(labels$point_id) else + labels$point_id + dup <- unique(ids[duplicated(ids)]) + if (length(dup) > 0L) { + stop("`point_id` must be unique; duplicated: ", + paste(utils::head(dup, 5), collapse = ", "), + if (length(dup) > 5L) ", ..." else "", ".", call. = FALSE) + } + if ("use" %in% names(labels)) { + # a blank cell (read.csv() reads one as "") counts as NA, as the + # contract allows + use <- as.character(labels$use) + bad_use <- setdiff(unique(use[!accuracy_is_missing(use)]), + c("accuracy", "training")) + if (length(bad_use) > 0L) { + stop("`use` must be \"accuracy\", \"training\" or NA; got: ", + paste(bad_use, collapse = ", "), ".", call. = FALSE) + } + } + + if (!is.null(strata)) { + accuracy_check_strata(strata) + lab_strata <- unique(accuracy_key(labels$stratum)) + known <- accuracy_key(strata$stratum) + unknown <- setdiff(lab_strata, known) + if (length(unknown) > 0L) { + stop("`labels` has stratum value(s) absent from `strata`: ", + paste(unknown, collapse = ", "), ".", call. = FALSE) + } + empty <- intersect(lab_strata, known[strata$n_cells == 0]) + if (length(empty) > 0L) { + stop("`labels` has points in stratum/strata with no cells in ", + "`strata`: ", paste(empty, collapse = ", "), ". The labels and ", + "the strata table come from different draws.", call. = FALSE) + } + uncovered <- setdiff(known[strata$n_cells > 0], lab_strata) + if (length(uncovered) > 0L) { + stop("Stratum/strata with cells but no labels: ", + paste(uncovered, collapse = ", "), ". The labels do not cover ", + "the population the weights describe.", call. = FALSE) + } + } + invisible(labels) +} + +# Shared by dft_accuracy_labels() and dft_accuracy_estimate() +accuracy_check_strata <- function(strata) { + if (!is.data.frame(strata)) { + stop("`strata` must be a data frame, such as `$strata` from ", + "dft_accuracy_sample().", call. = FALSE) + } + missing_cols <- setdiff(c("stratum", "n_cells"), names(strata)) + if (length(missing_cols) > 0L) { + stop("`strata` is missing required column(s): ", + paste(missing_cols, collapse = ", "), ".", call. = FALSE) + } + if (any(accuracy_is_missing(strata$stratum)) || + anyDuplicated(accuracy_key(strata$stratum))) { + stop("`strata$stratum` must be unique and non-missing.", call. = FALSE) + } + n_ok <- is.numeric(strata$n_cells) && !anyNA(strata$n_cells) && + all(strata$n_cells >= 0) + if (!n_ok) { + stop("`strata$n_cells` must be non-negative numbers.", call. = FALSE) + } + invisible(strata) +} + +# The one string form of a class or stratum value, used for every +# comparison. as.character() writes the double 100000 as "1e+05" and the +# integer as "100000", so a double column and an integer column holding the +# same code would not match; whole numbers are written out in full instead. +accuracy_key <- function(x) { + out <- as.character(x) + # a factor or character value can carry R's own rendering of a number: + # levels(factor(100000)) is "1e+05". Read such strings back as numbers so + # they meet the numeric column holding the same code. + num <- if (is.numeric(x)) as.numeric(x) else suppressWarnings(as.numeric(out)) + whole <- !is.na(num) & is.finite(num) & num == round(num) & abs(num) < 2^53 + out[whole] <- sprintf("%.0f", num[whole]) + out +} + +# NA, or a blank string (a CSV reader's empty cell) +accuracy_is_missing <- function(x) { + miss <- is.na(x) + if (is.character(x) || is.factor(x)) { + miss <- miss | !nzchar(trimws(as.character(x))) + } + miss +} diff --git a/R/dft_accuracy_sample.R b/R/dft_accuracy_sample.R new file mode 100644 index 0000000..ccfc06a --- /dev/null +++ b/R/dft_accuracy_sample.R @@ -0,0 +1,368 @@ +#' Draw a stratified random sample of points for accuracy assessment +#' +#' Draw random cells within each stratum of a raster -- the map classes, a +#' transition map, or strata the caller built, such as "change attributed to +#' fire / unattributed / stable" -- as the sample design behind +#' [dft_accuracy_estimate()], following Olofsson et al. (2014). +#' +#' Points, not patches: sampling patches over-weights large ones, and area is +#' the quantity being estimated. +#' +#' @param strata A single-layer [terra::SpatRaster] of integer stratum codes +#' in a projected CRS. `NA` cells are outside the population. A factor +#' raster (e.g. `$raster` from [dft_rast_transition()] or +#' [dft_rast_break_class()]) keeps its codes as `stratum` and its labels as +#' `stratum_label`. +#' @param n Sample size per stratum: a single number for equal allocation, or +#' a vector named by stratum code (e.g. from [dft_accuracy_size()]) covering +#' every stratum present. Each must be at least 2. +#' @param seed Integer seed. Required: the draw is part of the record, and a +#' sample nobody can redraw cannot be audited. +#' @param map Optional. The map being assessed, on the same grid as `strata`: +#' a single-layer SpatRaster (becomes `map_class`), a multi-layer SpatRaster +#' or a named list of single-layer SpatRasters (become `map_`, e.g. +#' `map_2017` for a series). Values are read as raw codes, so a transition +#' map gives its `from * 1000 + to` id. +#' +#' @return A list: +#' - `points` -- `sf` points at cell centres: `point_id`, `stratum`, +#' `stratum_label`, `cell`, and any `map_*` columns. The geometry is the +#' location; there are no coordinate columns to disagree with it. +#' - `strata` -- tibble: `stratum`, `stratum_label`, `n_cells` (`N_h`), `area` +#' (ha), `weight` (`N_h / N`), `n` (points drawn). This is the `strata` +#' argument [dft_accuracy_estimate()] takes. +#' - `design` -- list recording how to redraw it: seed, RNG kinds, allocation +#' requested, census strata, R and terra versions, and the grid's +#' dimensions, extent and CRS. +#' +#' @section Reproducible, and extensible from a pilot: +#' The draw uses only R's own random number generator and the raster's cell +#' order -- not [terra::spatSample()], whose output is not promised stable +#' across terra versions. Each stratum gets its own random stream, seeded from +#' `seed` and the stratum code, with the generator kinds pinned +#' (Mersenne-Twister, Inversion, Rejection), so: +#' - the same `seed`, `n` and raster give the same points on any machine; +#' - **raising `n` extends the sample**: the first 30 points of a stratum at +#' `n = 50` are the 30 points drawn at `n = 30`, so labels from a pilot +#' carry into the full sample (`point_id` is stable too); +#' - changing one stratum's `n`, or adding a stratum, leaves every other +#' stratum's points unchanged. +#' +#' The caller's random number state is restored afterwards. +#' +#' `cell` is a cell number on this grid. Cropping or extending the raster +#' renumbers cells, so draw from the grid the map is on and keep it. +#' +#' @section Small strata: +#' A stratum with no more cells than its allocation is taken whole -- a +#' census -- and a message names it. Transition strata with a handful of cells +#' are normal. [dft_accuracy_estimate()] applies the finite population +#' correction, so a census stratum contributes no variance. +#' +#' @section Memory: +#' The raster is read twice in row chunks of about ten million cells -- once +#' to count cells per stratum, once to find the drawn cells -- so a +#' floodplain-scale raster is never held in memory whole. +#' +#' @section Polygon strata: +#' Rasterise onto the map's grid first, in memory, and mask to the map's +#' footprint so the population is the mapped area: +#' `strata <- terra::mask(terra::rasterize(polys, map, field = "stratum"), map)`. +#' Do not rasterise straight to a file with an integer `datatype`: terra then +#' writes cells no polygon covers as 0 rather than `NA`, and 0 becomes a +#' stratum. +#' +#' @references +#' Olofsson, P., Foody, G.M., Herold, M., Stehman, S.V., Woodcock, C.E. and +#' Wulder, M.A. (2014). Good practices for estimating area and assessing +#' accuracy of land change. *Remote Sensing of Environment* 148, 42-57. +#' \doi{10.1016/j.rse.2014.02.015} +#' +#' @seealso [dft_accuracy_estimate()], [dft_accuracy_size()], +#' [dft_accuracy_labels()]. +#' +#' @export +#' @examples +#' map <- terra::rast(system.file("extdata", "example_2017.tif", +#' package = "drift")) +#' +#' # strata are the map classes; two strata here have only 2 cells, which are +#' # taken whole +#' s <- dft_accuracy_sample(map, n = 20, seed = 81, map = map) +#' s$strata +#' head(s$points) +#' +#' # a pilot of 20 per stratum extends to 40 without moving the first 20 +#' s40 <- dft_accuracy_sample(map, n = 40, seed = 81) +#' all(s$points$cell %in% s40$points$cell) +dft_accuracy_sample <- function(strata, n, seed, map = NULL) { + if (!inherits(strata, "SpatRaster")) { + stop("`strata` must be a SpatRaster.", call. = FALSE) + } + if (terra::nlyr(strata) != 1L) { + stop("`strata` must have one layer; it has ", terra::nlyr(strata), ".", + call. = FALSE) + } + dft_check_crs(strata, "dft_accuracy_sample") + if (missing(seed) || !is.numeric(seed) || length(seed) != 1L || + is.na(seed) || seed != trunc(seed)) { + stop("`seed` must be a single whole number; it is required so the ", + "sample can be redrawn.", call. = FALSE) + } + map_list <- accuracy_map_list(map, strata) + + # pass 1: cells per stratum, and a refusal of non-integer codes + counts <- accuracy_scan(strata, function(v, cell0, acc) { + v <- v[!is.na(v)] + if (length(v) == 0L) return(acc) + if (any(v != round(v))) { + stop("`strata` must hold integer stratum codes; found ", + v[v != round(v)][1], ".", call. = FALSE) + } + u <- unique(v) + new <- setdiff(u, acc$code) + acc$code <- c(acc$code, new) + acc$count <- c(acc$count, numeric(length(new))) + acc$count <- acc$count + tabulate(match(v, acc$code), length(acc$code)) + acc + }, list(code = numeric(0), count = numeric(0))) + if (length(counts$code) == 0L) { + stop("`strata` has no non-NA cells.", call. = FALSE) + } + o <- order(counts$code) + code <- counts$code[o] + n_cells <- counts$count[o] + + n_req <- accuracy_allocation(n, code) + census <- n_req >= n_cells + n_draw <- pmin(n_req, n_cells) + if (any(census)) { + message("Taking all cells (a census) of ", sum(census), + " stratum/strata with no more cells than allocated: ", + paste0(code[census], " (", n_cells[census], ")", collapse = ", "), + ".") + } + + ranks <- accuracy_draw(code, n_cells, n_draw, seed) + + # pass 2: resolve each stratum's within-stratum ranks to cell numbers + wanted <- lapply(ranks, sort) + found <- accuracy_scan(strata, function(v, cell0, acc) { + keep <- which(!is.na(v)) + vv <- v[keep] + in_chunk <- tabulate(match(vv, code), length(code)) + lo <- acc$seen + hi <- acc$seen + in_chunk + hit <- vapply(seq_along(code), function(h) { + w <- wanted[[h]] + any(w > lo[h] & w <= hi[h]) + }, logical(1)) + if (any(hit)) { + ord <- order(vv, method = "radix") + start <- match(code, vv[ord]) + for (h in which(hit)) { + w <- wanted[[h]] + local <- w[w > lo[h] & w <= hi[h]] - lo[h] + pos <- keep[ord[start[h] + local - 1L]] + acc$cells[[h]] <- c(acc$cells[[h]], cell0 + pos) + acc$rank[[h]] <- c(acc$rank[[h]], local + lo[h]) + } + } + acc$seen <- hi + acc + }, list(seen = numeric(length(code)), + cells = vector("list", length(code)), + rank = vector("list", length(code)))) + if (!identical(as.numeric(found$seen), as.numeric(n_cells))) { + stop("Internal error: the two passes over `strata` counted different ", + "cells. Was the raster modified during the draw?", call. = FALSE) + } + + labels <- accuracy_stratum_labels(strata, code) + pts <- lapply(seq_along(code), function(h) { + # back to draw order, so point_id k is the k-th point drawn in its stratum + cells <- found$cells[[h]][match(ranks[[h]], found$rank[[h]])] + data.frame( + point_id = sprintf("%s_%05d", accuracy_key(code[h]), + seq_along(cells)), + stratum = code[h], + stratum_label = labels[h], + cell = cells, + stringsAsFactors = FALSE + ) + }) + pts <- do.call(rbind, pts) + for (nm in names(map_list)) { + pts[[nm]] <- as.vector(terra::extract(map_list[[nm]], pts$cell, + raw = TRUE)[, 1]) + } + xy <- terra::xyFromCell(strata, pts$cell) + points <- sf::st_as_sf( + cbind(pts, x = xy[, 1], y = xy[, 2]), + coords = c("x", "y"), + crs = terra::crs(strata) + ) + + cell_area_ha <- prod(terra::res(strata)) * 1e-4 + strata_tbl <- tibble::tibble( + stratum = code, + stratum_label = labels, + n_cells = n_cells, + area = n_cells * cell_area_ha, + weight = n_cells / sum(n_cells), + n = n_draw + ) + + design <- list( + seed = seed, + rng_kind = c(kind = "Mersenne-Twister", normal.kind = "Inversion", + sample.kind = "Rejection"), + n_requested = stats::setNames(n_req, accuracy_key(code)), + census = code[census], + r_version = R.version.string, + terra_version = as.character(utils::packageVersion("terra")), + dims = c(nrow = terra::nrow(strata), ncol = terra::ncol(strata)), + extent = as.vector(terra::ext(strata)), + crs = terra::crs(strata) + ) + + list(points = points, strata = strata_tbl, design = design) +} + +# Read a single-layer raster in row chunks, folding `fun(values, cell0, acc)` +# over them. cell0 is the cell number before the chunk's first cell. +accuracy_scan <- function(r, fun, acc) { + nr <- terra::nrow(r) + nc <- terra::ncol(r) + rows <- getOption("drift.accuracy_rows_chunk", + max(1L, floor(1e7 / nc))) + terra::readStart(r) + on.exit(terra::readStop(r), add = TRUE) + row <- 1L + while (row <= nr) { + k <- min(rows, nr - row + 1L) + v <- terra::readValues(r, row = row, nrows = k) + acc <- fun(v, (row - 1) * nc, acc) + row <- row + k + } + acc +} + +# Resolve `n` to one allocation per stratum code +accuracy_allocation <- function(n, code) { + if (!is.numeric(n) || anyNA(n) || any(n != trunc(n))) { + stop("`n` must be whole numbers.", call. = FALSE) + } + key <- accuracy_key(code) + if (length(n) == 1L && is.null(names(n))) { + out <- rep(n, length(code)) + } else { + if (is.null(names(n)) || anyNA(names(n)) || anyDuplicated(names(n))) { + stop("A per-stratum `n` must be named by stratum code, uniquely.", + call. = FALSE) + } + # names are strings already: normalise numeric-looking ones through the + # same key, so a name written "1e+05" by setNames() still finds 100000 + nn <- names(n) + num <- suppressWarnings(as.numeric(nn)) + nn[!is.na(num)] <- accuracy_key(num[!is.na(num)]) + if (anyDuplicated(nn)) { + stop("A per-stratum `n` names one stratum twice.", call. = FALSE) + } + names(n) <- nn + absent <- setdiff(key, names(n)) + if (length(absent) > 0L) { + stop("`n` has no allocation for stratum/strata present in `strata`: ", + paste(absent, collapse = ", "), ". Every stratum with cells must ", + "be sampled, or the estimate cannot cover it.", call. = FALSE) + } + extra <- setdiff(names(n), key) + if (length(extra) > 0L) { + stop("`n` names stratum/strata with no cells in `strata`: ", + paste(extra, collapse = ", "), ".", call. = FALSE) + } + out <- unname(n[key]) + } + if (any(out < 2)) { + stop("Every stratum needs at least 2 points (its variance is undefined ", + "with 1).", call. = FALSE) + } + as.numeric(out) +} + +# One independent, pinned random stream per stratum. Restores the caller's +# RNG state, including its kinds, or removes .Random.seed if there was none. +accuracy_draw <- function(code, n_cells, n_draw, seed) { + env <- globalenv() + had_seed <- exists(".Random.seed", envir = env, inherits = FALSE) + if (had_seed) old_seed <- get(".Random.seed", envir = env, inherits = FALSE) + old_kind <- RNGkind() + on.exit({ + if (had_seed) { + assign(".Random.seed", old_seed, envir = env) # nolint: object_name_linter. R's own name + } else { + RNGkind(old_kind[1], old_kind[2], old_kind[3]) + if (exists(".Random.seed", envir = env, inherits = FALSE)) { + rm(".Random.seed", envir = env) + } + } + }, add = TRUE) + lapply(seq_along(code), function(h) { + # a census still comes from the stream: returning cell order here would + # renumber a stratum's points once a larger n tipped it into a census, + # and pilot labels joined by point_id would land on other cells + set.seed(accuracy_stream_seed(seed, code[h]), kind = "Mersenne-Twister", + normal.kind = "Inversion", sample.kind = "Rejection") + # useHash pinned: by default sample.int() switches algorithm on n and + # size, which would move a draw when a pilot's n grows. The non-hash path + # is a partial shuffle, so a larger size extends a smaller one; it holds + # N_h integers, a few MB at floodplain scale + sample.int(n_cells[h], n_draw[h], useHash = FALSE) + }) +} + +# A stratum's stream seed: a stable 32-bit hash of the seed and the code +accuracy_stream_seed <- function(seed, code) { + digest::digest2int(paste0("drift-accuracy:", accuracy_key(seed), ":", + accuracy_key(code))) +} + +accuracy_stratum_labels <- function(strata, code) { + if (!terra::is.factor(strata)[1]) return(rep(NA_character_, length(code))) + lv <- terra::levels(strata)[[1]] + as.character(lv[[2]][match(code, lv[[1]])]) +} + +# Normalise `map` to a named list of single-layer rasters on the strata grid +accuracy_map_list <- function(map, strata) { + if (is.null(map)) return(list()) + if (inherits(map, "SpatRaster")) { + if (terra::nlyr(map) == 1L) { + out <- list(map_class = map) + } else { + out <- terra::as.list(map) + names(out) <- paste0("map_", names(map)) + } + } else if (is.list(map) && length(map) > 0L && !is.null(names(map)) && + all(nzchar(names(map))) && + all(vapply(map, inherits, logical(1), "SpatRaster"))) { + out <- map + names(out) <- paste0("map_", names(map)) + } else { + stop("`map` must be a SpatRaster or a named list of SpatRasters.", + call. = FALSE) + } + if (anyDuplicated(names(out))) { + stop("`map` layer names must be unique.", call. = FALSE) + } + for (nm in names(out)) { + if (terra::nlyr(out[[nm]]) != 1L) { + stop("Each raster in a `map` list must have one layer.", call. = FALSE) + } + if (!isTRUE(terra::compareGeom(strata, out[[nm]], stopOnError = FALSE))) { + stop("`map` (", nm, ") is not on the `strata` grid; cell numbers would ", + "point at different ground.", call. = FALSE) + } + } + out +} diff --git a/R/dft_accuracy_size.R b/R/dft_accuracy_size.R new file mode 100644 index 0000000..75ee90c --- /dev/null +++ b/R/dft_accuracy_size.R @@ -0,0 +1,179 @@ +#' Size and allocate a stratified accuracy sample +#' +#' How many reference points a stratified sample needs to hit a target +#' standard error, and how to split them across strata -- Olofsson et al. +#' (2014), Eq. 13 and section 5.1. +#' +#' @param weights Numeric stratum weights (`N_h / N`), summing to 1, named by +#' stratum -- `$strata$weight` from [dft_accuracy_sample()], or from +#' [dft_accuracy_estimate()]'s `$stratum` table. +#' @param se_target The standard error to achieve for the estimate being +#' planned for, as a proportion (0.01 is one percentage point). For an +#' area, that is the class's share of total area. +#' @param s_h Per-stratum standard deviation of the indicator behind that +#' estimate, in the order of `weights` (or matched by name when both are +#' named). From a pilot, take `sd` from +#' [dft_accuracy_estimate()]`$stratum` for the target: `"agreement"` to plan +#' for overall accuracy, or a reference class to plan for its area. +#' @param ua Alternatively, anticipated user's accuracy per stratum. Valid +#' only when **the strata are the map classes**, where the SD of the +#' agreement indicator in stratum `i` is `sqrt(ua_i * (1 - ua_i))` +#' (Olofsson Eq. 13). Give `s_h` or `ua`, not both. +#' @param allocation How to split `n` across strata: +#' - `"equal"` -- `n / H` each; +#' - `"proportional"` -- `n * weight`; +#' - `"proportional_min"` (default) -- at least `n_min` in every stratum, +#' the remainder proportional to weight among the others. This is +#' Olofsson's recommendation for change maps, where the change strata are +#' rare and proportional allocation would give them a handful of points. +#' @param n_min Minimum points per stratum for `"proportional_min"`. Olofsson +#' suggests 50-100 per change stratum. Default 50. +#' +#' @return A list: `n`, the total from Eq. 13 (rounded up); `allocation`, the +#' per-stratum sizes named by stratum, ready for [dft_accuracy_sample()]'s +#' `n`; and the `s_h` used. Allocations are rounded, so they can sum to a +#' point or two either side of `n`. +#' +#' @details +#' Eq. 13 is `n = (sum(W_h * S_h) / SE)^2`. Its finite-population term is +#' dropped, which is safe when strata hold millions of cells and conservative +#' otherwise. +#' +#' The allocation changes which estimates are precise, not whether they are +#' unbiased: any allocation with at least two points per stratum gives +#' unbiased estimates through [dft_accuracy_estimate()]. Check the +#' anticipated standard errors of the estimates that matter -- the rare change +#' classes' areas, typically -- rather than overall accuracy alone. +#' +#' `s_h` of 0 (a pilot stratum where every point agreed) is legitimate and +#' makes that stratum contribute nothing to `n`; the minimum allocation still +#' samples it. A pilot that small understates the stratum's variance, so +#' treat a zero with suspicion. +#' +#' @references +#' Olofsson, P., Foody, G.M., Herold, M., Stehman, S.V., Woodcock, C.E. and +#' Wulder, M.A. (2014). Good practices for estimating area and assessing +#' accuracy of land change. *Remote Sensing of Environment* 148, 42-57. +#' \doi{10.1016/j.rse.2014.02.015} +#' +#' @seealso [dft_accuracy_sample()], [dft_accuracy_estimate()]. +#' +#' @export +#' @examples +#' # Olofsson et al. (2014) section 5.1.1: four strata that are the map +#' # classes, anticipated user's accuracies, and a target SE of 0.01 for +#' # overall accuracy -> n = 641 +#' w <- c(deforestation = 0.020, forest_gain = 0.015, +#' stable_forest = 0.320, stable_nonforest = 0.645) +#' dft_accuracy_size(w, se_target = 0.01, ua = c(0.70, 0.60, 0.90, 0.95), +#' allocation = "proportional") +#' +#' # From a pilot: size a full sample for the area of one class +#' map <- terra::rast(system.file("extdata", "example_2017.tif", +#' package = "drift")) +#' ref <- terra::rast(system.file("extdata", "example_2023.tif", +#' package = "drift")) +#' pilot <- dft_accuracy_sample(map, n = 20, seed = 1, map = map) +#' pts <- sf::st_drop_geometry(pilot$points) +#' pts$ref_class <- terra::values(ref)[pts$cell, 1] # stand-in reference +#' est <- dft_accuracy_estimate(pts, pilot$strata) +#' trees <- est$stratum[est$stratum$target == "2", ] # IO LULC 2 = Trees +#' plan <- dft_accuracy_size(stats::setNames(trees$weight, trees$stratum), +#' se_target = 0.02, s_h = trees$sd, n_min = 20) +#' plan$allocation +dft_accuracy_size <- function(weights, se_target, s_h = NULL, ua = NULL, + allocation = c("proportional_min", "equal", + "proportional"), + n_min = 50) { + allocation <- match.arg(allocation) + if (!is.numeric(weights) || length(weights) == 0L || anyNA(weights)) { + stop("`weights` must be a non-empty numeric vector without NA.", + call. = FALSE) + } + if (any(weights <= 0)) { + stop("Every weight must be positive; a stratum with no cells has ", + "nothing to sample. Drop it.", call. = FALSE) + } + if (abs(sum(weights) - 1) > 1e-6) { + stop("`weights` must sum to 1; they sum to ", signif(sum(weights), 6), + ".", call. = FALSE) + } + if (!is.numeric(se_target) || length(se_target) != 1L || + is.na(se_target) || se_target <= 0) { + stop("`se_target` must be a single positive number.", call. = FALSE) + } + if (is.null(s_h) == is.null(ua)) { + stop("Give exactly one of `s_h` or `ua`.", call. = FALSE) + } + if (!is.null(ua)) { + if (!is.numeric(ua) || anyNA(ua) || any(ua < 0 | ua > 1)) { + stop("`ua` must be proportions between 0 and 1.", call. = FALSE) + } + s_h <- sqrt(ua * (1 - ua)) + what <- "ua" + } else { + if (!is.numeric(s_h) || anyNA(s_h) || any(s_h < 0)) { + stop("`s_h` must be non-negative numbers without NA. A stratum with ", + "one pilot point has no SD; label another.", call. = FALSE) + } + what <- "s_h" + } + if (length(s_h) != length(weights)) { + stop("`", what, "` must have one value per stratum (", length(weights), + "); it has ", length(s_h), ".", call. = FALSE) + } + # match by name when both are named, so a differently ordered vector is not + # paired with the wrong weights + vals <- if (what == "ua") ua else s_h + if (!is.null(names(vals)) && !is.null(names(weights))) { + nv <- accuracy_key(names(vals)) + nw <- accuracy_key(names(weights)) + if (!setequal(nv, nw) || anyDuplicated(nv)) { + stop("`", what, "` and `weights` are both named but name different ", + "strata.", call. = FALSE) + } + s_h <- s_h[match(nw, nv)] + } + ws <- sum(weights * s_h) + if (ws == 0) { + stop("Every stratum has an SD of 0, so Eq. 13 asks for no sample at all. ", + "A pilot that agreed everywhere is too small to plan from.", + call. = FALSE) + } + n <- ceiling((ws / se_target)^2) + + nm <- names(weights) + alloc <- switch( + allocation, + equal = rep(round(n / length(weights)), length(weights)), + proportional = round(n * weights), + proportional_min = accuracy_alloc_min(n, weights, n_min) + ) + names(alloc) <- nm + names(s_h) <- nm + list(n = n, allocation = alloc, s_h = s_h) +} + +# At least n_min per stratum, the rest proportional among the others. Strata +# pushed below n_min by the reallocation join the minimum group in turn. +accuracy_alloc_min <- function(n, weights, n_min) { + if (!is.numeric(n_min) || length(n_min) != 1L || is.na(n_min) || + n_min < 2 || n_min != trunc(n_min)) { + stop("`n_min` must be a single whole number of at least 2.", + call. = FALSE) + } + h <- length(weights) + if (n_min * h > n) { + return(rep(n_min, h)) + } + small <- rep(FALSE, h) + repeat { + rest <- n - n_min * sum(small) + share <- rest * weights / sum(weights[!small]) + grow <- !small & share < n_min + if (!any(grow)) break + small <- small | grow + } + out <- ifelse(small, n_min, round(share)) + unname(out) +} diff --git a/man/dft_accuracy_estimate.Rd b/man/dft_accuracy_estimate.Rd new file mode 100644 index 0000000..e542628 --- /dev/null +++ b/man/dft_accuracy_estimate.Rd @@ -0,0 +1,137 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dft_accuracy_estimate.R +\name{dft_accuracy_estimate} +\alias{dft_accuracy_estimate} +\title{Accuracy and error-adjusted area from a stratified reference sample} +\usage{ +dft_accuracy_estimate(labels, strata, level = 0.95) +} +\arguments{ +\item{labels}{A data frame meeting the \code{\link[=dft_accuracy_labels]{dft_accuracy_labels()}} contract: +\code{point_id}, \code{stratum}, \code{map_class}, \code{ref_class}.} + +\item{strata}{The \verb{$strata} table from \code{\link[=dft_accuracy_sample]{dft_accuracy_sample()}}, or any data +frame with \code{stratum}, \code{n_cells} (the stratum's size in cells, \code{n_cells_h}) and +\code{area} (its mapped area in hectares).} + +\item{level}{Confidence level for the intervals. Default \code{0.95}.} +} +\value{ +A list: +\itemize{ +\item \code{matrix} -- tibble, long: \code{map_class}, \code{ref_class}, \code{proportion} (the +estimated share of total area), one row per pair of classes, zeros +included. +\item \code{accuracy} -- tibble: \code{measure} (\code{"overall"}, \code{"user"}, \code{"producer"}), +\code{class} (\code{NA} for overall), \code{estimate}, \code{se}, \code{lower}, \code{upper}. +\item \code{area} -- tibble, one row per class: \code{proportion}, \code{proportion_se}, +\code{area}, \code{area_se}, \code{lower}, \code{upper} (hectares). +\item \code{stratum} -- tibble, long, per stratum and \code{target}: \code{n}, \code{n_cells}, +\code{weight}, \code{mean} and \code{sd} of an indicator within the stratum. \code{target} is +\code{"agreement"} (map equals reference) or a reference class, written in +full as a string (\code{"100000"}, never \code{"1e+05"}). \code{sd} is what +\code{\link[=dft_accuracy_size]{dft_accuracy_size()}} takes to size a full sample from a pilot. +\item \code{level}, \code{area_total} (ha). +} +} +\description{ +Turn reference labels at stratified random sample points into the numbers +a change map should be reported with: an error matrix in estimated +proportions of area, overall / user's / producer's accuracy, and +\strong{error-adjusted area per class with a confidence interval} -- following +Olofsson et al. (2014) and Stehman (2014). +} +\details{ +Mapped area is a biased estimator whenever the map has errors, and for a +change map the bias is usually large: a single wrong label on either date +manufactures a transition. The error-adjusted area is the area the +\emph{reference} labels imply, weighted back to the population by each +stratum's size. +} +\section{Estimators}{ + +The estimates come from \code{\link[mapaccuracy:stehman2014]{mapaccuracy::stehman2014()}}, which implements +Stehman (2014)'s estimators for stratified random sampling. Those hold +\strong{whether or not the strata are the map classes} -- strata such as "change +attributed to fire", "change unattributed" and "stable" are fine. When the +strata are the map classes they give Olofsson et al. (2014)'s estimates, +and the package test suite reproduces that paper's worked example. +\itemize{ +\item The \strong{finite population correction} \code{(1 - n_h / n_cells_h)} is always +applied, so a stratum sampled in full (a census) contributes no variance. +Olofsson et al. omit it; at pixel-scale \code{n_cells_h} the difference is below +anything reported. +\item Intervals are Wald intervals, \verb{estimate +/- z * se} with +\code{z = qnorm(1 - (1 - level) / 2)}, and are \strong{not} truncated at 0 or 1: a +lower bound below 0 says the class is too rare for the sample to bound +away from zero, and truncating it would hide that. +\item \strong{A rare class hiding in a large stratum makes the interval too narrow +at small samples.} If 1\\% of a big stratum is really class \emph{j} and +the stratum gets 25 points, most samples see none of it, estimate that +stratum's variance for \emph{j} as 0, and report an interval that misses. On +the bundled tiles (98 of 7,127 "Trees" cells reference Water) Water's +95\\% interval covered 61\\% of the time at 25 points per stratum, 83\\% +at 75 and 89\\% at 150; the estimate itself stayed unbiased. Omission +hides in the large stable strata, so do not starve them. +\item When the strata are not the map classes, the error matrix's row totals +(the map-class shares) are \strong{estimated} from the sample, not the known +stratum weights. That surprises readers used to Olofsson's tables. +\item Producer's accuracy is \code{NA} for a class no reference label falls in. +\item Run time grows with roughly the 2.5th power of the number of classes: +about 13 s for 1,000 points over 80 transition classes. Collapse classes +you will not report before estimating. +} +} + +\section{Targets that are unions of classes}{ + +"Tree loss" is every Trees -> non-Trees transition; an "unattributed +residual" is loss outside any fire or harvest. Recode \code{map_class} and +\code{ref_class} to the target (say \code{"loss"} / \code{"other"}) and estimate again. The +standard error of a union is \strong{not} the sum of its members' standard +errors, because the members' estimates are correlated. +} + +\section{Training points}{ + +A row with \code{use == "training"} is refused. Points that trained a +classifier cannot also measure it -- the estimate would be optimistic by +construction. Keep the accuracy sample separate from the start; a subset of +accuracy points chosen for training by judgement also stops being a random +sample of its stratum. +} + +\examples{ +# Olofsson et al. (2014), Table 8: 640 labelled points in four strata that +# are the map classes, over 10 million 30 m pixels +cls <- c("deforestation", "forest_gain", "stable_forest", "stable_nonforest") +counts <- matrix(c(66, 0, 5, 4, 0, 55, 8, 12, 1, 0, 153, 11, 2, 1, 9, 313), + nrow = 4, byrow = TRUE) +idx <- which(counts > 0, arr.ind = TRUE) +map <- rep(cls[idx[, 1]], counts[idx]) +ref <- rep(cls[idx[, 2]], counts[idx]) +labels <- data.frame(point_id = seq_along(map), stratum = map, + map_class = map, ref_class = ref) +pixels <- c(200000, 150000, 3200000, 6450000) +strata <- data.frame(stratum = cls, n_cells = pixels, area = pixels * 0.09) + +res <- dft_accuracy_estimate(labels, strata) +res$area # deforestation: 21,158 ha +/- 6,157 -- mapped was 18,000 +res$accuracy +} +\references{ +Olofsson, P., Foody, G.M., Herold, M., Stehman, S.V., Woodcock, C.E. and +Wulder, M.A. (2014). Good practices for estimating area and assessing +accuracy of land change. \emph{Remote Sensing of Environment} 148, 42-57. +\doi{10.1016/j.rse.2014.02.015} + +Stehman, S.V. (2014). Estimating area and map accuracy for stratified +random sampling when the strata are different from the map classes. +\emph{International Journal of Remote Sensing} 35(13), 4923-4939. +\doi{10.1080/01431161.2014.930207} +} +\seealso{ +\code{\link[=dft_accuracy_sample]{dft_accuracy_sample()}} to draw the points, +\code{\link[=dft_accuracy_labels]{dft_accuracy_labels()}} for the label contract, +\code{\link[=dft_accuracy_size]{dft_accuracy_size()}} to size a sample from a pilot. +} diff --git a/man/dft_accuracy_labels.Rd b/man/dft_accuracy_labels.Rd new file mode 100644 index 0000000..a4d97b6 --- /dev/null +++ b/man/dft_accuracy_labels.Rd @@ -0,0 +1,89 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dft_accuracy_labels.R +\name{dft_accuracy_labels} +\alias{dft_accuracy_labels} +\title{Check a reference-label table against the accuracy-assessment contract} +\usage{ +dft_accuracy_labels(labels, strata = NULL) +} +\arguments{ +\item{labels}{A data frame with one row per labelled sample point.} + +\item{strata}{Optional. The \verb{$strata} table from \code{\link[=dft_accuracy_sample]{dft_accuracy_sample()}}, or +any data frame with columns \code{stratum} and \code{n_cells}. When given, the +labels are also checked for coverage of it.} +} +\value{ +\code{labels}, invisibly. Every fault is an error that names it. +} +\description{ +drift does not store reference labels; the caller does, in whatever review +tool they use. This is the contract such a table must meet before +\code{\link[=dft_accuracy_estimate]{dft_accuracy_estimate()}} will read it, so a review tool can be checked +against it directly. +} +\section{Columns}{ + +\tabular{lll}{ +\code{point_id} \tab required \tab unique, non-missing sample point id \cr +\code{stratum} \tab required \tab the stratum the point was drawn from \cr +\code{map_class} \tab required \tab the map's class at the point \cr +\code{ref_class} \tab required \tab the reference class the reviewer assigned \cr +\code{confidence} \tab optional \tab reviewer confidence, carried and not read \cr +\code{reviewer} \tab optional \tab who labelled it, carried and not read \cr +\code{use} \tab optional \tab \code{"accuracy"}, \code{"training"} or \code{NA} +} + +\code{stratum} and \code{map_class} come from the design, not from the reviewer: +\code{\link[=dft_accuracy_sample]{dft_accuracy_sample()}} writes both (\code{map_class} through its \code{map} +argument). Only \code{ref_class}, and the optional columns, are the reviewer's. +} + +\section{What is refused, and why}{ + +\itemize{ +\item \strong{A missing or blank \code{ref_class}} (\code{read.csv()} reads an empty cell as +\code{""}, not \code{NA}; both are refused). A point the reviewer could not label (cloud, no +imagery) is nonresponse, and dropping it changes the stratum's sample +size. That is a design decision, so the caller makes it explicitly -- +remove the rows and say so -- rather than having it happen silently here. +\item \strong{A duplicated \code{point_id}}, which would count one point twice. +\item \strong{A stratum absent from \code{strata}}, a stratum with no cells that +nevertheless has labels, or a stratum with cells but no labels: either +way the labels do not cover the population the weights describe, and no +estimator can repair that. +\item \strong{A \code{use} value other than \code{"accuracy"}, \code{"training"} or \code{NA}.} +} + +Filtering by \code{confidence} before estimating is possible and is also a +design change: the kept points are no longer a random sample of their +stratum when low confidence is not random (it rarely is -- edges, mixed +cells). Report it if you do it. +} + +\section{Change maps}{ + +For accuracy of a transition map, \code{map_class} is the transition id and +\code{ref_class} is composed from the reviewer's two endpoint labels with the +same scheme \code{\link[=dft_rast_transition]{dft_rast_transition()}} uses: \code{ref_from * 1000 + ref_to}. The +contract needs no endpoint columns; carry them alongside if useful. +} + +\examples{ +labels <- data.frame( + point_id = c("1_00001", "1_00002", "2_00001", "2_00002"), + stratum = c(1, 1, 2, 2), + map_class = c(1, 1, 2, 2), + ref_class = c(1, 2, 2, 2) +) +dft_accuracy_labels(labels) + +# a point nobody could label is refused, not dropped +labels$ref_class[2] <- NA +try(dft_accuracy_labels(labels)) +} +\seealso{ +\code{\link[=dft_accuracy_estimate]{dft_accuracy_estimate()}}, which calls this; +\code{\link[=dft_accuracy_sample]{dft_accuracy_sample()}}, which writes \code{point_id}, \code{stratum} and +\code{map_class}. +} diff --git a/man/dft_accuracy_sample.Rd b/man/dft_accuracy_sample.Rd new file mode 100644 index 0000000..b04dd75 --- /dev/null +++ b/man/dft_accuracy_sample.Rd @@ -0,0 +1,123 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dft_accuracy_sample.R +\name{dft_accuracy_sample} +\alias{dft_accuracy_sample} +\title{Draw a stratified random sample of points for accuracy assessment} +\usage{ +dft_accuracy_sample(strata, n, seed, map = NULL) +} +\arguments{ +\item{strata}{A single-layer \link[terra:SpatRaster-class]{terra::SpatRaster} of integer stratum codes +in a projected CRS. \code{NA} cells are outside the population. A factor +raster (e.g. \verb{$raster} from \code{\link[=dft_rast_transition]{dft_rast_transition()}} or +\code{\link[=dft_rast_break_class]{dft_rast_break_class()}}) keeps its codes as \code{stratum} and its labels as +\code{stratum_label}.} + +\item{n}{Sample size per stratum: a single number for equal allocation, or +a vector named by stratum code (e.g. from \code{\link[=dft_accuracy_size]{dft_accuracy_size()}}) covering +every stratum present. Each must be at least 2.} + +\item{seed}{Integer seed. Required: the draw is part of the record, and a +sample nobody can redraw cannot be audited.} + +\item{map}{Optional. The map being assessed, on the same grid as \code{strata}: +a single-layer SpatRaster (becomes \code{map_class}), a multi-layer SpatRaster +or a named list of single-layer SpatRasters (become \verb{map_}, e.g. +\code{map_2017} for a series). Values are read as raw codes, so a transition +map gives its \code{from * 1000 + to} id.} +} +\value{ +A list: +\itemize{ +\item \code{points} -- \code{sf} points at cell centres: \code{point_id}, \code{stratum}, +\code{stratum_label}, \code{cell}, and any \verb{map_*} columns. The geometry is the +location; there are no coordinate columns to disagree with it. +\item \code{strata} -- tibble: \code{stratum}, \code{stratum_label}, \code{n_cells} (\code{N_h}), \code{area} +(ha), \code{weight} (\code{N_h / N}), \code{n} (points drawn). This is the \code{strata} +argument \code{\link[=dft_accuracy_estimate]{dft_accuracy_estimate()}} takes. +\item \code{design} -- list recording how to redraw it: seed, RNG kinds, allocation +requested, census strata, R and terra versions, and the grid's +dimensions, extent and CRS. +} +} +\description{ +Draw random cells within each stratum of a raster -- the map classes, a +transition map, or strata the caller built, such as "change attributed to +fire / unattributed / stable" -- as the sample design behind +\code{\link[=dft_accuracy_estimate]{dft_accuracy_estimate()}}, following Olofsson et al. (2014). +} +\details{ +Points, not patches: sampling patches over-weights large ones, and area is +the quantity being estimated. +} +\section{Reproducible, and extensible from a pilot}{ + +The draw uses only R's own random number generator and the raster's cell +order -- not \code{\link[terra:sample]{terra::spatSample()}}, whose output is not promised stable +across terra versions. Each stratum gets its own random stream, seeded from +\code{seed} and the stratum code, with the generator kinds pinned +(Mersenne-Twister, Inversion, Rejection), so: +\itemize{ +\item the same \code{seed}, \code{n} and raster give the same points on any machine; +\item \strong{raising \code{n} extends the sample}: the first 30 points of a stratum at +\code{n = 50} are the 30 points drawn at \code{n = 30}, so labels from a pilot +carry into the full sample (\code{point_id} is stable too); +\item changing one stratum's \code{n}, or adding a stratum, leaves every other +stratum's points unchanged. +} + +The caller's random number state is restored afterwards. + +\code{cell} is a cell number on this grid. Cropping or extending the raster +renumbers cells, so draw from the grid the map is on and keep it. +} + +\section{Small strata}{ + +A stratum with no more cells than its allocation is taken whole -- a +census -- and a message names it. Transition strata with a handful of cells +are normal. \code{\link[=dft_accuracy_estimate]{dft_accuracy_estimate()}} applies the finite population +correction, so a census stratum contributes no variance. +} + +\section{Memory}{ + +The raster is read twice in row chunks of about ten million cells -- once +to count cells per stratum, once to find the drawn cells -- so a +floodplain-scale raster is never held in memory whole. +} + +\section{Polygon strata}{ + +Rasterise onto the map's grid first, in memory, and mask to the map's +footprint so the population is the mapped area: +\code{strata <- terra::mask(terra::rasterize(polys, map, field = "stratum"), map)}. +Do not rasterise straight to a file with an integer \code{datatype}: terra then +writes cells no polygon covers as 0 rather than \code{NA}, and 0 becomes a +stratum. +} + +\examples{ +map <- terra::rast(system.file("extdata", "example_2017.tif", + package = "drift")) + +# strata are the map classes; two strata here have only 2 cells, which are +# taken whole +s <- dft_accuracy_sample(map, n = 20, seed = 81, map = map) +s$strata +head(s$points) + +# a pilot of 20 per stratum extends to 40 without moving the first 20 +s40 <- dft_accuracy_sample(map, n = 40, seed = 81) +all(s$points$cell \%in\% s40$points$cell) +} +\references{ +Olofsson, P., Foody, G.M., Herold, M., Stehman, S.V., Woodcock, C.E. and +Wulder, M.A. (2014). Good practices for estimating area and assessing +accuracy of land change. \emph{Remote Sensing of Environment} 148, 42-57. +\doi{10.1016/j.rse.2014.02.015} +} +\seealso{ +\code{\link[=dft_accuracy_estimate]{dft_accuracy_estimate()}}, \code{\link[=dft_accuracy_size]{dft_accuracy_size()}}, +\code{\link[=dft_accuracy_labels]{dft_accuracy_labels()}}. +} diff --git a/man/dft_accuracy_size.Rd b/man/dft_accuracy_size.Rd new file mode 100644 index 0000000..8e2515d --- /dev/null +++ b/man/dft_accuracy_size.Rd @@ -0,0 +1,107 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dft_accuracy_size.R +\name{dft_accuracy_size} +\alias{dft_accuracy_size} +\title{Size and allocate a stratified accuracy sample} +\usage{ +dft_accuracy_size( + weights, + se_target, + s_h = NULL, + ua = NULL, + allocation = c("proportional_min", "equal", "proportional"), + n_min = 50 +) +} +\arguments{ +\item{weights}{Numeric stratum weights (\code{N_h / N}), summing to 1, named by +stratum -- \verb{$strata$weight} from \code{\link[=dft_accuracy_sample]{dft_accuracy_sample()}}, or from +\code{\link[=dft_accuracy_estimate]{dft_accuracy_estimate()}}'s \verb{$stratum} table.} + +\item{se_target}{The standard error to achieve for the estimate being +planned for, as a proportion (0.01 is one percentage point). For an +area, that is the class's share of total area.} + +\item{s_h}{Per-stratum standard deviation of the indicator behind that +estimate, in the order of \code{weights} (or matched by name when both are +named). From a pilot, take \code{sd} from +\code{\link[=dft_accuracy_estimate]{dft_accuracy_estimate()}}\verb{$stratum} for the target: \code{"agreement"} to plan +for overall accuracy, or a reference class to plan for its area.} + +\item{ua}{Alternatively, anticipated user's accuracy per stratum. Valid +only when \strong{the strata are the map classes}, where the SD of the +agreement indicator in stratum \code{i} is \code{sqrt(ua_i * (1 - ua_i))} +(Olofsson Eq. 13). Give \code{s_h} or \code{ua}, not both.} + +\item{allocation}{How to split \code{n} across strata: +\itemize{ +\item \code{"equal"} -- \code{n / H} each; +\item \code{"proportional"} -- \code{n * weight}; +\item \code{"proportional_min"} (default) -- at least \code{n_min} in every stratum, +the remainder proportional to weight among the others. This is +Olofsson's recommendation for change maps, where the change strata are +rare and proportional allocation would give them a handful of points. +}} + +\item{n_min}{Minimum points per stratum for \code{"proportional_min"}. Olofsson +suggests 50-100 per change stratum. Default 50.} +} +\value{ +A list: \code{n}, the total from Eq. 13 (rounded up); \code{allocation}, the +per-stratum sizes named by stratum, ready for \code{\link[=dft_accuracy_sample]{dft_accuracy_sample()}}'s +\code{n}; and the \code{s_h} used. Allocations are rounded, so they can sum to a +point or two either side of \code{n}. +} +\description{ +How many reference points a stratified sample needs to hit a target +standard error, and how to split them across strata -- Olofsson et al. +(2014), Eq. 13 and section 5.1. +} +\details{ +Eq. 13 is \code{n = (sum(W_h * S_h) / SE)^2}. Its finite-population term is +dropped, which is safe when strata hold millions of cells and conservative +otherwise. + +The allocation changes which estimates are precise, not whether they are +unbiased: any allocation with at least two points per stratum gives +unbiased estimates through \code{\link[=dft_accuracy_estimate]{dft_accuracy_estimate()}}. Check the +anticipated standard errors of the estimates that matter -- the rare change +classes' areas, typically -- rather than overall accuracy alone. + +\code{s_h} of 0 (a pilot stratum where every point agreed) is legitimate and +makes that stratum contribute nothing to \code{n}; the minimum allocation still +samples it. A pilot that small understates the stratum's variance, so +treat a zero with suspicion. +} +\examples{ +# Olofsson et al. (2014) section 5.1.1: four strata that are the map +# classes, anticipated user's accuracies, and a target SE of 0.01 for +# overall accuracy -> n = 641 +w <- c(deforestation = 0.020, forest_gain = 0.015, + stable_forest = 0.320, stable_nonforest = 0.645) +dft_accuracy_size(w, se_target = 0.01, ua = c(0.70, 0.60, 0.90, 0.95), + allocation = "proportional") + +# From a pilot: size a full sample for the area of one class +map <- terra::rast(system.file("extdata", "example_2017.tif", + package = "drift")) +ref <- terra::rast(system.file("extdata", "example_2023.tif", + package = "drift")) +pilot <- dft_accuracy_sample(map, n = 20, seed = 1, map = map) +pts <- sf::st_drop_geometry(pilot$points) +pts$ref_class <- terra::values(ref)[pts$cell, 1] # stand-in reference +est <- dft_accuracy_estimate(pts, pilot$strata) +trees <- est$stratum[est$stratum$target == "2", ] # IO LULC 2 = Trees +plan <- dft_accuracy_size(stats::setNames(trees$weight, trees$stratum), + se_target = 0.02, s_h = trees$sd, n_min = 20) +plan$allocation +} +\references{ +Olofsson, P., Foody, G.M., Herold, M., Stehman, S.V., Woodcock, C.E. and +Wulder, M.A. (2014). Good practices for estimating area and assessing +accuracy of land change. \emph{Remote Sensing of Environment} 148, 42-57. +\doi{10.1016/j.rse.2014.02.015} +} +\seealso{ +\code{\link[=dft_accuracy_sample]{dft_accuracy_sample()}}, \code{\link[=dft_accuracy_estimate]{dft_accuracy_estimate()}}. +} diff --git a/planning/active/findings.md b/planning/active/findings.md index ec54691..c5681e1 100644 --- a/planning/active/findings.md +++ b/planning/active/findings.md @@ -77,7 +77,35 @@ The values are transcribed into `tests/testthat/helper-accuracy.R` with page and - It only *warns* on a stratum with one observation. - It applies the FPC in eq. 25 and eq. 28. +## Phase 2 measurements (2026-09-29) + +- **Restore-the-bug:** passing `rep(mean(N_h))` as `Nh_strata` inside the wrapper turns 9 of the 50 estimator assertions red (areas, half-widths, Table 9, PA, OA, the union test). Restored. +- **`stehman2014()` runtime** at 1,000 points and 8 strata: 20 classes 0.33 s, 40 classes 2.0 s, 80 classes 13.0 s. That is about k^2.5, from its `classes²` indicator columns and two `aggregate()` calls. It is workable for a one-off estimate, so nothing was filed upstream. The roxygen notes it. +- **testthat 3e ignores `scale =`**, so `tolerance` was relative. The pins use an `expect_within()` absolute helper. + +## Phase 3: sampler measurements (2026-09-29) + +**BULK scale** (`classified_2017.tif`, 14651 × 11552 = 169,248,352 cells, 4,108,972 non-NA; 64 GB machine; RSS sampled every 1–2 s): + +| step | time | peak RSS | +|---|---|---| +| one chunked pass (count) | 1.7 s | | +| `dft_accuracy_sample(r17, n = 100)` (2 passes + draw) | 3.3 s | 0.95 GiB (baseline after load 0.27 GiB) | +| classify 2017+2023 → `dft_rast_transition()` → sample the factor transition map (63 strata, n = 30, 10 censuses) | 1.5 + 3.4 + 3.0 s | 4.63 GiB, dominated by the in-memory transition raster | +| `dft_accuracy_estimate()` on those 1,677 points, 63 classes (perfect labels) | 10.8 s | | + +Logs: scratchpad `bulk/bulk_run{2,3}.log`, `bulk/rss*.log` (not committed). + +**Census oracle (bundled tiles, 300 draws).** Map-class strata at n = 25: the estimate is unbiased (|bias z| < 1), but Water's 95% CI covered 61% (SE ratio 0.76). 98 of the 7,127 map-Trees cells (1.4%) are reference Water, so most draws see none, and that stratum's variance is estimated as 0. Coverage was 83% at n = 75 and 89% at n = 150 (SE ratio 1.04). This is the Wald interval's small-sample weakness, not an estimator bug, and it is documented on `dft_accuracy_estimate()`. With changed/stable strata (not the map classes) at n = 60: |bias z| ≤ 0.83, SE ratio 1.00–1.07, coverage 0.91–0.96. + +**`sample.int(useHash = TRUE)` refuses size > n/2**, so the sampler pins `useHash = FALSE`. The golden draw did not move, because the two paths gave identical draws at the tested sizes. + +**Found on the way:** `dft_rast_classify()` mutates its input raster in place (`set.cats()`), so a test that later stacked the raw tile saw `class_name` layers. Filed as #89; not fixed here. + ## Errors Encountered | Error | Resolution | |-------|------------| +| `sample.int(useHash = TRUE)`: "This algorithm is for size <= n/2" | Pin `useHash = FALSE` (a partial shuffle that keeps its prefix) | +| RSS sampler watched a subshell and wrote `rss.log` into the repo: `cd X && prog &` backgrounds the whole list | `cd` first, then `nohup prog &` alone, with the log at an absolute path | +| `expect_equal(..., tolerance, scale = 1)` warned "Unused arguments (scale = 1)" and compared relatively | `expect_within()` helper with an absolute bound | diff --git a/planning/active/progress.md b/planning/active/progress.md index a0e5725..e177e5f 100644 --- a/planning/active/progress.md +++ b/planning/active/progress.md @@ -8,3 +8,5 @@ - Next: Phase 2 estimator code and tests; published-value pins wait on the PDFs (Phase 1) - Decision: the estimator wraps `mapaccuracy::stehman2014()` (user, 2026-09-28) instead of reimplementing it; plan revised - Phase 1 done: Olofsson 2014 read from the local Zotero PDF; values in `helper-accuracy.R`; three printed values contradict the paper's own equations (findings). Issue body revised. Stehman 2014 is not needed (mapaccuracy tests it) +- Phases 2–4 code and tests written; code-check rounds 1–2 found 6 defects (fixed), round 3 running; BULK scale run 3.3 s / 0.95 GiB for the sampler; filed #89 (classify mutates input) +- Code-check: 3 rounds. R1: 4 findings (1 bug). R2: 2 (1 bug, 1 inside the R1 fix). R3: 6 (2 bugs, 1 inside the R2 fix); ended by enumerating 29 identity sites, and every fix was mutation-checked diff --git a/planning/active/review-round1.md b/planning/active/review-round1.md new file mode 100644 index 0000000..a57c8eb --- /dev/null +++ b/planning/active/review-round1.md @@ -0,0 +1,50 @@ +# Code-check review, round 1 (#81, phase 2 staged diff) + +Scope: R/dft_accuracy_estimate.R, R/dft_accuracy_labels.R, tests/testthat/test-dft_accuracy_estimate.R, +tests/testthat/helper-accuracy.R, DESCRIPTION, NAMESPACE, man/. Read mapaccuracy::stehman2014() (0.1.2, CRAN) +and its .check_labels / .check_length. The suite passes on a scratch copy: `[ FAIL 0 | WARN 0 | SKIP 0 | PASS 50 ]`. +The Olofsson Table 8/9 and section 5.2 transcriptions were checked against the PDF (pp. 13-14) and match. +An independent Olofsson-style computation (12 strata, strata = map classes, so the s1..s12 ids and the +s1-vs-s10 regex are exercised) agrees with the wrapper's area proportions to 1e-17. The SE ratio is +0.997-0.999, which is the FPC. + +## Findings + +- **[severity: bug]** R/dft_accuracy_estimate.R:156 and :180-181. The class set comes from + `c(labels$map_class, labels$ref_class)`, but the labels passed to stehman2014() come from + `as.character()` on each column separately (:157-158). These are two normalisations of the same labels. + When exactly one column is a factor, `c()` swaps that factor for its integer codes. That happens with + numeric + factor, and with factor + character, where `c.factor` falls back to `unlist`. The codes then + become phantom classes. Stehman still gets the right strings, so no error is raised, and the output + carries extra classes with area 0, zero rows in the matrix, NA accuracies, and extra stratum-table + targets. Reproduced: `map_class = c(1001, 2002)` numeric with `ref_class` a factor gives an `area` table + of classes 1, 2, 1001, 2002, and classes 1 and 2 read as "0 ha, SE 0". With IO LULC codes (1, 2, 4, 5, + 7, 8, 11) the phantom codes 3 and 6 look exactly like real classes that have zero area. A factor + `ref_class` is plausible from a review-tool export. The fix is to build the class set from + `c(map_chr, ref_chr)`, and to restore numeric type only when both columns are numeric. + +- **[severity: fragile]** tests/testthat/test-dft_accuracy_estimate.R:131-132. The assertion + `expect_false(isTRUE(all.equal(unname(rows[olofsson_classes]), olofsson_pixels / 1e7)))` cannot fail. + `rows` comes from `tapply()`, so it is a 1-d array, and `unname()` keeps the `dim`. `all.equal()` then + reports "target is array, current is numeric" whatever the values are. Measured: on + `tapply(1:4, letters[1:4], sum)` against an identical numeric vector, `isTRUE(all.equal(...))` is FALSE. + The claim under test does hold (the estimated shares are 0.117 / 0.117 / 0.259 / 0.508 against W_i + 0.020 / 0.015 / 0.320 / 0.645), but the test does not guard it. Wrap `rows[...]` in `as.vector()`. + +- **[severity: fragile]** tests/testthat/helper-accuracy.R:104. `expect_within()`'s length guard compares + `length(d)`, and `d` has already been recycled by the subtraction. So an object shorter than `expected` + that divides into it passes. Measured: `expect_within(c(1, 2), c(1, 2, 1, 2), 0.1)` passes. Its practical + reach with the current distinct published values is nil, because a short object cannot match four + distinct numbers, but the guard does not check what it says it checks. Compare + `length(act$val) == length(expected)`. + +- **[severity: fragile]** R/dft_accuracy_estimate.R:127 with R/dft_accuracy_labels.R:113-120. Labels whose + stratum has `n_cells == 0` pass `dft_accuracy_labels()`: the stratum is known, and the zero-cell stratum + is excluded from "uncovered". Line 127 then drops the stratum before the `n_h > n_cells_h` check, so the + guard written for this case ("More labels than cells") never sees it. `s_lab` becomes NA and the call + fails inside mapaccuracy with `Arguments should include only: s1, s2, s3, s4`, which names internal ids + the caller never passed. The failure is loud, not silent, but the message does not locate the fault. + Count `n_h` against the unfiltered strata, or refuse labels in zero-cell strata in `dft_accuracy_labels()`. + +No security issues. DESCRIPTION/NAMESPACE are consistent: mapaccuracy is in Imports, and rlang and tibble, +used by the helper and the tests, are already there. diff --git a/planning/active/review-round2.md b/planning/active/review-round2.md new file mode 100644 index 0000000..f9f725d --- /dev/null +++ b/planning/active/review-round2.md @@ -0,0 +1,41 @@ +# Code review, round 2 (#81 staged diff: dft_accuracy_estimate / dft_accuracy_labels) + +All four round-1 fixes hold up. Staged-only copy: `devtools::test(filter = "accuracy")` gives +55 pass / 0 fail, and `devtools::document()` leaves man/ and NAMESPACE unchanged. I also checked +stratum order against an independent Stehman (2014) computation: 12 strata, names unsorted, the +strata table in reverse order, label rows shuffled. Area proportions, their SEs, OA, UA and the +stratum table's n / n_cells all match. So the s1..sH remap and mapaccuracy's aggregate() and +table() ordering agree with each other. + +## Findings + +- **[bug]** R/dft_accuracy_labels.R:83-95. A blank `ref_class` gets past the nonresponse + refusal and is estimated as a class named `""`. The check is `is.na()`, but `read.csv()` + and `data.table::fread()` both read an empty cell in a character column as `""`, not `NA`. + Only readr turns it into `NA`. So a reviewer's unlabelled point in a CSV passes the check + the contract says refuses it, and it goes into mapaccuracy as a real reference class. + Reproduced with 6 points in 2 strata, one blank: + - OA comes out 0.667. The blank counts as a disagreement. + - Class `x`'s area is 0.5 of the total, and the missing 1/6 of area goes nowhere. + - The `""` row's `area`, `proportion`, `user` and `producer` are all `NA`, because + `est$area[""]` indexes by name and `x[""]` is always `NA`. So the proportions no longer + sum to 1, and no error is raised. + + These are silently wrong numbers, on exactly the path the "refused, not dropped" section + documents. The same gap applies to `point_id`, `stratum` and `map_class`, though a blank + stratum at least fails loudly as "absent from `strata`". Fix: treat `!nzchar(trimws(x))` + as missing, alongside `is.na(x)`, in character and factor columns. + +- **[fragile]** R/dft_accuracy_estimate.R:159-160 (and :138-139 / dft_accuracy_labels.R:113-114 + for stratum keys). Fix 1 normalises each column with `as.character()`, but `as.character()` + depends on the type. A double with 5+ trailing zeros prints in scientific notation, an + integer never does: `as.character(100000)` is `"1e+05"` and `as.character(100000L)` is + `"100000"`. If `map_class` is double (a raster extract) and `ref_class` is integer (e.g. + `ref_from * 1000L + ref_to`), class 100000 splits into two classes. Reproduced with 6 + points: every point agrees and OA comes out 0.5, and `$area` has two rows with class + `100000`, one of them with zero area. Codes like this only occur when `to = 0` and `from` + is a multiple of 100, so real data should rarely hit it. It is still the round-1 defect + one axis over: two derivations of one class that disagree. For stratum keys, a CSV + round trip (double) against an integer `strata$stratum` fails loudly as "absent from + `strata`", not silently. Fix: for numeric columns, `format(x, scientific = FALSE, + trim = TRUE, digits = 15)` (or convert both to double) before `as.character()`. diff --git a/planning/active/review-round3.md b/planning/active/review-round3.md new file mode 100644 index 0000000..8b1e8f5 --- /dev/null +++ b/planning/active/review-round3.md @@ -0,0 +1,151 @@ +# Code review round 3 (#81): sampler, sizer, and the key mechanism + +Reviewer worked in a scratch copy of the repo; the repo tree was not modified. Every +finding below was reproduced with a probe script, and the results are quoted. + +## Mechanism + +**One identity is derived more than once, each time independently, and the code is +correct only while those derivations happen to agree.** A class, stratum, point or +allocation identity gets re-derived at each use site from whatever representation +reaches it: `as.character()` of a double, the integer codes `c()` returns for a factor, +a factor's levels, `names()` written by `setNames()`, a position in a parallel vector, +or a draw order. Nothing joins on one canonical key. R1 (factor codes vs strings), R2 +(`""` vs `NA`) and R2b (`"1e+05"` vs `"100000"`) are all two derivations that disagree +on an input the author did not picture. `accuracy_key()` gives a canonical form for one +representation, numeric. Factor, character and positional identities still go around +it, and so does one new kind of identity: **a point's id is derived from its draw order, +and the census branch derives draw order differently.** + +Every place in the four R files where an identity is derived, compared or aligned: + +| # | Site | Sound? | +|---|------|--------| +| 1 | labels.R:86 `accuracy_is_missing()` on the required columns | sound | +| 2 | labels.R:99 `point_id` uniqueness through `accuracy_key` | sound | +| 3 | labels.R:107 `use` compared raw against `"accuracy"`/`"training"`; blank `""` is not NA | **fragile** (F4, loud) | +| 4 | labels.R:117-130 label strata vs `strata$stratum`, both keyed; `known[strata$n_cells == 0]` is positional within one frame | sound | +| 5 | labels.R:151-152 `accuracy_check_strata` uniqueness via key | sound | +| 6 | labels.R:167-174 `accuracy_key()` factor branch returns levels verbatim, and `factor(100000)` has level `"1e+05"` | **bug** (F2) | +| 7 | estimate.R:147-149 `n_h` via `factor(lab_key, levels = key)` | sound | +| 8 | estimate.R:152/160 `key[...]` positional in the same filtered frame | sound | +| 9 | estimate.R:169-171 class set from the keys of both columns | sound, except input from #6 | +| 10 | estimate.R:172/195 `classes_numeric` -> `as.numeric(classes)` (keys are full-digit, so the round trip is exact) | sound | +| 11 | estimate.R:178-180 `sid[match(lab_key, key)]`; stehman2014 anchors its stratum regex (`^s1$`, checked in mapaccuracy 0.1.2 source) | sound | +| 12 | estimate.R:197-221 `est$matrix`/`UA`/`PA`/`area` indexed by `classes`; stehman2014 names these by `order` = `classes` | sound | +| 13 | estimate.R:239-243 `accuracy_class_order` (numeric or radix, locale-free) | sound | +| 14 | estimate.R:251-272 stratum table, all keyed; `strata$stratum[h]` positional in the same frame as `key` | sound | +| 15 | sample.R:115-128 pass-1 codes via `setdiff`/`match` on exact whole doubles | sound | +| 16 | sample.R:132-136 codes sorted once and shared by allocation, draw, pass 2 and labels | sound | +| 17 | sample.R:153/162 pass-2 `match(vv, code)`, `match(code, vv[ord])`; `order(method = "radix")` is stable, so the k-th in sorted order is the k-th in cell order. Verified exact against `which(v == h)[rank]` at rows 1/3/600 | sound | +| 18 | sample.R:176 pass-1 vs pass-2 totals guard | sound | +| 19 | sample.R:184 `match(ranks, found$rank)` (integer vs double whole numbers) | sound | +| 20 | sample.R:186 `point_id` = key of code + **draw position** | **bug** (F1) | +| 21 | sample.R:256-284 allocation names normalised through key; `n[key]` | sound | +| 22 | sample.R:311 census returns `seq_len(N)`, which is cell order, not the stream's draw order | **bug** (F1) | +| 23 | sample.R:323-326 stream seed from keys of seed and code | sound | +| 24 | sample.R:331 labels via `match(code, lv[[1]])`; `levels()` returns ID + active category | sound | +| 25 | sample.R:195-198 map values by cell after `compareGeom`; `raw = TRUE` gives codes (factor test passes) | sound | +| 26 | sample.R:220 `design$n_requested` names written by `setNames()` (`"1e+05"`); re-fed as `n`, they are normalised by #21 | sound | +| 27 | size.R:294-307 `weights` vs `s_h`/`ua` aligned **by position**; names on `s_h`/`ua` are ignored, then overwritten at :323 | **fragile** (F3) | +| 28 | size.R:315-322 `names(alloc) <- names(weights)` -> sampler's #21 | sound | +| 29 | estimate.R:245-247 -> size.R:297 census stratum sd is `NA` at n = 1 | **fragile** (F5, loud) | + +Also checked and sound: + +- **RNG save/restore:** + - existing seed: kind and seed restored, including `sample.kind = "Rounding"` and L'Ecuyer; + - absent seed with a non-default kind: kind restored and `.Random.seed` left absent. + + All probed. +- **Chunk-size independence:** the ranks are drawn before either pass, and the resolver is + deterministic. +- **terra-version stability:** only `readValues` row-major order is relied on. +- **Pilot extension for non-census strata:** `sample.int(useHash = FALSE)` has the partial- + shuffle prefix property under Rejection sampling. +- **Size allocation arithmetic:** the `proportional_min` loop terminates, including when + every stratum reaches the floor. + +## Findings + +- **[severity: bug]** R/dft_accuracy_sample.R:311 (with :184-187). **A census stratum + breaks the documented pilot-extension and `point_id`-stability contract.** A census + returns `seq_len(N_h)`, which is cell order. A non-census draw returns the stream's + order. A stratum with `n_pilot < N_h <= n_full` is drawn from the stream in the pilot + and taken as a census in the full sample, so its `point_id`s are re-assigned. The + roxygen promises that "labels from a pilot carry into the full sample (`point_id` is + stable too)". A caller who joins pilot labels by `point_id` therefore puts reference + labels on the wrong cells, and the error is silent. + + The case is common: transition strata of 20-50 cells, with a pilot of 20 and a full + sample of 50 from `dft_accuracy_size()`. Probe: 40-cell stratum, seed 81. + - n = 30 gives cells `14 2 26 35 16 ...`; n = 50 gives `1 2 3 4 5 ...`. + - **28 of 30 pilot `point_id`s point at a different cell in the full sample.** + + The test "a larger n extends the pilot" cannot reach this: no stratum on r17 has + 30 < N <= 50. Its counts are 941/7127/2/998/55/2/3186. + + Fix: drop the early return, and draw `sample.int(N_h, min(n, N_h), useHash = FALSE)` + for every stratum. With size = N it is a full permutation, and it extends any smaller + draw (verified: `sample.int(40, 40)[1:30] == sample.int(40, 30)` under the same seed). + The golden test pins stratum 2 and stratum 1's first ids, which are not census strata, + so it is unaffected. The package is unreleased, so no stored sample moves. + +- **[severity: bug]** R/dft_accuracy_labels.R:168. **The factor branch of `accuracy_key()` + lets the R2b split back in through a factor column.** `factor()` builds its levels with + `as.character()`, so `factor(100000)` has the level `"1e+05"`, while a numeric + `ref_class` of 100000 keys to `"100000"`. Probe: `map_class = + factor(rep(c(100000, 2), each = 4))` with a numeric `ref_class`. + - The result has classes `2`, `1e+05` and `100000`. + - User's accuracy for `1e+05` is 0, and for `100000` it is `NA`. + - **Overall accuracy is 0.375; the truth is 0.75.** There is no error or warning. + + For `stratum` the same split is refused loudly (unknown stratum). For `map_class` / + `ref_class` it is silent. Round codes of 1e5 and above arise only from custom schemes, + such as a transition to code 0, so this is a narrow input. It is the same mechanism, + though, on the branch R2b's fix did not normalise. Fix: key a factor's levels through + the same numeric normalisation (`accuracy_key` of `as.numeric(levels)` where every level + parses as a number). The alternative is to normalise any character or factor value that + parses as a whole number. + +- **[severity: fragile]** R/dft_accuracy_size.R:294-307, 323. `s_h` and `ua` are aligned + to `weights` by position. When the caller passes named `s_h`/`ua` in a different order + from `weights`, the names are ignored, n from Eq. 13 is silently wrong, and line 323 + then overwrites the names so that the returned `s_h` looks aligned. Where both carry + names, match them, or refuse a mismatch. + +- **[severity: fragile]** R/dft_accuracy_labels.R:107. The contract permits `use` to be + `NA`, but `read.csv()` of a partly filled `use` column gives `""`. The table is then + refused with "got: ." (probed). This is R2's blank-vs-NA mechanism reaching the one + column `accuracy_is_missing()` was not applied to. The failure is loud, not silent, + but a label table that meets the contract is rejected. Fix: treat + `accuracy_is_missing(labels$use)` as NA in both the check and estimate.R:128. + +- **[severity: fragile]** R/dft_accuracy_estimate.R:245-247 -> R/dft_accuracy_size.R:297. + A census stratum of one cell is legitimate: the sampler emits it, and the estimator + accepts it. It gets `sd = NA` in `$stratum`. Fed to `dft_accuracy_size()` as the + documented pilot path does, it is refused with "label another", which cannot be done + because the stratum has one cell (probed: 1-cell stratum, `sd` NA, sizer errors). Its + true S_h contribution is 0. Fix: set `sd` to 0 where `n == n_cells`. + +- **[severity: fragile]** tests/testthat/test-dft_accuracy_sample.R:60-66. The block is + labelled "brute force: the k-th cell of a stratum in cell order is which(v == h)[k]". + It asserts only `p$cell %in% cells_h`, so a resolver that mapped a rank to the wrong + cell *within* the right stratum would pass. The resolver is in fact correct: checked + exactly against `which(v == code)[ranks]` at rows 1/3/600. The test simply cannot fail + for the defect it names. + +## Disposition (parent session, 2026-09-29) + +All five non-sound rows are fixed, each with a test, and the mutation check turned every test red with its fix undone (scratch copy): +- #3: a blank `use` is NA. Test: "a blank `use` cell is NA". +- #6: `accuracy_key()` reads numeric-looking strings and factor levels back as numbers. Mutation: 2 failures. +- #20/#22: a census draws from the stream. Test: "a stratum tipped into a census … keeps its pilot ids". Mutation: 1 failure. +- #27: `s_h`/`ua` are matched to `weights` by keyed name. Test: "named s_h is matched to weights by name". +- #29: a census stratum reports sd 0. Mutation: 2 failures. + +Two sound rows were also routed through `accuracy_key()`, so every site now uses one key: #26 (`$design$n_requested` names) and the new name match in the sizer. + +#2 changed: `point_id` is now compared raw, not keyed, so "01" and "1" stay distinct ids. + +**Termination:** the round-3 enumeration covers every identity derivation in the four files (29 sites), and none now sits outside `accuracy_key()` or a same-frame positional index. That ends the loop. diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index 3da7cf9..ba28928 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -65,18 +65,18 @@ The package always applies the FPC `(1 − n_h/N_h)`, so a census stratum contri - [x] `findings.md`: the estimator equations with numbers, and a check of the Olofsson example by hand arithmetic (deforestation 21,158 ha is reproducible from the row counts; confirm against the PDF) ### Phase 2: Estimator (tests first) -- [ ] `mapaccuracy` in Imports (DESCRIPTION) -- [ ] `test-dft_accuracy_estimate.R`, published: Olofsson 2014 Tables 8–9 through `dft_accuracy_estimate()`, with **absolute** tolerance at the published precision. The forest-gain and stable-non-forest PA CIs pin the Eq. 7 values (±0.254, ±0.018), not the printed ±0.23 / ±0.01, and cite the discrepancy (findings). Add the Stehman 2014 example once its PDF is in hand. The fixture expands counts to per-point rows with base `rep()` -- [ ] Must-fail: run the **estimator** with equal `n_cells` on the Olofsson counts and assert it differs from the published values -- [ ] Census oracle (independent truth, no PDF needed): map = 2017 and "reference" = 2023 on the bundled tile, with the true error matrix and areas from `terra::crosstab`. Run about 500 stratified draws under map-class strata and under a changed/stable split; check bias ≈ 0, empirical SD ≈ mean SE, and CI coverage ≈ level. Skippable if slow -- [ ] Perfect labels (`ref = map`) give UA = PA = OA = 1, SE = 0, and adjusted area equal to mapped area. Recoding to a 2-class union gives an SE that is not the sum of the SEs -- [ ] Contract refusals: training rows, NA `ref_class`, duplicate ids, an unknown stratum, a stratum with no labels, `n_h = 1`. A reference-only class appears in the matrix, and PA is NA at `p̂_·j = 0` -- [ ] Measure `stehman2014()` runtime at #93 scale (about 1,000 points × about 80 transition classes; it builds `classes²` indicator columns). File upstream and report if it is impractical -- [ ] `R/dft_accuracy_estimate.R` + `R/dft_accuracy_labels.R`. Freeze the `$strata` and `$stratum` shapes here -- [ ] Restore-the-bug check: pass unweighted `N_h` inside the wrapper and confirm the published-value tests go red +- [x] `mapaccuracy` in Imports (DESCRIPTION) +- [x] `test-dft_accuracy_estimate.R`, published: Olofsson 2014 Tables 8–9 through `dft_accuracy_estimate()`, with **absolute** tolerance at the published precision. The forest-gain and stable-non-forest PA CIs pin the Eq. 7 values (±0.254, ±0.018), not the printed ±0.23 / ±0.01, and cite the discrepancy (findings). Add the Stehman 2014 example once its PDF is in hand. The fixture expands counts to per-point rows with base `rep()` +- [x] Must-fail: run the **estimator** with equal `n_cells` on the Olofsson counts and assert it differs from the published values +- [x] Census oracle (independent truth, no PDF needed): map = 2017 and "reference" = 2023 on the bundled tile, with the true error matrix and areas from `terra::crosstab`. Run about 500 stratified draws under map-class strata and under a changed/stable split; check bias ≈ 0, empirical SD ≈ mean SE, and CI coverage ≈ level. Skippable if slow +- [x] Perfect labels (`ref = map`) give UA = PA = OA = 1, SE = 0, and adjusted area equal to mapped area. Recoding to a 2-class union gives an SE that is not the sum of the SEs +- [x] Contract refusals: training rows, NA `ref_class`, duplicate ids, an unknown stratum, a stratum with no labels, `n_h = 1`. A reference-only class appears in the matrix, and PA is NA at `p̂_·j = 0` +- [x] Measure `stehman2014()` runtime at #93 scale (about 1,000 points × about 80 transition classes; it builds `classes²` indicator columns). File upstream and report if it is impractical +- [x] `R/dft_accuracy_estimate.R` + `R/dft_accuracy_labels.R`. Freeze the `$strata` and `$stratum` shapes here +- [x] Restore-the-bug check: pass unweighted `N_h` inside the wrapper and confirm the published-value tests go red ### Phase 3: Sampler (tests first) -- [ ] `test-dft_accuracy_sample.R`: +- [x] `test-dft_accuracy_sample.R`: - golden `point_id`s and cells for a seeded draw on `example_2017.tif`, with an allocation the tile can satisfy (classes 4 and 9 have 2 cells) - chunk invariance: identical results at row chunks 1, 7, 50 and full, matching brute-force `which()` - pilot extension: the first 30 per stratum at n = 30 are identical at n = 50, and adding a stratum leaves the others unchanged @@ -88,13 +88,13 @@ The package always applies the FPC `(1 − n_h/N_h)`, so a census stratum contri - refusals: lonlat, non-integer, multi-layer, and a named allocation that omits a stratum - NA cells are never drawn - `map =` extraction, including a grid-mismatch refusal -- [ ] `R/dft_accuracy_sample.R` -- [ ] Sampler → estimator integration: draw, fake labels from a reference raster, estimate (covered by the census oracle once both exist) -- [ ] Scale test on BULK: `classified_2017.tif` and its `dft_rast_transition()` factor output (the #93 input), with an RSS sampler. Record pass-1 and pass-2 time and peak RSS in the PR body +- [x] `R/dft_accuracy_sample.R` +- [x] Sampler → estimator integration: draw, fake labels from a reference raster, estimate (covered by the census oracle once both exist) +- [x] Scale test on BULK: `classified_2017.tif` and its `dft_rast_transition()` factor output (the #93 input), with an RSS sampler. Record pass-1 and pass-2 time and peak RSS in the PR body ### Phase 4: Sizing -- [ ] `test-dft_accuracy_size.R`: reproduce Olofsson §5.1.1's sample-size example (n and allocation); the `s_h` form from a pilot `$stratum` agrees with the `ua` form when strata = map classes; edge cases (UA = 1, a zero weight) -- [ ] `R/dft_accuracy_size.R` +- [x] `test-dft_accuracy_size.R`: reproduce Olofsson §5.1.1's sample-size example (n and allocation); the `s_h` form from a pilot `$stratum` agrees with the `ua` form when strata = map classes; edge cases (UA = 1, a zero weight) +- [x] `R/dft_accuracy_size.R` ### Phase 5: Docs and release - [ ] A runnable `@examples` block on every function (the estimator example uses Olofsson's published counts; the sampler uses the bundled tile with a satisfiable allocation) @@ -108,6 +108,6 @@ No vignette. One made with fabricated labels would illustrate a number nobody me ## Validation - [ ] Tests pass -- [ ] `/code-check` clean on each commit +- [x] `/code-check` clean on each commit - [ ] PWF checkboxes match landed work - [ ] `/planning-archive` on completion diff --git a/tests/testthat/helper-accuracy.R b/tests/testthat/helper-accuracy.R index b44b672..a9dceef 100644 --- a/tests/testthat/helper-accuracy.R +++ b/tests/testthat/helper-accuracy.R @@ -14,10 +14,10 @@ olofsson_classes <- c("deforestation", "forest_gain", "stable_forest", # Table 8 (p. 55): sample counts n_ij, rows = map (= stratum), cols = reference olofsson_counts <- matrix( - c(66, 0, 5, 4, - 0, 55, 8, 12, - 1, 0, 153, 11, - 2, 1, 9, 313), + c(66, 0, 5, 4, + 0, 55, 8, 12, + 1, 0, 153, 11, + 2, 1, 9, 313), nrow = 4, byrow = TRUE, dimnames = list(olofsson_classes, olofsson_classes) ) @@ -92,3 +92,18 @@ olofsson_strata <- function(n_cells = olofsson_pixels) { weight = unname(n_cells) / sum(n_cells) ) } + +# testthat's `tolerance` is RELATIVE (and 3e ignores `scale`): `tolerance = 1` +# on a 21,158 ha area accepts anything within 100%, while the default demands +# agreement a published figure rounded to the ha can never give. The +# published values need an absolute bound at their printed precision. +expect_within <- function(object, expected, tol) { + act <- testthat::quasi_label(rlang::enquo(object), arg = "object") + d <- abs(as.numeric(act$val) - as.numeric(expected)) + testthat::expect( + length(act$val) == length(expected) && !anyNA(d) && all(d <= tol), + sprintf("%s differs from the expected values by up to %g (> %g).", + act$lab, suppressWarnings(max(d, na.rm = TRUE)), tol) + ) + invisible(act$val) +} diff --git a/tests/testthat/test-dft_accuracy_estimate.R b/tests/testthat/test-dft_accuracy_estimate.R new file mode 100644 index 0000000..185d528 --- /dev/null +++ b/tests/testthat/test-dft_accuracy_estimate.R @@ -0,0 +1,278 @@ +# Published values are in helper-accuracy.R, each cited to Olofsson et al. +# (2014)'s page and table. Tolerances are ABSOLUTE and set at the published +# precision: 0.005 for proportions printed to 2 dp, 1 ha for areas printed to +# the ha, and 1.5 ha for area half-widths, which Stehman's finite population +# correction and exact z (vs the paper's 1.96, no FPC) move by up to 1.1 ha +# (findings.md). + +ol <- dft_accuracy_estimate(olofsson_labels(), olofsson_strata()) +by_class <- function(tbl, measure) { + x <- tbl[tbl$measure == measure, ] + x[match(olofsson_classes, x$class), ] +} + +test_that("error-adjusted areas reproduce Olofsson et al. (2014) section 5.2.2", { + a <- ol$area[match(olofsson_classes, ol$area$class), ] + expect_within(a$area, olofsson_area, 1) + expect_within(a$area - a$lower, olofsson_area_hw, 1.5) + expect_within(a$upper - a$area, olofsson_area_hw, 1.5) + expect_equal(ol$area_total, 900000) +}) + +test_that("the error matrix reproduces Olofsson Table 9 to 4 dp", { + got <- matrix(NA_real_, 4, 4, dimnames = dimnames(olofsson_p)) + got[cbind(ol$matrix$map_class, ol$matrix$ref_class)] <- ol$matrix$proportion + expect_within(got, olofsson_p, 5e-5) + expect_equal(sum(ol$matrix$proportion), 1) + # zeros are zeros, not NA (mapaccuracy returns empty cells as NA) + expect_false(anyNA(ol$matrix$proportion)) +}) + +test_that("accuracies reproduce Olofsson section 5.2.1 to 2 dp", { + u <- by_class(ol$accuracy, "user") + p <- by_class(ol$accuracy, "producer") + o <- ol$accuracy[ol$accuracy$measure == "overall", ] + expect_within(u$estimate, olofsson_user, 0.005) + expect_within(u$upper - u$estimate, olofsson_user_hw, 0.005) + expect_within(p$estimate, olofsson_prod, 0.005) + # two of these are Eq. (7) values, not the printed ones (helper-accuracy.R) + expect_within(p$upper - p$estimate, olofsson_prod_hw, 0.005) + expect_within(o$estimate, olofsson_overall, 0.005) + expect_within(o$upper - o$estimate, olofsson_overall_hw, 0.005) + expect_true(is.na(o$class)) +}) + +test_that("the stratum weights are load-bearing: equal weights give the wrong answer", { + # the must-fail: the same labels through the estimator with every stratum the + # same size, which is what an unweighted confusion matrix assumes + flat <- dft_accuracy_estimate(olofsson_labels(), + olofsson_strata(rep(2.5e6, 4))) + a <- flat$area[match(olofsson_classes, flat$area$class), ] + expect_gt(abs(a$area[1] - olofsson_area[1]), 100000) # vs 21,158 ha + p <- by_class(flat$accuracy, "producer") + expect_gt(max(abs(p$estimate - olofsson_prod)), 0.05) +}) + +test_that("perfect labels give accuracy 1, SE 0 and area equal to mapped area", { + lab <- olofsson_labels() + lab$ref_class <- lab$map_class + res <- dft_accuracy_estimate(lab, olofsson_strata()) + expect_equal(res$accuracy$estimate, rep(1, 9)) + expect_equal(res$accuracy$se, rep(0, 9)) + a <- res$area[match(olofsson_classes, res$area$class), ] + expect_equal(a$area, olofsson_pixels * 0.09, ignore_attr = TRUE) + expect_equal(a$area_se, rep(0, 4)) +}) + +test_that("a union's SE comes from recoding, and is not the sum of its members'", { + lab <- olofsson_labels() + change <- c("deforestation", "forest_gain") + lab_u <- lab + lab_u$map_class <- ifelse(lab$map_class %in% change, "change", "stable") + lab_u$ref_class <- ifelse(lab$ref_class %in% change, "change", "stable") + u <- dft_accuracy_estimate(lab_u, olofsson_strata()) + a <- ol$area[match(change, ol$area$class), ] + ch <- u$area[u$area$class == "change", ] + # the area is additive ... + expect_equal(ch$area, sum(a$area)) + # ... the standard error is not + expect_false(isTRUE(all.equal(ch$area_se, sum(a$area_se)))) + expect_lt(ch$area_se, sum(a$area_se)) +}) + +test_that("the stratum table carries per-stratum means and SDs for sizing", { + st <- ol$stratum + expect_setequal(unique(st$target), c("agreement", olofsson_classes)) + ag <- st[st$target == "agreement", ] + ag <- ag[match(olofsson_classes, ag$stratum), ] + expect_equal(ag$n, c(75, 75, 165, 325)) + expect_equal(ag$mean, unname(diag(olofsson_counts)) / c(75, 75, 165, 325)) + p <- ag$mean + expect_equal(ag$sd, sqrt(p * (1 - p) * c(75, 75, 165, 325) / + (c(75, 75, 165, 325) - 1))) + expect_equal(sum(ag$weight), 1) +}) + +test_that("numeric classes are ordered numerically and returned numeric", { + lab <- data.frame(point_id = 1:8, stratum = rep(c(2, 10), each = 4), + map_class = rep(c(2, 10), each = 4), + ref_class = c(2, 2, 2, 10, 10, 10, 10, 2)) + st <- data.frame(stratum = c(2, 10), n_cells = c(100, 300), area = c(1, 3)) + res <- dft_accuracy_estimate(lab, st) + expect_identical(res$area$class, c(2, 10)) + expect_type(res$matrix$map_class, "double") +}) + +test_that("a reference-only class enters the matrix and producer's accuracy is NA off the reference", { + lab <- olofsson_labels() + lab$ref_class[1] <- "water" # a class the map never uses + res <- dft_accuracy_estimate(lab, olofsson_strata()) + expect_true("water" %in% res$area$class) + expect_true("water" %in% res$matrix$ref_class) + # no point is labelled forest_gain in the reference -> producer's NA + lab2 <- olofsson_labels() + lab2$ref_class[lab2$ref_class == "forest_gain"] <- "stable_forest" + res2 <- dft_accuracy_estimate(lab2, olofsson_strata()) + pa <- res2$accuracy[res2$accuracy$measure == "producer" & + res2$accuracy$class == "forest_gain", ] + expect_true(is.na(pa$estimate)) +}) + +test_that("strata that are not the map classes are estimated, with estimated row totals", { + # two strata that cut across all four map classes + lab <- olofsson_labels() + lab$stratum <- ifelse(seq_len(nrow(lab)) %% 2 == 0, "a", "b") + st <- data.frame(stratum = c("a", "b"), n_cells = c(4e6, 6e6), + area = c(4e6, 6e6) * 0.09) + res <- dft_accuracy_estimate(lab, st) + rows <- tapply(res$matrix$proportion, res$matrix$map_class, sum) + expect_equal(sum(rows), 1) + # the map-class shares are estimated from the sample, not the Olofsson W_i + shares <- as.vector(rows[olofsson_classes]) + expect_gt(max(abs(shares - olofsson_pixels / 1e7)), 0.05) + expect_true(all(is.finite(res$area$area_se))) +}) + +test_that("a census stratum contributes no variance and is not refused", { + lab <- data.frame(point_id = 1:7, + stratum = c(1, 1, 1, 1, 1, 2, 2), + map_class = c(1, 1, 1, 1, 1, 2, 2), + ref_class = c(1, 1, 2, 1, 1, 2, 1)) + # stratum 2 has exactly two cells, both labelled + st <- data.frame(stratum = c(1, 2), n_cells = c(1000, 2), area = c(10, 0.02)) + res <- dft_accuracy_estimate(lab, st) + s2 <- res$stratum[res$stratum$stratum == 2 & res$stratum$target == "agreement", ] + expect_equal(s2$n, 2) + expect_true(all(is.finite(res$area$area_se))) + # a one-point census is allowed; the warning mapaccuracy raises is muffled + lab1 <- lab[-7, ] + st1 <- data.frame(stratum = c(1, 2), n_cells = c(1000, 1), area = c(10, 0.01)) + expect_no_warning(dft_accuracy_estimate(lab1, st1)) +}) + +test_that("the design refusals name what is wrong", { + lab <- olofsson_labels() + st <- olofsson_strata() + + lab_t <- lab + lab_t$use <- "accuracy" + lab_t$use[1:3] <- "training" + expect_error(dft_accuracy_estimate(lab_t, st), "3 row\\(s\\) have `use == \"training\"`") + + lab_na <- lab + lab_na$ref_class[5] <- NA + expect_error(dft_accuracy_estimate(lab_na, st), "`ref_class` has 1 missing.*nonresponse") + + lab_dup <- lab + lab_dup$point_id[2] <- lab_dup$point_id[1] + expect_error(dft_accuracy_estimate(lab_dup, st), "must be unique; duplicated: p0001") + + lab_unk <- lab + lab_unk$stratum[1] <- "mystery" + expect_error(dft_accuracy_estimate(lab_unk, st), "absent from `strata`: mystery") + + st_extra <- rbind(st, tibble::tibble(stratum = "water", n_cells = 10, + area = 0.9, weight = 0)) + expect_error(dft_accuracy_estimate(lab, st_extra), "cells but no labels: water") + + first_gain <- match("forest_gain", lab$stratum) + lab_one <- lab[lab$stratum != "forest_gain" | seq_len(nrow(lab)) == first_gain, ] + expect_error(dft_accuracy_estimate(lab_one, st), "single labelled point: forest_gain") + + expect_error(dft_accuracy_estimate(lab, st[, c("stratum", "n_cells")]), + "needs an `area` column") + st_grid <- st + st_grid$area[1] <- st_grid$area[1] * 2 + expect_error(dft_accuracy_estimate(lab, st_grid), "describe different grids") + expect_error(dft_accuracy_estimate(lab, st, level = 95), "between 0 and 1") + + lab_use <- lab + lab_use$use <- "test" + expect_error(dft_accuracy_estimate(lab_use, st), "got: test") +}) + +test_that("a factor in one class column adds no phantom classes", { + # c() of a factor and a numeric falls back to the factor's integer codes + lab <- data.frame(point_id = 1:6, stratum = c(1, 1, 1, 2, 2, 2), + map_class = c(1001, 1001, 2002, 2002, 2002, 1001), + ref_class = factor(c(1001, 2002, 2002, 2002, 1001, 1001))) + st <- data.frame(stratum = 1:2, n_cells = c(100, 100), area = c(1, 1)) + res <- dft_accuracy_estimate(lab, st) + expect_identical(res$area$class, c("1001", "2002")) + lab$ref_class <- as.character(lab$ref_class) + lab$map_class <- factor(lab$map_class) + expect_identical(dft_accuracy_estimate(lab, st)$area$class, c("1001", "2002")) +}) + +test_that("a blank reference label is nonresponse, refused like NA", { + # read.csv() and fread() read an empty cell as "", not NA + lab <- olofsson_labels() + lab$ref_class[3] <- "" + expect_error(dft_accuracy_estimate(lab, olofsson_strata()), + "`ref_class` has 1 missing or blank") + lab$ref_class[3] <- " " + expect_error(dft_accuracy_estimate(lab, olofsson_strata()), "missing or blank") +}) + +test_that("a double and an integer column holding the same code are one class", { + # as.character(100000) is "1e+05"; as.character(100000L) is "100000" + lab <- data.frame(point_id = 1:6, stratum = c(1L, 1L, 1L, 2L, 2L, 2L), + map_class = c(100000, 100000, 2002, 2002, 2002, 100000), + ref_class = c(100000L, 100000L, 2002L, 2002L, 2002L, 100000L)) + st <- data.frame(stratum = c(1, 2), n_cells = c(100, 100), area = c(1, 1)) + res <- dft_accuracy_estimate(lab, st) + expect_identical(res$area$class, c(2002, 100000)) + expect_equal(res$accuracy$estimate[res$accuracy$measure == "overall"], 1) + expect_true("100000" %in% res$stratum$target) +}) + +test_that("a factor of numeric codes meets a numeric column as one class", { + # levels(factor(100000)) is "1e+05" + lab <- data.frame(point_id = 1:8, stratum = rep(1:2, each = 4), + map_class = factor(rep(c(100000, 2002), each = 4)), + ref_class = c(100000, 100000, 100000, 2002, + 2002, 2002, 2002, 100000)) + st <- data.frame(stratum = 1:2, n_cells = c(100, 100), area = c(1, 1)) + res <- dft_accuracy_estimate(lab, st) + expect_identical(res$area$class, c("2002", "100000")) + expect_equal(res$accuracy$estimate[res$accuracy$measure == "overall"], 0.75) +}) + +test_that("a blank `use` cell is NA, as the contract allows", { + lab <- olofsson_labels() + lab$use <- "" + lab$use[1] <- "accuracy" + expect_no_error(dft_accuracy_estimate(lab, olofsson_strata())) +}) + +test_that("a census stratum reports sd 0, so the sizer can use it", { + lab <- data.frame(point_id = 1:6, stratum = c(1, 1, 1, 1, 1, 2), + map_class = c(1, 1, 1, 1, 1, 2), + ref_class = c(1, 1, 2, 1, 1, 2)) + st <- data.frame(stratum = c(1, 2), n_cells = c(1000, 1), area = c(10, 0.01)) + res <- dft_accuracy_estimate(lab, st) + s2 <- res$stratum[res$stratum$stratum == 2, ] + expect_true(all(s2$sd == 0)) + ag <- res$stratum[res$stratum$target == "agreement", ] + expect_no_error(dft_accuracy_size(stats::setNames(ag$weight, ag$stratum), + se_target = 0.05, s_h = ag$sd, n_min = 2)) +}) + +test_that("labels in a stratum with no cells are refused by name", { + lab <- olofsson_labels() + st <- olofsson_strata() + st$n_cells[st$stratum == "forest_gain"] <- 0 + st$area[st$stratum == "forest_gain"] <- 0 + expect_error(dft_accuracy_estimate(lab, st), "no cells in `strata`: forest_gain") +}) + +test_that("expect_within() refuses a short vector that recycles", { + expect_failure(expect_within(c(1, 2), c(1, 2, 1, 2), 0.1)) + expect_success(expect_within(c(1, 2), c(1.05, 2), 0.1)) +}) + +test_that("level sets the interval width through qnorm", { + r90 <- dft_accuracy_estimate(olofsson_labels(), olofsson_strata(), level = 0.90) + hw <- r90$area$upper - r90$area$area + expect_equal(hw, stats::qnorm(0.95) * r90$area$area_se) +}) diff --git a/tests/testthat/test-dft_accuracy_sample.R b/tests/testthat/test-dft_accuracy_sample.R new file mode 100644 index 0000000..77b27ae --- /dev/null +++ b/tests/testthat/test-dft_accuracy_sample.R @@ -0,0 +1,294 @@ +tile <- function(yr = 2017) { + terra::rast(system.file("extdata", paste0("example_", yr, ".tif"), + package = "drift")) +} +r17 <- tile(2017) +# classes 4 and 9 have 2 cells each on this tile, so n = 5 censuses them +alloc <- 5 +draw <- function(...) suppressMessages(dft_accuracy_sample(...)) + +test_that("a seeded draw is pinned: same seed, same points, on any terra", { + s <- draw(r17, n = alloc, seed = 81) + # Golden values. If this fails, the draw moved: a change to the RNG + # handling, the stream seeds, cell order, or the resolver -- every stored + # sample.gpkg drawn before the change no longer redraws. Do not re-pin + # without saying so in NEWS. + expect_identical( + s$points$cell[s$points$stratum == 2], + c(33234, 46009, 26228, 28183, 27211) + ) + expect_identical(utils::head(s$points$point_id, 3), + c("1_00001", "1_00002", "1_00003")) + expect_identical(draw(r17, n = alloc, seed = 81)$points$cell, s$points$cell) + expect_false(identical(draw(r17, n = alloc, seed = 82)$points$cell, + s$points$cell)) +}) + +test_that("strata carry counts, areas and weights that match the raster", { + s <- draw(r17, n = alloc, seed = 1) + tab <- table(terra::values(r17)[, 1]) + expect_identical(s$strata$stratum, as.numeric(names(tab))) + expect_equal(s$strata$n_cells, as.vector(tab)) + expect_equal(sum(s$strata$weight), 1) + summ <- dft_rast_summarize(r17, source = "io-lulc", unit = "ha") + expect_equal(sum(s$strata$area), sum(summ$area)) + expect_equal(s$strata$n, pmin(alloc, as.vector(tab))) + expect_true(all(is.na(s$strata$stratum_label))) +}) + +test_that("points fall in their stratum, never on NA, one per cell", { + s <- draw(r17, n = 30, seed = 3) + v <- terra::values(r17)[s$points$cell, 1] + expect_false(anyNA(v)) + expect_equal(v, s$points$stratum) + expect_false(anyDuplicated(s$points$cell) > 0) + expect_false(anyDuplicated(s$points$point_id) > 0) + xy <- sf::st_coordinates(s$points) + expect_equal(terra::cellFromXY(r17, xy), s$points$cell) + expect_identical(sf::st_crs(s$points)$wkt, sf::st_crs(terra::crs(r17))$wkt) +}) + +test_that("the chunked resolver matches brute force at every chunk size", { + # the tile fits one default chunk, so force the block-boundary carry + ref <- draw(r17, n = 30, seed = 7) + for (rows in c(1L, 7L, 50L, 314L)) { + withr::local_options(drift.accuracy_rows_chunk = rows) + got <- draw(r17, n = 30, seed = 7) + expect_identical(got$points$cell, ref$points$cell, info = paste("rows", rows)) + expect_identical(got$strata, ref$strata, info = paste("rows", rows)) + } + # brute force: the k-th cell of a stratum in cell order is which(v == h)[k], + # so the drawn ranks, resolved by which(), must give exactly these cells + v <- terra::values(r17)[, 1] + ranks <- drift:::accuracy_draw(ref$strata$stratum, ref$strata$n_cells, + ref$strata$n, seed = 7) + for (i in seq_along(ranks)) { + h <- ref$strata$stratum[i] + expect_identical(ref$points$cell[ref$points$stratum == h], + as.numeric(which(v == h)[ranks[[i]]]), + info = paste("stratum", h)) + } +}) + +test_that("a larger n extends the pilot, and other strata do not move", { + s30 <- draw(r17, n = 30, seed = 11) + s50 <- draw(r17, n = 50, seed = 11) + for (h in s30$strata$stratum) { + a <- s30$points[s30$points$stratum == h, ] + b <- s50$points[s50$points$stratum == h, ] + k <- nrow(a) + expect_identical(b$cell[seq_len(k)], a$cell, info = paste("stratum", h)) + expect_identical(b$point_id[seq_len(k)], a$point_id, info = paste("stratum", h)) + } + # raising one stratum's n leaves the others' draws alone + n1 <- stats::setNames(rep(30, nrow(s30$strata)), s30$strata$stratum) + n1["11"] <- 60 + s_mix <- draw(r17, n = n1, seed = 11) + keep <- s30$points$stratum != 11 + expect_identical(s_mix$points$cell[s_mix$points$stratum != 11], + s30$points$cell[keep]) +}) + +test_that("a stratum tipped into a census by a larger n keeps its pilot ids", { + # a 40-cell stratum: drawn at n = 30, taken whole at n = 50 + r <- terra::rast(nrows = 10, ncols = 10, xmin = 0, xmax = 100, ymin = 0, + ymax = 100, crs = "EPSG:32609", + vals = rep(c(1, 2), times = c(40, 60))) + p30 <- draw(r, n = 30, seed = 81) + p50 <- draw(r, n = 50, seed = 81) + a <- p30$points[p30$points$stratum == 1, ] + b <- p50$points[p50$points$stratum == 1, ] + expect_equal(nrow(b), 40) + expect_identical(b$cell[match(a$point_id, b$point_id)], a$cell) +}) + +test_that("the caller's RNG state is restored, or left absent", { + set.seed(123) + before <- .Random.seed + draw(r17, n = alloc, seed = 1) + expect_identical(.Random.seed, before) + + withr::local_seed(1) # restored at the end of this test + rm(".Random.seed", envir = globalenv()) + draw(r17, n = alloc, seed = 1) + expect_false(exists(".Random.seed", envir = globalenv(), inherits = FALSE)) + + set.seed(5, kind = "L'Ecuyer-CMRG") + kind_before <- RNGkind() + seed_before <- .Random.seed + s1 <- draw(r17, n = alloc, seed = 1) + expect_identical(RNGkind(), kind_before) + expect_identical(.Random.seed, seed_before) + # and the caller's generator kind does not change the draw + suppressWarnings(RNGkind("Mersenne-Twister", "Inversion", "Rounding")) + expect_identical(suppressWarnings(draw(r17, n = alloc, seed = 1))$points$cell, + s1$points$cell) + RNGkind("default", "default", "default") +}) + +test_that("a stratum no larger than its allocation is taken whole, with a message", { + expect_message(s <- dft_accuracy_sample(r17, n = alloc, seed = 1), + "census.*4 \\(2\\), 9 \\(2\\)") + expect_equal(s$strata$n[s$strata$stratum %in% c(4, 9)], c(2, 2)) + expect_identical(s$design$census, c(4, 9)) +}) + +test_that("a factor transition raster keeps codes as stratum and labels alongside", { + # fresh tiles: dft_rast_classify() mutates its input in place (#89) + cl <- dft_rast_classify(list("2017" = tile(2017), "2023" = tile(2023)), + source = "io-lulc") + tr <- dft_rast_transition(cl, from = "2017", to = "2023")$raster + s <- draw(tr, n = 4, seed = 2, map = tr) + lv <- terra::levels(tr)[[1]] + expect_true(all(s$strata$stratum %in% lv[[1]])) + expect_identical(s$strata$stratum_label, + as.character(lv[[2]][match(s$strata$stratum, lv[[1]])])) + expect_true("Water -> Water" %in% s$strata$stratum_label) + # map values are raw codes, so a transition map's class is its id + expect_identical(s$points$map_class, s$points$stratum) + expect_equal(sum(s$strata$n_cells), sum(!is.na(terra::values(tr)))) +}) + +test_that("map values are read at the sampled cells, as one column or a series", { + r23 <- tile(2023) + s <- draw(r17, n = 10, seed = 4, map = list(`2017` = r17, `2023` = r23)) + expect_equal(s$points$map_2017, terra::values(r17)[s$points$cell, 1]) + expect_equal(s$points$map_2023, terra::values(r23)[s$points$cell, 1]) + stack <- c(r17, r23) + # both tiles' layers are called "data", which would collide + expect_error(draw(r17, n = 10, seed = 4, map = stack), "names must be unique") + names(stack) <- c("y2017", "y2023") + s2 <- draw(r17, n = 10, seed = 4, map = stack) + expect_equal(s2$points$map_y2023, s$points$map_2023) + shifted <- terra::shift(r23, dx = 10) + expect_error(draw(r17, n = 10, seed = 4, map = shifted), "not on the `strata` grid") +}) + +test_that("draws within a stratum are uniform over its cells", { + # a 2-class raster: stratum 1 is 1000 cells; draw 100 each time and count + # how often each cell is picked over many seeds + r <- terra::rast(nrows = 40, ncols = 50, xmin = 0, xmax = 500, ymin = 0, + ymax = 400, crs = "EPSG:32609", vals = rep(1:2, each = 1000)) + hits <- integer(1000) + for (sd in 1:200) { + s <- draw(r, n = 100, seed = sd) + c1 <- s$points$cell[s$points$stratum == 1] + hits <- hits + tabulate(c1, 1000) + } + # expected 20 per cell; a chi-square on 999 df + p <- stats::pchisq(sum((hits - 20)^2 / 20), df = 999, lower.tail = FALSE) + expect_gt(p, 0.001) + # and the first and last cells are reachable + expect_gt(hits[1], 0) + expect_gt(hits[1000], 0) +}) + +test_that("an allocation named the way setNames() writes a number still matches", { + big <- terra::rast(nrows = 10, ncols = 10, xmin = 0, xmax = 100, ymin = 0, + ymax = 100, crs = "EPSG:32609", + vals = rep(c(2, 100000), each = 50)) + n <- stats::setNames(c(5, 5), c(2, 100000)) # names "2", "1e+05" + expect_identical(names(n)[2], "1e+05") + s <- draw(big, n = n, seed = 1) + expect_equal(s$strata$n, c(5, 5)) + expect_true(all(startsWith(s$points$point_id[s$points$stratum == 100000], + "100000_"))) +}) + +test_that("the design record carries what a redraw needs", { + s <- draw(r17, n = alloc, seed = 81) + d <- s$design + expect_identical(d$seed, 81) + expect_identical(unname(d$dims), c(terra::nrow(r17), terra::ncol(r17))) + expect_identical(d$terra_version, as.character(utils::packageVersion("terra"))) + expect_identical(unname(d$rng_kind), c("Mersenne-Twister", "Inversion", "Rejection")) +}) + +test_that("bad inputs are refused, naming the fault", { + expect_error(draw(r17, n = alloc), "`seed` must be a single whole number") + expect_error(draw(r17, n = alloc, seed = 1.5), "`seed` must be") + expect_error(draw(c(r17, r17), n = alloc, seed = 1), "one layer; it has 2") + ll <- terra::project(r17, "EPSG:4326", method = "near") + expect_error(draw(ll, n = alloc, seed = 1), "projected CRS") + fl <- r17 + 0.5 + expect_error(draw(fl, n = alloc, seed = 1), "integer stratum codes") + expect_error(draw(r17, n = 1, seed = 1), "at least 2 points") + expect_error(draw(r17, n = c(`1` = 5, `2` = 5), seed = 1), + "no allocation for stratum/strata present in `strata`: 4, 5, 7, 9, 11") + full <- stats::setNames(rep(5, 7), c(1, 2, 4, 5, 7, 9, 11)) + expect_error(draw(r17, n = c(full, `3` = 5), seed = 1), "no cells in `strata`: 3") + empty <- terra::rast(r17) + terra::values(empty) <- NA + expect_error(draw(empty, n = alloc, seed = 1), "no non-NA cells") + expect_error(draw(r17, n = alloc, seed = 1, map = "x"), "SpatRaster or a named list") +}) + +test_that("the sample feeds the estimator end to end", { + # reference = 2023 read at the points, map = 2017 + r23 <- tile(2023) + s <- draw(r17, n = 30, seed = 9, map = r17) + pts <- sf::st_drop_geometry(s$points) + pts$ref_class <- terra::values(r23)[pts$cell, 1] + pts <- pts[!is.na(pts$ref_class), ] + res <- dft_accuracy_estimate(pts, s$strata) + expect_equal(sum(res$area$area), sum(s$strata$area)) +}) + +# Census oracle: the tiles are small enough to know the truth. Take 2017 as +# the map and 2023 as "reference", draw many stratified samples, and check the +# estimator against the population -- unbiased, SEs that match the spread of +# the estimates, and intervals that cover at close to the nominal rate. The +# truth comes from every cell, not from the code under test. +census_oracle <- function(strata, reps = 300, n = 25) { + ref <- tile(2023) + sv <- terra::values(strata)[, 1] + rv <- terra::values(ref)[, 1] + pop <- !is.na(sv) + truth <- table(rv[pop]) / sum(pop) + big <- names(truth)[truth > 0.02] # Wald intervals need some mass + est <- se <- matrix(NA_real_, reps, length(big), dimnames = list(NULL, big)) + for (i in seq_len(reps)) { + s <- draw(strata, n = n, seed = i, map = r17) + pts <- sf::st_drop_geometry(s$points) + pts$ref_class <- rv[pts$cell] + a <- dft_accuracy_estimate(pts, s$strata)$area + k <- match(as.numeric(big), a$class) + est[i, ] <- ifelse(is.na(k), 0, a$proportion[k]) + se[i, ] <- ifelse(is.na(k), 0, a$proportion_se[k]) + } + z <- stats::qnorm(0.975) + tr <- as.vector(truth[big]) + list( + bias_z = (colMeans(est) - tr) / (apply(est, 2, stats::sd) / sqrt(reps)), + se_ratio = colMeans(se) / apply(est, 2, stats::sd), + coverage = colMeans(abs(est - rep(tr, each = reps)) <= z * se) + ) +} + +test_that("census oracle: strata are the map classes", { + skip_on_cran() + # 98 of the 7,127 map-Trees cells are reference Water (1.4%). At n = 25 a + # draw expects 0.34 of them, most draws see none, the stratum's variance is + # estimated as 0, and Water's interval covered 61% (SE ratio 0.76) while the + # estimate stayed unbiased. Coverage was 83% at n = 75 and 89% at n = 150: + # the Wald interval's small-sample weakness, documented on the estimator. + o <- census_oracle(r17, n = 150) + expect_true(all(abs(o$bias_z) < 4), info = paste(round(o$bias_z, 2), collapse = " ")) + expect_true(all(o$se_ratio > 0.8 & o$se_ratio < 1.25), + info = paste(round(o$se_ratio, 3), collapse = " ")) + expect_true(all(o$coverage > 0.85 & o$coverage < 0.99), + info = paste(round(o$coverage, 3), collapse = " ")) +}) + +test_that("census oracle: strata are not the map classes (changed / stable)", { + skip_on_cran() + # strata from a different pair of years than the map-vs-reference one, so + # they cut across the map classes + chg <- terra::ifel(r17 != tile(2019), 1L, 0L) + o <- census_oracle(chg, n = 60) + expect_true(all(abs(o$bias_z) < 4), info = paste(round(o$bias_z, 2), collapse = " ")) + expect_true(all(o$se_ratio > 0.8 & o$se_ratio < 1.25), + info = paste(round(o$se_ratio, 3), collapse = " ")) + expect_true(all(o$coverage > 0.88 & o$coverage < 0.99), + info = paste(round(o$coverage, 3), collapse = " ")) +}) diff --git a/tests/testthat/test-dft_accuracy_size.R b/tests/testthat/test-dft_accuracy_size.R new file mode 100644 index 0000000..a000b89 --- /dev/null +++ b/tests/testthat/test-dft_accuracy_size.R @@ -0,0 +1,85 @@ +w_ol <- stats::setNames(olofsson_plan_weights, olofsson_classes) + +test_that("Eq. 13 reproduces Olofsson section 5.1.1: n = 641", { + plan <- dft_accuracy_size(w_ol, se_target = 0.01, ua = olofsson_plan_ua) + expect_identical(plan$n, olofsson_plan_n) + # the S_i column of Table 5 + expect_equal(unname(plan$s_h), c(0.458, 0.490, 0.300, 0.218), tolerance = 1e-3) +}) + +test_that("equal and proportional allocations reproduce Table 5", { + eq <- dft_accuracy_size(w_ol, 0.01, ua = olofsson_plan_ua, allocation = "equal") + expect_equal(unname(eq$allocation), olofsson_plan_equal) + pr <- dft_accuracy_size(w_ol, 0.01, ua = olofsson_plan_ua, + allocation = "proportional") + expect_equal(unname(pr$allocation), olofsson_plan_prop) + expect_identical(names(pr$allocation), olofsson_classes) +}) + +test_that("proportional_min floors rare strata and splits the rest by weight", { + # Olofsson's Alloc1 rule: 100 per change stratum, remainder proportional + # among the stable ones. The paper prints 149/292 for the stable pair, which + # its own rule does not give (441 * 0.32/0.965 = 146.2); this asserts the rule + pm <- dft_accuracy_size(w_ol, 0.01, ua = olofsson_plan_ua, n_min = 100) + expect_equal(unname(pm$allocation), c(100, 100, 146, 295)) + # a stratum pushed under the floor by the reallocation joins the floor + w3 <- c(a = 0.01, b = 0.09, c = 0.90) + p3 <- dft_accuracy_size(w3, se_target = 0.05, s_h = c(0.5, 0.5, 0.5), + n_min = 20) + expect_identical(p3$n, 100) + expect_equal(unname(p3$allocation[c("a", "b")]), c(20, 20)) + expect_equal(unname(p3$allocation["c"]), 60) + # a floor bigger than the sample takes the floor everywhere + p4 <- dft_accuracy_size(w3, se_target = 0.2, s_h = c(0.5, 0.5, 0.5), n_min = 20) + expect_equal(unname(p4$allocation), c(20, 20, 20)) +}) + +test_that("the s_h form agrees with the ua form when strata are the map classes", { + s <- sqrt(olofsson_plan_ua * (1 - olofsson_plan_ua)) + a <- dft_accuracy_size(w_ol, 0.01, ua = olofsson_plan_ua) + b <- dft_accuracy_size(w_ol, 0.01, s_h = s) + expect_identical(a$n, b$n) + expect_identical(a$allocation, b$allocation) +}) + +test_that("named s_h is matched to weights by name, not position", { + s <- stats::setNames(sqrt(olofsson_plan_ua * (1 - olofsson_plan_ua)), + olofsson_classes) + a <- dft_accuracy_size(w_ol, 0.01, s_h = s) + b <- dft_accuracy_size(w_ol, 0.01, s_h = rev(s)) + expect_identical(a$n, b$n) + expect_identical(a$s_h, b$s_h) + bad <- stats::setNames(s, c("a", "b", "c", "d")) + expect_error(dft_accuracy_size(w_ol, 0.01, s_h = bad), "name different strata") +}) + +test_that("a pilot's stratum table feeds the sizer", { + est <- dft_accuracy_estimate(olofsson_labels(), olofsson_strata()) + ag <- est$stratum[est$stratum$target == "agreement", ] + plan <- dft_accuracy_size(stats::setNames(ag$weight, ag$stratum), + se_target = 0.01, s_h = ag$sd) + expect_gt(plan$n, 0) + expect_identical(names(plan$allocation), ag$stratum) +}) + +test_that("degenerate and malformed inputs are refused", { + expect_error(dft_accuracy_size(w_ol, 0.01, ua = rep(1, 4)), "SD of 0") + expect_error(dft_accuracy_size(c(a = 0, b = 1), 0.01, s_h = c(0.1, 0.1)), + "must be positive") + expect_error(dft_accuracy_size(c(a = 0.5, b = 0.6), 0.01, s_h = c(0.1, 0.1)), + "sum to 1") + expect_error(dft_accuracy_size(w_ol, 0.01), "exactly one of") + expect_error(dft_accuracy_size(w_ol, 0.01, s_h = rep(.1, 4), ua = rep(.9, 4)), + "exactly one of") + expect_error(dft_accuracy_size(w_ol, 0.01, ua = c(0.9, 1.2, 0.9, 0.9)), + "between 0 and 1") + expect_error(dft_accuracy_size(w_ol, 0.01, s_h = c(0.1, NA, 0.1, 0.1)), + "one pilot point") + expect_error(dft_accuracy_size(w_ol, 0.01, s_h = c(0.1, 0.1)), "one value per stratum") + expect_error(dft_accuracy_size(w_ol, 0, ua = olofsson_plan_ua), "positive") + expect_error(dft_accuracy_size(w_ol, 0.01, ua = olofsson_plan_ua, n_min = 1), + "at least 2") + # a zero SD in one stratum is allowed and contributes nothing to n + z <- dft_accuracy_size(w_ol, 0.01, s_h = c(0, 0.49, 0.3, 0.218)) + expect_lt(z$n, olofsson_plan_n) +}) From f9b3b0341eb3d1b89e6af8b627fed74d8452d198 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Tue, 29 Sep 2026 07:18:42 -0700 Subject: [PATCH 07/10] NEWS, CLAUDE.md pipeline block and tile size for #81 Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- CLAUDE.md | 17 ++++++++++++++++- NEWS.md | 13 +++++++++++++ planning/active/progress.md | 1 + planning/active/task_plan.md | 12 ++++++------ 4 files changed, 36 insertions(+), 7 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 985ea01..9a97f66 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -85,6 +85,21 @@ dft_map_interactive(classified, aoi = aoi, rgb = tc) # or rgb = chips <- lapply(seq_len(nrow(buf)), \(i) dft_stac_composite(buf[i, ], years = 2023)) ``` +Accuracy and error-adjusted area (#81; Olofsson et al. 2014, Stehman 2014). Mapped area is +biased wherever the map is wrong. These correct it from reference labels at stratified random points: + +```r +s <- dft_accuracy_sample(strata, n = 30, seed = 81, map = transition) # strata need not be map classes +# ... reviewer adds ref_class per point (drift stores no labels; contract: dft_accuracy_labels()) +est <- dft_accuracy_estimate(labels, s$strata) # wraps mapaccuracy::stehman2014(); $area in ha with CIs +plan <- dft_accuracy_size(weights, se_target = 0.01, s_h = pilot_sd) # size the full sample from a pilot +``` + +Same seed + larger `n` **extends** the draw (per-stratum streams), so pilot labels carry over. For +change maps `ref_class = ref_from * 1000 + ref_to`. Union targets (all tree loss) are recode-then-estimate: +the SE of a union is not the sum of SEs. Wald CIs run narrow for a rare class hiding in a big stratum +at small n (measured on the tile, see the estimator docs), so do not starve the stable strata. + ## Key Patterns - **Dual-mode maps:** `dft_map_interactive()` uses `addRasterImage()` for local SpatRasters, `addTiles()` via titiler for remote COGs @@ -118,7 +133,7 @@ devtools::install() # needed before rendering vignettes ### Scale-test raster functions on the BULK floodplain before the PR -The bundled Neexdzii Kwa tile is 600 x 600 cells and cannot reach memory or runtime failure modes. A real floodplain-scale pair is two downloads away — the BULK watershed group item on stac-floodplains-bc, already clipped to the floodplain, ~1.2 MB each: +The bundled Neexdzii Kwa tile is 314 x 326 cells and cannot reach memory or runtime failure modes. A real floodplain-scale pair is two downloads away — the BULK watershed group item on stac-floodplains-bc, already clipped to the floodplain, ~1.2 MB each: - `https://stac-floodplains-bc.s3.us-west-2.amazonaws.com/bulk_co_ff04/classified_2017.tif` - `https://stac-floodplains-bc.s3.us-west-2.amazonaws.com/bulk_co_ff04/classified_2023.tif` diff --git a/NEWS.md b/NEWS.md index 656cab3..1c58308 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,16 @@ +# drift 0.19.0 + +- **Accuracy and error-adjusted area (#81).** drift reported mapped area only. Mapped area is biased wherever the map is wrong, and badly biased for change maps, because a single wrong label on either date manufactures a transition. Four new functions implement the standard remedy (Olofsson et al. 2014): a stratified random sample of points, reference labels, and estimators that weight the labels back to the population. + - `dft_accuracy_sample()` draws the points from any strata raster: the map classes, a transition map, or strata such as "change attributed / unattributed / stable". It records each stratum's size (the weights) and, through `map =`, the map's class at every point. + - `dft_accuracy_estimate()` returns the error matrix in proportions of area, overall / user's / producer's accuracy, and error-adjusted area in hectares with confidence intervals. It wraps [mapaccuracy](https://cran.r-project.org/package=mapaccuracy)'s implementation of Stehman (2014), which holds whether or not the strata are the map classes. + - `dft_accuracy_labels()` is the contract a reference-label table must meet, so any review tool can feed the estimator. drift stores no labels. + - `dft_accuracy_size()` sizes and allocates a sample for a target standard error (Olofsson Eq. 13). It takes per-stratum SDs from a pilot, so it works for strata that are not the map classes. +- **It reproduces the paper, and the paper has typos.** Olofsson et al. 2014's worked example (Tables 8–9) comes back to its printed precision: deforestation 21,158 ± 6,157 ha against the mapped 18,000 ha. Three printed values contradict the paper's own equations: two producer's-accuracy half-widths (forest gain ±0.23, where Eq. 7 gives ±0.254; stable non-forest ±0.01, where it gives ±0.018) and the SE of deforestation area (34,097 px, where its own margin implies 34,907). The tests pin the equation values. A must-fail test runs the same labels with equal stratum weights, which is what an unweighted confusion matrix assumes: deforestation comes out at 200,748 ha instead of 21,158. +- **An independent check against a census.** The bundled tiles are small enough to know the truth, with 2017 as the map and 2023 as the reference. Over 300 draws the estimates are unbiased, with strata equal to and different from the map classes. The intervals cover close to 95%, except in one case the docs now warn about: a rare class hiding in a large stratum makes the interval too narrow at small samples. Water inside the Trees stratum covered 61% at 25 points per stratum, 83% at 75 and 89% at 150. +- **Draws are reproducible and extend from a pilot.** Each stratum draws from its own seeded stream with pinned RNG kinds, not `terra::spatSample()`, whose output terra does not promise to keep stable. The same seed with a larger `n` returns the pilot's points first, so pilot labels carry into the full sample. Changing one stratum's allocation leaves the others' points alone. A golden draw is pinned in the tests. A stratum no larger than its allocation is taken whole. +- **Floodplain scale.** The raster is read in chunks of about ten million cells, twice. On BULK (169M cells, 4.1M of them data) a draw takes 3.3 s and peaks at 0.95 GiB. Estimating 63 transition classes from 1,677 points takes 11 s: the estimator's run time grows with roughly the 2.5th power of the class count. +- **Found on the way:** `dft_rast_classify()` modifies the raster it is given (#89). + # drift 0.18.0 - **Dated reference imagery (#79).** The Esri and Google basemaps say what a place is but not when it changed. New `dft_stac_composite()` returns a cloud-masked median composite of any Sentinel-2 band roles for each year's run of calendar months, for example true colour for June–July in 2017 and in 2023. Values are surface reflectance, not display bytes, so a composite is data and can later be classifier input. It shares its read path with `dft_stac_cube()`. The cache entry is written as a Cloud Optimized GeoTIFF, so it can be served through titiler unchanged. floodplains#93 is the first consumer: a stratified accuracy assessment of IO LULC, reviewed at sample points across whole floodplains. diff --git a/planning/active/progress.md b/planning/active/progress.md index e177e5f..0ddd6d2 100644 --- a/planning/active/progress.md +++ b/planning/active/progress.md @@ -10,3 +10,4 @@ - Phase 1 done: Olofsson 2014 read from the local Zotero PDF; values in `helper-accuracy.R`; three printed values contradict the paper's own equations (findings). Issue body revised. Stehman 2014 is not needed (mapaccuracy tests it) - Phases 2–4 code and tests written; code-check rounds 1–2 found 6 defects (fixed), round 3 running; BULK scale run 3.3 s / 0.95 GiB for the sampler; filed #89 (classify mutates input) - Code-check: 3 rounds. R1: 4 findings (1 bug). R2: 2 (1 bug, 1 inside the R1 fix). R3: 6 (2 bugs, 1 inside the R2 fix); ended by enumerating 29 identity sites, and every fix was mutation-checked +- Phase 5: examples run; NEWS 0.19.0 drafted; CLAUDE.md accuracy block added and tile size corrected (314 x 326); floodplains#93 body updated with the delivered API. Full suite 1289 pass / 15 skip; R CMD check 0 errors, 1 WARNING (non-ASCII in R/dft_stac_fetch.R, already on main, untouched here), 1 NOTE (future timestamps) diff --git a/planning/active/task_plan.md b/planning/active/task_plan.md index ba28928..326a400 100644 --- a/planning/active/task_plan.md +++ b/planning/active/task_plan.md @@ -97,17 +97,17 @@ The package always applies the FPC `(1 − n_h/N_h)`, so a census stratum contri - [x] `R/dft_accuracy_size.R` ### Phase 5: Docs and release -- [ ] A runnable `@examples` block on every function (the estimator example uses Olofsson's published counts; the sampler uses the bundled tile with a satisfiable allocation) -- [ ] `devtools::document()`, `lintr::lint_package()`, `pkgdown::check_pkgdown()`. `_pkgdown.yml` has no `reference:` index, so the check cannot catch an omission; it stays that way (out of scope) +- [x] A runnable `@examples` block on every function (the estimator example uses Olofsson's published counts; the sampler uses the bundled tile with a satisfiable allocation) +- [x] `devtools::document()`, `lintr::lint_package()`, `pkgdown::check_pkgdown()`. `_pkgdown.yml` has no `reference:` index, so the check cannot catch an omission; it stays that way (out of scope) - [ ] NEWS.md, then version 0.19.0 as the final commit -- [ ] CLAUDE.md Core Pipeline: add the accuracy block, and correct the bundled tile's size (314 x 326, not 600 x 600) -- [ ] Update the floodplains#93 body's "What lives where" table: the `map =` argument, the `ref_class` composition recipe, the pilot-extension rule, and per-stratum sizing +- [x] CLAUDE.md Core Pipeline: add the accuracy block, and correct the bundled tile's size (314 x 326, not 600 x 600) +- [x] Update the floodplains#93 body's "What lives where" table: the `map =` argument, the `ref_class` composition recipe, the pilot-extension rule, and per-stratum sizing No vignette. One made with fabricated labels would illustrate a number nobody measured. The worked example belongs in floodplains#93 once real labels exist. ## Validation -- [ ] Tests pass +- [x] Tests pass - [x] `/code-check` clean on each commit -- [ ] PWF checkboxes match landed work +- [x] PWF checkboxes match landed work - [ ] `/planning-archive` on completion From ba399a8692db56bc566d7f2384b4cef016cb25a3 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Tue, 29 Sep 2026 07:19:13 -0700 Subject: [PATCH 08/10] Archive planning files for issue #81 Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- planning/active/.gitkeep | 0 .../README.md | 29 +++++++++++++++++++ .../findings.md | 0 .../progress.md | 0 .../review-plan.md | 0 .../review-round1.md | 0 .../review-round2.md | 0 .../review-round3.md | 0 .../task_plan.md | 2 +- 9 files changed, 30 insertions(+), 1 deletion(-) create mode 100644 planning/active/.gitkeep create mode 100644 planning/archive/2026-09-issue-81-accuracy-assessment/README.md rename planning/{active => archive/2026-09-issue-81-accuracy-assessment}/findings.md (100%) rename planning/{active => archive/2026-09-issue-81-accuracy-assessment}/progress.md (100%) rename planning/{active => archive/2026-09-issue-81-accuracy-assessment}/review-plan.md (100%) rename planning/{active => archive/2026-09-issue-81-accuracy-assessment}/review-round1.md (100%) rename planning/{active => archive/2026-09-issue-81-accuracy-assessment}/review-round2.md (100%) rename planning/{active => archive/2026-09-issue-81-accuracy-assessment}/review-round3.md (100%) rename planning/{active => archive/2026-09-issue-81-accuracy-assessment}/task_plan.md (99%) diff --git a/planning/active/.gitkeep b/planning/active/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/planning/archive/2026-09-issue-81-accuracy-assessment/README.md b/planning/archive/2026-09-issue-81-accuracy-assessment/README.md new file mode 100644 index 0000000..a4270c6 --- /dev/null +++ b/planning/archive/2026-09-issue-81-accuracy-assessment/README.md @@ -0,0 +1,29 @@ +## Outcome + +Added the `dft_accuracy_*` family (#81): stratified random point sampling from any strata raster, a reference-label contract, error-adjusted area and accuracy with confidence intervals, and sample sizing from a pilot, following Olofsson et al. 2014 and Stehman 2014. + +- **Estimator:** wraps `mapaccuracy::stehman2014()`. A soul session surfaced that CRAN package, and the user chose it over writing a second implementation. drift owns the sampler, the contract and the train guard, CIs in ha, and the per-stratum table the sizer needs. +- **Reference values:** the tests reproduce Olofsson Tables 8–9 through `dft_accuracy_estimate()`. Three printed values in the paper contradict its own equations; the tests pin the equation values and record the printed ones (`tests/testthat/helper-accuracy.R`). +- **Plan review:** it caught that `terra::freq()` reports labels while `readValues()` reports codes on factor strata, so the sampler's first design could not have joined its two passes. +- **Code-check:** three rounds, twice finding a defect inside the previous fix. All shared one mechanism: the same class/stratum/point identity derived separately at each use site. The loop ended by enumerating all 29 identity sites and routing every one through `accuracy_key()`, with each fix mutation-checked. +- **Found on the way:** `dft_rast_classify()` mutates its input raster in place (#89). + +## Measurement + +- **Olofsson 2014 through the wrapper:** areas match to the ha (21,157.8 / 11,686.2 / 285,769.9 / 581,386.2), and Table 9 matches to 4 dp. Area half-widths are within 1.1 ha; Stehman's FPC and exact z against the paper's 1.96 without FPC account for the difference. With equal stratum weights, deforestation comes out at 200,748 ha instead of 21,158. +- **Census oracle** (bundled tiles, 2017 map vs 2023 reference, 300 draws): + - Unbiased throughout: |bias z| < 2. + - Map-class strata at n = 25: Water's 95% CI covered **61%** (SE ratio 0.76), because 98 of the 7,127 map-Trees cells (1.4%) are reference Water and most draws see none. Coverage rose to 83% at n = 75 and 89% at n = 150. This is the Wald interval's small-sample weakness, now documented on the estimator, not an estimator bug. + - Changed/stable strata at n = 60: coverage 0.91–0.96. +- **`stehman2014()` runtime** at 1,000 points: 0.33 s / 2.0 s / 13.0 s at 20 / 40 / 80 classes, about k^2.5. Not filed upstream; documented. +- **BULK scale** (169,248,352 cells, 4.1M non-NA, 64 GB machine): + - `dft_accuracy_sample()`: 3.3 s (1.7 s per chunked pass), peak RSS 0.95 GiB. + - The 63-stratum transition map with n = 30 and 10 censuses: 3.0 s. The pipeline peaks at 4.63 GiB, dominated by the in-memory transition raster. + - The estimator on those 1,677 points: 10.8 s. +- **Suite:** 1289 pass / 15 skip. R CMD check: 0 errors. It has 1 WARNING, non-ASCII in `R/dft_stac_fetch.R`, which is already on main and untouched here. + +## Evidence + +Review records: `review-plan.md`, `review-round*.md` in this directory. The BULK run logs were in the session scratchpad (`bulk/bulk_run*.log`, `bulk/rss*.log`) and are not committed; the numbers above are their summary. + +Closed by: PR for branch `81-accuracy-assessment-for-change-maps-stra` diff --git a/planning/active/findings.md b/planning/archive/2026-09-issue-81-accuracy-assessment/findings.md similarity index 100% rename from planning/active/findings.md rename to planning/archive/2026-09-issue-81-accuracy-assessment/findings.md diff --git a/planning/active/progress.md b/planning/archive/2026-09-issue-81-accuracy-assessment/progress.md similarity index 100% rename from planning/active/progress.md rename to planning/archive/2026-09-issue-81-accuracy-assessment/progress.md diff --git a/planning/active/review-plan.md b/planning/archive/2026-09-issue-81-accuracy-assessment/review-plan.md similarity index 100% rename from planning/active/review-plan.md rename to planning/archive/2026-09-issue-81-accuracy-assessment/review-plan.md diff --git a/planning/active/review-round1.md b/planning/archive/2026-09-issue-81-accuracy-assessment/review-round1.md similarity index 100% rename from planning/active/review-round1.md rename to planning/archive/2026-09-issue-81-accuracy-assessment/review-round1.md diff --git a/planning/active/review-round2.md b/planning/archive/2026-09-issue-81-accuracy-assessment/review-round2.md similarity index 100% rename from planning/active/review-round2.md rename to planning/archive/2026-09-issue-81-accuracy-assessment/review-round2.md diff --git a/planning/active/review-round3.md b/planning/archive/2026-09-issue-81-accuracy-assessment/review-round3.md similarity index 100% rename from planning/active/review-round3.md rename to planning/archive/2026-09-issue-81-accuracy-assessment/review-round3.md diff --git a/planning/active/task_plan.md b/planning/archive/2026-09-issue-81-accuracy-assessment/task_plan.md similarity index 99% rename from planning/active/task_plan.md rename to planning/archive/2026-09-issue-81-accuracy-assessment/task_plan.md index 326a400..1bb6483 100644 --- a/planning/active/task_plan.md +++ b/planning/archive/2026-09-issue-81-accuracy-assessment/task_plan.md @@ -110,4 +110,4 @@ No vignette. One made with fabricated labels would illustrate a number nobody me - [x] Tests pass - [x] `/code-check` clean on each commit - [x] PWF checkboxes match landed work -- [ ] `/planning-archive` on completion +- [x] `/planning-archive` on completion From ce567683829a79c75660ff613646a05f25d9d7ae Mon Sep 17 00:00:00 2001 From: almac2022 Date: Tue, 29 Sep 2026 07:19:45 -0700 Subject: [PATCH 09/10] Release v0.19.0 (#81) Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- DESCRIPTION | 4 ++-- .../archive/2026-09-issue-81-accuracy-assessment/task_plan.md | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index dd3cfcb..c08e9b7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: drift Title: Detecting Riparian and Inland Floodplain Transitions -Version: 0.18.0 -Date: 2026-09-06 +Version: 0.19.0 +Date: 2026-09-29 Authors@R: c( person("Allan", "Irvine", , "al@newgraphenvironment.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-3495-2128")), diff --git a/planning/archive/2026-09-issue-81-accuracy-assessment/task_plan.md b/planning/archive/2026-09-issue-81-accuracy-assessment/task_plan.md index 1bb6483..c0cfcbf 100644 --- a/planning/archive/2026-09-issue-81-accuracy-assessment/task_plan.md +++ b/planning/archive/2026-09-issue-81-accuracy-assessment/task_plan.md @@ -99,7 +99,7 @@ The package always applies the FPC `(1 − n_h/N_h)`, so a census stratum contri ### Phase 5: Docs and release - [x] A runnable `@examples` block on every function (the estimator example uses Olofsson's published counts; the sampler uses the bundled tile with a satisfiable allocation) - [x] `devtools::document()`, `lintr::lint_package()`, `pkgdown::check_pkgdown()`. `_pkgdown.yml` has no `reference:` index, so the check cannot catch an omission; it stays that way (out of scope) -- [ ] NEWS.md, then version 0.19.0 as the final commit +- [x] NEWS.md, then version 0.19.0 as the final commit - [x] CLAUDE.md Core Pipeline: add the accuracy block, and correct the bundled tile's size (314 x 326, not 600 x 600) - [x] Update the floodplains#93 body's "What lives where" table: the `map =` argument, the `ref_class` composition recipe, the pilot-extension rule, and per-stratum sizing From 0a360bb857f367e2678c788b81110fd2354f44d2 Mon Sep 17 00:00:00 2001 From: almac2022 Date: Tue, 29 Sep 2026 07:20:16 -0700 Subject: [PATCH 10/10] Qualify bare #93 references as floodplains#93 in the #81 archive Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01PRhUJsuKABLfpBGktPoiBN --- .../2026-09-issue-81-accuracy-assessment/review-plan.md | 2 +- .../archive/2026-09-issue-81-accuracy-assessment/task_plan.md | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/planning/archive/2026-09-issue-81-accuracy-assessment/review-plan.md b/planning/archive/2026-09-issue-81-accuracy-assessment/review-plan.md index 18afe39..ff0fd3c 100644 --- a/planning/archive/2026-09-issue-81-accuracy-assessment/review-plan.md +++ b/planning/archive/2026-09-issue-81-accuracy-assessment/review-plan.md @@ -5,7 +5,7 @@ Read-only reviewer, so its findings came back as reply text; I transcribed them | id | finding | disposition | |---|---|---| | B1 | `terra::freq()` returns labels on factor strata while `readValues()` returns codes; `freq()` also rounds floats. Count and resolve steps could not join | **Confirmed by probe** (`Water -> Water` vs `1001`). Two passes on one chunked reader; `stratum_label` from `levels()` | -| B2 | Sizing through `ua` assumes strata = map classes; #93 needs per-stratum SD for an area target | Primary form takes `s_h` from the estimator's `$stratum`; `ua` kept as a convenience | +| B2 | Sizing through `ua` assumes strata = map classes; floodplains#93 needs per-stratum SD for an area target | Primary form takes `s_h` from the estimator's `$stratum`; `ua` kept as a convenience | | B3 | FPC undecided; `n_h > N_h` refusal breaks on tiny strata (the tile has 2- and 3-cell strata) | A stratum with `n_h ≥ N_h` becomes a census, with a message. FPC is now fixed on by `mapaccuracy`, so there is no `fpc` argument (superseded 2026-09-28) | | G1 | The sampler does not emit `map_class` | `map =` argument, extracted at each cell | | G2 | One `blocks()` chunk on the tile, so the carry logic is untested | Explicit row chunks, capped at about 1e7 cells; chunk-invariance test | diff --git a/planning/archive/2026-09-issue-81-accuracy-assessment/task_plan.md b/planning/archive/2026-09-issue-81-accuracy-assessment/task_plan.md index c0cfcbf..fed22e1 100644 --- a/planning/archive/2026-09-issue-81-accuracy-assessment/task_plan.md +++ b/planning/archive/2026-09-issue-81-accuracy-assessment/task_plan.md @@ -71,7 +71,7 @@ The package always applies the FPC `(1 − n_h/N_h)`, so a census stratum contri - [x] Census oracle (independent truth, no PDF needed): map = 2017 and "reference" = 2023 on the bundled tile, with the true error matrix and areas from `terra::crosstab`. Run about 500 stratified draws under map-class strata and under a changed/stable split; check bias ≈ 0, empirical SD ≈ mean SE, and CI coverage ≈ level. Skippable if slow - [x] Perfect labels (`ref = map`) give UA = PA = OA = 1, SE = 0, and adjusted area equal to mapped area. Recoding to a 2-class union gives an SE that is not the sum of the SEs - [x] Contract refusals: training rows, NA `ref_class`, duplicate ids, an unknown stratum, a stratum with no labels, `n_h = 1`. A reference-only class appears in the matrix, and PA is NA at `p̂_·j = 0` -- [x] Measure `stehman2014()` runtime at #93 scale (about 1,000 points × about 80 transition classes; it builds `classes²` indicator columns). File upstream and report if it is impractical +- [x] Measure `stehman2014()` runtime at floodplains#93 scale (about 1,000 points × about 80 transition classes; it builds `classes²` indicator columns). File upstream and report if it is impractical - [x] `R/dft_accuracy_estimate.R` + `R/dft_accuracy_labels.R`. Freeze the `$strata` and `$stratum` shapes here - [x] Restore-the-bug check: pass unweighted `N_h` inside the wrapper and confirm the published-value tests go red @@ -90,7 +90,7 @@ The package always applies the FPC `(1 − n_h/N_h)`, so a census stratum contri - `map =` extraction, including a grid-mismatch refusal - [x] `R/dft_accuracy_sample.R` - [x] Sampler → estimator integration: draw, fake labels from a reference raster, estimate (covered by the census oracle once both exist) -- [x] Scale test on BULK: `classified_2017.tif` and its `dft_rast_transition()` factor output (the #93 input), with an RSS sampler. Record pass-1 and pass-2 time and peak RSS in the PR body +- [x] Scale test on BULK: `classified_2017.tif` and its `dft_rast_transition()` factor output (the floodplains#93 input), with an RSS sampler. Record pass-1 and pass-2 time and peak RSS in the PR body ### Phase 4: Sizing - [x] `test-dft_accuracy_size.R`: reproduce Olofsson §5.1.1's sample-size example (n and allocation); the `s_h` form from a pilot `$stratum` agrees with the `ua` form when strata = map classes; edge cases (UA = 1, a zero weight)