"""Deforestation detection in Python with Hansen/UMD Global Forest Change (v1.13, 2001-2025),
correct pixel areas, a Landsat/Sentinel-2 before-and-after from the Planetary Computer, an NBR
check against the GFC loss year and a MODIS burned-area cross-check, for central Rondonia, Brazil.

TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/deforestation-detection-python/
Data: Hansen/UMD/Google/USGS/NASA (CC BY 4.0); Copernicus Sentinel-2, USGS Landsat and NASA MODIS
via the Microsoft Planetary Computer.
"""
import sys
import warnings
from importlib.metadata import version
from pathlib import Path

import numpy as np
import pandas as pd
import rasterio
from rasterio.enums import Resampling
from rasterio.transform import from_origin
from rasterio.vrt import WarpedVRT
from rasterio.windows import from_bounds
from pyproj import Geod
from shapely.geometry import box, shape
import pystac_client
import planetary_computer
import matplotlib
matplotlib.use("Agg")   # render to files only; no GUI backend
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, BoundaryNorm

warnings.filterwarnings("ignore")
plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9})
DATA = Path("data"); DATA.mkdir(exist_ok=True)
OUT = Path("output"); OUT.mkdir(exist_ok=True)

VERSION = "GFC-2025-v1.13"
GFC = f"https://storage.googleapis.com/earthenginepartners-hansen/{VERSION}/Hansen_{VERSION}_{{layer}}_10S_070W.tif"
AOI = (-64.0, -11.0, -63.0, -10.0)   # lon/lat: central Rondonia around Buritis and Campo Novo de Rondonia
RES = 0.00025                        # GFC pixel size in degrees (0.9 arc-seconds)
FOREST_CANOPY = 30                   # % canopy cover in 2000 that counts as forest (the GFW default)
STAC_URL = "https://planetarycomputer.microsoft.com/api/stac/v1"


def section(title):
    print(f"\n=== {title} ===")


section("versions")
print(f"python {sys.version.split()[0]} | rasterio {version('rasterio')} | numpy {version('numpy')} | "
      f"pyproj {version('pyproj')} | pystac-client {version('pystac-client')} | "
      f"planetary-computer {version('planetary-computer')}")

# ---------------------------------------------------------------- Step 1: read GFC layers for the AOI
section("Global Forest Change layers")


def read_gfc(layer):
    """Window of one 10x10 degree GFC tile, cached locally as a tiled GeoTIFF."""
    path = DATA / f"{VERSION}_{layer}_rondonia.tif"
    if not path.exists():
        with rasterio.open("/vsicurl/" + GFC.format(layer=layer)) as src:
            win = from_bounds(*AOI, transform=src.transform)
            arr = src.read(1, window=win)
            prof = src.profile | {"width": arr.shape[1], "height": arr.shape[0], "transform": src.window_transform(win),
                                  "tiled": True, "blockxsize": 512, "blockysize": 512, "compress": "deflate"}
        with rasterio.open(path, "w", **prof) as dst:
            dst.write(arr, 1)
    with rasterio.open(path) as src:
        return src.read(1), src.transform


treecover, transform = read_gfc("treecover2000")
lossyear, _ = read_gfc("lossyear")
lossyear = lossyear.astype(np.int16)                 # uint8 overflows when we add 2000
datamask, _ = read_gfc("datamask")
H, W = treecover.shape
print(f"{VERSION}: {W} x {H} pixels of {RES} deg over lon {AOI[0]}..{AOI[2]}, lat {AOI[1]}..{AOI[3]}")
print(f"lossyear values: {lossyear.min()}..{lossyear.max()} (0 = no loss, 1..25 = 2001..2025)")
print(f"datamask: land {(datamask == 1).mean():.2%}, water {(datamask == 2).mean():.2%}, no data {(datamask == 0).mean():.2%}")

# ---------------------------------------------------------------- Step 2: true pixel areas
section("pixel area")
geod = Geod(ellps="WGS84")
lats = transform.f + transform.e * np.arange(H + 1)            # row edges, north to south
row_area = np.array([abs(geod.polygon_area_perimeter([0, RES, RES, 0], [lats[i], lats[i], lats[i + 1], lats[i + 1]])[0])
                     for i in range(H)])                         # m2 of one pixel in each row
w_m = geod.inv(0, -10.5, RES, -10.5)[2]
h_m = geod.inv(0, -10.5 + RES / 2, 0, -10.5 - RES / 2)[2]
print(f"at 10.5 S one pixel is {w_m:.2f} m wide x {h_m:.2f} m tall = {row_area[H // 2]:.1f} m2 ({row_area[H // 2] / 1e4:.4f} ha)")
print(f"row areas range {row_area.min():.1f}..{row_area.max():.1f} m2 across the AOI")
for lat in (0, 10.5, 30, 50, 60, 70):
    a = abs(geod.polygon_area_perimeter([0, RES, RES, 0], [lat, lat, lat + RES, lat + RES])[0])
    print(f"  pixel area at {lat:4.1f} deg: {a:7.1f} m2  ({a / 900:.2f} x a 30 m square)")
area_ha = np.repeat(row_area[:, None] / 1e4, W, axis=1).astype(np.float32)   # hectares per pixel


def ha(mask):
    return float(area_ha[mask].sum())


land = datamask == 1
forest2000 = land & (treecover >= FOREST_CANOPY)
print(f"AOI land area {ha(land):,.0f} ha; forest (>= {FOREST_CANOPY}% canopy) in 2000 {ha(forest2000):,.0f} ha "
      f"({ha(forest2000) / ha(land):.1%} of land)")
naive = forest2000.sum() * 0.09
print(f"pixel count x 0.09 ha would give {naive:,.0f} ha, {naive / ha(forest2000) - 1:+.1%} too much")

# ---------------------------------------------------------------- Step 3: loss by year
section("forest loss by year (canopy >= 30% in 2000)")
years = np.arange(2001, 2026)
loss = forest2000 & (lossyear > 0)
by_year = np.bincount(lossyear[forest2000].ravel(), weights=area_ha[forest2000].ravel(), minlength=26)[1:26]
table = pd.DataFrame({"year": years, "loss_ha": by_year})
table["pct_of_2000_forest"] = 100 * table.loss_ha / ha(forest2000)
table["moving_avg_3yr"] = table.loss_ha.rolling(3, center=True).mean()
print(table.to_string(index=False, float_format=lambda v: f"{v:,.1f}"))
table.to_csv(OUT / "loss_by_year.csv", index=False)
print(f"total 2001-2025: {ha(loss):,.0f} ha = {ha(loss) / ha(forest2000):.1%} of the 2000 forest; "
      f"peak {years[by_year.argmax()]} ({by_year.max():,.0f} ha)")
print(f"2001-2010 mean {by_year[:10].mean():,.0f} ha/yr; 2011-2020 mean {by_year[10:20].mean():,.0f} ha/yr; "
      f"2021-2025 mean {by_year[20:].mean():,.0f} ha/yr")
for thr in (10, 30, 50, 75):
    f = land & (treecover >= thr)
    print(f"canopy threshold {thr:2d}%: forest 2000 {ha(f):,.0f} ha, loss 2001-2025 {ha(f & (lossyear > 0)):,.0f} ha")
print(f"loss pixels outside the 30% forest mask: {ha(land & (lossyear > 0) & ~forest2000):,.0f} ha")

# figure 1: tree cover 2000 and loss year
ly_cmap = plt.get_cmap("Spectral_r", 25)
fig, axes = plt.subplots(1, 2, figsize=(12.5, 5.6), constrained_layout=True)
ext = (AOI[0], AOI[2], AOI[1], AOI[3])
im0 = axes[0].imshow(np.where(land, treecover, np.nan), cmap="YlGn", vmin=0, vmax=100, extent=ext,
                     interpolation="nearest")
fig.colorbar(im0, ax=axes[0], shrink=0.75, label="tree canopy cover in 2000 (%)")
axes[0].set_title("Tree cover 2000", fontsize=10)
axes[1].imshow(np.where(forest2000 & (lossyear == 0), 1, np.nan), cmap=ListedColormap(["#c7dfb7"]), extent=ext,
               interpolation="nearest")
im1 = axes[1].imshow(np.where(loss, lossyear + 2000, np.nan), cmap=ly_cmap, vmin=2000.5, vmax=2025.5, extent=ext,
                     interpolation="nearest")
fig.colorbar(im1, ax=axes[1], shrink=0.75, label="year of forest loss", ticks=[2001, 2005, 2010, 2015, 2020, 2025])
axes[1].set_title("Forest loss 2001-2025 (pale green = forest remaining)", fontsize=10)
for ax in axes:
    ax.set_xlabel("longitude"); ax.set_ylabel("latitude")
fig.savefig(OUT / "01-treecover-lossyear.png")
plt.close(fig)

# figure 2: annual loss
fig, ax = plt.subplots(figsize=(10, 4.2), constrained_layout=True)
ax.bar(years, table.loss_ha / 1000, color=[ly_cmap(i) for i in range(25)], edgecolor="none")
ax.plot(years, table.moving_avg_3yr / 1000, color="#333333", lw=1.6, label="3-year moving average")
ax.axvspan(2000.5, 2010.5, color="#999999", alpha=0.12, lw=0)
ax.text(2005.5, ax.get_ylim()[1] * 0.93, "2001-2010: original method,\nnot yet reprocessed", ha="center", va="top",
        fontsize=8, color="#555555")
ax.text(2018, ax.get_ylim()[1] * 0.93, "2011 onward: reprocessed with Landsat 8 (2013+) and 9 (2022+)",
        ha="center", va="top", fontsize=8, color="#555555")
ax.set_ylabel("forest loss (thousand ha)")
ax.set_xlim(2000.3, 2025.7)
ax.legend(frameon=False, loc="center right")
ax.set_title(f"Annual tree cover loss, central Rondonia (canopy >= {FOREST_CANOPY}% in 2000), GFC v1.13", fontsize=10)
fig.savefig(OUT / "02-loss-by-year.png")
plt.close(fig)

# figure 3: pixel area against latitude
lat_axis = np.arange(0, 80.1, 0.5)
pix = np.array([abs(geod.polygon_area_perimeter([0, RES, RES, 0], [l, l, l + RES, l + RES])[0]) for l in lat_axis])
fig, ax = plt.subplots(figsize=(7.5, 3.8), constrained_layout=True)
ax.plot(lat_axis, pix, color="#3288bd", lw=2, label="true area of a 0.00025 deg GFC pixel")
ax.axhline(900, color="#d53e4f", lw=1.4, ls="--", label="'30 m pixel' assumption (900 m2)")
ax.axvspan(10, 11, color="#fdae61", alpha=0.5, lw=0, label="this study area")
ax.set_xlabel("latitude (degrees north or south)"); ax.set_ylabel("pixel area (m2)")
ax.set_xlim(0, 80); ax.set_ylim(0, 1000)
ax.legend(frameon=False, loc="lower left")
fig.savefig(OUT / "03-pixel-area-latitude.png")
plt.close(fig)

# ---------------------------------------------------------------- Step 4: pick a zoom area of recent loss
section("zoom area")
Z = 480                                              # 0.12 degrees, about 13 km
recent = (forest2000 & (lossyear >= 18) & (lossyear <= 24)).astype(np.float32)
blocks = recent[: H // Z * Z, : W // Z * Z].reshape(H // Z, Z, W // Z, Z).sum(axis=(1, 3))
bi, bj = np.unravel_index(blocks.argmax(), blocks.shape)
r0, c0 = bi * Z, bj * Z
zt = transform * rasterio.Affine.translation(c0, r0)
zb = (zt.c, zt.f - Z * RES, zt.c + Z * RES, zt.f)
zsl = (slice(r0, r0 + Z), slice(c0, c0 + Z))
zly, zf, zarea = lossyear[zsl], forest2000[zsl], area_ha[zsl]
print(f"zoom: lon {zb[0]:.3f}..{zb[2]:.3f}, lat {zb[1]:.3f}..{zb[3]:.3f} ({Z} x {Z} GFC pixels)")
print(f"forest 2000 {zarea[zf].sum():,.0f} ha; loss 2001-2017 {zarea[zf & (zly >= 1) & (zly <= 17)].sum():,.0f} ha; "
      f"loss 2018-2024 {zarea[zf & (zly >= 18) & (zly <= 24)].sum():,.0f} ha; loss 2025 {zarea[zf & (zly == 25)].sum():,.0f} ha")

catalog = pystac_client.Client.open(STAC_URL, modifier=planetary_computer.sign_inplace)


def read_to_grid(href, bounds, shape, resampling=Resampling.average, dtype="float32"):
    """Read a COG band straight onto our lat/lon grid (GFC pixels) with GDAL's warper."""
    t = from_origin(bounds[0], bounds[3], RES, RES)
    with rasterio.open(href) as src, WarpedVRT(src, crs="EPSG:4326", transform=t, width=shape[1], height=shape[0],
                                              resampling=resampling) as vrt:
        return vrt.read(1).astype(dtype)


def s2_reflectance(item, band, bounds, shape):
    dn = read_to_grid(item.assets[band].href, bounds, shape)
    offset = 1000 if item.properties.get("s2:processing_baseline", "00.00") >= "04.00" else 0
    return np.where(dn == 0, np.nan, (dn - offset) / 10000)


def best_s2(year, bounds, months="06-01/08-31", n=3):
    items = catalog.search(collections=["sentinel-2-l2a"], bbox=bounds, datetime=f"{year}-{months.replace('/', f'/{year}-')}",
                           query={"eo:cloud_cover": {"lt": 20}}).item_collection()
    items = sorted((it for it in items if shape(it.geometry).contains(box(*bounds))),
                   key=lambda it: it.properties["eo:cloud_cover"])
    return items[:n]


def nbr_composite(items, bounds, shape):
    """Median NBR of the clearest scenes, with cloud, shadow, cirrus and nodata masked by SCL."""
    stack = []
    for it in items:
        scl = read_to_grid(it.assets["SCL"].href, bounds, shape, Resampling.nearest)
        nir, swir = s2_reflectance(it, "B8A", bounds, shape), s2_reflectance(it, "B12", bounds, shape)
        nbr = (nir - swir) / (nir + swir)
        nbr[np.isin(scl, [0, 1, 3, 8, 9, 10])] = np.nan
        stack.append(nbr)
    return np.nanmedian(np.stack(stack), axis=0)


def stretch(rgb, lo=2, hi=98):
    p1, p2 = np.nanpercentile(rgb, [lo, hi])
    return np.clip((rgb - p1) / (p2 - p1), 0, 1)


# ---------------------------------------------------------------- Step 5: before and after imagery
section("before/after imagery")
landsat = catalog.search(collections=["landsat-c2-l2"], bbox=zb, datetime="2000-06-01/2000-08-31",
                         query={"platform": {"eq": "landsat-7"}, "eo:cloud_cover": {"lt": 5}}).item_collection()
ls = sorted(landsat, key=lambda it: (it.properties["eo:cloud_cover"], it.datetime))[0]
print(f"Landsat 7 before: {ls.id} cloud {ls.properties['eo:cloud_cover']}%")
ls_rgb = np.dstack([read_to_grid(ls.assets[b].href, zb, (Z, Z)) * 0.0000275 - 0.2 for b in ("red", "green", "blue")])
s2_before, s2_after = best_s2(2017, zb, n=1)[0], best_s2(2025, zb, n=1)[0]
rgbs = {}
for label, it in [("before", s2_before), ("after", s2_after)]:
    print(f"Sentinel-2 {label}: {it.id} cloud {it.properties['eo:cloud_cover']:.1f}% "
          f"baseline {it.properties.get('s2:processing_baseline')}")
    rgbs[label] = np.dstack([s2_reflectance(it, b, zb, (Z, Z)) for b in ("B04", "B03", "B02")])
zext = (zb[0], zb[2], zb[1], zb[3])
fig, axes = plt.subplots(1, 4, figsize=(16, 4.5), constrained_layout=True)
panels = [(stretch(ls_rgb), f"Landsat 7 true colour, {ls.datetime:%d %b %Y}"),
          (stretch(rgbs["before"]), f"Sentinel-2 true colour, {s2_before.datetime:%d %b %Y}"),
          (stretch(rgbs["after"]), f"Sentinel-2 true colour, {s2_after.datetime:%d %b %Y}")]
for ax, (img, title) in zip(axes, panels):
    ax.imshow(img, extent=zext); ax.set_title(title, fontsize=10)
axes[3].imshow(np.where(zf & (zly == 0), 1, np.nan), cmap=ListedColormap(["#c7dfb7"]), extent=zext, interpolation="nearest")
im = axes[3].imshow(np.where(zf & (zly > 0), zly + 2000, np.nan), cmap=ly_cmap, vmin=2000.5, vmax=2025.5, extent=zext,
                    interpolation="nearest")
fig.colorbar(im, ax=axes[3], shrink=0.8, label="GFC loss year")
axes[3].set_title("GFC v1.13 loss year", fontsize=10)
for ax in axes:
    ax.set_xticks([]); ax.set_yticks([])
fig.savefig(OUT / "04-before-after.png")
plt.close(fig)

# ---------------------------------------------------------------- Step 6: NBR check against the GFC loss year
section("Sentinel-2 dry-season NBR, 2017-2025")
nbr_years = {}
for yr in range(2017, 2026):
    items = best_s2(yr, zb)
    nbr_years[yr] = nbr_composite(items, zb, (Z, Z))
    print(f"{yr}: {len(items)} scenes ({', '.join(f'{it.datetime:%d %b}' for it in items)}), "
          f"valid {np.isfinite(nbr_years[yr]).mean():.1%}")
stable = zf & (zly == 0)
traj = {"stable forest": [np.nanmedian(nbr_years[y][stable]) for y in nbr_years]}
for ly in range(18, 26):
    m = zf & (zly == ly)
    if m.sum() > 200:
        traj[f"lost in {2000 + ly}"] = [np.nanmedian(nbr_years[y][m]) for y in nbr_years]
print("median NBR by year (rows: GFC class, columns: 2017..2025)")
for k, v in traj.items():
    print(f"  {k:14s} " + " ".join(f"{x:5.2f}" for x in v))

# year of the largest year-on-year NBR drop vs the GFC year, pixel by pixel (loss 2018-2025)
yrs = np.array(list(nbr_years))
cube = np.stack([nbr_years[y] for y in yrs])
drops = cube[:-1] - cube[1:]                         # positive = NBR fell from year t to t+1
valid = np.isfinite(drops).all(axis=0)
drop_year = yrs[1:][np.nanargmax(np.where(np.isfinite(drops), drops, -9), axis=0)]
biggest = np.nanmax(np.where(np.isfinite(drops), drops, -9), axis=0)
sel = zf & (zly >= 18) & valid
diff = drop_year[sel] - (zly[sel] + 2000)
print(f"loss pixels 2018-2025 with a full NBR record: {sel.sum():,}")
print(f"largest annual NBR drop in the GFC year: {(diff == 0).mean():.1%}; one year later: {(diff == 1).mean():.1%}; "
      f"one year earlier: {(diff == -1).mean():.1%}; further off: {(abs(diff) > 1).mean():.1%}")
print(f"largest drop >= 0.2 for {(biggest[sel] >= 0.2).mean():.1%} of GFC loss pixels vs "
      f"{(biggest[stable & valid] >= 0.2).mean():.1%} of stable-forest pixels")
dnbr = nbr_years[2017] - nbr_years[2025]
print(f"dNBR 2017->2025 median: stable forest {np.nanmedian(dnbr[stable]):.3f}, "
      f"lost 2018-2025 {np.nanmedian(dnbr[zf & (zly >= 18)]):.3f}, "
      f"lost 2001-2017 {np.nanmedian(dnbr[zf & (zly >= 1) & (zly <= 17)]):.3f}")

fig, axes = plt.subplots(1, 2, figsize=(13, 4.6), constrained_layout=True, gridspec_kw={"width_ratios": [1, 1.25]})
im = axes[0].imshow(dnbr, cmap="RdYlGn_r", vmin=-0.4, vmax=0.8, extent=zext, interpolation="nearest")
fig.colorbar(im, ax=axes[0], shrink=0.85, label="dNBR = NBR 2017 - NBR 2025")
axes[0].contour(np.flipud(zf & (zly >= 18)).astype(float), levels=[0.5], colors="#222222", linewidths=0.5,
                extent=zext)
axes[0].set_title("Sentinel-2 dNBR, dry season 2017 to 2025\n(black outline: GFC loss 2018-2025)", fontsize=10)
axes[0].set_xticks([]); axes[0].set_yticks([])
cols = lambda i: plt.get_cmap("plasma")(0.85 * i / max(1, len(traj) - 2))
for i, (k, v) in enumerate(traj.items()):
    style = dict(color="#1a9850", lw=2.5) if k == "stable forest" else dict(color=cols(i - 1), lw=1.4)
    axes[1].plot(yrs, v, marker="o", ms=3, label=k, **style)
axes[1].set_xlabel("dry season (June-August)"); axes[1].set_ylabel("median NBR")
axes[1].set_title("NBR of pixels grouped by their GFC loss year", fontsize=10)
axes[1].legend(frameon=False, fontsize=8, ncol=2, loc="lower left")
fig.savefig(OUT / "05-nbr-check.png")
plt.close(fig)

# ---------------------------------------------------------------- Step 7: how much of the loss burned? (MODIS MCD64A1)
section("MODIS burned area vs GFC loss")
burn_rows = []
for yr in range(2019, 2026):
    burned = np.zeros((H, W), dtype=bool)
    items = catalog.search(collections=["modis-64A1-061"], bbox=AOI, datetime=f"{yr}-01-01/{yr}-12-31").item_collection()
    for it in items:
        bd = read_to_grid(it.assets["Burn_Date"].href, (AOI[0], AOI[1], AOI[2], AOI[3]), (H, W), Resampling.nearest,
                          dtype="int16")
        burned |= bd > 0
    m = forest2000 & (lossyear == yr - 2000)
    months = len({it.properties["start_datetime"][:7] for it in items})
    burn_rows.append((yr, months, ha(m), ha(m & burned), ha(burned & land)))
burn = pd.DataFrame(burn_rows, columns=["year", "modis_months", "gfc_loss_ha", "loss_in_modis_burn_ha", "modis_burned_ha"])
burn["share_burned"] = (burn.loss_in_modis_burn_ha / burn.gfc_loss_ha).where(burn.modis_months == 12)
print(burn.to_string(index=False, float_format=lambda v: f"{v:,.2f}" if v < 1 else f"{v:,.0f}"))
burn.to_csv(OUT / "loss_vs_modis_burned.csv", index=False)

fig, ax = plt.subplots(figsize=(9, 3.9), constrained_layout=True)
ax.bar(burn.year, burn.gfc_loss_ha / 1000, color="#abdda4", label="GFC forest loss")
ax.bar(burn.year, burn.loss_in_modis_burn_ha / 1000, color="#d53e4f", label="of which inside a MODIS burn scar that year")
for _, r in burn.iterrows():
    label = f"{r.share_burned:.0%}" if r.modis_months == 12 else f"{r.modis_months:.0f} months\nof MODIS"
    ax.text(r.year, r.gfc_loss_ha / 1000 + 0.15, label, ha="center", fontsize=8)
ax.set_ylim(0, burn.gfc_loss_ha.max() / 1000 * 1.25)
ax.set_ylabel("thousand ha"); ax.legend(frameon=False, loc="upper left", ncol=2)
ax.set_title("How much of each year's GFC loss burned in the same year (MODIS MCD64A1, 500 m)", fontsize=10)
fig.savefig(OUT / "06-loss-vs-fire.png")
plt.close(fig)
full = burn[burn.modis_months == 12]
print(f"complete years {full.year.min()}-{full.year.max()}: {full.loss_in_modis_burn_ha.sum() / full.gfc_loss_ha.sum():.1%} "
      f"of GFC loss sat inside a same-year MODIS burn scar")
print("\nwrote figures 01-06, output/loss_by_year.csv and output/loss_vs_modis_burned.csv")