{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fire"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Fire tools provide per-cell raster functions for burn severity mapping, fire behavior modeling, and drought indexing.\n",
    "\n",
    "- [dNBR](#dNBR): Differenced Normalized Burn Ratio (pre minus post NBR).\n",
    "- [RdNBR](#RdNBR): Relative dNBR, normalized by pre-fire vegetation density.\n",
    "- [Burn Severity Classification](#Burn-Severity-Classification): USGS 7-class severity from dNBR.\n",
    "- [Fireline Intensity](#Fireline-Intensity): Byram's fireline intensity (kW/m).\n",
    "- [Flame Length](#Flame-Length): Flame length from intensity (m).\n",
    "- [Rate of Spread](#Rate-of-Spread): Simplified Rothermel with Anderson 13 fuel models (m/min).\n",
    "- [KBDI](#KBDI): Keetch-Byram Drought Index, single time-step update."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Importing packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import xarray as xr\n",
    "\n",
    "from datashader.transfer_functions import shade, stack, Images\n",
    "\n",
    "from xrspatial.fire import (\n",
    "    dnbr, rdnbr, burn_severity_class,\n",
    "    fireline_intensity, flame_length,\n",
    "    rate_of_spread, kbdi,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate synthetic data\n",
    "\n",
    "We create a 200x200 landscape with a simulated burn scar. Pre-fire NBR is higher where vegetation is denser; after the fire, NBR drops inside an elliptical burn perimeter.\n",
    "\n",
    "In a real workflow you would compute NBR from satellite imagery using `xrspatial.multispectral.nbr`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "H, W = 200, 200\n",
    "rng = np.random.default_rng(42)\n",
    "\n",
    "ys = np.linspace(H - 1, 0, H)\n",
    "xs = np.linspace(0, W - 1, W)\n",
    "\n",
    "def make_da(data, name):\n",
    "    return xr.DataArray(data.astype(np.float32), dims=['y', 'x'],\n",
    "                        coords={'y': ys, 'x': xs}, name=name)\n",
    "\n",
    "yy, xx = np.meshgrid(np.linspace(0, 1, H), np.linspace(0, 1, W), indexing='ij')\n",
    "veg = 0.3 + 0.3 * np.sin(2 * np.pi * yy) * np.cos(np.pi * xx)\n",
    "pre_nbr = np.clip(veg + rng.normal(0, 0.03, (H, W)), 0.05, 0.85)\n",
    "\n",
    "dist = np.sqrt(((yy - 0.5) / 0.25) ** 2 + ((xx - 0.5) / 0.35) ** 2)\n",
    "burn_mask = dist < 1.0\n",
    "burn_intensity = np.clip(1.0 - dist, 0, 1)\n",
    "\n",
    "post_nbr = pre_nbr.copy()\n",
    "post_nbr[burn_mask] -= burn_intensity[burn_mask] * (0.3 + rng.uniform(0, 0.3, burn_mask.sum()))\n",
    "post_nbr = np.clip(post_nbr, -0.5, 0.85)\n",
    "\n",
    "pre_nbr_agg = make_da(pre_nbr, 'pre_nbr')\n",
    "post_nbr_agg = make_da(post_nbr, 'post_nbr')\n",
    "\n",
    "pre_img = shade(pre_nbr_agg, cmap=['brown', 'yellow', 'green'], how='linear')\n",
    "pre_img.name = 'Pre-fire NBR'\n",
    "post_img = shade(post_nbr_agg, cmap=['brown', 'yellow', 'green'], how='linear')\n",
    "post_img.name = 'Post-fire NBR'\n",
    "imgs = Images(pre_img, post_img)\n",
    "imgs.num_cols = 2\n",
    "imgs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### dNBR\n",
    "\n",
    "The differenced Normalized Burn Ratio is `pre_NBR - post_NBR`. Positive values mean vegetation loss; negative values mean regrowth. USGS and BAER teams use dNBR as input to the severity classification thresholds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dnbr_agg = dnbr(pre_nbr_agg, post_nbr_agg)\n",
    "\n",
    "print(f\"dNBR range: {float(dnbr_agg.min()):.3f} to {float(dnbr_agg.max()):.3f}\")\n",
    "shade(dnbr_agg, cmap=['green', 'lightyellow', 'orange', 'red', 'darkred'], how='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RdNBR\n",
    "\n",
    "Relative dNBR normalizes severity by pre-fire vegetation density: `dNBR / sqrt(abs(pre_NBR / 1000))`. This lets you compare burn severity across vegetation types. Pixels where pre-fire NBR is near zero are set to NaN."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rdnbr_agg = rdnbr(dnbr_agg, pre_nbr_agg)\n",
    "\n",
    "print(f\"RdNBR range: {float(np.nanmin(rdnbr_agg.data)):.3f} to \"\n",
    "      f\"{float(np.nanmax(rdnbr_agg.data)):.3f}\")\n",
    "shade(rdnbr_agg, cmap=['green', 'lightyellow', 'orange', 'red', 'darkred'], how='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Burn Severity Classification\n",
    "\n",
    "`burn_severity_class` bins dNBR into the standard USGS 7-class scheme (int8 output, 0 = nodata). This function accepts Datasets via `@supports_dataset`.\n",
    "\n",
    "| Class | Label | dNBR range |\n",
    "|-------|-------|------------|\n",
    "| 1 | Enhanced regrowth (high) | < -0.251 |\n",
    "| 2 | Enhanced regrowth (low) | -0.251 to -0.101 |\n",
    "| 3 | Unburned | -0.101 to 0.099 |\n",
    "| 4 | Low severity | 0.099 to 0.269 |\n",
    "| 5 | Moderate-low severity | 0.269 to 0.439 |\n",
    "| 6 | Moderate-high severity | 0.439 to 0.659 |\n",
    "| 7 | High severity | >= 0.659 |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "severity = burn_severity_class(dnbr_agg)\n",
    "\n",
    "severity_float = severity.astype(np.float32)\n",
    "severity_float.values = np.where(severity_float.values == 0, np.nan, severity_float.values)\n",
    "shade(severity_float,\n",
    "      cmap=['darkgreen', 'green', 'lightgreen', 'yellow', 'orange', 'red', 'darkred'],\n",
    "      how='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fireline Intensity\n",
    "\n",
    "Byram's fireline intensity: `I = H * w * R` where *H* is heat content (kJ/kg), *w* is fuel consumed (kg/m²), and *R* is spread rate (m/s). Output is kW/m. Fires below ~350 kW/m can be attacked by hand crews; above ~4,000 kW/m they typically need indirect attack or aerial resources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fuel = make_da((veg * 3.0 + rng.uniform(0, 0.5, (H, W))).astype(np.float32), 'fuel')\n",
    "spread = make_da((0.02 + 0.03 * rng.uniform(0, 1, (H, W))).astype(np.float32), 'spread')\n",
    "\n",
    "intensity_agg = fireline_intensity(fuel, spread, heat_content=18000)\n",
    "\n",
    "print(f\"Intensity range: {float(intensity_agg.min()):.1f} to {float(intensity_agg.max()):.1f} kW/m\")\n",
    "shade(intensity_agg, cmap=['lightyellow', 'orange', 'red', 'darkred'], how='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flame Length\n",
    "\n",
    "Flame length from fireline intensity: `L = 0.0775 * I^0.46`. Zero or negative intensity gives zero flame length. Accepts Datasets via `@supports_dataset`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fl_agg = flame_length(intensity_agg)\n",
    "\n",
    "print(f\"Flame length range: {float(fl_agg.min()):.2f} to {float(fl_agg.max()):.2f} m\")\n",
    "shade(fl_agg, cmap=['lightyellow', 'orange', 'red'], how='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rate of Spread\n",
    "\n",
    "`rate_of_spread` uses a simplified Rothermel (1972) model with the Anderson 13 fuel model table. Inputs are slope (degrees), mid-flame wind speed (km/h), and dead fuel moisture (fraction 0-1). The `fuel_model` parameter (1-13) selects fuel bed properties.\n",
    "\n",
    "Below, slope increases from bottom to top and wind increases from left to right, so spread rate is highest in the top-right corner."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope_agg = make_da((5.0 + 20.0 * yy).astype(np.float32), 'slope')\n",
    "wind_agg = make_da((5.0 + 15.0 * xx).astype(np.float32), 'wind')\n",
    "moisture_agg = make_da(np.full((H, W), 0.06, dtype=np.float32), 'moisture')\n",
    "\n",
    "ros_agg = rate_of_spread(slope_agg, wind_agg, moisture_agg, fuel_model=1)\n",
    "\n",
    "print(f\"Rate of spread: {float(ros_agg.min()):.2f} to {float(ros_agg.max()):.2f} m/min\")\n",
    "shade(ros_agg, cmap=['lightyellow', 'orange', 'red', 'darkred'], how='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparing fuel models with the same inputs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for fm, name in [(1, 'Short grass'), (3, 'Tall grass'), (4, 'Chaparral'), (8, 'Timber litter')]:\n",
    "    r = rate_of_spread(slope_agg, wind_agg, moisture_agg, fuel_model=fm)\n",
    "    print(f\"  Model {fm:2d} ({name:15s}): {float(r.min()):8.2f} to {float(r.max()):8.2f} m/min\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### KBDI\n",
    "\n",
    "The Keetch-Byram Drought Index tracks cumulative soil moisture deficit (0-800 mm). It gets updated daily from max temperature (Celsius) and precipitation (mm). `annual_precip` is a scalar for mean annual rainfall.\n",
    "\n",
    "Below we start from KBDI = 300 (moderate drought), run 30 hot dry days, drop 40 mm of rain, then continue."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "current = make_da(np.full((H, W), 300.0, dtype=np.float32), 'kbdi')\n",
    "hot = make_da(np.full((H, W), 35.0, dtype=np.float32), 'temp')\n",
    "no_rain = make_da(np.zeros((H, W), dtype=np.float32), 'precip')\n",
    "\n",
    "history = [float(current.mean())]\n",
    "for _ in range(30):\n",
    "    current = kbdi(current, hot, no_rain, annual_precip=1200.0)\n",
    "    history.append(float(current.mean()))\n",
    "\n",
    "rain = make_da(np.full((H, W), 40.0, dtype=np.float32), 'precip')\n",
    "current = kbdi(current, hot, rain, annual_precip=1200.0)\n",
    "history.append(float(current.mean()))\n",
    "\n",
    "for _ in range(5):\n",
    "    current = kbdi(current, hot, no_rain, annual_precip=1200.0)\n",
    "    history.append(float(current.mean()))\n",
    "\n",
    "print(f\"Day  0: {history[0]:.1f}\")\n",
    "print(f\"Day 30: {history[30]:.1f}  (pre-rain)\")\n",
    "print(f\"Day 31: {history[31]:.1f}  (post-rain)\")\n",
    "print(f\"Day 36: {history[-1]:.1f}  (5 days after rain)\")\n",
    "\n",
    "shade(current, cmap=['green', 'yellow', 'orange', 'red'], how='linear')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References\n",
    "\n",
    "- Key, C.H. and Benson, N.C. (2006). Landscape Assessment. In: *FIREMON*, USDA Forest Service Gen. Tech. Rep. RMRS-GTR-164-CD.\n",
    "- Rothermel, R.C. (1972). A mathematical model for predicting fire spread in wildland fuels. USDA Forest Service Res. Pap. INT-115.\n",
    "- Anderson, H.E. (1982). Aids to determining fuel models for estimating fire behavior. USDA Forest Service Gen. Tech. Rep. INT-122.\n",
    "- Keetch, J.J. and Byram, G.M. (1968). A drought index for forest fire control. USDA Forest Service Res. Pap. SE-38.\n",
    "- USGS Burn Severity Portal: https://burnseverity.cr.usgs.gov/"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}