README
======

Bayesian Weibull Modeling of Post-Fire Burn Severity (dNBR) — Intermountain West
----------------------------------------------------------------------------------

This repository contains the code and supporting tables used to fit, evaluate,
and apply a Bayesian statistical model that predicts the probability
distribution of differenced Normalized Burn Ratio (dNBR) as a function of
vegetation type, vegetation density (EVC), annual aridity (AHM), and
topography, across 126 vegetation types in the Intermountain West.

The full methodology, results, and discussion are described in the associated
manuscript/thesis (Fernandez Leger, Zorzetto, and Cadol; New Mexico Institute
of Mining and Technology). Please refer to that paper for the statistical model, priors, and
interpretation of results.

This README documents what is in this repo, what is expected to be added
locally (large raster/sample data), and how the pieces fit together.


1. Repository contents
-----------------------

A. Model fitting (HPC / Alpine cluster, SLURM)
  - run_model.py
      Loads sampled per-fire pixel data (dNBR + predictors), applies EVT
      remapping/adjustments, EVC/FVC unit harmonization, and fits a Bayesian
      Weibull model (PyMC, BlackJAX NUTS sampler) for a single vegetation
      type (--veg_code) and a single topographic predictor configuration
      (--model_type, index into the model_list defined in the script).
      Saves posterior traces as NetCDF.
  - job_configs.txt
      One (EVT_code, model_type) pair per line; consumed by submit_jobs.sh
      to run one SLURM array task per line via run_model.py.
  - submit_jobs.sh
      SLURM array submission script that reads job_configs.txt and launches
      run_model.py for each (veg_code, model_type) combination.

B. Model comparison (ELPD / LOO / WAIC ranking)
  - compare.py
      For each EVT code, loads all fitted NetCDF traces matching that EVT,
      computes ArviZ LOO-based model comparison across the candidate
      topographic-predictor models, and writes a per-EVT comparison CSV.
  - submit_compare.sh
      SLURM array wrapper for compare.py (one array task per EVT code).

C. Posterior extraction/trimming
  - extract_posterior.py
      Strips NetCDF trace files down to posterior (+ observed_data /
      sample_stats) groups only, to reduce file size before archiving or
      downloading traces from the cluster.
  - submit_extract.sh
      SLURM array wrapper for extract_posterior.py.

D. Landscape-scale prediction / simulation
  - landscape_prediction.ipynb
      Notebook that takes fitted model parameters (either full posterior
      NetCDF traces or the summarized CSV of posterior means) and simulates
      a raster of predicted dNBR (mean and credible interval bounds) across
      a landscape extent, using current LandFire EVT/EVC and topographic
      rasters. Also supports a "historic_diagnostic" mode that compares
      simulated dNBR/severity against an observed fire scar. All run
      parameters (paths, mode, model, AHM scenario, severity classification
      method) are set in the CONFIG dictionary at the top of the notebook —
      update paths there before running locally.

E. Reference / lookup tables / datasets
  - EVT_naming_code.csv
      LandFire Existing Vegetation Type (EVT) code -> human-readable name
      crosswalk.
  - evt_code_matching14_studyarea_complete.csv
      Crosswalk from the pre-2016 LandFire EVT classification (VALUE14) to
      the current (post-2016 remap) EVT classification (VALUE24), used to
      harmonize fires sampled before 2020.
  - evt_adjustments.csv
      Manual EVT code adjustments/consolidations (sample_values -> change)
      applied after remapping, to merge or correct specific EVT codes.
  - valid_vegtype_summary.csv
      Per-EVT summary table used throughout the pipeline: sample size,
      severity class fractions, EVT/fuel/physiognomic group labels,
      ecoregion, and the ranked topographic-predictor models (rank0-rank4)
      from the model comparison step. This is the master list of the 126
      EVTs included in the study and is a required input to
      landscape_prediction.ipynb.
  - evt_model_params_means.csv
      Posterior mean parameter values (Weibull shape/scale regression
      coefficients) per EVT for the selected topographic model
      (variability1000_HLI), for use in "csv_means" mode of
      landscape_prediction.ipynb (i.e., simulation without needing the full
      NetCDF posterior traces).
  - evt_lowprob_draws_mean_p025_p975.csv
      Per-EVT ordinal-logistic severity-classification parameters (fusion of
      the Weibull dNBR distribution with low/moderate/high severity class
      probabilities) at the posterior mean, 2.5th, and 97.5th percentile
      draws, used by the "ordlogit_fusion" severity classification option in
      landscape_prediction.ipynb.
  - Posterior trace NetCDF files (*_dnbr_sev__traces.nc)
      Full PyMC/ArviZ posterior output from run_model.py. Large; only needed
      if running landscape_prediction.ipynb in "netcdf" parameter_source
      mode instead of "csv_means".
  - dnbr_samples
      Per-fire CSVs of sampled burned pixels (dNBR, RdNBR, MTBS severity
      class, EVT, EVC, AHM, and topographic predictor values), plus MTBS
      burn-boundary shapefiles (used to extract each fire's dNBR offset).
      This is the raw training data consumed by run_model.py.

2. Data NOT included in this repository
----------------------------------------

The following inputs rasters are required to reproduce the full pipeline but are not
included here due to file size. Paths referenced in the scripts/notebook
(e.g., /projects/afernandez@xsede.org/..., /scratch/alpine/..., D:\...) will
need to be updated to wherever these are stored locally or on your own
cluster.

  - Topographic predictor rasters
      30 m DEM-derived rasters (Slope, HLI, TWI, VRM, TPI, Dissection Index,
      elevation variability at 300 m / 1 km / 2 km, etc.) used both in
      training-data sampling and in landscape_prediction.ipynb.
  - LandFire EVT / EVC (or FVC) rasters
      By year, used both for training-data sampling (historic) and for
      landscape-scale prediction (current conditions).


3. Typical workflow
---------------------

  1. Sample dNBR and predictors per fire -> dnbr_samples/*.csv (external to
     this repo).
  2. Fit Bayesian Weibull models per (EVT, topographic model) combination on
     the cluster: submit_jobs.sh -> run_model.py, reading job_configs.txt.
  3. Rank candidate topographic models per EVT: submit_compare.sh ->
     compare.py -> model_comparison_.csv, feeding the rank0-rank4
     columns in valid_vegtype_summary.csv.
  4. (Optional) Trim posterior NetCDF files for archiving/download:
     submit_extract.sh -> extract_posterior.py.
  5. Summarize posterior means across EVTs into evt_model_params_means.csv
     and severity-fusion parameters into
     evt_lowprob_draws_mean_p025_p975.csv (summary step not included as a
     standalone script here; see landscape_prediction.ipynb for how these
     tables are consumed).
  6. Run landscape-scale prediction/diagnostics: landscape_prediction.ipynb.


4. Environment
----------------

Model fitting and comparison were run on the CU Boulder Alpine HPC cluster
under a conda environment named mc_conda_env, with the following key
packages:
  - pymc (5.21.1), pymc.sampling.jax / blackjax (NUTS sampler)
  - arviz
  - numpy, pandas, scipy
  - fiona (MTBS burn-boundary shapefile reading)
  - nest_asyncio

landscape_prediction.ipynb additionally requires:
  - rasterio
  - matplotlib
  - (historic_diagnostic mode only) seaborn, scikit-learn

Note: run_model.py sets PYTENSOR_FLAGS / a custom compiledir to avoid shared
compiledir conflicts between concurrent SLURM array tasks on Alpine — adjust
this if running outside that environment.


5. Vegetation/EVT code notes
------------------------------

EVT codes follow the LandFire Existing Vegetation Type classification.
Because the classification changed in the 2016 LandFire remap, fires sampled
before 2020 are remapped from the pre-2016 (VALUE14) to the post-2016
(VALUE24) system using evt_code_matching14_studyarea_complete.csv, and
evt_adjustments.csv applies additional manual corrections/consolidations.
EVT_code = -9999 or 7944 and non-wildland classes (developed, agricultural,
ruderal, barren/non-vegetated) are excluded from model fitting.


6. Citation
-------------

If you use this code or the associated tables, please cite the manuscript
(citation to be added once published) or contact the corresponding author,
Daniel Cadol (daniel.cadol@nmt.edu), Department of Earth and Environmental
Science, New Mexico Institute of Mining and Technology.