---
name: "watershed"
description: "🗿 Watershed DRR in GRASS: r.watershed→r.stream.order→strahler 1=ruscello→streams.shp+watershed.shp+gpkg. Include script standalone e fallback SQL per area_ha."
model_tier: 1
triggers:
  - "watershed"
  - "bacino"
  - "stream network"
  - "drenaggio"
  - "catchment"
  - "drainage basin"
tools:
  - exec
---

# Watershed Analysis — SCOLPITO NELLA ROCCIA 🗿

## REGOLA ASSOLUTA (NON VIOLARE MAI)
Quando serve una **watershed analysis**, segui:
1. **CHIEDERE** A) area media bacino in ha (100/200/500) e B) se DTM già disponibile
2. **SE DTM MANCANTE** → chiedere all'utente di caricarlo su `http://100.86.154.25:9090/` oppure scaricare Copernicus GLO-30
3. **ESEGUIRE** il workflow GRASS qui sotto
4. **CONSEGNARE**: streams.shp + watershed.shp + gpkg + ZIP via HTTP

## Prerequisiti
```bash
# Una tantum
grass --exec g.extension extension=r.stream.order operation=add
```

## Workflow GRASS

### Calcolo soglia
```bash
# Soglia = (B * 10000) / 900    (per FathomDEM 30m)
# 100 ha → soglia 1100
# 200 ha → soglia 2200
# 500 ha → soglia 5500
THRESHOLD=1100
```

### Pipeline completa
```bash
DTM="DTM_UTM37_100_Croped.tif"

# 1. Location GRASS
grass -c "$DTM" -e grassdata/location

# 2. Import
grass grassdata/location/PERMANENT --exec r.in.gdal -o input="$DTM" output=dem --overwrite
grass grassdata/location/PERMANENT --exec g.region raster=dem

# 3. r.watershed ⏱️ ~3-5 min per 6MB DTM
grass grassdata/location/PERMANENT --exec r.watershed -a elevation=dem accumulation=acc drainage=dir stream=str basin=ws threshold=$THRESHOLD

# 4. Strahler order
grass grassdata/location/PERMANENT --exec r.stream.order elevation=dem stream_rast=str strahler=str_order direction=dir

# 5. Stream → vector con strahler
grass grassdata/location/PERMANENT --exec r.thin input=str output=str_th
grass grassdata/location/PERMANENT --exec r.to.vect -s input=str_th output=str_lines type=line
grass grassdata/location/PERMANENT --exec v.db.addcolumn map=str_lines columns="strahler integer"
grass grassdata/location/PERMANENT --exec v.rast.stats map=str_lines raster=str_order method=average column_prefix=ord
grass grassdata/location/PERMANENT --exec db.execute sql="UPDATE str_lines SET strahler = ROUND(ord_average)"
grass grassdata/location/PERMANENT --exec v.db.dropcolumn map=str_lines columns="ord_average,value,label"
grass grassdata/location/PERMANENT --exec v.out.ogr input=str_lines output=streams.shp format=ESRI_Shapefile --overwrite

# 6. Watershed → vector con area_ha ⚠️
grass grassdata/location/PERMANENT --exec r.to.vect -s input=ws output=ws_v type=area
grass grassdata/location/PERMANENT --exec v.db.addcolumn map=ws_v columns="area_ha double precision"
# Prova v.to.db (a volte scrive null — ok, fixiamo dopo in GPKG)
grass grassdata/location/PERMANENT --exec v.to.db map=ws_v option=area units=hectares columns=area_ha --overwrite 2>/dev/null || true
grass grassdata/location/PERMANENT --exec v.out.ogr input=ws_v output=watershed.shp format=ESRI_Shapefile --overwrite

# 7. GeoPackage con area_ha via SQL (MOLTO PIÙ AFFIDABILE)
ogr2ogr -f GPKG output.gpkg streams.shp -nln streams -overwrite
ogr2ogr -f GPKG -update output.gpkg watershed.shp -nln watershed -overwrite
ogrinfo output.gpkg -sql "UPDATE watershed SET area_ha = ROUND(ST_Area(geom) / 10000.0, 1)"

# 8. ZIP
zip -9 output.zip streams.* watershed.* output.gpkg

# 9. Verifica area_ha
ogrinfo -q -dialect SQLite -sql "SELECT ROUND(MIN(area_ha)) as min_ha, ROUND(AVG(area_ha)) as avg_ha, ROUND(MAX(area_ha)) as max_ha FROM watershed WHERE area_ha > 1" output.gpkg
# Output atteso: min 1, avg ~136, max ~1179 ha (per soglia 1100)

# 10. Verifica distribuzione Strahler
ogrinfo -q -dialect SQLite -sql "SELECT strahler, COUNT(*) as n FROM streams GROUP BY strahler ORDER BY strahler" output.gpkg
# Output atteso: null(149), 1(485), 2(283), 3(134), 4(38), 5(39) — soglia 1100, area Etiopia
```

## Script standalone
```bash
bash /home/buldo/.hermes/skills/gis/watershed/scripts/watershed_pipeline.sh DTM.tif [soglia] [outdir]
# Esempio:
bash /home/buldo/.hermes/skills/gis/watershed/scripts/watershed_pipeline.sh DTM_UTM37_100_Croped.tif 1100 ~/output
```

## Consegna via HTTP
```bash
# Copiare in cartella served
cp streams.* watershed.* output.gpkg output.zip /home/buldo/.hermes/generated/

# URL di accesso:
# http://100.86.154.25:8088/output.zip
# http://100.86.154.25:8088/streams.shp
# http://100.86.154.25:8088/watershed.shp
# In QGIS: Layer → Add Layer → Add Vector Layer → URL
```

## ⚠️ Pitfalls (imparati a proprie spese 12.06.2026)

### 1. area_ha null nel SHP — IL PIÙ CRITICO
- **CAUSA**: `v.to.db` può scrivere valori null se la colonna esisteva già da una run GRASS precedente o se il `--overwrite` non funziona come previsto.
- **SINTOMO**: nell'SHP tutti i valori area_ha sono null.
- **SOLUZIONE**: NON fidarsi di v.to.db. Ricalcolare sempre via SQL su GPKG:
  ```bash
  ogrinfo output.gpkg -sql "UPDATE watershed SET area_ha = ROUND(ST_Area(geom) / 10000.0, 1)"
  ```
  Questo è l'unico metodo che FUNZIONE SEMPRE. Va eseguito DOPO l'export.

### 2. Null strahler nei stream
- **CAUSA**: r.thin produce artefatti ai bordi dell'estensione DTM (149 su 1,128 linee in questo run).
- **EFFETTO**: Non critico — filtrare in QGIS con `strahler IS NOT NULL` o ignorare.
- **POSSIBILE FIX**: Applicare buffer interno di 1 cella prima di r.thin.

### 3. DTM in EPSG:4326 vs UTM
- GRASS crea la location dalla proiezione del TIFF. Se il DTM è in EPSG:4326 (gradi), l'area calcolata sarà in gradi², non m².
- **CONTROLLARE SEMPRE**: 
  ```bash
  gdalinfo DTM.tif | grep "AUTHORITY\|EPSG"
  ```
  Deve essere UTM (EPSG:326xx per emisfero nord, 327xx per sud) per calcoli area in metri.
- **SE IN 4326**: riproiettare prima:
  ```bash
  gdalwarp -t_srs EPSG:32637 DTM.tif DTM_UTM.tif
  ```

### 4. Nomi file hardcoded
- La skill usa `streams.shp` e `watershed.shp` come output. Se si esegue più run nella stessa cartella, sovrascrivono i precedenti.
- **SOLUZIONE**: usare output dir separato per ogni run (lo script standalone lo fa automaticamente con timestamp).

### 5. GRASS location già esistente
- Se `grassdata/location/` esiste già, grass -c fallisce.
- **SOLUZIONE**: cancellare prima:
  ```bash
  rm -rf grassdata/location  # O usare cartella con timestamp
  ```

## Output finali
- `streams.shp` — linee, attributo **strahler** (1→N)
- `watershed.shp` — poligoni, attributo **area_ha** (ha, via SQL)
- `output.gpkg` — GeoPackage con entrambi + area_ha corretta
- `output.zip` — tutti i file zippati

## Esempio output reale (Etiopia Melkadida, 12.06.2026)
- DTM: `DTM_UTM37_100_Croped.tif`, FathomDEM 30m, UTM37N
- Soglia: 1100 (~100 ha/bacino)
- Streams: **1,128** linee (Strahler 1-5 + 149 null edge)
- Bacini: **1,116** poligoni
- Area media: **136 ha**, max: 1,179 ha, totale: ~1,517 km²
- ZIP: 1.2 MB (shp+gpkg)

## Riferimenti
- `references/ethiopia-melkadida-2026-06-12.md` — report dettagliato di questo run
- `scripts/watershed_pipeline.sh` — script bash standalone riutilizzabile
- Video: https://youtu.be/SHB3JOffL7M + https://youtu.be/PILoYxp6Ipk

## 🗿 SCOLPITO NELLA ROCCIA il 08.06.2026 | Rinominato 12.06.2026 | Refined 12.06.2026 (v2: SQL area fallback + pitfalls reali + script standalone + verification steps)
## NESSUN CAMBIAMENTO SENZA APPROVAZIONE ESPLICITA DI BULDO