# Real demo: forest-change screening in Đắk Glong, Central Highlands (Viet Nam)

Every number here was produced by the scripts in `/home/yitec/redstone_forest/data/scripts/` from public open data on 2026-10-18. No number is invented. Where a figure depends on an assumption, the assumption is stated next to it. **All machine-readable numbers are in `demo_summary.json`.**

## 1. Area of interest (AOI)

- **Location:** 10 × 10 km window, 100% inside **Đắk Glong district, Đắk Nông province**. These are pre-1/7/2025 boundaries; after the 2025 merger the area belongs to Lâm Đồng province.
- **Bounds:** lon 107.8737–107.9664 E, lat 11.9544–12.0456 N (UTM 48N 813000–823000 E, 1323240–1333240 N). Elevation is 535–781 m (Copernicus DEM).
- **Why this window:** I scanned Hansen GFC v1.13 over the 5 Central Highlands provinces in ~10 km blocks. This was the highest-loss cluster in Đắk Nông for 2019–2024 (`hansen_central_highlands.json`, `hotspots_10km`).
- **What the landscape is:** a forest-to-agriculture mosaic, not intact forest. Hansen marks 6,855 ha of loss in 2001–2025 out of 9,491 ha with treecover2000 > 30%. Loss peaked in 2016 (758 ha), and the remaining forest core sits in the north-centre of the window.

## 2. Data (provenance)

| Layer | Source | Detail |
|---|---|---|
| Tree-cover loss and treecover2000 | Hansen/UMD/Google/USGS/NASA **GFC-2025-v1.13** (loss to 2025). https://storage.googleapis.com/earthenginepartners-hansen/GFC-2025-v1.13/download.html | Tile 20N_100E, 30 m. Resampled nearest-neighbour onto the 10 m grid. |
| Sentinel-2 L2A, before | `S2A_48PZU_20190128_0_L2A`, 2019-01-28 03:24 UTC | Scene cloud 0.50%; AOI cloud/shadow/nodata (SCL) 0.004%. Processing baseline 02.11. |
| Sentinel-2 L2A, after | `S2B_48PZU_20260126_0_L2A`, 2026-01-26 03:24 UTC | Scene cloud 13.85%; AOI cloud/shadow/nodata 0.003%. Baseline 05.11. |
| STAC | Element84 Earth Search v1, https://earth-search.aws.element84.com/v1, collection `sentinel-2-l2a` | Windowed COG reads of B02, B03, B04, B08 (10 m), B11, B12 and SCL (20 m → 10 m bilinear/nearest). |
| DEM | Copernicus DEM GLO-30, public on AWS (`copernicus-dem-30m`): tiles N11_E107 and N12_E107 | Bilinear to 256×256 (39 m cells). |
| Boundaries | geoBoundaries gbOpen VNM ADM1/ADM2 | Pre-2025 provinces and districts. |
| Viet Nam forest area | Decision 1106/QĐ-BNNMT (31/3/2026): **14,971,553 ha** at 31/12/2025 | https://luatvietnam.vn/nong-nghiep/quyet-dinh-1106-qd-bnnmt-2026-cong-bo-hien-trang-rung-toan-quoc-nam-2025-431538-d1.html |
| LLM prices | https://platform.claude.com/docs/en/about-claude/pricing (fetched 2026-10-18) | Image tokens = ⌈w/28⌉·⌈h/28⌉ per https://platform.claude.com/docs/en/build-with-claude/vision |

**Scene choice:** both scenes are dry-season and only 2 days of year apart (28 Jan and 26 Jan), which minimises phenology differences.

**Reflectance offset:** for Earth Search items with `earthsearch:boa_offset_applied = true`, the DNs are already harmonised. The STAC `raster:bands` offset of −0.1 must therefore **not** be applied again, because doing so gives negative reflectance. I verified this empirically and handle it in `03_extract.py`.

## 3. Method

1. **Hansen statistics** (`01`, `02`, `04`): annual loss in ha, defined as Hansen loss on pixels with treecover2000 > 30% (the GFW default).
   - Computed for the AOI and for each of the 5 Central Highlands provinces (Kon Tum, Gia Lai, Đắk Lắk, Đắk Nông, Lâm Đồng).
   - Province areas use a per-row cos(lat) pixel-area approximation, accurate to about 0.5%.
   - **Central Highlands total:** 28,423 ha (2019), 28,310 (2020), 19,525 (2021), 15,836 (2022), 13,575 (2023), 14,346 (2024), 12,239 (2025).
   - **AOI:** 367, 236, 348, 208, 205, 118, 103 ha for 2019–2025, totalling 1,582 ha.
2. **Imagery** (`04`): true colour with one common percentile stretch for both dates so they are comparable, SWIR false colour, NDVI, NBR, ΔNDVI and dNBR maps.
3. **System-1 screener** (`05`): LightGBM, 400 trees, 63 leaves, 32 per-pixel features.
   - Features: 6 bands × 2 dates, NDVI/NBR/NDMI for both dates, their deltas, and 5×5 (50 m) neighbourhood means/std.
   - **Domain:** the "forest at start of 2019" mask, i.e. treecover2000 > 30% and no Hansen loss 2001–2018. That leaves 421,766 pixels, 37.5% of them positive.
   - **Labels (weak):** 1 = Hansen loss year 2019–2025; 0 = no Hansen loss 2001–2025.
   - **Validation:** spatial 2-fold split (train west half → test east half, and the reverse), with a 300 m buffer strip between halves. Every metric and every alert pixel is **out-of-fold**. The decision threshold is tuned on an inner spatial validation strip inside the training half only.
4. **Alerts** (`06`): connected components (8-connected) of out-of-fold p ≥ 0.41 inside the forest-2019 mask, after a 2×2 morphological opening, keeping area > 0.5 ha.
   - Each alert gets: area, centroid, mean probability (confidence), mean ΔNDVI/dNBR, the fraction overlapping Hansen 2019–2025 loss, and context (pre-2019 loss and stable forest in a 300 m ring, which hints at mosaic or plantation versus a forest frontier).
   - **Tiers:** High = ≥ 5 ha AND confidence ≥ 0.75 AND ΔNDVI ≤ −0.15. Medium = ≥ 2 ha AND confidence ≥ 0.6. Everything else is Low.
   - **Ranking:** by area × confidence × max(0, −ΔNDVI).
5. **Cost model** (`07`): measured CPU throughput extrapolated to 14.97 M ha, compared with three LLM scenarios.

## 4. Results

**Model metrics** (agreement with Hansen, out-of-fold):

| Fold | Threshold | Precision | Recall | F1 | ROC-AUC | PR-AUC | ΔNDVI-rule baseline F1 |
|---|---|---|---|---|---|---|---|
| train W → test E | 0.41 (tuned) | 0.593 | 0.743 | **0.659** | 0.820 | 0.680 | 0.491 |
| train W → test E | 0.50 | 0.632 | 0.661 | 0.646 | | | |
| train E → test W | 0.28 (tuned) | 0.645 | 0.844 | **0.731** | 0.829 | 0.771 | 0.533 |
| train E → test W | 0.50 | 0.752 | 0.605 | 0.670 | | | |

- **Confusion matrix** (W→E, threshold 0.41): TN 106,967, FP 35,350, FN 17,809, TP 51,453.
- **Without the forest mask** (all valid pixels in the east half), precision falls to 0.24 (recall 0.74, F1 0.36), because cropland cycles look like loss. A good forest mask matters as much as the model.
- **Speed** (Intel Xeon E5-2686 v4, 36 threads, shared host under load), end-to-end feature computation plus prediction:
  - 36 threads: **337,000 px/s** (420,000 px/s prediction only).
  - 1 thread: 17,900 px/s.
- **Alerts:** 564 alerts > 0.5 ha (1,316 ha in total): **8 High, 142 Medium, 414 Low**.
  - The top alert, DG-001, is 54.4 ha at 12.0274 N, 107.9370 E, confidence 0.76, ΔNDVI −0.195. 72% of it overlaps Hansen loss, dominant year 2019.

**Cost comparison: screening all 14.97 M ha (1.497 billion 10 m pixels) once**

- **System 1 (this model):** 1.23 h on a 36-thread CPU, about **US$1.8** of compute at US$1.428/h (AWS c7i.8xlarge on-demand, per instances.vantage.sh). This **excludes** downloading and reading the imagery, which dominates at national scale (≈50 s per 1 M-pixel window, sequential, in this demo) and must be parallelised.
- **LLM, per-pixel text** (same job, **assumption:** 100 input + 2 output tokens per pixel, so 150 billion input tokens):
  - Claude Opus 5.5: ≈ US$659k (US$329k with the Batch API).
  - Claude Sonnet 5.5: ≈ US$329k (US$165k batch).
  - Claude Haiku 5.5: ≈ US$16.5k (US$8.2k batch).
  - This also takes weeks under typical rate limits (account-dependent, not computed).
- **LLM, looking at image tiles** (about 3,000 tiles of 1000×1000 px, before and after): only about US$53 (Opus 5.5) or US$1.3 (Haiku 5.5). But it is **a weaker, different job**:
  - The model sees an RGB rendering only, with no NIR/SWIR values.
  - Locations are approximate.
  - There is no calibrated per-pixel probability, and results are not reproducible pixel-for-pixel.
  - I report this scenario so nobody can say we hid the cheap LLM option.
- **Proposed System 1 → System 2:** the screener reduces 1,000,000 pixels to 150 High+Medium alerts (a 6,667× reduction).
  - An LLM reviewing each alert dossier costs about US$0.02 with Opus 5.5 (**assumption:** 1,900 input and 600 output tokens: four 256×256 chips plus context). That is about US$2.9 for this AOI.
  - A naive national extrapolation from this unusually active AOI gives about 530k alerts and about US$10k. This is an **upper bound**, since most of Viet Nam's forest changes far less than this frontier.
- **The real argument:** System 1 is what makes per-pixel, every-revisit (5-day) screening at 10 m possible and auditable. The LLM adds value where judgement, context, legal reasoning and documents are needed, on a few hundred cases rather than billions of pixels.

## 5. Limitations (read before presenting)

1. **Weak labels:** the metrics measure *agreement with Hansen GFC*, not ground truth. Hansen has its own errors: it misses selective logging and degradation, and it has a commission/omission rate that varies by region. No field data or high-resolution verification was used.
2. **Plantation harvest is counted as "loss":** Hansen loss is any stand-replacement disturbance, so it includes acacia/eucalyptus rotation harvest (5–7 year cycles) and coffee/rubber/pepper clearing and replanting. About 33% of Viet Nam's forest is planted (4.89 M ha of 14.97 M ha, Decision 1106/QĐ-BNNMT).
   - **Neither Hansen nor this model can tell illegal clearing from legal harvest.** Many alerts here are probably agricultural or plantation dynamics.
   - That distinction needs the forest-status map (hiện trạng rừng), forest-owner records and harvest permits, which is exactly where the System-2 / knowledge-graph layer comes in.
3. **Label/date mismatch:** Hansen 2019 loss that happened before 28 Jan 2019, and early-2019 loss that regrew by Jan 2026, are noisy labels. Hansen 30 m pixels were also resampled to 10 m.
4. **One AOI, two dates:** results come from one 100 km² mosaic landscape with 37.5% positive rate in the domain. Precision/recall will differ in intact forest (lower base rate means lower precision at the same threshold). The model is not trained for degradation, fire, or cloudy-season imagery. The national cost extrapolation is an order of magnitude estimate.
5. **Speed:** measured on a shared host. It excludes imagery I/O, cloud compositing and training (training took 12–76 s per fold on about 200k pixels).
6. **Boundaries:** geoBoundaries ADM1/ADM2 reflect pre-2025 provinces and districts and are not official cadastral boundaries.
7. **FIRMS active-fire data was not included:** the FIRMS servers did not respond from this host, and the APIs need a key (`firms_summary.json`).
8. **Token assumptions** (100 tokens per pixel row; 1,900/600 per alert) were not measured with `count_tokens` because no API key was available. Treat them as ±30%.

## 6. Files

- **Summary:** `demo_summary.json` (all deck numbers), `cost_model.json`, `model_metrics.json`, `hansen_aoi_loss_by_year.json`, `hansen_central_highlands.json` (provinces and 10 km hotspots), `hansen_hotspots.json` (tile-wide scan, includes Laos/Cambodia), `alerts_summary.json`, `alerts_top20.md`, `alerts.geojson` (all 564, full properties), `firms_summary.json`.
- **Images (1024 px PNG):** `truecolor_before.png`, `truecolor_after.png`, `falsecolor_swir_{before,after}.png`, `ndvi_{before,after}.png`, `nbr_{before,after}.png`, `change_dndvi.png` (brown = NDVI drop), `change_dnbr.png`, `hansen_lossyear_map.png` (+ `_legend.png`; grey = loss 2001–2018, yellow→red = 2019→2025), `model_probability_map.png`, `model_agreement_map.png` (yellow = model and Hansen agree on loss, magenta = model-only, cyan = Hansen-only), `alerts_overlay.png`, `chart_aoi_loss_by_year.png`, `chart_central_highlands_loss_by_year.png`.
- **`web/`:**
  - Images: `truecolor_before.jpg` and `truecolor_after.jpg` (1024 px, ≤ 400 KB), `change_dndvi.jpg`, `alerts_overlay.jpg`.
  - Data: `loss_by_year.json` (AOI, Central Highlands, per province), `alerts.geojson` (High+Medium, simplified), `alerts_points.json` (all 564 centroids).
  - DEM: `dem_256.json` (Float32 values plus bounds; row 0 = north), `dem_256_16bit.png` (height = min + v/65535·(max−min), min 535.4 m, max 780.8 m), `dem_256_float32.bin` (little-endian raw).
  - The true-colour JPGs cover exactly the same 10×10 km footprint as the DEM, so they can be draped directly.

## 7. Reproduce

```
cd /home/yitec/redstone_forest/data
.venv/bin/python -I scripts/01_hansen_scan.py   # optional tile-wide scan
.venv/bin/python -I scripts/02_hansen_stats.py
.venv/bin/python -I scripts/03_extract.py
.venv/bin/python -I scripts/04_visuals.py
.venv/bin/python -I scripts/05_model.py
.venv/bin/python -I scripts/06_alerts.py
.venv/bin/python -I scripts/07_summary.py
```

Raw inputs are in `data/raw/{hansen,boundaries}`. Sentinel-2 and the DEM are streamed by windowed reads, and intermediates are in `data/work/`.

## 8. v2 additions after the independent technical critique (`research/critique_tech.md`)

All v1 files, fields and numbers above are **unchanged**. Everything below is additive.

### 8.1 Naming and scores (E)
- **The "alerts" are retrospective change areas** between 2019-01-28 and 2026-01-26 (Vietnamese: *vùng biến động 2019→2026*), not real-time alerts.
  - Dominant Hansen loss year per alert: 2019: 115, 2020: 86, 2021: 129, 2022: 66, 2023: 66, 2024: 45, 2025: 37, none: 20.
  - All 8 High alerts are 2019–2021 events.
- **New fields on every alert** (in `alerts.geojson`, `alerts_v2.geojson`, `web/alerts.geojson`, `web/alerts_v2.geojson`):
  - `risk_score_uncalibrated`, which equals the existing `confidence` (kept for the slide builder). It is the raw mean LightGBM score, **not a probability**.
  - `observed_between: "2019-01-28/2026-01-26"`.
  - `alert_kind`, `tier_rule_note` (tiers = *quy tắc minh họa, chưa hiệu chỉnh*, an illustrative rule, not calibrated), `treecover2000_mean`, `context_class`.

### 8.2 Landscape context (A): a heuristic from Hansen history only, not a legal or land-use determination
**Rule:**
- `mosaic_rotation_likely` if ≥ 40% of the 300 m ring was lost in 2001–2018;
- else `forest_edge` if ≥ 60% of the ring is stable forest AND the polygon's mean treecover2000 is ≥ 70;
- else `mixed`.

| context_class | alerts | area (ha) | High / Medium / Low |
|---|---|---|---|
| mosaic_rotation_likely | 443 | 1,075.4 | 8 / 115 / 320 |
| mixed | 111 | 226.6 | 0 / 24 / 87 |
| forest_edge | 10 | 13.6 | 0 / 3 / 7 |

**82% of the alert area lies in already-converted mosaic.** The screener mostly finds land-use dynamics (rotation, re-clearing, crops). Telling lawful harvest from conversion needs the forest-status map, parcel records and harvest files, which is System 2's job.

**Files:** `alerts_context.png` and `web/alerts_context.jpg` (blue = mosaic, aqua = mixed, orange = forest edge, dark-green tint = Hansen stable forest), `alerts_context_legend.png` (with legend), `web/alerts_context_legend.json`, `alerts_summary.json` → `context_classes`.

### 8.3 Honest metrics pack (C): `metrics_pack.json`, `web/metrics_charts.json`, `chart_pr_curves.png`, `chart_reliability.png`
- **Both test halves:** LightGBM ROC-AUC 0.820 / 0.829 and PR-AUC 0.680 / 0.771. The halves have different prevalence (0.33 vs 0.43).
  - Pooled out-of-fold: ROC-AUC 0.821, PR-AUC 0.714.
  - Baselines on the same pixels, pooled: −ΔNDVI ROC-AUC 0.633 / PR-AUC 0.541; ΔNBR (before − after) 0.608 / 0.555.
- **PR curves:** 50 points each, recall/precision/threshold.
- **Base-rate table** (same TPR/FPR, assumed to transfer, which is optimistic). Precision at the tuned threshold, W→E (TPR 0.743, FPR 0.248):

| prevalence | 37.5% | 10% | 5% | 1% | 0.1% |
|---|---|---|---|---|---|
| precision | 0.64 | 0.25 | 0.14 | 0.03 | 0.003 |

  - E→W (TPR 0.844, FPR 0.344) gives 0.60 / 0.21 / 0.11 / 0.02 / 0.0025.
  - **At a realistic ~1% loss prevalence, about 33–40 of every 41 flags would be false.** The demo shows the mechanics in a hotspot; it is **not a deployable screener**. Operational systems add time-series confirmation, a minimum mapping unit, forest masks and multi-sensor agreement.
- **Ablation** (re-trained, same folds, 200 rounds for every group), ROC-AUC W→E / E→W:
  - before-only 0.650 / 0.688 (no post-event information, so this is a pure site prior: the model also learns which land tends to get cleared);
  - after-only 0.730 / 0.726;
  - delta-only 0.788 / 0.792;
  - all pixel features 0.803 / 0.819;
  - full model with 5×5 context 0.819 / 0.830.
- **Calibration** (10 bins), against Hansen at this AOI's 37.5% prevalence only:
  - Raw score: ECE 0.018, Brier 0.166.
  - Isotonic cross-fitted between halves: ECE 0.037, Brier 0.175. It is *worse*, because calibration fitted at one prevalence (0.43) does not transfer to another (0.33).
  - Lesson: re-calibrate on field verdicts at the operating prevalence before saying "% likely".
- **Speed:** 337k px/s on 36 threads, so ≈ 74 min for one national pass (1.497 billion px), compute only. This excludes imagery I/O. The host was shared, and load reached ~110 on 36 threads during the later runs.

### 8.4 Chips for human photo-interpretation (D): `web/chips/`, `web/chips.json`
- 30 alerts: the top 20 plus 10 random Low-tier (seed 20261019).
- For each: `<id>_before.jpg`, `<id>_after.jpg` (true colour, same stretch on both dates, 256 px) and `<id>_{before,after}_swir.jpg` (B11/B08/B04). A yellow line marks the polygon. Footprints are ≥ 640 m.
- `chips.json` has empty `human_label` fields with options: natural_forest_clearing / plantation_harvest / crop_or_seasonal / no_change_visible / unclear.
- **No human labels exist yet.** Only report them if someone actually labels the chips, and state n.

### 8.5 One-year run 2025-01-26 → 2026-01-26 (B): `oneyear_2025_2026/`
- **Scenes and domain:** before `S2C_48PZU_20250126_0_L2A` (scene cloud 12.4%, AOI cloud 0.000%) and the same after scene as v1, at the same day of year. Labels are Hansen loss year 2025 (1) against no loss 2001–2025 (0). The domain is forest at the start of 2025: 273,816 px, of which 3.75% are positive. Hansen 2025 loss in the AOI is 102.5 ha (canopy > 30%).
- **Temporal transfer** (the 2019→2026 model applied unchanged to 2025→2026, east half): ROC-AUC 0.860, PR-AUC 0.140.
  - At its original threshold 0.41: precision 0.157, recall 0.674.
  - The −ΔNDVI baseline gets ROC-AUC 0.726 / PR-AUC 0.176.
  - The ranking transfers reasonably well; the threshold does not.
- **Re-trained one-year model** (same spatial folds), tuned thresholds 0.18 / 0.16:
  - W→E: precision 0.41, recall 0.54, F1 0.47, ROC-AUC 0.93, PR-AUC 0.43.
  - E→W: precision 0.39, recall 0.46, F1 0.43, ROC-AUC 0.90, PR-AUC 0.39.
  - Base-rate precision at 1% prevalence is 0.13–0.17, better than the 7-year model's 0.02–0.03, but still mostly false flags.
- **Change areas > 0.5 ha:** 21, totalling 20.6 ha (High 1, Medium 6, Low 14; one-year tiers: High ≥ 2 ha, score ≥ 0.6, ΔNDVI ≤ −0.15; Medium ≥ 1 ha, score ≥ 0.5).
  - All 21 are `mosaic_rotation_likely`.
  - They cover only about 20% of Hansen's 102.5 ha, because one year of loss is small and fragmented and much of it falls below 0.5 ha after the 2×2 opening. Object-level recall is low.
- **Caveats:** Hansen 2025 loss that happened before 26 Jan 2025 is mislabelled as change, and loss in Jan 2026 is unlabelled.
- **Files:** `oneyear_metrics.json`, `alerts_2025_2026.geojson`, `truecolor_2025-01-26.png`, `alerts_2025_2026_overlay.png`, and `web/` versions (JPG ≤ 400 KB, GeoJSON).

### 8.6 New data caveat: Earth Search offset metadata is unreliable
`S2C_48PZU_20250126_0_L2A` reports `earthsearch:boa_offset_applied = false`, yet its B02 DNs go down to about 150. That is impossible for raw baseline ≥ 04.00 data, which carry +1000.

`10_one_year.py` therefore decides the offset **empirically, once per item, from B02** (offset present only if the 0.5th-percentile DN is ≥ 1000). A per-band test is wrong, because NIR/SWIR exceed 1000 DN anyway; I hit that bug and fixed it.

The 2019 and 2026 scenes used in v1 pass the same empirical test, so v1 is unaffected. After the fix, the 2025 vs 2026 medians agree: B02 0.030 vs 0.033, B08 0.315 vs 0.326, B12 0.115 vs 0.115.

### 8.7 Scripts added
- `08_metrics_pack.py`
- `09_chips.py`
- `10_one_year.py`
- `06_alerts.py` (extended additively)
- `07_summary.py` (adds `metrics_pack`, `one_year_2025_2026`, `chips` and context keys to `demo_summary.json`)