# RiverBankErosionAccretion v4

QGIS plugin that quantifies riverbank erosion and accretion directly
from satellite-derived waterbody polygons.

## Differentiators vs. DSAS / manual-bankline QGIS tools

| Gap in existing tools | This plugin |
|---|---|
| Require manually digitised banklines | Accepts satellite-derived waterbody **polygons** directly |
| Report linear rate *or* area, never both | Produces **km² areal and m/yr linear** rates from the same inputs |
| Assume single-thread channels | **Regime-aware** processing with auto-classification |
| No input-data uncertainty quantification | **JRC occurrence-threshold sensitivity** module |

## Architecture — three modules

```
INPUT: Year1 waterbody polygon (.shp), Year2 waterbody polygon (.shp),
       projected CRS, output folder

Module 1 - Grid-based change detection        (ALWAYS RUNS)
Module 2 - Transect-based rates                (optional, auto-gated)
Module 3 - Uncertainty quantification          (optional)

OUTPUT: structured folder tree of shapefiles + CSVs
```

### Module 1 — grid-based change detection (always runs)

Reprojects both epochs to the chosen CRS, computes a **dynamic riparian
buffer** (`max(channel_width x 2, 100 m)`, where `channel_width =
union(Year1).area / union(Year1).length`), runs connected-component
labelling + a Minimum Mapping Unit area filter on both epochs, builds a
30 m fishnet clipped to the riparian corridor, differences it per cell,
aggregates along-reach into user-selectable stretches (90 m / 3 km /
6 km / 9 km / custom) using a Voronoi-skeleton centerline for true
along-channel `From_m`/`To_m` and sinuosity, ranks Top-N hotspots
(default 15, 1-30), and computes a braiding index (BI) and 5-class
regime (`REGULATED_SINGLE_THREAD`, `MEANDERING`, `BRAIDED`,
`HYPER_BRAIDED`, `DELTAIC`) per stretch.

**Erosion/accretion convention** (do not reverse): EROSION = new water
in Year2 not present in Year1 (bank retreat) = `difference(Year2,
Year1)`. ACCRETION = water lost between Year1 and Year2 (bank advance)
= `difference(Year1, Year2)`.

### Module 2 — transect-based rates (optional)

Generates perpendicular transects along the Module 1 centerline and
measures left/right bank retreat rates (m/yr) plus a grid-transect
Consistency Index (CI). **Auto-disabled** when any stretch has BI ≥ 3
(hyper-braided) or a channel cluster count > 3 (deltaic); shows a
MODERATE-reliability warning for 2.0 ≤ BI < 3.0.

### Module 3 — uncertainty quantification (optional)

Re-runs Module 1 at up to three JRC Global Surface Water occurrence
thresholds (70% / 85% / 100%) and reports per-cell/per-stretch
min/max/range/mean erosion & accretion envelopes. If only the standard
100% pair is supplied, this degrades gracefully to a single-threshold
best estimate. An optional multi-epoch sub-analysis computes a rate
acceleration metric `RA = (rate_interval2 - rate_interval1) /
rate_interval1`.

## Module files

| File | Role |
|---|---|
| `riverbank_erosion_accretion_v2.py` | QGIS plugin registration + pipeline orchestration (only file that touches PyQGIS layer objects) |
| `riverbank_erosion_accretion_v2_dialog.py` | 5-section GUI (inputs, Module 1/2/3 params, run/log) |
| `change_analysis.py` | Module 1 engine |
| `transect_module.py` | Module 2 engine (also supplies the centerline Module 1 uses for stretch ordering) |
| `uncertainty_module.py` | Module 3 engine |
| `regime_classifier.py` | 5-class regime taxonomy |
| `input_validator.py` | Pre-run parameter checks + post-run output sanity checks |
| `output_manager.py` | Folder tree, shapefile/CSV schema, metadata.json, auto-symbology |
| `validation.py` | Separate: JRC-reference accuracy validation for satellite water extraction (not part of the 3-module pipeline) |

All logic except the GUI and `output_manager.apply_symbology()` is
pure geopandas/shapely/numpy, independent of PyQGIS, so it can be
exercised and unit-tested outside a running QGIS instance.

## Output structure

```
output_folder/
├── Grid/            Erosion_Grids.shp, Accretion_Grids.shp, Unchanged_Grids.shp, Fishnet_Study_Area.shp
├── Stretch/          Stretch_Summary.shp/.csv, Top15_Erosion.shp, Top15_Accretion.shp
├── Classification/   Reach_Regime.shp
├── Transect/         Centerline.shp, Transects.shp, Transect_Rates.shp/.csv, BI_per_stretch.csv   (Module 2 only)
├── Uncertainty/      Uncertainty_Envelope.shp, Uncertainty_Stretch.csv, Consistency_Index.csv,
│                     Rate_Acceleration.csv, Interval_Comparison.shp                                (Module 3 only)
└── Summary/          Analysis_Summary.csv, Analysis_Log.txt, metadata.json
```

Shapefile (.dbf) field names are capped at 10 characters, so
shapefiles use abbreviated field names (e.g. `Eros_km2`); CSV exports
use the full names (e.g. `Erosion_km2`). See the `*_SHP_RENAME` /
`*_CSV_RENAME` maps in `output_manager.py`.

## Known limitation — very large corridors

Module 1's per-cell differencing uses bulk `geopandas.overlay` operations
(not a per-cell Python loop), verified to handle corridors up to
~2.1M grid cells (e.g. a 150 km single-thread reach) in a few minutes.

On an extreme corridor — a very wide, heavily braided reach spanning a
large fraction of its own bounding box, >2M cells with input polygons
split into 700+ disjoint parts — a single whole-grid `overlay()` call
for the final grid-clip step did not complete in a 25-minute session
budget. This is fixed: `compute_per_grid_change` now routes the clip
through `_overlay_intersection_tiled` (grids above
`GRID_TILE_THRESHOLD` = 150,000 cells) in `change_analysis.py`, which
splits the fishnet into geographic tiles (~`TARGET_CELLS_PER_TILE` =
5,000 cells each) and, per tile, first truncates candidate raw
difference polygons to that tile's bounding box (padded by one cell
size, so cells whose 30 m box straddles a tile edge still get their
full, correct intersection area) before overlaying against just that
tile's cells. Truncating first is what makes this fast: without it, a
single differencing sliver that follows the whole braided corridor
gets re-intersected at full geometric complexity against every tile it
touches, which is *slower* than no tiling at all.

Verified on the Brahmaputra/Majuli hyper-braided stress test (~2.13M
cells, corridor 1831.9 km²): the grid-clip step (all three of
erosion/accretion/unchanged) now completes in ~328 s (~5.5 min, down
from not completing in 25 min), and a full `run_module1()` call
completes in ~1175 s (~19.6 min) — the remainder beyond the clip step
is stretch-aggregation (`dissolve`) and braiding-index computation over
2.13M cells, which was not part of this fix's scope and is a candidate
for future optimization if faster whole-Module-1 times are needed.
Results match hand-verified expectations (617.125 km² erosion / 259.532
km² accretion; Module 2 correctly auto-disabled, BI up to 18.46). The
Narmada and Palasbari baselines are unchanged by this refactor (Narmada:
82.546 / 6.795 km², exact match); Palasbari incidentally got much faster
too (106 s vs. the previous 793 s), since its heavily braided geometry
hit the same complexity trap on a smaller scale.

Note: Module 3 (optional uncertainty quantification) re-runs
`run_module1()` once per supplied JRC occurrence threshold, so enabling
it multiplies total time accordingly on a corridor this size — expected
behavior of that optional feature, not a regression from this fix.

## Requirements

- QGIS 3.28+, Python 3.8+
- `geopandas`, `numpy`, `shapely`, `pandas`
- Optional (legacy satellite-fetch tab): `pystac-client`, `rasterio`

## How to use

1. Open the plugin from the QGIS toolbar.
2. Select Year 1 / Year 2 waterbody polygon layers, set years, pick a
   **projected** CRS, and choose an output folder.
3. Adjust Module 1 parameters (stretch distance, MMU, Top-N) if needed.
4. Optionally enable Module 2 (transects) — greyed out automatically
   for braided/deltaic channels.
5. Optionally enable Module 3 (uncertainty) and supply 70%/85% JRC
   threshold polygons if available.
6. Click **Start Analysis**.

## Credits

Developed by **Dr. Srija Roy**
