init commit, needs final approval
This commit is contained in:
BastaMasta
+34
@@ -0,0 +1,34 @@
|
||||
# python / uv
|
||||
.venv/
|
||||
venv/
|
||||
__pycache__/
|
||||
*.py[cod]
|
||||
*.egg-info/
|
||||
.python-version
|
||||
|
||||
# downloaded + generated data (tens of GB, never commit)
|
||||
data/**
|
||||
|
||||
# but keep the empty gauge folders so people know where the xlsx files go
|
||||
!data/gauge/
|
||||
!data/gauge/wse/
|
||||
!data/gauge/discharge/
|
||||
!data/gauge/wse/.gitkeep
|
||||
!data/gauge/discharge/.gitkeep
|
||||
!data/gauge/wse/README.txt
|
||||
!data/gauge/discharge/README.txt
|
||||
|
||||
# credentials - never commit these
|
||||
.netrc
|
||||
_netrc
|
||||
*.env
|
||||
.env
|
||||
|
||||
# editor / os junk
|
||||
.vscode/
|
||||
.idea/
|
||||
.DS_Store
|
||||
Thumbs.db
|
||||
|
||||
# matplotlib font cache etc
|
||||
.cache/
|
||||
@@ -0,0 +1,70 @@
|
||||
# swot river discharge pipeline
|
||||
|
||||
## how it works (short version)
|
||||
|
||||
- SWOT L2 HR RiverSP version D (`SWOT_L2_HR_RiverSP_D`), node + reach shapefiles from nasa earthdata
|
||||
- ALOS World 3D 30m DSM via GEE (or a local geotiff if you dont have GEE)
|
||||
- gauge WSE + discharge from india-WRIS
|
||||
- match each gauge to the nearest SWOT node (haversine), match times within ±3h
|
||||
- depth = WSE − bed, A = width × depth, R = A/P, then **Q = (1/n) · A · R^(2/3) · S^(1/2)**
|
||||
- compare against gauge discharge, get R²/RMSE/NSE/bias + plots
|
||||
|
||||
## setup
|
||||
|
||||
```
|
||||
# Install UV
|
||||
|
||||
# Windows:
|
||||
powershell -ExecutionPolicy ByPass -c "irm https://astral.sh/uv/install.ps1 | iex"
|
||||
|
||||
#linux/macOS:
|
||||
curl -LsSf https://astral.sh/uv/install.sh | sh
|
||||
|
||||
```
|
||||
|
||||
**nasa earthdata** - first run will ask for username/password once and save them. also approve the "PO.DAAC" app on the earthdata site if downloads fail.
|
||||
|
||||
**google earth engine** - register at https://code.earthengine.google.com/register, then `uv run earthengine authenticate` and put your project id in `config.yaml` under `gee.project`.
|
||||
|
||||
**gauge data** - put the WRIS xlsx exports here:
|
||||
|
||||
```
|
||||
data/gauge/wse/ River Water Level_<Station>.xlsx
|
||||
data/gauge/discharge/ River Water Discharge_<Station>.xlsx
|
||||
```
|
||||
|
||||
and add each station's lat/lon to `stations.csv`. station name has to match the filename part after the underscore ("(seasonal)" is ignored). get the exact coordinates from the station page on WRIS - if they're wrong the pipeline matches the wrong river and you get garbage (ask me how i know).
|
||||
|
||||
## running it
|
||||
|
||||
```
|
||||
uv run python run_pipeline.py # all 8 steps
|
||||
uv run python run_pipeline.py --steps 3-8 # skip the downloads
|
||||
uv run python run_pipeline.py --steps 6,7,8 # just recompute discharge + export
|
||||
```
|
||||
|
||||
| # | name | what it does | output |
|
||||
| --- | --------- | ------------------------------------------------------------------------------ | -------------------------------------------- |
|
||||
| 1 | download | SWOT granules from earthdata (~50gb for 2 yrs, skips already-downloaded files) | `data/swot/` |
|
||||
| 2 | process | merge shapefiles, clip to basin, clean (slow, rerun only if AOI/dates change) | `data/processed/nodes_clean.csv` etc |
|
||||
| 3 | gauge | parse the WRIS xlsx files (picks the MSL series, converts IST→UTC) | `gauge_wse.csv`, `gauge_discharge.csv` |
|
||||
| 4 | wse | node matching + WSE validation | `matched_wse.csv`, `outputs/wse_metrics.csv` |
|
||||
| 5 | dem | bed elevation (min DSM within 300m buffer) | `node_bed_elevation.csv` |
|
||||
| 6 | discharge | manning's equation | `swot_discharge.csv` |
|
||||
| 7 | validate | SWOT vs gauge discharge | `outputs/discharge_metrics.csv`, plots |
|
||||
| 8 | export | everything into one excel workbook for analysis | `outputs/<project>_swot_analysis.xlsx` |
|
||||
|
||||
## tuning knobs (config.yaml)
|
||||
|
||||
- `hydraulics.mannings_n` - global roughness, default 0.035
|
||||
- `hydraulics.mannings_n_per_station` - override n per station, this is the main thing to play with. if a station is consistently overestimating by ~4x, try ~4x the n (it scales Q linearly). the n used shows up in the outputs so you always know which run used what
|
||||
- `hydraulics.baseline_depth_min_m` / `baseline_depth_max_m` - the DEM is unreliable at the channel (bank pixels, voids) so the implied low-water depth gets clamped into this band, default [0.5, 2.0]
|
||||
- `quality.max_node_q`, `quality.swot_wse_mad_k`, `matching.wse_outlier_mad_k` - outlier filtering. SWOT throws the occasional ±15m garbage observation, these catch it
|
||||
- `matching.*` - distance limit + time windows
|
||||
- `hydraulics.scope: all` - compute discharge at every node in the basin, not just gauge-matched ones (step 5 gets slower)
|
||||
|
||||
## things to know before running the pipeline
|
||||
|
||||
- the "bed elevation" from the DEM is really the water surface the satellite saw at acquisition time, so depths are relative, not true bathymetry. thats why I created a clamp
|
||||
- SWOT WSE is EGM2008, gauges are local MSL/datum - expect a constant offset (plots auto-shift when its > 1m, metrics report the raw bias)
|
||||
- dry season: SWOT sees pooled water behind weirs and manning's cant tell a standing pool from flowing water, so theres a discharge floor where the gauge says zero. like estimating traffic from how full the parking lot is - utter garbage data
|
||||
+87
@@ -0,0 +1,87 @@
|
||||
# ============================================================
|
||||
# SWOT River Discharge Pipeline - Configuration
|
||||
# Study: River Discharge Estimation Using SWOT Satellite
|
||||
# Observations (Krishna Basin methodology, NRSC/ISRO)
|
||||
# ============================================================
|
||||
|
||||
project_name: krishna
|
||||
|
||||
# Area of interest (bounding box, WGS84). Default = Krishna Basin.
|
||||
aoi:
|
||||
lon_min: 73.0
|
||||
lat_min: 13.0
|
||||
lon_max: 81.5
|
||||
lat_max: 19.5
|
||||
|
||||
# Analysis period
|
||||
time:
|
||||
start: "2023-07-01"
|
||||
end: "2025-09-01"
|
||||
|
||||
# SWOT Level-2 HR River Single-Pass Vector product, Version D
|
||||
swot:
|
||||
short_name: SWOT_L2_HR_RiverSP_D
|
||||
continent_code: AS # Asia granules only
|
||||
|
||||
# ALOS World 3D 30m DSM via Google Earth Engine
|
||||
gee:
|
||||
project: earth-engine-river
|
||||
dem_collection: JAXA/ALOS/AW3D30/V3_2
|
||||
band: DSM
|
||||
scale_m: 30
|
||||
buffer_m: 300 # take the MINIMUM DSM within this radius of the node
|
||||
# (point samples often hit banks/bridges, not the channel)
|
||||
|
||||
# Optional fallback: local DEM GeoTIFF (used if GEE unavailable / project null)
|
||||
dem:
|
||||
local_tif: null # e.g. data/dem/alos_krishna.tif
|
||||
|
||||
# SWOT observation quality screening
|
||||
quality:
|
||||
max_node_q: 2 # keep node obs with node_q <= this (0 best, 3 worst);
|
||||
# applied in step 2 when the node_q column is present
|
||||
swot_wse_mad_k: 5 # per-node robust outlier filter on SWOT WSE before
|
||||
# discharge computation: drop obs further than
|
||||
# k * MAD from the node median (0 disables)
|
||||
|
||||
# Gauge data (user places India-WRIS/CWC xlsx exports here)
|
||||
gauge:
|
||||
wse_dir: data/gauge/wse
|
||||
discharge_dir: data/gauge/discharge
|
||||
stations_csv: stations.csv
|
||||
|
||||
# Spatial + temporal matching
|
||||
matching:
|
||||
max_distance_km: 10.0 # max gauge-to-node distance (Haversine)
|
||||
wse_time_window_hours: 3 # +/- window for WSE validation
|
||||
discharge_time_window_hours: 24 # +/- window for discharge validation
|
||||
wse_outlier_mad_k: 7 # drop matched pairs whose (SWOT - gauge)
|
||||
# residual is > k * MAD from the station
|
||||
# median residual (0 disables)
|
||||
|
||||
# Hydraulic computation (Manning's equation)
|
||||
hydraulics:
|
||||
mannings_n: 0.035 # default roughness coefficient (natural channel)
|
||||
mannings_n_per_station: # optional per-station overrides for experimentation;
|
||||
# stations not listed use mannings_n above
|
||||
# Keesara: 0.12
|
||||
# Kurundwad: 0.045
|
||||
min_depth_m: 0.10 # discard non-physical depths below this
|
||||
scope: matched # "matched" = gauge-matched nodes only, "all" = every node
|
||||
bed_mode: auto # "auto": clamp the low-water baseline depth implied by
|
||||
# the DEM (P5 of SWOT WSE minus bed) into the band below -
|
||||
# the ALOS DSM is unreliable at the channel in both
|
||||
# directions; "dem": use the DEM value as-is
|
||||
baseline_depth_min_m: 0.5
|
||||
baseline_depth_max_m: 2.0
|
||||
slope_source: median # "median" = reach median over all passes (stable);
|
||||
# "pass" = same-pass slope (noisy, CV often > 100%)
|
||||
width_source: median # "median" = node median width; "pass" = per-pass
|
||||
|
||||
# Folder layout (relative to pipeline root)
|
||||
paths:
|
||||
swot_raw: data/swot/raw
|
||||
swot_nodes: data/swot/nodes
|
||||
swot_reaches: data/swot/reaches
|
||||
processed: data/processed
|
||||
outputs: outputs
|
||||
@@ -0,0 +1,26 @@
|
||||
PUT THE DISCHARGE XLSX FILES HERE
|
||||
==================================
|
||||
|
||||
what goes in this folder: the "River Water Discharge" excel exports
|
||||
from india-WRIS (indiawris.gov.in), one file per gauge station.
|
||||
|
||||
filenames MUST look like this (this is what WRIS gives you by default,
|
||||
so just don't rename them):
|
||||
|
||||
River Water Discharge_Arjunwad (seasonal).xlsx
|
||||
River Water Discharge_Kurundwad.xlsx
|
||||
|
||||
the part after the underscore is the station name. it has to EXACTLY
|
||||
match the "station" column in stations.csv at the project root
|
||||
(the "(seasonal)" part is ignored) AND match the matching water level
|
||||
file in ../wse/ - same station, same spelling.
|
||||
|
||||
discharge data is used for VALIDATION (comparing the satellite-derived
|
||||
discharge against what was actually measured). a station without a
|
||||
discharge file still gets its water level validated, it just won't
|
||||
show up in the discharge metrics.
|
||||
|
||||
do NOT put anything else in here. the pipeline only reads .xlsx/.xls
|
||||
files from this folder. this README is ignored.
|
||||
|
||||
then run: uv run python run_pipeline.py --steps 3-8
|
||||
@@ -0,0 +1,27 @@
|
||||
PUT THE WATER LEVEL XLSX FILES HERE
|
||||
====================================
|
||||
|
||||
what goes in this folder: the "River Water Level" excel exports from
|
||||
india-WRIS (indiawris.gov.in), one file per gauge station.
|
||||
|
||||
filenames MUST look like this (this is what WRIS gives you by default,
|
||||
so just don't rename them):
|
||||
|
||||
River Water Level_Arjunwad (seasonal).xlsx
|
||||
River Water Level_Kurundwad.xlsx
|
||||
|
||||
the part after the underscore is the station name. it has to EXACTLY
|
||||
match the "station" column in stations.csv at the project root
|
||||
(the "(seasonal)" part is ignored, don't worry about it).
|
||||
|
||||
for every file you add here:
|
||||
1. also add the station's lat/lon to stations.csv (get the exact
|
||||
coordinates from the station's page on WRIS - wrong coordinates =
|
||||
the pipeline matches the wrong river = garbage results)
|
||||
2. if you have discharge data for it too, put that file in
|
||||
../discharge/
|
||||
|
||||
do NOT put anything else in here. the pipeline only reads .xlsx/.xls
|
||||
files from this folder. this README is ignored.
|
||||
|
||||
then run: uv run python run_pipeline.py --steps 3-8
|
||||
@@ -0,0 +1,6 @@
|
||||
def main():
|
||||
print("Hello from swot-flow-pipeline!")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -0,0 +1,7 @@
|
||||
station,n,r2,rmse,nse,bias
|
||||
Arjunwad,38,0.9377272963911925,194.76574090329413,0.9243355373639429,58.52270121758602
|
||||
Dameracherla,28,0.2701885692936038,118.35985971046465,-0.07823247511048304,61.388205407483746
|
||||
Huvinhedigi,48,0.5852912736923251,135.39496292469656,0.36830792328092954,77.96814914653531
|
||||
Keesara,49,0.640828095229739,295.0904056862152,-12.647848628217174,271.5735139765754
|
||||
Kurundwad,32,0.6273400159711809,447.7583964500478,0.5349628652228668,167.32312943263042
|
||||
Wadenepally,51,0.5613120947826807,855.5593600465282,-8.367856786512865,785.4747002371549
|
||||
|
Binary file not shown.
Binary file not shown.
|
After Width: | Height: | Size: 108 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 104 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 104 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 121 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 101 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 104 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 117 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 103 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 112 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 100 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 112 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 96 KiB |
@@ -0,0 +1,7 @@
|
||||
station,n,r2,rmse,nse,bias
|
||||
Arjunwad,36,0.9484101295300449,1.070147937088865,0.8504611422678376,0.8638955555555552
|
||||
Dameracherla,23,0.01387788681088579,16.828991764372887,-78.67803667119695,16.700217391304353
|
||||
Huvinhedigi,53,0.9944665724968776,0.9900936623103742,0.7641690886715478,0.9782883018867918
|
||||
Keesara,44,0.9430136816999417,0.48048752582899845,0.7479876970288146,0.4223281818181816
|
||||
Kurundwad,26,0.9887656982062827,0.9898095884931875,0.9092813528685897,0.9217849999999911
|
||||
Wadenepally,42,0.9944175574354176,0.5269873038544758,0.9457369142085821,0.4938671428571427
|
||||
|
@@ -0,0 +1,18 @@
|
||||
[project]
|
||||
name = "swot-flow-pipeline"
|
||||
version = "0.1.0"
|
||||
description = "Add your description here"
|
||||
readme = "README.md"
|
||||
requires-python = ">=3.12"
|
||||
dependencies = [
|
||||
"earthaccess>=0.18.0",
|
||||
"earthengine-api>=1.7.33",
|
||||
"geopandas>=1.1.4",
|
||||
"matplotlib>=3.11.0",
|
||||
"numpy>=2.5.1",
|
||||
"openpyxl>=3.1.5",
|
||||
"pandas>=3.0.3",
|
||||
"pyyaml>=6.0.3",
|
||||
"rasterio>=1.5.0",
|
||||
"tqdm>=4.68.3",
|
||||
]
|
||||
@@ -0,0 +1,93 @@
|
||||
#!/usr/bin/env python3
|
||||
"""river discharge estimation pipeline - complete controller.
|
||||
|
||||
Steps:
|
||||
1 download - SWOT L2 HR RiverSP_D node/reach shapefiles (NASA Earthdata)
|
||||
2 process - merge, clip to AOI, extract parameters, clean
|
||||
3 gauge - load CWC/WRIS gauge xlsx (user-placed in data/gauge/)
|
||||
4 wse - Haversine node matching + WSE validation (+/-3 h)
|
||||
5 dem - bed elevation from ALOS AW3D30 (GEE) or local GeoTIFF
|
||||
6 discharge - depth/area/hydraulic radius + Manning's equation
|
||||
7 validate - SWOT vs gauge discharge (R^2, RMSE, NSE, plots)
|
||||
8 export - all inputs/intermediates/results into one Excel workbook
|
||||
|
||||
Usage:
|
||||
python run_pipeline.py # run everything
|
||||
python run_pipeline.py --steps 2-4 # a range
|
||||
python run_pipeline.py --steps gauge,wse,discharge,validate
|
||||
python run_pipeline.py --config other.yaml
|
||||
"""
|
||||
|
||||
from __future__ import annotations
|
||||
|
||||
import argparse
|
||||
import sys
|
||||
import time
|
||||
|
||||
from swotflow.common import load_config, log, setup_logging
|
||||
|
||||
STEPS = {
|
||||
1: ("download", "swotflow.step1_download_swot"),
|
||||
2: ("process", "swotflow.step2_process_swot"),
|
||||
3: ("gauge", "swotflow.step3_load_gauge"),
|
||||
4: ("wse", "swotflow.step4_match_validate_wse"),
|
||||
5: ("dem", "swotflow.step5_bed_elevation"),
|
||||
6: ("discharge", "swotflow.step6_compute_discharge"),
|
||||
7: ("validate", "swotflow.step7_validate_discharge"),
|
||||
8: ("export", "swotflow.step8_export_excel"),
|
||||
}
|
||||
NAME_TO_NUM = {name: num for num, (name, _) in STEPS.items()}
|
||||
|
||||
|
||||
def parse_steps(spec: str) -> list[int]:
|
||||
if not spec:
|
||||
return list(STEPS)
|
||||
nums: set[int] = set()
|
||||
for part in spec.split(","):
|
||||
part = part.strip().lower()
|
||||
if "-" in part and part.replace("-", "").isdigit():
|
||||
a, b = part.split("-")
|
||||
nums.update(range(int(a), int(b) + 1))
|
||||
elif part.isdigit():
|
||||
nums.add(int(part))
|
||||
elif part in NAME_TO_NUM:
|
||||
nums.add(NAME_TO_NUM[part])
|
||||
else:
|
||||
sys.exit(f"Unknown step: {part!r}. Valid: {list(NAME_TO_NUM)}")
|
||||
bad = nums - set(STEPS)
|
||||
if bad:
|
||||
sys.exit(f"Invalid step numbers: {sorted(bad)}")
|
||||
return sorted(nums)
|
||||
|
||||
|
||||
def main() -> None:
|
||||
ap = argparse.ArgumentParser(
|
||||
description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter
|
||||
)
|
||||
ap.add_argument("--config", default="config.yaml")
|
||||
ap.add_argument("--steps", default="", help="e.g. '1-7', '2,4', 'gauge,wse'")
|
||||
args = ap.parse_args()
|
||||
|
||||
setup_logging()
|
||||
cfg = load_config(args.config)
|
||||
steps = parse_steps(args.steps)
|
||||
|
||||
log.info(
|
||||
"Pipeline: %s | steps: %s",
|
||||
cfg.get("project_name"),
|
||||
[STEPS[s][0] for s in steps],
|
||||
)
|
||||
for s in steps:
|
||||
name, module_name = STEPS[s]
|
||||
log.info("=" * 60)
|
||||
log.info("STEP %d: %s", s, name.upper())
|
||||
log.info("=" * 60)
|
||||
t0 = time.time()
|
||||
module = __import__(module_name, fromlist=["run"])
|
||||
module.run(cfg)
|
||||
log.info("Step %d done in %.1f s", s, time.time() - t0)
|
||||
log.info("Pipeline finished. Results in %s", cfg["_root"] / cfg["paths"]["outputs"])
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
@@ -0,0 +1,13 @@
|
||||
# CWC gauge stations: name must match the part after "_" in the xlsx filenames
|
||||
# e.g. "River Water Level_Arjunwad (seasonal).xlsx" -> station "Arjunwad"
|
||||
# Coordinates below are from the project presentation (nearest-node locations).
|
||||
# Dameracherla is APPROXIMATE (CWC site on the Musi river near the Wadapally/
|
||||
# Krishna confluence) - replace with the exact lat/lon from the station page
|
||||
# on India-WRIS (indiawris.gov.in).
|
||||
station,lat,lon
|
||||
Arjunwad,16.782515,74.632709
|
||||
Dameracherla,16.793,79.677
|
||||
Huvinhedigi,16.491348,76.923209
|
||||
Keesara,16.715175,80.316223
|
||||
Kurundwad,16.684692,74.602998
|
||||
Wadenepally,16.793817,80.069633
|
||||
|
Can't render this file because it contains an unexpected character in line 1 and column 54.
|
@@ -0,0 +1,2 @@
|
||||
"""SWOT-based river discharge estimation pipeline (Krishna Basin methodology)."""
|
||||
__version__ = "1.0.0"
|
||||
@@ -0,0 +1,90 @@
|
||||
"""Shared helpers: config loading, paths, Haversine, metrics."""
|
||||
from __future__ import annotations
|
||||
|
||||
import logging
|
||||
import math
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
import yaml
|
||||
|
||||
log = logging.getLogger("swotflow")
|
||||
|
||||
SWOT_FILL_VALUES = (-999999999999.0, -99999999999.0, -9999.0)
|
||||
|
||||
|
||||
def setup_logging() -> None:
|
||||
logging.basicConfig(
|
||||
level=logging.INFO,
|
||||
format="%(asctime)s | %(levelname)-7s | %(message)s",
|
||||
datefmt="%H:%M:%S",
|
||||
)
|
||||
|
||||
|
||||
def load_config(path: str | Path = "config.yaml") -> dict:
|
||||
path = Path(path)
|
||||
with open(path, encoding="utf-8") as f:
|
||||
cfg = yaml.safe_load(f)
|
||||
cfg["_root"] = path.resolve().parent
|
||||
return cfg
|
||||
|
||||
|
||||
def p(cfg: dict, *keys: str) -> Path:
|
||||
"""Resolve a path from config relative to the pipeline root, creating it."""
|
||||
node = cfg
|
||||
for k in keys:
|
||||
node = node[k]
|
||||
out = cfg["_root"] / node
|
||||
out.mkdir(parents=True, exist_ok=True)
|
||||
return out
|
||||
|
||||
|
||||
def haversine_km(lat1, lon1, lat2, lon2):
|
||||
"""Haversine great-circle distance in km. Accepts scalars or numpy arrays."""
|
||||
R = 6371.0
|
||||
lat1, lon1, lat2, lon2 = map(np.radians, (lat1, lon1, lat2, lon2))
|
||||
dlat = lat2 - lat1
|
||||
dlon = lon2 - lon1
|
||||
a = np.sin(dlat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2
|
||||
return 2 * R * np.arcsin(np.sqrt(a))
|
||||
|
||||
|
||||
def clean_fill(series: pd.Series) -> pd.Series:
|
||||
"""Replace SWOT fill values with NaN."""
|
||||
s = pd.to_numeric(series, errors="coerce")
|
||||
for fv in SWOT_FILL_VALUES:
|
||||
s = s.mask(np.isclose(s, fv), np.nan)
|
||||
return s
|
||||
|
||||
|
||||
def regression_metrics(obs: np.ndarray, est: np.ndarray) -> dict:
|
||||
"""R^2 (Pearson^2), RMSE, NSE, bias between observed and estimated arrays."""
|
||||
obs = np.asarray(obs, float)
|
||||
est = np.asarray(est, float)
|
||||
mask = np.isfinite(obs) & np.isfinite(est)
|
||||
obs, est = obs[mask], est[mask]
|
||||
n = len(obs)
|
||||
out = {"n": n, "r2": np.nan, "rmse": np.nan, "nse": np.nan, "bias": np.nan}
|
||||
if n < 2:
|
||||
return out
|
||||
if np.std(obs) > 0 and np.std(est) > 0:
|
||||
out["r2"] = float(np.corrcoef(obs, est)[0, 1] ** 2)
|
||||
out["rmse"] = float(math.sqrt(np.mean((est - obs) ** 2)))
|
||||
denom = np.sum((obs - obs.mean()) ** 2)
|
||||
if denom > 0:
|
||||
out["nse"] = float(1 - np.sum((est - obs) ** 2) / denom)
|
||||
out["bias"] = float(np.mean(est - obs))
|
||||
return out
|
||||
|
||||
|
||||
def nearest_in_time(target: pd.Timestamp, times: pd.Series, values: pd.Series,
|
||||
window_hours: float):
|
||||
"""Value whose timestamp is nearest to `target` within +/- window. None if none."""
|
||||
dt = (times - target).abs()
|
||||
if dt.empty:
|
||||
return None, None
|
||||
i = dt.idxmin()
|
||||
if dt.loc[i] <= pd.Timedelta(hours=window_hours):
|
||||
return values.loc[i], times.loc[i]
|
||||
return None, None
|
||||
@@ -0,0 +1,84 @@
|
||||
"""Step 1 - Download SWOT L2 HR RiverSP (Version D) node + reach shapefiles.
|
||||
|
||||
Source: NASA Earthdata (PO.DAAC), via the `earthaccess` library.
|
||||
Requires an Earthdata login: set EARTHDATA_USERNAME / EARTHDATA_PASSWORD
|
||||
environment variables, or let earthaccess prompt once and persist to ~/.netrc.
|
||||
|
||||
Granules are zipped shapefiles, one per pass/continent, e.g.
|
||||
SWOT_L2_HR_RiverSP_Node_..._AS_...zip (points - WSE, width, lat/lon)
|
||||
SWOT_L2_HR_RiverSP_Reach_..._AS_...zip (polylines - slope)
|
||||
They are unzipped into separate node/ and reach/ folders.
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
import zipfile
|
||||
from pathlib import Path
|
||||
|
||||
from tqdm import tqdm
|
||||
|
||||
from .common import log, p
|
||||
|
||||
|
||||
def run(cfg: dict) -> None:
|
||||
import earthaccess
|
||||
|
||||
raw_dir = p(cfg, "paths", "swot_raw")
|
||||
node_dir = p(cfg, "paths", "swot_nodes")
|
||||
reach_dir = p(cfg, "paths", "swot_reaches")
|
||||
aoi = cfg["aoi"]
|
||||
cc = cfg["swot"].get("continent_code")
|
||||
|
||||
log.info("Logging in to NASA Earthdata ...")
|
||||
earthaccess.login(persist=True)
|
||||
|
||||
log.info("Searching %s granules %s -> %s ...", cfg["swot"]["short_name"],
|
||||
cfg["time"]["start"], cfg["time"]["end"])
|
||||
results = earthaccess.search_data(
|
||||
short_name=cfg["swot"]["short_name"],
|
||||
temporal=(cfg["time"]["start"], cfg["time"]["end"]),
|
||||
bounding_box=(aoi["lon_min"], aoi["lat_min"], aoi["lon_max"], aoi["lat_max"]),
|
||||
)
|
||||
if cc:
|
||||
results = [g for g in results if f"_{cc}_" in _granule_name(g)]
|
||||
log.info("Found %d granules (Node + Reach).", len(results))
|
||||
if not results:
|
||||
log.warning("Nothing found - check AOI/date range/short_name.")
|
||||
return
|
||||
|
||||
# Download in batches with an overall progress bar (earthaccess shows
|
||||
# per-file bars; this one tracks total progress across all granules).
|
||||
files = []
|
||||
batch = 100
|
||||
with tqdm(total=len(results), desc="Downloading SWOT granules",
|
||||
unit="granule") as bar:
|
||||
for i in range(0, len(results), batch):
|
||||
files += earthaccess.download(results[i:i + batch], str(raw_dir))
|
||||
bar.update(len(results[i:i + batch]))
|
||||
log.info("Downloaded %d files to %s", len(files), raw_dir)
|
||||
|
||||
n_node = n_reach = 0
|
||||
for f in map(Path, tqdm(files, desc="Extracting granules", unit="zip")):
|
||||
if f.suffix.lower() != ".zip":
|
||||
continue
|
||||
name = f.name
|
||||
if "RiverSP_Node" in name:
|
||||
dest, n_node = node_dir, n_node + 1
|
||||
elif "RiverSP_Reach" in name:
|
||||
dest, n_reach = reach_dir, n_reach + 1
|
||||
else:
|
||||
continue
|
||||
if (dest / f"{f.stem}.shp").exists(): # already extracted on a previous run
|
||||
continue
|
||||
try:
|
||||
with zipfile.ZipFile(f) as z:
|
||||
z.extractall(dest)
|
||||
except zipfile.BadZipFile:
|
||||
log.warning("Bad zip skipped: %s", name)
|
||||
log.info("Extracted %d node and %d reach granules.", n_node, n_reach)
|
||||
|
||||
|
||||
def _granule_name(granule) -> str:
|
||||
try:
|
||||
return granule["umm"]["GranuleUR"]
|
||||
except Exception:
|
||||
return str(granule)
|
||||
@@ -0,0 +1,95 @@
|
||||
"""Step 2 - Merge SWOT node & reach shapefiles, clip to AOI, extract & clean.
|
||||
|
||||
Per the project methodology:
|
||||
- merge all node granules / all reach granules
|
||||
- clip to the study-area bounding box
|
||||
- extract: node_id, reach_id, time, lat, lon, wse, width (nodes)
|
||||
reach_id, time, slope (reaches)
|
||||
- clean: remove fill values / NaN / negative or unrealistic values / duplicates
|
||||
|
||||
Outputs (CSV, in paths.processed):
|
||||
nodes_clean.csv, reaches_clean.csv
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
from pathlib import Path
|
||||
|
||||
import pandas as pd
|
||||
from tqdm import tqdm
|
||||
|
||||
from .common import clean_fill, log, p
|
||||
|
||||
NODE_COLS = ["node_id", "reach_id", "time_str", "lat", "lon", "wse", "width",
|
||||
"node_q"]
|
||||
REACH_COLS = ["reach_id", "time_str", "slope", "reach_q"]
|
||||
|
||||
|
||||
def _read_shapefiles(folder: Path, bbox: tuple, columns: list[str]) -> pd.DataFrame:
|
||||
import geopandas as gpd
|
||||
|
||||
shps = sorted(folder.rglob("*.shp"))
|
||||
if not shps:
|
||||
raise FileNotFoundError(f"No shapefiles in {folder} - run step 1 first.")
|
||||
log.info("Reading %d shapefiles from %s ...", len(shps), folder.name)
|
||||
frames = []
|
||||
for shp in tqdm(shps, desc=f"Merging {folder.name}", unit="file"):
|
||||
try:
|
||||
gdf = gpd.read_file(shp, bbox=bbox)
|
||||
except Exception as e: # corrupt granule should not kill the run
|
||||
log.warning("Skipping %s (%s)", shp.name, e)
|
||||
continue
|
||||
if gdf.empty:
|
||||
continue
|
||||
keep = [c for c in columns if c in gdf.columns]
|
||||
frames.append(pd.DataFrame(gdf[keep]))
|
||||
if not frames:
|
||||
return pd.DataFrame(columns=columns)
|
||||
return pd.concat(frames, ignore_index=True)
|
||||
|
||||
|
||||
def run(cfg: dict) -> None:
|
||||
aoi = cfg["aoi"]
|
||||
bbox = (aoi["lon_min"], aoi["lat_min"], aoi["lon_max"], aoi["lat_max"])
|
||||
out_dir = p(cfg, "paths", "processed")
|
||||
|
||||
# ---------------- nodes ----------------
|
||||
nodes = _read_shapefiles(p(cfg, "paths", "swot_nodes"), bbox, NODE_COLS)
|
||||
log.info("Merged nodes: %d records", len(nodes))
|
||||
|
||||
for c in ("lat", "lon", "wse", "width"):
|
||||
nodes[c] = clean_fill(nodes[c])
|
||||
nodes["time"] = pd.to_datetime(nodes["time_str"], errors="coerce", utc=True)
|
||||
nodes = nodes.dropna(subset=["lat", "lon", "wse", "width", "time"])
|
||||
nodes = nodes[(nodes["width"] > 0) & (nodes["wse"] > -100) & (nodes["wse"] < 9000)]
|
||||
if "node_q" in nodes.columns:
|
||||
max_q = cfg.get("quality", {}).get("max_node_q", 2)
|
||||
before = len(nodes)
|
||||
nodes = nodes[pd.to_numeric(nodes["node_q"], errors="coerce") <= max_q]
|
||||
log.info("Quality filter node_q <= %d: dropped %d of %d records.",
|
||||
max_q, before - len(nodes), before)
|
||||
nodes = nodes.drop_duplicates(subset=["node_id", "time_str"])
|
||||
nodes["node_id"] = nodes["node_id"].astype(str)
|
||||
nodes["reach_id"] = nodes["reach_id"].astype(str)
|
||||
log.info("Clean nodes: %d records, %d unique nodes",
|
||||
len(nodes), nodes["node_id"].nunique())
|
||||
nodes.drop(columns=["time_str"]).to_csv(out_dir / "nodes_clean.csv", index=False)
|
||||
|
||||
# ---------------- reaches ----------------
|
||||
reaches = _read_shapefiles(p(cfg, "paths", "swot_reaches"), bbox, REACH_COLS)
|
||||
log.info("Merged reaches: %d records", len(reaches))
|
||||
|
||||
reaches["slope"] = clean_fill(reaches["slope"])
|
||||
reaches["time"] = pd.to_datetime(reaches["time_str"], errors="coerce", utc=True)
|
||||
reaches = reaches.dropna(subset=["slope", "time"])
|
||||
reaches = reaches[reaches["slope"] > 0] # need positive slope for Manning
|
||||
if "reach_q" in reaches.columns:
|
||||
max_q = cfg.get("quality", {}).get("max_node_q", 2)
|
||||
before = len(reaches)
|
||||
reaches = reaches[pd.to_numeric(reaches["reach_q"], errors="coerce") <= max_q]
|
||||
log.info("Quality filter reach_q <= %d: dropped %d of %d records.",
|
||||
max_q, before - len(reaches), before)
|
||||
reaches = reaches.drop_duplicates(subset=["reach_id", "time_str"])
|
||||
reaches["reach_id"] = reaches["reach_id"].astype(str)
|
||||
log.info("Clean reaches: %d records, %d unique reaches",
|
||||
len(reaches), reaches["reach_id"].nunique())
|
||||
reaches.drop(columns=["time_str"]).to_csv(out_dir / "reaches_clean.csv", index=False)
|
||||
@@ -0,0 +1,97 @@
|
||||
"""Step 3 - Load CWC / India-WRIS gauge xlsx files placed by the user.
|
||||
|
||||
Expected folders (see config.yaml):
|
||||
data/gauge/wse/ River Water Level_<Station>.xlsx
|
||||
data/gauge/discharge/ River Water Discharge_<Station>.xlsx
|
||||
|
||||
File format (WRIS export): first row is a header row -
|
||||
Data Type Code | Data Type Description | Data Time | Data Value | Unit
|
||||
|
||||
Station name = filename part after the first "_", minus extension and any
|
||||
trailing "(seasonal)" tag. Must match `station` in stations.csv.
|
||||
|
||||
Outputs: gauge_wse.csv, gauge_discharge.csv (station, time, value)
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
import re
|
||||
from pathlib import Path
|
||||
|
||||
import pandas as pd
|
||||
from tqdm import tqdm
|
||||
|
||||
from .common import log, p
|
||||
|
||||
|
||||
def _station_from_filename(path: Path) -> str:
|
||||
stem = path.stem
|
||||
name = stem.split("_", 1)[1] if "_" in stem else stem
|
||||
name = re.sub(r"\s*\(.*?\)\s*$", "", name) # drop "(seasonal)" etc.
|
||||
return name.strip()
|
||||
|
||||
|
||||
def _load_folder(folder: Path, label: str) -> pd.DataFrame:
|
||||
files = sorted(list(folder.glob("*.xlsx")) + list(folder.glob("*.xls")))
|
||||
if not files:
|
||||
log.warning("No %s xlsx files found in %s", label, folder)
|
||||
return pd.DataFrame(columns=["station", "time", "value"])
|
||||
frames = []
|
||||
for f in tqdm(files, desc=f"Loading {label} xlsx", unit="file"):
|
||||
try:
|
||||
df = pd.read_excel(f, header=None)
|
||||
except Exception as e:
|
||||
log.warning("Could not read %s (%s)", f.name, e)
|
||||
continue
|
||||
# locate the header row containing "Data Time"
|
||||
hdr = None
|
||||
for i in range(min(10, len(df))):
|
||||
if df.iloc[i].astype(str).str.contains("Data Time", case=False).any():
|
||||
hdr = i
|
||||
break
|
||||
if hdr is None:
|
||||
log.warning("No 'Data Time' header in %s - skipped", f.name)
|
||||
continue
|
||||
df.columns = df.iloc[hdr]
|
||||
df = df.iloc[hdr + 1:]
|
||||
tcol = next(c for c in df.columns if "time" in str(c).lower())
|
||||
vcol = next(c for c in df.columns if "value" in str(c).lower())
|
||||
# WRIS files can hold several series (e.g. WL in MSL *and* local gauge
|
||||
# datum). Prefer the MSL series; otherwise keep the most common one.
|
||||
dcol = next((c for c in df.columns if "description" in str(c).lower()), None)
|
||||
if dcol is not None and df[dcol].notna().any():
|
||||
desc = df[dcol].astype(str)
|
||||
if desc.str.contains("MSL", case=False).any():
|
||||
df = df[desc.str.contains("MSL", case=False)]
|
||||
elif desc.nunique() > 1:
|
||||
df = df[desc == desc.mode().iloc[0]]
|
||||
out = pd.DataFrame({
|
||||
"station": _station_from_filename(f),
|
||||
"time": pd.to_datetime(df[tcol], errors="coerce"),
|
||||
"value": pd.to_numeric(df[vcol], errors="coerce"),
|
||||
}).dropna(subset=["time", "value"])
|
||||
frames.append(out)
|
||||
log.info(" %-45s %6d records (%s)", f.name, len(out), label)
|
||||
if not frames:
|
||||
return pd.DataFrame(columns=["station", "time", "value"])
|
||||
all_df = pd.concat(frames, ignore_index=True)
|
||||
# WRIS times are IST (UTC+05:30); convert to UTC to match SWOT times
|
||||
all_df["time"] = (all_df["time"].dt.tz_localize("Asia/Kolkata", ambiguous="NaT",
|
||||
nonexistent="NaT")
|
||||
.dt.tz_convert("UTC"))
|
||||
all_df = all_df.dropna(subset=["time"]).drop_duplicates(["station", "time"])
|
||||
return all_df
|
||||
|
||||
|
||||
def run(cfg: dict) -> None:
|
||||
out_dir = p(cfg, "paths", "processed")
|
||||
|
||||
wse = _load_folder(p(cfg, "gauge", "wse_dir"), "water level")
|
||||
q = _load_folder(p(cfg, "gauge", "discharge_dir"), "discharge")
|
||||
|
||||
wse = wse[wse["value"] > -100] # sanity
|
||||
q = q[q["value"] >= 0]
|
||||
|
||||
wse.to_csv(out_dir / "gauge_wse.csv", index=False)
|
||||
q.to_csv(out_dir / "gauge_discharge.csv", index=False)
|
||||
log.info("Gauge WSE: %d records / %d stations; Discharge: %d records / %d stations",
|
||||
len(wse), wse["station"].nunique(), len(q), q["station"].nunique())
|
||||
@@ -0,0 +1,154 @@
|
||||
"""Step 4 - Spatial + temporal matching, WSE validation.
|
||||
|
||||
Per methodology:
|
||||
- for each gauge station, find the nearest SWOT node (Haversine distance),
|
||||
within matching.max_distance_km
|
||||
- for every SWOT observation at that node, find the gauge WSE recorded
|
||||
within +/- matching.wse_time_window_hours (default 3 h)
|
||||
- validate SWOT WSE against gauge WSE (R^2, RMSE, bias) and plot
|
||||
|
||||
Outputs:
|
||||
processed/station_nodes.csv - station -> nearest node/reach mapping
|
||||
processed/matched_wse.csv - paired SWOT vs gauge WSE
|
||||
outputs/wse_metrics.csv, outputs/plots/wse_<station>.png
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
from tqdm import tqdm
|
||||
|
||||
from .common import haversine_km, log, nearest_in_time, p, regression_metrics
|
||||
|
||||
|
||||
def _load_stations(cfg) -> pd.DataFrame:
|
||||
path = cfg["_root"] / cfg["gauge"]["stations_csv"]
|
||||
df = pd.read_csv(path, comment="#")
|
||||
df["station"] = df["station"].astype(str).str.strip()
|
||||
return df
|
||||
|
||||
|
||||
def nearest_nodes(cfg, nodes: pd.DataFrame, stations: pd.DataFrame) -> pd.DataFrame:
|
||||
"""Map each gauge station to its nearest SWOT node."""
|
||||
uniq = (nodes.groupby("node_id")
|
||||
.agg(lat=("lat", "median"), lon=("lon", "median"),
|
||||
reach_id=("reach_id", "first"))
|
||||
.reset_index())
|
||||
rows = []
|
||||
for _, st in stations.iterrows():
|
||||
d = haversine_km(st["lat"], st["lon"], uniq["lat"].values, uniq["lon"].values)
|
||||
i = int(np.argmin(d))
|
||||
row = {"station": st["station"], "gauge_lat": st["lat"], "gauge_lon": st["lon"],
|
||||
"node_id": uniq.iloc[i]["node_id"], "reach_id": uniq.iloc[i]["reach_id"],
|
||||
"node_lat": uniq.iloc[i]["lat"], "node_lon": uniq.iloc[i]["lon"],
|
||||
"distance_km": float(d[i])}
|
||||
if row["distance_km"] > cfg["matching"]["max_distance_km"]:
|
||||
log.warning("%s: nearest node is %.1f km away (> %.1f km) - excluded",
|
||||
st["station"], row["distance_km"],
|
||||
cfg["matching"]["max_distance_km"])
|
||||
continue
|
||||
rows.append(row)
|
||||
log.info("%s -> node %s (%.2f km, reach %s)",
|
||||
st["station"], row["node_id"], row["distance_km"], row["reach_id"])
|
||||
return pd.DataFrame(rows)
|
||||
|
||||
|
||||
def run(cfg: dict) -> None:
|
||||
import matplotlib
|
||||
matplotlib.use("Agg")
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
proc = p(cfg, "paths", "processed")
|
||||
out = p(cfg, "paths", "outputs")
|
||||
plots = out / "plots"
|
||||
plots.mkdir(exist_ok=True)
|
||||
|
||||
nodes = pd.read_csv(proc / "nodes_clean.csv", parse_dates=["time"],
|
||||
dtype={"node_id": str, "reach_id": str})
|
||||
gauge = pd.read_csv(proc / "gauge_wse.csv", parse_dates=["time"])
|
||||
stations = _load_stations(cfg)
|
||||
|
||||
mapping = nearest_nodes(cfg, nodes, stations)
|
||||
mapping.to_csv(proc / "station_nodes.csv", index=False)
|
||||
if mapping.empty:
|
||||
log.warning("No station-node matches; check stations.csv and AOI.")
|
||||
return
|
||||
|
||||
window = cfg["matching"]["wse_time_window_hours"]
|
||||
pairs = []
|
||||
for _, m in mapping.iterrows():
|
||||
sn = nodes[nodes["node_id"] == m["node_id"]].sort_values("time")
|
||||
gw = gauge[gauge["station"] == m["station"]].reset_index(drop=True)
|
||||
if gw.empty:
|
||||
log.warning("%s: no gauge WSE file loaded - skipping WSE validation",
|
||||
m["station"])
|
||||
continue
|
||||
for _, obs in tqdm(sn.iterrows(), total=len(sn), unit="obs", leave=False,
|
||||
desc=f"Matching WSE - {m['station']}"):
|
||||
val, t = nearest_in_time(obs["time"], gw["time"], gw["value"], window)
|
||||
if val is not None:
|
||||
pairs.append({"station": m["station"], "node_id": m["node_id"],
|
||||
"reach_id": m["reach_id"], "time_swot": obs["time"],
|
||||
"time_gauge": t, "swot_wse": obs["wse"],
|
||||
"gauge_wse": val, "width": obs["width"]})
|
||||
matched = pd.DataFrame(pairs)
|
||||
log.info("Matched WSE pairs: %d", len(matched))
|
||||
|
||||
# Robust outlier rejection: flag pairs whose (SWOT - gauge) residual is
|
||||
# far from the station's median residual. Removes bad satellite obs
|
||||
# (dark water / layover) without touching the seasonal signal.
|
||||
k = cfg["matching"].get("wse_outlier_mad_k", 0)
|
||||
if not matched.empty:
|
||||
matched["outlier"] = False
|
||||
if k:
|
||||
for st, g in matched.groupby("station"):
|
||||
resid = g["swot_wse"] - g["gauge_wse"]
|
||||
med = resid.median()
|
||||
mad = max((resid - med).abs().median(), 0.10) # floor: 10 cm
|
||||
bad = (resid - med).abs() > k * mad
|
||||
matched.loc[g.index[bad], "outlier"] = True
|
||||
if bad.any():
|
||||
log.info("%s: flagged %d/%d WSE outliers (>%d*MAD).",
|
||||
st, int(bad.sum()), len(g), k)
|
||||
# datum / location sanity check
|
||||
if abs(med) > 20:
|
||||
log.warning("%s: median SWOT-gauge offset is %.1f m - the "
|
||||
"station coordinates or the gauge datum are "
|
||||
"almost certainly wrong. Check stations.csv "
|
||||
"and the WRIS export.", st, med)
|
||||
matched.to_csv(proc / "matched_wse.csv", index=False)
|
||||
|
||||
metrics = []
|
||||
for st, g in matched[~matched["outlier"]].groupby("station"):
|
||||
mt = regression_metrics(g["gauge_wse"], g["swot_wse"])
|
||||
mt["station"] = st
|
||||
metrics.append(mt)
|
||||
|
||||
# If SWOT and the gauge use different vertical datums (constant offset,
|
||||
# e.g. EGM2008 vs a local gauge zero), shift SWOT by the median offset
|
||||
# in the PLOTS so shapes are comparable. Metrics keep the raw bias.
|
||||
off = float((g["swot_wse"] - g["gauge_wse"]).median())
|
||||
shift = off if abs(off) > 1.0 else 0.0
|
||||
swot_lbl = f"SWOT − {shift:.1f} m (datum shift)" if shift else "SWOT"
|
||||
swot_plot = g["swot_wse"] - shift
|
||||
|
||||
fig, ax = plt.subplots(1, 2, figsize=(11, 4.2))
|
||||
ax[0].scatter(g["gauge_wse"], swot_plot, s=18, alpha=0.7)
|
||||
lims = [min(g["gauge_wse"].min(), swot_plot.min()),
|
||||
max(g["gauge_wse"].max(), swot_plot.max())]
|
||||
ax[0].plot(lims, lims, "k--", lw=1)
|
||||
ax[0].set_xlabel("Gauge WSE (m)")
|
||||
ax[0].set_ylabel(swot_lbl + " (m)")
|
||||
ax[0].set_title(f"{st} R²={mt['r2']:.3f} n={mt['n']}")
|
||||
gg = g.sort_values("time_swot")
|
||||
ax[1].plot(gg["time_swot"], gg["gauge_wse"], "o-", ms=4, label="Gauge")
|
||||
ax[1].plot(gg["time_swot"], gg["swot_wse"] - shift, "s--", ms=4,
|
||||
label=swot_lbl)
|
||||
ax[1].set_ylabel("WSE (m)"); ax[1].legend(); ax[1].set_title("Time series")
|
||||
fig.autofmt_xdate(); fig.tight_layout()
|
||||
fig.savefig(plots / f"wse_{st}.png", dpi=150); plt.close(fig)
|
||||
|
||||
if metrics:
|
||||
mdf = pd.DataFrame(metrics)[["station", "n", "r2", "rmse", "nse", "bias"]]
|
||||
mdf.to_csv(out / "wse_metrics.csv", index=False)
|
||||
log.info("WSE validation metrics:\n%s", mdf.to_string(index=False))
|
||||
@@ -0,0 +1,119 @@
|
||||
"""Step 5 - Bed elevation from ALOS World 3D 30m (AW3D30 v3.2, DSM band).
|
||||
|
||||
Primary source: Google Earth Engine (JAXA/ALOS/AW3D30/V3_2), as in the study.
|
||||
Requires one-time auth: earthengine authenticate (and gee.project in config)
|
||||
Fallback: a local DEM GeoTIFF (dem.local_tif) sampled with rasterio.
|
||||
|
||||
Sampling points:
|
||||
hydraulics.scope = "matched" -> only gauge-matched nodes (station_nodes.csv)
|
||||
hydraulics.scope = "all" -> every unique node in the AOI
|
||||
|
||||
Output: processed/node_bed_elevation.csv (node_id, lat, lon, bed_elev)
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
import pandas as pd
|
||||
from tqdm import tqdm
|
||||
|
||||
from .common import log, p
|
||||
|
||||
|
||||
def _points_to_sample(cfg, proc) -> pd.DataFrame:
|
||||
scope = cfg["hydraulics"].get("scope", "matched")
|
||||
if scope == "all":
|
||||
nodes = pd.read_csv(proc / "nodes_clean.csv", dtype={"node_id": str})
|
||||
pts = (nodes.groupby("node_id")
|
||||
.agg(lat=("lat", "median"), lon=("lon", "median"))
|
||||
.reset_index())
|
||||
else:
|
||||
m = pd.read_csv(proc / "station_nodes.csv", dtype={"node_id": str})
|
||||
pts = m.rename(columns={"node_lat": "lat", "node_lon": "lon"})[
|
||||
["node_id", "lat", "lon"]].drop_duplicates("node_id")
|
||||
log.info("Sampling bed elevation at %d node locations (scope=%s).", len(pts), scope)
|
||||
return pts
|
||||
|
||||
|
||||
def _sample_gee(cfg, pts: pd.DataFrame) -> pd.DataFrame:
|
||||
import ee
|
||||
|
||||
project = cfg["gee"].get("project")
|
||||
ee.Initialize(project=project)
|
||||
dem = (ee.ImageCollection(cfg["gee"]["dem_collection"])
|
||||
.select(cfg["gee"]["band"]).mosaic())
|
||||
scale = cfg["gee"].get("scale_m", 30)
|
||||
buffer_m = cfg["gee"].get("buffer_m", 0)
|
||||
|
||||
out = []
|
||||
chunk = 500
|
||||
for i in tqdm(range(0, len(pts), chunk), desc="GEE DEM sampling", unit="batch"):
|
||||
sub = pts.iloc[i:i + chunk]
|
||||
feats = []
|
||||
for r in sub.itertuples():
|
||||
geom = ee.Geometry.Point([r.lon, r.lat])
|
||||
if buffer_m:
|
||||
geom = geom.buffer(buffer_m) # min within buffer ~ channel floor
|
||||
feats.append(ee.Feature(geom, {"node_id": r.node_id}))
|
||||
sampled = dem.reduceRegions(collection=ee.FeatureCollection(feats),
|
||||
reducer=ee.Reducer.min(),
|
||||
scale=scale).getInfo()
|
||||
for f in sampled["features"]:
|
||||
props = f["properties"]
|
||||
out.append({"node_id": str(props["node_id"]),
|
||||
"bed_elev": props.get("min")})
|
||||
return pts.merge(pd.DataFrame(out), on="node_id", how="left")
|
||||
|
||||
|
||||
def _sample_local(cfg, pts: pd.DataFrame) -> pd.DataFrame:
|
||||
import rasterio
|
||||
|
||||
import numpy as np
|
||||
|
||||
tif = cfg["_root"] / cfg["dem"]["local_tif"]
|
||||
buffer_m = cfg["gee"].get("buffer_m", 0)
|
||||
log.info("Sampling local DEM %s (min within %dm) ...", tif, buffer_m)
|
||||
vals = []
|
||||
with rasterio.open(tif) as src:
|
||||
# buffer radius in pixels (approx; DEM assumed ~30 m in degrees or meters)
|
||||
px = max(int(round(buffer_m / 30)), 0)
|
||||
for lon, lat in zip(pts["lon"], pts["lat"]):
|
||||
row, col = src.index(lon, lat)
|
||||
win = rasterio.windows.Window(col - px, row - px, 2 * px + 1, 2 * px + 1)
|
||||
block = src.read(1, window=win, boundless=True,
|
||||
fill_value=src.nodata if src.nodata is not None else 0)
|
||||
block = block.astype(float)
|
||||
if src.nodata is not None:
|
||||
block[block == src.nodata] = np.nan
|
||||
vals.append(np.nanmin(block) if np.isfinite(block).any() else np.nan)
|
||||
pts = pts.copy()
|
||||
pts["bed_elev"] = vals
|
||||
return pts
|
||||
|
||||
|
||||
def run(cfg: dict) -> None:
|
||||
proc = p(cfg, "paths", "processed")
|
||||
pts = _points_to_sample(cfg, proc)
|
||||
|
||||
if cfg["gee"].get("project"):
|
||||
try:
|
||||
result = _sample_gee(cfg, pts)
|
||||
except Exception as e:
|
||||
log.warning("GEE sampling failed (%s); trying local DEM fallback.", e)
|
||||
result = None
|
||||
else:
|
||||
log.info("gee.project not set - skipping GEE.")
|
||||
result = None
|
||||
|
||||
if result is None:
|
||||
if cfg["dem"].get("local_tif"):
|
||||
result = _sample_local(cfg, pts)
|
||||
else:
|
||||
raise RuntimeError(
|
||||
"No DEM source available. Set gee.project (after "
|
||||
"`earthengine authenticate`) or dem.local_tif in config.yaml.")
|
||||
|
||||
result["bed_elev"] = pd.to_numeric(result["bed_elev"], errors="coerce")
|
||||
n_bad = result["bed_elev"].isna().sum()
|
||||
if n_bad:
|
||||
log.warning("%d nodes have no DEM value.", n_bad)
|
||||
result.to_csv(proc / "node_bed_elevation.csv", index=False)
|
||||
log.info("Bed elevation written for %d nodes.", result["bed_elev"].notna().sum())
|
||||
@@ -0,0 +1,138 @@
|
||||
"""Step 6 - Hydraulic computation: discharge via Manning's equation.
|
||||
|
||||
Per methodology:
|
||||
depth D = WSE - bed elevation (SWOT WSE, ALOS DEM bed)
|
||||
area A = width x depth (rectangular section)
|
||||
wetted perimeter P = width + 2*depth
|
||||
hydraulic radius R = A / P (~ depth for wide rivers)
|
||||
slope S : from the SWOT reach record of the SAME pass (same time),
|
||||
falling back to the reach's median slope
|
||||
discharge Q = (1/n) * A * R^(2/3) * S^(1/2) [m^3/s]
|
||||
|
||||
Scope (config hydraulics.scope):
|
||||
"matched" -> full SWOT time series at gauge-matched nodes
|
||||
"all" -> every node observation in the AOI
|
||||
|
||||
Output: processed/swot_discharge.csv
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
from .common import log, p
|
||||
|
||||
|
||||
def run(cfg: dict) -> None:
|
||||
proc = p(cfg, "paths", "processed")
|
||||
n_manning = float(cfg["hydraulics"]["mannings_n"])
|
||||
min_depth = float(cfg["hydraulics"]["min_depth_m"])
|
||||
scope = cfg["hydraulics"].get("scope", "matched")
|
||||
|
||||
nodes = pd.read_csv(proc / "nodes_clean.csv", parse_dates=["time"],
|
||||
dtype={"node_id": str, "reach_id": str})
|
||||
reaches = pd.read_csv(proc / "reaches_clean.csv", parse_dates=["time"],
|
||||
dtype={"reach_id": str})
|
||||
bed = pd.read_csv(proc / "node_bed_elevation.csv", dtype={"node_id": str})[
|
||||
["node_id", "bed_elev"]]
|
||||
bed["bed_elev"] = bed["bed_elev"].astype(float) # GEE DSM can come back int
|
||||
|
||||
if scope == "matched":
|
||||
mapping = pd.read_csv(proc / "station_nodes.csv", dtype={"node_id": str})
|
||||
nodes = nodes[nodes["node_id"].isin(mapping["node_id"])]
|
||||
nodes = nodes.merge(mapping[["node_id", "station"]], on="node_id", how="left")
|
||||
else:
|
||||
nodes = nodes.assign(station=pd.NA)
|
||||
log.info("Computing discharge for %d SWOT observations (scope=%s, n=%.3f).",
|
||||
len(nodes), scope, n_manning)
|
||||
|
||||
# Per-node robust outlier filter on SWOT WSE (bad satellite obs would
|
||||
# otherwise become absurd depths). Uses only SWOT data, so no circularity.
|
||||
k = cfg.get("quality", {}).get("swot_wse_mad_k", 0)
|
||||
if k:
|
||||
med = nodes.groupby("node_id")["wse"].transform("median")
|
||||
mad = ((nodes["wse"] - med).abs()
|
||||
.groupby(nodes["node_id"]).transform("median").clip(lower=0.10))
|
||||
before = len(nodes)
|
||||
nodes = nodes[(nodes["wse"] - med).abs() <= k * mad]
|
||||
log.info("SWOT WSE outlier filter (%d*MAD): dropped %d of %d obs.",
|
||||
k, before - len(nodes), before)
|
||||
|
||||
df = nodes.merge(bed, on="node_id", how="left").dropna(subset=["bed_elev"])
|
||||
|
||||
# Bed sanity ("auto" mode): the ALOS DSM is unreliable at the channel -
|
||||
# it can sit ABOVE the low-water surface (bank pixels -> all depths
|
||||
# negative) or far BELOW it (voids/pits -> huge permanent depths and a
|
||||
# discharge floor even when the gauge reads zero). Clamp the implied
|
||||
# low-water baseline depth (P5(WSE) - bed) into a plausible band.
|
||||
if cfg["hydraulics"].get("bed_mode", "auto") == "auto":
|
||||
lo = float(cfg["hydraulics"].get("baseline_depth_min_m", 0.5))
|
||||
hi = float(cfg["hydraulics"].get("baseline_depth_max_m", 2.0))
|
||||
p5 = df.groupby("node_id")["wse"].transform(lambda s: s.quantile(0.05))
|
||||
baseline = p5 - df["bed_elev"]
|
||||
clamped = baseline.clip(lo, hi)
|
||||
changed = ~np.isclose(baseline, clamped)
|
||||
if changed.any():
|
||||
fixed = df.loc[changed, "node_id"].unique()
|
||||
log.warning("DEM baseline depth outside [%.1f, %.1f] m at %d "
|
||||
"node(s) %s - clamping effective bed there.",
|
||||
lo, hi, len(fixed), list(fixed))
|
||||
df["bed_elev"] = p5 - clamped
|
||||
|
||||
# Slope: per-pass SWOT slopes are very noisy (CV can exceed 100%), so the
|
||||
# default is the reach's median slope over all passes ("median").
|
||||
# "pass" uses the same-pass value with median fallback.
|
||||
slope_source = cfg["hydraulics"].get("slope_source", "median")
|
||||
med_slope = reaches.groupby("reach_id")["slope"].median().rename("slope_median")
|
||||
df = df.merge(med_slope, on="reach_id", how="left")
|
||||
if slope_source == "pass":
|
||||
df = df.sort_values("time")
|
||||
reaches = reaches.sort_values("time")
|
||||
df = pd.merge_asof(df, reaches.rename(columns={"slope": "slope_pass"})[
|
||||
["reach_id", "time", "slope_pass"]],
|
||||
on="time", by="reach_id",
|
||||
tolerance=pd.Timedelta("30min"), direction="nearest")
|
||||
df["slope"] = df["slope_pass"].fillna(df["slope_median"])
|
||||
else:
|
||||
df["slope"] = df["slope_median"]
|
||||
df = df.dropna(subset=["slope"])
|
||||
df = df[df["slope"] > 0]
|
||||
|
||||
# Width: same story - per-pass widths are noisy and barely correlate with
|
||||
# flow; default is the node's median width ("median" vs "pass").
|
||||
if cfg["hydraulics"].get("width_source", "median") == "median":
|
||||
df["width"] = df.groupby("node_id")["width"].transform("median")
|
||||
|
||||
# hydraulic parameters
|
||||
df["depth"] = df["wse"] - df["bed_elev"]
|
||||
n_before = len(df)
|
||||
df = df[df["depth"] >= min_depth]
|
||||
log.info("Dropped %d observations with depth < %.2f m.", n_before - len(df), min_depth)
|
||||
|
||||
# Manning's n: global default, with optional per-station overrides so
|
||||
# researchers can experiment with different roughness values per site.
|
||||
n_over = cfg["hydraulics"].get("mannings_n_per_station") or {}
|
||||
df["mannings_n"] = df["station"].map(n_over).fillna(n_manning) \
|
||||
if n_over else n_manning
|
||||
if n_over:
|
||||
used = {s: n for s, n in n_over.items()
|
||||
if s in set(df["station"].dropna())}
|
||||
log.info("Per-station Manning's n overrides: %s (default %.3f)",
|
||||
used, n_manning)
|
||||
|
||||
df["area"] = df["width"] * df["depth"]
|
||||
df["wetted_perimeter"] = df["width"] + 2 * df["depth"]
|
||||
df["hydraulic_radius"] = df["area"] / df["wetted_perimeter"]
|
||||
df["discharge"] = ((1.0 / df["mannings_n"]) * df["area"]
|
||||
* np.power(df["hydraulic_radius"], 2.0 / 3.0)
|
||||
* np.sqrt(df["slope"]))
|
||||
|
||||
cols = ["station", "node_id", "reach_id", "time", "lat", "lon", "wse",
|
||||
"bed_elev", "depth", "width", "area", "wetted_perimeter",
|
||||
"hydraulic_radius", "slope", "mannings_n", "discharge"]
|
||||
df[cols].to_csv(proc / "swot_discharge.csv", index=False)
|
||||
log.info("Discharge computed for %d observations -> swot_discharge.csv", len(df))
|
||||
if len(df):
|
||||
log.info("Discharge range: %.1f - %.1f m3/s (median %.1f).",
|
||||
df["discharge"].min(), df["discharge"].max(),
|
||||
df["discharge"].median())
|
||||
@@ -0,0 +1,81 @@
|
||||
"""Step 7 - Validate SWOT-derived discharge against CWC gauge discharge.
|
||||
|
||||
Each SWOT discharge estimate is paired with the nearest gauge discharge
|
||||
observation within +/- matching.discharge_time_window_hours (default 24 h,
|
||||
since CWC discharge is typically daily/weekly at 08:00 IST).
|
||||
|
||||
Outputs:
|
||||
processed/matched_discharge.csv
|
||||
outputs/discharge_metrics.csv
|
||||
outputs/plots/discharge_<station>.png
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
import pandas as pd
|
||||
from tqdm import tqdm
|
||||
|
||||
from .common import log, nearest_in_time, p, regression_metrics
|
||||
|
||||
|
||||
def run(cfg: dict) -> None:
|
||||
import matplotlib
|
||||
matplotlib.use("Agg")
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
proc = p(cfg, "paths", "processed")
|
||||
out = p(cfg, "paths", "outputs")
|
||||
plots = out / "plots"
|
||||
plots.mkdir(exist_ok=True)
|
||||
|
||||
swot = pd.read_csv(proc / "swot_discharge.csv", parse_dates=["time"])
|
||||
gauge = pd.read_csv(proc / "gauge_discharge.csv", parse_dates=["time"])
|
||||
swot = swot.dropna(subset=["station"])
|
||||
if swot.empty:
|
||||
log.warning("No station-tagged SWOT discharge (scope=all?). "
|
||||
"Validation needs scope=matched estimates.")
|
||||
return
|
||||
|
||||
window = cfg["matching"]["discharge_time_window_hours"]
|
||||
pairs = []
|
||||
for st, g in swot.groupby("station"):
|
||||
gq = gauge[gauge["station"] == st].reset_index(drop=True)
|
||||
if gq.empty:
|
||||
log.warning("%s: no gauge discharge file - skipped.", st)
|
||||
continue
|
||||
for _, obs in tqdm(g.iterrows(), total=len(g), unit="obs", leave=False,
|
||||
desc=f"Matching Q - {st}"):
|
||||
val, t = nearest_in_time(obs["time"], gq["time"], gq["value"], window)
|
||||
if val is not None:
|
||||
pairs.append({"station": st, "time_swot": obs["time"],
|
||||
"time_gauge": t, "swot_q": obs["discharge"],
|
||||
"gauge_q": val, "depth": obs["depth"],
|
||||
"width": obs["width"], "slope": obs["slope"]})
|
||||
matched = pd.DataFrame(pairs)
|
||||
matched.to_csv(proc / "matched_discharge.csv", index=False)
|
||||
log.info("Matched discharge pairs: %d", len(matched))
|
||||
if matched.empty:
|
||||
return
|
||||
|
||||
metrics = []
|
||||
for st, g in matched.groupby("station"):
|
||||
mt = regression_metrics(g["gauge_q"], g["swot_q"])
|
||||
mt["station"] = st
|
||||
metrics.append(mt)
|
||||
|
||||
fig, ax = plt.subplots(1, 2, figsize=(11, 4.2))
|
||||
ax[0].scatter(g["gauge_q"], g["swot_q"], s=20, alpha=0.7, color="tab:green")
|
||||
lims = [0, max(g["gauge_q"].max(), g["swot_q"].max()) * 1.05]
|
||||
ax[0].plot(lims, lims, "k--", lw=1)
|
||||
ax[0].set_xlabel("Gauge discharge (m³/s)")
|
||||
ax[0].set_ylabel("SWOT discharge (m³/s)")
|
||||
ax[0].set_title(f"{st} R²={mt['r2']:.3f} n={mt['n']}")
|
||||
gg = g.sort_values("time_swot")
|
||||
ax[1].plot(gg["time_swot"], gg["gauge_q"], "o-", ms=4, label="Gauge")
|
||||
ax[1].plot(gg["time_swot"], gg["swot_q"], "s--", ms=4, label="SWOT (Manning)")
|
||||
ax[1].set_ylabel("Q (m³/s)"); ax[1].legend(); ax[1].set_title("Time series")
|
||||
fig.autofmt_xdate(); fig.tight_layout()
|
||||
fig.savefig(plots / f"discharge_{st}.png", dpi=150); plt.close(fig)
|
||||
|
||||
mdf = pd.DataFrame(metrics)[["station", "n", "r2", "rmse", "nse", "bias"]]
|
||||
mdf.to_csv(out / "discharge_metrics.csv", index=False)
|
||||
log.info("Discharge validation metrics:\n%s", mdf.to_string(index=False))
|
||||
@@ -0,0 +1,260 @@
|
||||
"""Step 8 - Export all inputs, intermediates and results to one Excel workbook.
|
||||
|
||||
Creates outputs/<project_name>_swot_analysis.xlsx with these sheets:
|
||||
|
||||
README what each sheet contains + the run configuration
|
||||
Stations gauge stations and their matched SWOT node/reach
|
||||
Bed_Elevation DEM bed, low-water WSE (P5) and the effective bed used
|
||||
SWOT_Node_Obs cleaned SWOT observations at the matched nodes (input)
|
||||
Reach_Slopes SWOT reach slope observations + per-reach median
|
||||
Gauge_WSE in-situ water levels for the listed stations (input)
|
||||
Gauge_Discharge in-situ discharge for the listed stations (input)
|
||||
WSE_Validation SWOT vs gauge WSE pairs, residual as a live formula
|
||||
Discharge_Results every discharge estimate with all hydraulic intermediates
|
||||
Discharge_Validation SWOT vs gauge discharge pairs, residual + % error formulas
|
||||
Metrics R^2 / RMSE / NSE / bias per station for WSE and discharge
|
||||
|
||||
All times are UTC. Residual columns are Excel formulas so the workbook stays
|
||||
consistent if values are edited.
|
||||
"""
|
||||
from __future__ import annotations
|
||||
|
||||
from datetime import datetime, timezone
|
||||
|
||||
import pandas as pd
|
||||
from openpyxl.styles import Alignment, Font, PatternFill
|
||||
from openpyxl.utils import get_column_letter
|
||||
|
||||
from .common import log, p
|
||||
|
||||
FONT = "Arial"
|
||||
HDR_FILL = PatternFill("solid", start_color="1F4E79")
|
||||
HDR_FONT = Font(name=FONT, bold=True, color="FFFFFF")
|
||||
|
||||
|
||||
def _read(proc, name, **kw):
|
||||
path = proc / name
|
||||
if not path.exists():
|
||||
log.warning("%s not found - sheet will be skipped.", name)
|
||||
return None
|
||||
return pd.read_csv(path, **kw)
|
||||
|
||||
|
||||
def _strip_tz(df):
|
||||
for c in df.columns:
|
||||
# object-dtype time columns (mixed formats defeat parse_dates)
|
||||
if df[c].dtype == object and "time" in str(c).lower():
|
||||
parsed = pd.to_datetime(df[c], errors="coerce", utc=True,
|
||||
format="mixed")
|
||||
if parsed.notna().mean() > 0.9:
|
||||
df[c] = parsed
|
||||
if pd.api.types.is_datetime64_any_dtype(df[c]):
|
||||
try:
|
||||
df[c] = df[c].dt.tz_convert("UTC").dt.tz_localize(None)
|
||||
except TypeError:
|
||||
pass
|
||||
return df
|
||||
|
||||
|
||||
def _write_sheet(writer, df, name, formulas=None, num_fmt=None):
|
||||
"""Write df; `formulas` = {col_name: template with {r}=excel row}."""
|
||||
df = _strip_tz(df.copy())
|
||||
df.to_excel(writer, sheet_name=name, index=False)
|
||||
ws = writer.sheets[name]
|
||||
cols = list(df.columns)
|
||||
if formulas:
|
||||
for col, template in formulas.items():
|
||||
ci = cols.index(col) + 1
|
||||
for r in range(2, len(df) + 2):
|
||||
ws.cell(row=r, column=ci).value = template.format(r=r)
|
||||
for ci, col in enumerate(cols, 1):
|
||||
cell = ws.cell(row=1, column=ci)
|
||||
cell.font = HDR_FONT
|
||||
cell.fill = HDR_FILL
|
||||
cell.alignment = Alignment(horizontal="center")
|
||||
width = max(len(str(col)) + 2, 12)
|
||||
if "time" in str(col).lower():
|
||||
width = 20
|
||||
ws.column_dimensions[get_column_letter(ci)].width = width
|
||||
if num_fmt and col in num_fmt:
|
||||
for r in range(2, len(df) + 2):
|
||||
ws.cell(row=r, column=ci).number_format = num_fmt[col]
|
||||
for row in ws.iter_rows(min_row=2):
|
||||
for cell in row:
|
||||
if cell.font.name != FONT:
|
||||
cell.font = Font(name=FONT)
|
||||
ws.freeze_panes = "A2"
|
||||
log.info(" sheet %-22s %6d rows", name, len(df))
|
||||
|
||||
|
||||
def _col(df, name):
|
||||
return get_column_letter(list(df.columns).index(name) + 1)
|
||||
|
||||
|
||||
def run(cfg: dict) -> None:
|
||||
proc = p(cfg, "paths", "processed")
|
||||
out = p(cfg, "paths", "outputs")
|
||||
xlsx = out / f"{cfg.get('project_name', 'swot')}_swot_analysis.xlsx"
|
||||
|
||||
stations = pd.read_csv(cfg["_root"] / cfg["gauge"]["stations_csv"], comment="#")
|
||||
st_names = set(stations["station"].astype(str).str.strip())
|
||||
mapping = _read(proc, "station_nodes.csv", dtype={"node_id": str, "reach_id": str})
|
||||
nodes = _read(proc, "nodes_clean.csv", parse_dates=["time"],
|
||||
dtype={"node_id": str, "reach_id": str})
|
||||
reaches = _read(proc, "reaches_clean.csv", parse_dates=["time"],
|
||||
dtype={"reach_id": str})
|
||||
bed = _read(proc, "node_bed_elevation.csv", dtype={"node_id": str})
|
||||
gwse = _read(proc, "gauge_wse.csv", parse_dates=["time"])
|
||||
gq = _read(proc, "gauge_discharge.csv", parse_dates=["time"])
|
||||
mwse = _read(proc, "matched_wse.csv", parse_dates=["time_swot", "time_gauge"],
|
||||
dtype={"node_id": str, "reach_id": str})
|
||||
disq = _read(proc, "swot_discharge.csv", parse_dates=["time"],
|
||||
dtype={"node_id": str, "reach_id": str})
|
||||
mq = _read(proc, "matched_discharge.csv", parse_dates=["time_swot", "time_gauge"])
|
||||
wm = _read(out, "wse_metrics.csv")
|
||||
qm = _read(out, "discharge_metrics.csv")
|
||||
|
||||
log.info("Writing %s ...", xlsx.name)
|
||||
with pd.ExcelWriter(xlsx, engine="openpyxl") as writer:
|
||||
# README ------------------------------------------------------------
|
||||
h = cfg["hydraulics"]
|
||||
m = cfg["matching"]
|
||||
readme = pd.DataFrame({
|
||||
"Item": ["Generated", "Project", "Period", "Manning's n",
|
||||
"Baseline depth clamp (m)", "Slope / width source",
|
||||
"WSE match window (h)", "Discharge match window (h)",
|
||||
"Max gauge-node distance (km)", "Times", "", "Sheets:",
|
||||
"Stations", "Bed_Elevation", "SWOT_Node_Obs",
|
||||
"Reach_Slopes", "Gauge_WSE", "Gauge_Discharge",
|
||||
"WSE_Validation", "Discharge_Results",
|
||||
"Discharge_Validation", "Metrics"],
|
||||
"Value": [datetime.now(timezone.utc).strftime("%Y-%m-%d %H:%M UTC"),
|
||||
cfg.get("project_name", ""),
|
||||
f"{cfg['time']['start']} to {cfg['time']['end']}",
|
||||
h["mannings_n"],
|
||||
f"[{h.get('baseline_depth_min_m', 0.5)}, "
|
||||
f"{h.get('baseline_depth_max_m', 2.0)}]",
|
||||
f"{h.get('slope_source', 'median')} / "
|
||||
f"{h.get('width_source', 'median')}",
|
||||
m["wse_time_window_hours"],
|
||||
m["discharge_time_window_hours"],
|
||||
m["max_distance_km"], "all times UTC", "", "",
|
||||
"gauge stations + matched SWOT node/reach + distance",
|
||||
"DEM bed vs low-water WSE (P5) and effective bed used",
|
||||
"cleaned SWOT node observations at matched nodes (input)",
|
||||
"SWOT reach slope observations (input)",
|
||||
"in-situ water level series (input)",
|
||||
"in-situ discharge series (input)",
|
||||
"SWOT vs gauge WSE pairs with residual",
|
||||
"all discharge estimates with hydraulic intermediates",
|
||||
"SWOT vs gauge discharge pairs with residual and % error",
|
||||
"R2/RMSE/NSE/bias per station"]})
|
||||
_write_sheet(writer, readme, "README")
|
||||
writer.sheets["README"].column_dimensions["A"].width = 32
|
||||
writer.sheets["README"].column_dimensions["B"].width = 60
|
||||
|
||||
# Stations ------------------------------------------------------------
|
||||
if mapping is not None:
|
||||
stx = stations.merge(mapping.drop(columns=["gauge_lat", "gauge_lon"]),
|
||||
on="station", how="left")
|
||||
_write_sheet(writer, stx, "Stations",
|
||||
num_fmt={"distance_km": "0.000"})
|
||||
|
||||
# Bed elevation ---------------------------------------------------
|
||||
if bed is not None and nodes is not None:
|
||||
p5 = (nodes.groupby("node_id")["wse"].quantile(0.05)
|
||||
.rename("wse_p5_low_water"))
|
||||
bx = bed.merge(p5, on="node_id", how="left")
|
||||
if disq is not None and "bed_elev" in disq.columns:
|
||||
eff = (disq.groupby("node_id")["bed_elev"].first()
|
||||
.rename("effective_bed_used"))
|
||||
bx = bx.merge(eff, on="node_id", how="left")
|
||||
bx["dem_bed_clamped"] = ~pd.Series(
|
||||
bx["bed_elev"] == bx["effective_bed_used"]).fillna(False)
|
||||
bx = bx.rename(columns={"bed_elev": "dem_bed_elev"})
|
||||
_write_sheet(writer, bx, "Bed_Elevation",
|
||||
num_fmt={c: "0.00" for c in
|
||||
("dem_bed_elev", "wse_p5_low_water",
|
||||
"effective_bed_used")})
|
||||
|
||||
# SWOT node observations (inputs, matched nodes only) --------------
|
||||
if nodes is not None and mapping is not None:
|
||||
sub = nodes[nodes["node_id"].isin(mapping["node_id"])].merge(
|
||||
mapping[["node_id", "station"]], on="node_id", how="left")
|
||||
sub = sub.sort_values(["station", "time"])
|
||||
first = ["station", "node_id", "reach_id", "time"]
|
||||
sub = sub[first + [c for c in sub.columns if c not in first]]
|
||||
_write_sheet(writer, sub, "SWOT_Node_Obs",
|
||||
num_fmt={"wse": "0.000", "width": "0.0",
|
||||
"lat": "0.000000", "lon": "0.000000"})
|
||||
|
||||
# Reach slopes ------------------------------------------------------
|
||||
if reaches is not None and mapping is not None:
|
||||
rsub = reaches[reaches["reach_id"].isin(mapping["reach_id"])].copy()
|
||||
rsub["reach_median_slope"] = rsub.groupby("reach_id")["slope"] \
|
||||
.transform("median")
|
||||
_write_sheet(writer, rsub.sort_values(["reach_id", "time"]),
|
||||
"Reach_Slopes",
|
||||
num_fmt={"slope": "0.00000000",
|
||||
"reach_median_slope": "0.00000000"})
|
||||
|
||||
# Gauge inputs ------------------------------------------------------
|
||||
for df, name in ((gwse, "Gauge_WSE"), (gq, "Gauge_Discharge")):
|
||||
if df is not None:
|
||||
_write_sheet(writer,
|
||||
df[df["station"].isin(st_names)]
|
||||
.sort_values(["station", "time"]),
|
||||
name, num_fmt={"value": "0.000"})
|
||||
|
||||
# WSE validation ----------------------------------------------------
|
||||
if mwse is not None and len(mwse):
|
||||
mwse = mwse.copy()
|
||||
mwse["residual_m"] = 0.0
|
||||
cs, cg = _col(mwse, "swot_wse"), _col(mwse, "gauge_wse")
|
||||
_write_sheet(writer, mwse.sort_values(["station", "time_swot"]),
|
||||
"WSE_Validation",
|
||||
formulas={"residual_m": f"={cs}{{r}}-{cg}{{r}}"},
|
||||
num_fmt={"swot_wse": "0.000", "gauge_wse": "0.000",
|
||||
"residual_m": "0.000"})
|
||||
|
||||
# Discharge results (all intermediates) -----------------------------
|
||||
if disq is not None:
|
||||
_write_sheet(writer, disq.sort_values(["station", "time"]),
|
||||
"Discharge_Results",
|
||||
num_fmt={"wse": "0.000", "bed_elev": "0.00",
|
||||
"depth": "0.00", "width": "0.0",
|
||||
"area": "0.0", "wetted_perimeter": "0.0",
|
||||
"hydraulic_radius": "0.00",
|
||||
"slope": "0.00000000", "mannings_n": "0.000",
|
||||
"discharge": "0.0"})
|
||||
|
||||
# Discharge validation ----------------------------------------------
|
||||
if mq is not None and len(mq):
|
||||
mq = mq.copy()
|
||||
mq["residual_m3s"] = 0.0
|
||||
mq["pct_error"] = 0.0
|
||||
cs, cg = _col(mq, "swot_q"), _col(mq, "gauge_q")
|
||||
_write_sheet(writer, mq.sort_values(["station", "time_swot"]),
|
||||
"Discharge_Validation",
|
||||
formulas={
|
||||
"residual_m3s": f"={cs}{{r}}-{cg}{{r}}",
|
||||
"pct_error": (f"=IF({cg}{{r}}=0,\"\","
|
||||
f"({cs}{{r}}-{cg}{{r}})/{cg}{{r}})")},
|
||||
num_fmt={"swot_q": "0.0", "gauge_q": "0.0",
|
||||
"residual_m3s": "0.0", "pct_error": "0.0%"})
|
||||
|
||||
# Metrics -------------------------------------------------------
|
||||
frames = []
|
||||
if wm is not None:
|
||||
frames.append(wm.assign(variable="WSE"))
|
||||
if qm is not None:
|
||||
frames.append(qm.assign(variable="Discharge"))
|
||||
if frames:
|
||||
allm = pd.concat(frames, ignore_index=True)
|
||||
allm = allm[["variable"] + [c for c in allm.columns
|
||||
if c != "variable"]]
|
||||
_write_sheet(writer, allm, "Metrics",
|
||||
num_fmt={"r2": "0.000", "rmse": "0.00",
|
||||
"nse": "0.000", "bias": "0.00"})
|
||||
|
||||
log.info("Excel export complete: %s", xlsx)
|
||||
Reference in New Issue
Block a user