# Deviations, implementation choices and extensions

Each entry states what the paper did (or did not specify), what this replication
does, and why. "Gap n" refers to the numbered list of unspecified details in
`nanopubs/drafts/00_paper_summary.md`. The code that implements each choice is
named so the entry can be checked against it.

## Scope

| # | Paper | This replication | Code |
|---|---|---|---|
| S1 | Pantropical (±25°, tropical/subtropical moist, dry, coniferous forest biomes) | **Colombia only** (GADM 4.1 level 0) ∩ RESOLVE 2017 biomes 1–3. Costa Rica is used only for a step-1 cross-check. | `02b`, `03a` |
| S2 | Authors' code "available on request" | Independent code written from the Methods and SI; the authors' code was not requested or used. Step 1 reuses the authors' published **outputs** (Zenodo 10.5281/zenodo.7428804). | all |

## Step 1 — reproduction from the published map

| # | Paper | This replication | Code |
|---|---|---|---|
| R1 | Area "calculated in Mollweide projection", country sums with GADM (2022) | Primary: **exact WGS84 ellipsoidal area** of each 0.00025° pixel. Sensitivity: sphere of radius *a* (identical to PROJ's WGS84 Mollweide) and nominal 0.09 ha per pixel. GADM 4.1, pixel-centre rule. | `03a` |
| R2 | Continuous product stored as integer % (gap 14) | Integer storage and sub-1 % NoData reported as bounds, not corrected. | `03a` |
| R3 | "Available for restoration" mask undefined (gap 13) | Candidate masks reported (continuous-valid, binary-valid, and binary-valid within our study domain). | `03a`, `03` |

## Step 2 — independent random-forest replication

| # | Paper (gap) | This replication | Why | Code |
|---|---|---|---|---|
| M1 | R `randomForest`; hyperparameters not reported (gap 1) | `sklearn.ensemble.RandomForestClassifier` with the R defaults for classification: `n_estimators=500`, `max_features=floor(sqrt(p))` (3 for 10 variables), `min_samples_leaf=1`, bootstrap n with replacement, Gini. Probability = `predict_proba` (equals R vote fraction for pure leaves). Seed 20261003. | Closest open equivalent of the default R model | `03`, `03c` |
| M2 | Categorical land cover and biome as R factors | Integer-coded (sklearn random forests have no native categorical splits) | Library limitation; affects which splits are reachable, not the information available | `03` |
| M3 | Sample sizes: 6 M pool, 1 M final model, 10 × 500 k for selection, 4.87 M validation (gap 9) | Scaled by Colombia's share of Fagan regrowth patches (188,921 / 4.78 M = 3.95 %): final model 39,523 records, selection fits 19,762, pool + validation ≈ 429 k (split 6 : 4.87), balanced by class. Validation points never used in training. | "Proportional to the paper" | `02b`, `03` |
| M4 | Class-1 points sampled in regrowth polygons with `st_sample`; class 0 "ellipsoid-aware" in the domain | Bernoulli sampling of **30 m pixels**: class 1 = pixel centre inside a Fagan `regrowth` polygon; class 0 = pixel in the class-0 domain; inclusion probability ∝ exact pixel area. Points are pixel centres. | Covariates are per pixel; equivalent to point sampling followed by pixel lookup | `02b` |
| M5 | Fagan label version not stated (gap 18) | GFW S3 GeoPackage `pantropical_tree_plantation_expansion_2000_2012.gpkg` (ETag recorded in `data/raw/sources_predictors.json`); **all confidence ranks** of `regrowth` are class 1; `plantation` and `open` polygons are excluded from class 0. | The paper excludes "forestry" polygons; `open` (low-confidence) polygons are also tree gain and not clearly class 0 | `01b`, `02b` |
| M6 | "Tree cover in 2000" threshold not reported (gap 2) | Forest = Hansen GFC v1.13 `treecover2000` ≥ **30 %** | Common GFW default | `02b` |
| M7 | 2018 tree cover construction not reported (gap 2). Peer review file, Response 35: the authors used "the 2018 layer" from the Hansen GFC **v1.6** download page; v1.6 has no tree-cover-2018 layer, only treecover2000/gain/lossyear | Forest 2018 = forest 2000 AND NOT `lossyear` 2001–2018, from GFC **v1.13** (gain not added). Changed 2026-10-04 from an earlier "(forest 2000 OR gain) AND NOT loss" | Usual GFW convention; GFW advises against combining gain and loss. Version v1.13 vs v1.6: loss 2011+ was reprocessed in later releases | `02b` |
| M8 | Forest density: "1 km² area" / "1 km radius" / "1-km circular buffer" (gap 3) | Mean forest fraction in a **1 km-radius** disk on the 30 m grid (ellipsoidal metres per degree at the tile centre) | Two of the three descriptions say radius | `02b` |
| M9 | Distance to forest: method not reported | Euclidean distance on the 30 m grid (ellipsoidal metres per degree at tile centre), computed in 1° tiles with a 0.25° margin and **truncated at 25 km** | Bounded memory; regrowth is almost always within a few km of forest | `02b` |
| M10 | Land cover: ESA CCI simplified 31 → 11 classes, mapping not listed; training year implied 2000 (gap 4) | Our 11-class mapping (`data/clean/landcover_classes.csv`); 2000 for training, 2015 for prediction | — | `02` |
| M11 | Prediction: "2015 ESA CCI combined with 2018 tree cover"; NoData for 2018 forest, open water, urban, rock/bare (gap 11) | Prediction domain = study domain − forest 2018 − ESA 2015 classes 190, 200–202, 210. Sparse vegetation (150–153) is predicted (excluded only from class-0 training). "Combined" is read as: land-cover value from ESA 2015, forest mask from 2018 tree cover. | Literal reading of the Methods | `02b` |
| M12 | Resampling of 250 m / 300 m / 1 km layers to 30 m not reported (gap 12) | **Nearest neighbour** at the 30 m pixel centre for all coarse layers (categorical and continuous) | No information is invented; same as nearest resampling | `02b` |
| M13 | Soil: "SoilGrids250m", OCDENS / PHIHOX, top 30 cm depth-weighted (gap 8) | SoilGrids **v2017** (2017-03-10 release) OCDENS and PHIHOX, depths 0/5/15/30 cm, trapezoidal mean over 0–30 cm | Matches the paper's variable names (decision by Anne) | `02` |
| M14 | Bioclim PCA: `prcomp` on 1 M random land points; loadings, centring and scaling not published (Supp. Table 1 gives only SDs) | Re-derived: 1,000,000 points uniform by area on the sphere, kept where WorldClim is valid (global land incl. Antarctica); both `prcomp` variants fitted; the one whose component SDs match Supp. Table 1 is used (`data/clean/bioclim_pca_summary.csv`). **Result: centre-only (R default, `scale.=FALSE`)** — cumulative variance 0.791 / 0.940 / 0.974 / 0.991 / 0.996 vs paper 0.780 / 0.948 / 0.976 / 0.990 / 0.994; SDs 893 / 387 / 184 / 131 / 76 vs paper 780 / 362 / 147 / 103 / 60 (~15–25 % higher, consistent with a different land sample, e.g. Antarctica or the sampling scheme). Centre + scale gives SDs 3.1 / 2.1 / 1.4 / 1.1 / 0.8 and is clearly **not** what the paper used, so the brief's "centred + scaled" was not followed. | Loadings unavailable; the data decide the variant | `02` |
| M15 | Biome predictor (Dinerstein 2017) | RESOLVE Ecoregions 2017 `BIOME_NUM`, rasterised at 30 m (pixel centre) | — | `02`, `02b` |
| M16 | Areas in Mollweide | Exact WGS84 pixel area (primary); Mollweide-sphere and 0.09 ha as sensitivity columns. Prediction areas are estimated from a **1-in-100 systematic sample** of 30 m pixels (every 10th row and column), not from a full 30 m prediction. Domain areas and the authors' map within our domain are exact 30 m sums. | Full 30 m RF prediction over ~1 G pixels was not necessary for area totals; systematic sampling error is reported | `02b`, `03` |
| M17 | No Colombia accuracy in the paper (gap 19) | Compared with the global (87.9 % / 87.8 %) and Neotropics (86.5 %) figures | — | `03` |

## Extension (not a deviation): probability calibration

The paper sums uncalibrated RF vote fractions trained on 50/50 samples ("expected
area"). Those scores are not probabilities at the true regrowth prevalence. We
additionally report (`03_analysis.py`):

- prevalence π = exact area of Fagan regrowth pixels / exact area of the sampling
  domain (class-0 ∪ class-1 masks), from the 30 m tile sums of `02b`;
- prior-shift correction p′ = π·p / (π·p + (1 − π)(1 − p)) applied to the scores,
  with expected area Σ p′ × area and area with p′ > 0.5;
- cross-check: isotonic calibration fitted on a random half of the validation set
  with case weights that reweight the balanced set to prevalence π (equivalent in
  expectation to drawing at prevalence π); never used in training;
- Brier score and reliability diagram on the other (evaluation) half, weighted
  the same way.

## Step 3 — robustness

| # | Choice | Code |
|---|---|---|
| B1 | (a) Land cover from ESA CCI **1992** (earliest year) and **1999** replaces 2000 in training; samples unchanged (class-0 domain still defined with ESA 2000). A check of Evaluator 1's leakage concern, not a correction: the Fagan regrowth labels start after 2000 and 2000 is the Hansen forest baseline, so 2000 land cover describes the starting state. The predictors that do overlap the outcome period are tested in (c). | `03c` |
| B2 | (b) 5-fold CV on the training pool: stratified random folds vs `GroupKFold` grouped by HEALPix NESTED cell (healpix-geo, `ellipsoid="WGS84"`) at depths 6, 7, 8; each fold trains on a 39,523-record draw. Distance-to-training analysis: 3-D chord distance on the WGS84 ellipsoid, 0.5 km bins to 8 km | `03c` |
| B3 | (c) Selection: 10 fits on 19,762 balanced records; importance = permutation mean decrease in accuracy on a separate 10,000-record balanced hold-out from the pool (R uses OOB permutation); rank sum; forward addition evaluated on the validation set; selected size = smallest k within 0.1 percentage points of the best. Candidates: final 10 + PC5, slope (SRTMGL1), cropland density (ESA 2000, 5 km radius), distance to urban (ESA 2000); "with" adds NPP (MOD17A3 C5.5 mean 2000–2015, NTSG), burned area (GlobFire, fraction of 204 months burned 2001–2017 at 15″), road density (GRIP4 5′). | `03c` |
| B3a | Burned area source: GlobFire monthly `FinalArea` polygons (final perimeter of each fire; `ActiveArea` = daily spread, not used), rasterised at 15″. If any GlobFire month cannot be fetched after 8 attempts, the whole series switches to MODIS MCD64A1 v061 (`Burn Date` > 0, 6 sinusoidal tiles, nearest to 15″) — the product GlobFire is built from. The source actually used is recorded in `data/raw/burned_source.json` and in `burned.nc` attributes. | `01b`, `02` |
| B3b | Road density: GRIP4 `grip4_total_dens_m_km2.asc` (5′); the zip also contains `grip4_area_land_km2.asc`, which an early version read by mistake (caught in the smoke test, fixed before the full run). | `02` |
| B4 | Not included: 10 of the 12 SoilGrids properties, distance to water (Kummu 2011 raster not distributed), protected areas (WDPA Aug 2020 not publicly archived), socioeconomic variables | — |
| B5 | (d) `HistGradientBoostingClassifier`, sklearn defaults, native categorical splits for land cover and biome. Optional rule (`optional_variants`), run after (a)–(c). | `03c` (`VARIANTS=de`) |
| B6 | (e) **SENSITIVITY check, not a fix:** bioclim PC1–PC4 from CHELSA v2.1 1981–2010 (`https://os.zhdk.cloud.switch.ch/chelsav2/GLOBAL/climatologies/1981-2010/bio/CHELSA_bio{1..19}_1981-2010_V.2.1.tif`, files dated 2024-09-13) instead of WorldClim v2.1 1970–2000. Same PCA procedure and the same 1,000,000 sample points (WorldClim pixel centres mapped to the nearest CHELSA pixel; 998,736 valid — CHELSA ends at 84°N), same variant (centre only). The CHELSA files are strip-organised GeoTIFFs (not COGs), so the global read streams every row once over `/vsicurl/` (~5.8 GB transfer, nothing stored). Values converted with the files' scale/offset. CHELSA's period overlaps the 2000–2012 outcome period, so it is not a temporally cleaner predictor. Random and HEALPix-blocked CV are reported for it. Optional rule. | `02`, `03c` (`VARIANTS=de`) |

## Testing

Every rule was smoke-tested end to end (`snakemake --config smoke=1`: 1° tiles in
-75…-73°E, 4…6°N; 3,000 points per class; 100 trees; outputs in
`data/clean_smoke/`, `results/smoke/`, `figures/smoke/`) before the full run. The
smoke test also exercised the MCD64A1 fallback (2 months, tile h10v08) and the
Earthdata `~/.netrc` login from a detached (`setsid`) process.

## Diagnostic 3 — transferability (training scope)

**Question.** Our step-2 model (Colombia only) gives an uncalibrated expected
area of 4.0 Mha against 10.3 Mha for the authors' pantropical map in the same
domain. On Colombian non-regrowth validation points, the authors' map scores
p > 0.5 on 52.6 % (mean 0.52) and ours on 6–8 % (mean 0.17–0.19). Leading
hypothesis: training scope. Diagnostic 3 trains **the same model specification**
on a Neotropical sample and tests transfer to Colombia (`02c`, `03d`). This is a
diagnostic extension, not part of the paper's method.

| # | Choice | Why | Code |
|---|---|---|---|
| D1 | **Region:** RESOLVE 2017 biomes 1–3, latitude −25…25°, longitude −120…−30° (the Americas). No realm filter; the two Nearctic pine-oak ecoregions of northern Mexico (biome 3) inside ±25° are included. | Neotropical part of the paper's domain, read literally | `02c` §1 |
| D2 | **Sub-regions:** explicit ecoregion → sub-region table in `02c` (Mesoamerica, Caribbean, Andes incl. inter-Andean dry valleys, Amazon incl. Guiana Shield, Atlantic forest incl. Brazilian Atlantic dry forests, Cerrado/Chaco edge = Chiquitano, Mato Grosso dry forests, Caatinga, Maranhão Babaçu). A seventh stratum, **"Chocó-Pacific and northern lowlands"** (Chocó-Darién, western Ecuador, Tumbes-Piura, Magdalena-Urabá, Catatumbo, Sinú, Maracaibo, Lara-Falcón, Apure-Villavicencio, Orinoco Delta), holds the ecoregions that fit none of the six named ones. Caatinga is biome 2 in RESOLVE 2017, so it is in the domain. Cerrado and Chaco proper are savanna biomes, outside the domain. | Every domain pixel must belong to a stratum; the names are taken from RESOLVE, not invented | `02c` §1 |
| D3 | **Cells:** HEALPix depth 5, NESTED, `healpix_geo` with `ellipsoid="WGS84"` (41,509 km² each, ~200 km). They are described on a 0.05° lattice by domain area, dominant sub-region and dominant biome. A cell is **eligible** if ≥ 25 % of its area is domain: 313 of the 457 cells that touch the domain. Stratum = dominant sub-region × dominant biome (14 strata with eligible cells). **48 cells** are drawn at random (seed 20261003): one per stratum, the remaining 34 allocated proportional to stratum domain area (largest-remainder, capped by eligible cells). Result: Amazon/1: 26; Cerrado-Chaco edge/2: 5; Atlantic forest/1: 3; Andes/1: 3; Mesoamerica/1: 2; one each for Cerrado-Chaco/1, Mesoamerica/2, Mesoamerica/3, Chocó-Pacific/1, Chocó-Pacific/2, Atlantic forest/2, Andes/2, Caribbean/1, Caribbean/2. Four cells touch Colombia (7486, 7613, 7650, 7657). Cell list with counts: `data/clean/neotropics_cells.parquet`; map: `figures/diag3_training_cells.png`. | Proportional allocation keeps the sample close to self-weighting by area, like the paper's area-uniform class-0 sampling. The minimum of one cell per stratum guarantees that every sub-region and biome is represented. A square-root allocation, which would weight small strata more, was not used. | `02c` §2 |
| D4 | **Labels and pixels** as in `02b` (30 m Hansen lattice, pixel-centre rule). Class 1 = domain ∩ Fagan `regrowth`. Class 0 = domain − forest 2000 − ESA CCI 2000 {150–153, 190, 200–202, 210} − all Fagan polygons. Fagan polygons are read remotely (`/vsicurl/`, every layer, per cell bbox + 0.1°). The ESA CCI 2000 exclusion is applied **at the sampled points**, not on the 30 m mask. This is equivalent, because the Bernoulli draw is independent per pixel, so dropping excluded pixels afterwards is a thinning of the same draw. | Avoids a wall-to-wall ESA mosaic; same definition | `02c` §5, §7 |
| D5 | **Sampling:** Bernoulli per 30 m pixel, inclusion ∝ exact WGS84 pixel area, with **per-cell rates**. Class 1: 3 × cap / regrowth pixels in the cell (regrowth area = geodesic area of the Fagan regrowth polygons whose representative point is in the cell). Class 0: 30 × cap / domain pixels in the cell. Then **cap = 1,000 points per class and cell** (a 1.3 × cap pre-cap before the remote soil reads), drop points with a missing predictor, and balance the classes by random reduction of the larger class. Split 75 % training pool / 25 % validation. Target: ≤ 48,000 per class. Cells with little regrowth contribute fewer class-1 points, so the realised total is reported in `neotropics_meta.json`. | Cap so that regrowth hotspots (Mesoamerica, Atlantic forest) do not dominate. Training draws of 39,523 (as step 2) leave headroom for the "excluding Colombia" pool. | `02c` §5, §11 |
| D6 | **Predictors** (2000 only; the Neotropical sample is used only for training and CV). Forest = tree cover 2000 ≥ 30 %. Forest density (1 km disk) and distance to forest (EDT, capped at 25 km) are computed on 1° tiles with a 0.25° margin, **the same code as `02b`**. Forest 2018 (02b, 2026-10-04: forest 2000 AND NOT loss 2001–2018, no gain) is not needed for the Neotropical sample. All 2018 inputs in `03d` come from the Colombian `samples.parquet` / `pred_grid.parquet` of `02b`. Soil: SoilGrids v2017 OCDENS/PHIHOX sl1–sl4, trapezoidal 0–30 cm, read remotely at the points. Bioclim: WorldClim v2.1 at the points with the **existing** PCA loadings (`data/clean/bioclim_pca_loadings.csv`, not refitted), PC1–PC4. Land cover: ESA CCI 2000 at the points, the same 31 → 11 table (`data/clean/landcover_classes.csv`). Biome: RESOLVE raster at 30 m. All coarse layers use nearest neighbour. | Same definitions as step 2 | `02c` §5–9 |
| D7 | **Check:** points inside the Colombian window are compared with the `02` layers (`data/clean/*.nc`). In the smoke test, 379 of 379 points had identical pc1–pc4, ocdens, phihox and lc_2000. The full-run check is in `neotropics_meta.json` (`consistency_check`). | Proves that the remote point reads reproduce the Colombian predictors | `02c` §10 |
| D8 | **Disk:** no wall-to-wall Hansen, SoilGrids or WorldClim. Hansen tree-cover 2000 tiles (10°) are downloaded per group (md5 from `x-goog-hash`) into `data/raw/hansen_tmp_neotropics/` and deleted once no remaining group needs them. The download aborts if free space would drop below 10 GB. The Colombian tiles in `data/raw/hansen/` are reused and kept. The 30 m lossyear and gain layers are not downloaded (not needed). | Disk constraint (~22 GB free) | `02c` §4, §6 |
| D9 | **Test 1** (`03d`). Three models with the step-2 specification (`make_rf`: 500 trees, `max_features=3`, `min_samples_leaf=1`, bootstrap, seed 20261003) and 39,523 balanced records: (i) Colombia, the exact step-2 draw; (ii) the Neotropical pool without points in GADM Colombia; (iii) the full Neotropical pool. Neotropical points on the same 30 m pixel as a Colombian validation point are removed. Reported: accuracy on the Colombian validation set (2000 inputs); mean p and share p > 0.5 on Colombian non-regrowth validation points (2000 and 2018 inputs); uncalibrated expected area and p > 0.5 area on the Colombian 1-in-100 grid (2018 inputs, exact WGS84 areas, as in `03`); Pearson r with the authors' map per pixel and at HEALPix depth 8 (cells > 50 % in the domain, as in `03`). The authors' row uses `auth_pct / 100` at the same pixels: NoData counts as 0 in areas and is excluded from point metrics, and its accuracy is computed on the validation points where the map has a value. | — | `03d` |
| D10 | **Test 2** (`03d`). 5-fold CV on the Neotropical training pool (incl. Colombia): stratified random folds vs `GroupKFold` by HEALPix cell (WGS84) at depths 6, 5, 4 (~100, ~200, ~400 km). Each fold trains on min(39,523, training-fold size) records. Reported: mean ± sd accuracy, balanced accuracy (blocked folds are not class-balanced), and mean p / share p > 0.5 on held-out non-regrowth points. Depth-5 groups coincide with the sampled cells. | Same scheme as step 3(b), at coarser block sizes | `03d` |

- **Diagnostic 3, unclassified Fagan patches:** 276 Neotropical gain polygons have a null `pred3class` (unclassified in Fagan et al. 2022, e.g. missing Sentinel-1 data). They are excluded from both classes. Colombia has none, so step 2/3 are unaffected.

## Area of applicability (AOA) — extension

Diagnostic extension, not part of the paper's method (`03e_area_of_applicability.py`).
It is a Python port of CAST `trainDI()` / `aoa()` (Meyer & Pebesma 2021, MEE,
doi:10.1111/2041-210X.13650). The source read was github.com/HannaMeyer/CAST at
commit `779e13f`, files `R/trainDI.R`, `R/aoa.R`, `R/aoa-helpers.R` and
`R/caret-helpers.R`, plus caret's `rf` `varImp` module.

| # | CAST | This port | Why |
|---|---|---|---|
| A1 | Factors become one dummy per training level (`caret::dummyVars`, not full rank). A level that is not in the training data scores 0 on every dummy. Each dummy gets the **full** weight of its variable. | We do the same for `lc_2000` (`lc_2015` at prediction) and `biome`. They are integer-coded in our data, so they are declared categorical explicitly. The importance weight is not split across the dummies. | Follows CAST exactly. Note that `scale()` gives rare levels (e.g. LC class 10, with 3 training records) very large scaled values, so those points have high DI. This is CAST behaviour. |
| A2 | `scale()`: training mean and sample sd. Weights multiply the scaled columns. Euclidean distance. | Same (`ddof=1`). Zero-variance columns are dropped with a message (CAST would stop). This can only happen in the smoke test. | — |
| A3 | Weights = `caret::varImp(model, scale = FALSE)`. For an `rf` classifier this is the mean over classes of the per-class `randomForest::importance()` columns: OOB permutation mean decrease in accuracy, divided by its standard error (randomForest's default `scale = TRUE`). Negative weights are set to 0. | sklearn has no OOB permutation importance. We use the permutation decrease in **balanced accuracy** (= mean over the two classes of the per-class accuracy decrease), not divided by a standard error, on a balanced hold-out of 10,000 records drawn from each model's training pool outside its training draw. We use 5 repeats and set negative values to 0. | This is the closest available equivalent. The SE scaling of randomForest is not reproduced, so the relative weights can differ. |
| A4 | Training DI: distance to the nearest training point in the CV-training part of the fold that holds the point out, divided by the mean, over training points, of each point's mean distance to all other training points. Folds come from the caret model. | Our models are not fitted by CV. Folds are 5-fold `GroupKFold` by HEALPix cell, depth **6** (~100 km; NESTED, `healpix_geo`, `ellipsoid="WGS84"`), on each model's own training draw. The same depth is used in step 3(b) and diagnostic 3. Sensitivity thresholds use depth 7, random 5-fold, and no CV (`useCV = FALSE`). | CAST recommends spatial folds when the model is assessed with spatial CV. |
| A5 | Threshold = min(Q3 + 1.5 IQR, max) of the training DI (`.di_threshold`, R `quantile` type 7). New data are inside the AOA if DI ≤ threshold. | Same (numpy `linear` quantiles = type 7). | — |
| A6 | Nearest neighbour by brute force (`FNN`, chunked). | Exact `scipy.spatial.cKDTree`, in batches of 500,000 query points. The mean pairwise distance is computed exactly in chunks. | Same result, faster. |
| A7 | — | Applied to the step-2 Colombia model and to the diagnostic-3 Neotropical models (without and with Colombia), all refitted exactly as in `03` / `03d`. New data: the 1-in-100 Colombian grid (2018/2015 inputs, exact WGS84 areas) and the Colombian validation points (2000 inputs; 2018 inputs for the non-regrowth comparison). The authors' map is only **split by our AOA**; their AOA cannot be computed without their training data. | — |
| A8 | Meyer & Pebesma relate DI to CV error. | Accuracy is reported by DI decile for (i) the Colombian validation points (for the Neotropical models, these are the transfer points) and (ii) CV predictions on each training draw with the same depth-6 folds: one fold model per fold, with the same settings and seed `SEED + fold`. Accuracy per class and **balanced accuracy** are reported next to accuracy. The share of regrowth points falls steeply with DI (97 % in the lowest decile, 16 % in the highest, for the Colombia model), so plain accuracy by DI mixes a class effect with a dissimilarity effect. | Class confounding found in the first full run |

## MapBiomas labels — independent regrowth labels and experiments E1–E4 (extension)

This is a diagnostic extension, not part of the paper's method (`01c`, `02d`, `03f`). Labels come from **MapBiomas Colombia Collection 3** annual land cover, 1985–2024 (CC BY 4.0; cite "Proyecto MapBiomas Colombia – Colección 3"). The data are public COGs in `gs://mapbiomas-public/initiatives/colombia/collection_3/coverage/`.

**Sources for the rule.** The rule is taken from *Apéndice – Módulo de pérdida de vegetación y vegetación secundaria en Colección 3*, v1.0, October 2025. Its live URL (`colombia.mapbiomas.org/wp-content/uploads/sites/3/2025/10/Apendice_modulo_perdida_de_vegetacion.pdf`) returns 404. The Wayback Machine snapshot of 2026-05-14 is stored in `data/raw/mapbiomas_docs/`. Two other documents were also checked: the C3 general appendix (`Apendice-1-Colombia-Coleccion-3.0.pdf`) and the C3 legend code sheet.

- **No ready-made product could be used.** MapBiomas shows Colombian secondary-vegetation and age maps on its platform, but the public bucket has no such layer under `initiatives/colombia/` (all 684 objects were listed on 2026-10-04). Only coverage, transitions (C1), fire and statistics are there, so the rule is re-implemented.
- **No Colombian accuracy figures are published.** The C3 appendix (§6) describes the validation design but gives no numbers. The "Análisis de exactitud" page describes MapBiomas Brazil. The vegetation-loss module says that its validation protocol is "in development".

| # | MapBiomas source | This replication | Why |
|---|---|---|---|
| MB1 | Module Table 1 groups. Anthropic: pastures, agriculture, forest plantation, palm oil, agriculture/pasture mosaic, urban infrastructure, mining, aquaculture. Natural: forest, mangrove, flooded forest, flooded non-forest formation, herbaceous formation, other non-forest formation. Not included: hypersaline tidal flat, other non-vegetated areas, river/lake/ocean, glacier, not observed. | C3 codes. Anthropic: 9, 14, 21, 24, 30, 31, 35, 74 (banana), 75 (solar). Natural: 1, 3, 5, 6, 49 (wooded sand vegetation), 10, 11, 12, 13, 50, 81, 82. Other: 0 (NoData), 22, 23, 25, 26, 27, 29 (rocky outcrop, listed under "Otro" in Table 2), 32, 33, 34, 68. | Codes 49, 50, 81, 82, 74 and 75 are not named in Table 1. They are assigned by their legend parent (forest formation, natural non-forest formation, agricultural area); solar parks go with infrastructure. Aquaculture (31) is anthropic, as in Table 1, although the legend places it under water. |
| MB2 | §2.2: loss in t = Natural in t−1 and t−2 (maps of the previous steps), then Anthropic in t and t+1. Regeneration in t = Anthropic in t−1 and t−2, then Natural in t, t+1 and t+2 (input). Changes that do not persist are reverted to the previous year's class. The analysis starts in 1987. A pixel that is "Other" in any year is class Other. | Implemented exactly, vectorised over points (`02d.trajectories`). Primary vs secondary loss follows the pixel's status. | — |
| MB3 | Table 2 says "three, five or seven years" of persistence and refers to biome-specific thresholds in a Table 3 that is **not in the document**. The §2.4 post-processing (cumulative layers; a 1–3 ha spatial filter on 2024 loss only) does not affect 2000–2020. | The persistence stated in §2.1–2.2: 2 years before and 2 after for loss, 2 before and 3 after for regeneration. The label windows below add the conservative condition that secondary vegetation persists to the end of the window. | Table 3 is unavailable. The shorter loss rule ends secondary status more readily, which keeps class 1 conservative. |
| MB4 | MapBiomas "regeneration" is any Anthropic → Natural change, including to grassland. | **Regrowth label (ours)** for window (S, O, P). Class 1 = Anthropic in S; a regeneration event in t ∈ [S+1, O]; secondary vegetation without loss from t to P; raw class forest (1, 3, 5, 6, 49) in every year O…P. Class 0 = never Other; non-forest and not urban/mining/aquaculture/solar in S; never forest in S…P; no regeneration event in S+1…P. Everything else is excluded. Windows: E1 (2000, 2012, 2016), the closest analogue of Fagan 2000–2012 natural regrowth, with persistence sensitivity (2000, 2012, 2014) and (2000, 2012, 2018); forward test (2000, 2008, 2012) and (2012, 2020, 2024). Changed on 2026-10-04 (Anne) from (2000, 2008, 2010) / (2010, 2018, 2020): the forward test now uses the same 4-year forest persistence as E1 (2 years admitted more ephemeral regrowth) and equal 8-year onset windows, and period-2 predictors are measured in 2012, when the period-1 outcomes are complete. | Fagan's class is tree regrowth, not any natural vegetation. Class 0 mirrors step 2: non-forest in 2000, with the analogue of the ESA urban/bare exclusions. The "never forest" condition makes it cleaner than step 2's class 0. |
| MB5 | MapBiomas grid 0.000269494585°, origin (−83.94379, 14.01480). | **Not resampled.** Each point (the 1-in-100 grid of `02b`, plus the step-2 sample points) takes the code of the MapBiomas pixel that contains it. The two grids are not aligned (Hansen 0.00025°), so a point and its MapBiomas pixel overlap partly. | Avoids resampling a categorical series; the label noise this adds is reported, not corrected. |
| MB6 | — | MapBiomas-labelled training and validation points are the `02b` systematic grid points (self-weighting by area), split pixel-wise 6 : 4.87 as in step 2. Balanced draws of 39,523 records train the models; balanced validation uses all class 1 plus as many random class-0 points. Prevalence = class-1 area / labelled area (exact WGS84 areas). CV = 5-fold random and HEALPix d6/d7 (WGS84) on a balanced pool of 151,980, as in step 3(b). | No new 30 m predictor extraction is needed for the 2000 inputs. |
| MB7 | — | Forward-test predictors for 2012: forest 2012 = Hansen v1.13 tree cover 2000 ≥ 30 % and no loss 2001–2012, with forest density and distance computed with the `02b` code and tiles (the 2000 values recomputed in the same pass match `02b` exactly); land cover = ESA CCI v2.0.7 **2012** (CEDA, same 11 classes). Calibration of the forward predictions uses only the period-1 prevalence (an oracle π₂ is shown for reference). | Same definitions as step 2, shifted to 2012. |
| MB8 | — | E3 history predictors from raw MapBiomas codes 1985–1999: years since last forest before 2000 (16 = not forest in 1985–1999) and the number of Anthropic years. This is an **extension**. | — |

## Sampling design (03h), added 2026-10-09 (extension)

Prompted by C. Cloud's reanalysis of Williams et al. in Brazil (2026, Matters Arising manuscript and
code: https://github.com/CannonCloud/tropical-regeneration-reanalysis). Not part of the paper's method. The step-2
model, predictors, validation set and prediction grid are unchanged; only the training sample differs.

| # | Topic | Choice | Why |
|---|---|---|---|
| S1 | Training prevalence | Total fixed at 39,523 records; regrowth share 50 % (the step-2 draw), 20 %, 5 % and π = 1.644 %. Each class drawn at random from the step-2 pool (`random_state = SEED`). Areas on the 1-in-100 grid (step-2 reference), uncalibrated and with the prior shift from the training share to π. | Cloud's prevalence sweep. Same grid as step 2, so the 50 % row reproduces 3.989 Mha. |
| S2 | Distance to regrowth | Haversine distance (R = 6371.0088 km, `BallTree`) from each point to the nearest regrowth point of the 1-in-100 grid (300 m spacing): Fagan class 1 for S3, MapBiomas 2000→2012 class 1 for S5. | The full 30 m regrowth mask is not needed at a 3 km threshold; the grid spacing (median regrowth-to-grid distance 0.34 km) is small against it. |
| S3 | Paired design | Regrowth as in step 2; non-regrowth 80 % drawn at random from pool points within 3 km of regrowth, 20 % from the whole pool. | Cloud's annulus (3 km) and 80/20 split. Random draw within the annulus rather than one-to-one nearest-neighbour matching. |
| S4 | Evaluation | (i) Step-2 random validation set; (ii) *local* validation: all validation regrowth + an equal random number of validation non-regrowth within 3 km; (iii) permutation importance (AUC drop, 3 repeats) on 20,000 balanced validation points, climate = bioclim PC1–4 + biome; (iv) areas uncalibrated, prior-shifted, and isotonic-calibrated on the random validation set weighted to π; (v) top-5 % overlap of the two designs on the prediction grid. | (ii) asks whether the model separates nearby sites; (iv) isotonic calibration does not assume pure prior shift, which paired sampling violates. |
| S5 | Forward test | As `03f` E2 (period 1 training pool, period 2 evaluation, 2012 inputs); the global row reproduces E2. Precision at top 1/5/10 % = share of the highest-ranked labelled points that regrew in 2012–2024. Isotonic calibration fitted on period-1 validation points weighted to π₁. | Only period-1 information is used for calibration. With MapBiomas labels, 82 % of non-regrowth pool points are already within 3 km of regrowth, so this is a weak contrast between designs. |

## Full-resolution prediction and archive (05a, 05b), added 2026-10-04

| # | Topic | Choice | Why |
|---|---|---|---|
| F1 | Area of our model | Predicted on **every** 30 m pixel of the prediction domain (05a), same model and predictors as step 2; areas are the reference from now on | The step-2 1-in-100 systematic sample over-estimated expected areas: our model 3.989 vs 3.929 Mha (+1.5 %), the authors' map 10.351 vs 10.281 Mha (+0.7 %). Predictions are identical at all 2,456,636 sample pixels (max difference 3e-8), so the gap is a property of the fixed sampling lattice (possibly interference with the 300 m ESA CCI grid), not of the model. Its stated standard errors (~0.06 %) assumed independent points and were too optimistic. Accuracies are unaffected (they use the validation points). |
| F2 | HEALPix aggregation | `healpix_resample.ConservativeResampler` (v2026.7.0), WGS84: each pixel binned by its centre, value × exact WGS84 pixel area; identical to our own healpix-geo sums (tested) | Library implementation of the conservative method already used in steps 1–3. |
| F3 | Finest archived depth | **15** (~200 m, ~29 M cells over Colombia); coarser levels 14 → 8 are exact parent sums | Depth 17 (~49 m) would be ~460 M cells (~18 GB for ten layers). Binning at 15 equals binning at 17 and summing children (pixel centres lie in exactly one cell). |
| F4 | Sample-based layers | MapBiomas labels and the area of applicability exist only at the sample points; archived at depth ≤ 12 (MapBiomas, with `n_points`) and 8 (AOA), with a `support` attribute | A depth-12 cell holds ~30 sample points; finer levels would be noise. |
| F5 | Block sizes for spatial CV | HEALPix depths 4–8 were chosen to span ~25–400 km, not from a measured autocorrelation range (variogram of residuals, Roberts et al. 2017) | Open: a variogram-based choice would justify the block size. |
| F6 | Archive metadata | GRID4EARTH layout; `dggs` convention attributes and CF grid-mapping variable from `healpix_connector.dggs_zarr` (v0.1.0); CF variable attributes; read back and decoded with `xdggs` | Interoperability (FAIR4RS I1). |
| F7 | Blocking grid for spatial CV (03g) | Step 3(b) repeated with HEALPix on the sphere, latitude/longitude squares of equal area, and the squares shifted half a block | WGS84 vs sphere: ≤ 0.3 percentage points; HEALPix vs lat/lon squares: ≤ 0.2 points at ~25–100 km. Block placement matters more than grid type at ~100 km (shifted squares 0.835 vs 0.854; fold SD 0.024–0.043). Step 3(b) numbers reproduced exactly. |