Skip to content

Updating workflows for reservoirs and rivers - #54

Merged
kzalite merged 21 commits into
DHI:mainfrom
KittelC:main
Jul 17, 2026
Merged

Updating workflows for reservoirs and rivers#54
kzalite merged 21 commits into
DHI:mainfrom
KittelC:main

Conversation

@KittelC

@KittelC KittelC commented Jul 13, 2026

Copy link
Copy Markdown
Contributor

1. Kalman filter correctness + performance

Files: basic_filters.py, timeseries.py

  • Closed-form scalar update in _update(): for the n=1 case (
    only case used in this codebase), replaced the general matrix
    pinv-based Kalman update with a direct precision-weighted-average
    formula. Verified mathematically identical via 2000-trial fuzz test.
  • Epoch grouping fix: kalman() was grouping by exact time stamp, not calendar day, causing visible
    within-day convergence artifacts on dense passes. Fixed to group by
    date.floor("D"). Confirmed on real data: one case went from std
    0.432m (wrong, from hundreds of sequential single-point updates) to a
    single clean value matching the raw data's actual range.

2. svr_radial_max_iter / svr_linear_max_iter separation

Files: basic_filters.py, timeseries.py

  • svr_linear capped to 5000 iterations for speed.
  • Capping svr_radial's (RBF kernel, high C) max_iter was tried once
    and found actively unsafe: output quality was non-monotonic in
    max_iter (5000→10/5760 kept, 50000→827/5760, only unbounded→correct
    2154/5760). Reverted to -1 (unbounded) as the correct default.
  • Review focus: confirm svr_radial_max_iter and
    svr_linear_max_iter are genuinely separate everywhere, consistently
    named.

3. Extraction skip-if-exists + Kalman grouping performance

Files: basic_filters.py, flows.py, project.py

  • kalman() performance improved by grouping once
    upfront into a dict.
  • Dominant time cost: _extract_reservoirs_timeseries
    re-extracting every raw file on every run regardless of whether output
    already existed. Added overwrite_extraction (default False) to
    skip already-extracted targets, implemented separately for ICESat-2,
    Sentinel-3/6, and SWOT.
  • Bug caught in the process: overwrite_extraction was read via
    getattr(..., False) in flows.py but nothing ever set it from
    config — silently always False. Fixed in project.py for both
    reservoirs and rivers.

4. Bias correction tuning + config wiring

Files: flows.py, project.py

  • Implemented bias correction between missions using platform and relative orbit as
    stable key for Sentinel-3 - Sentinel-6 still uses "pass"
    as orbit_key
    , flagged as unverified whether it has the same
    instability or its own relative_orbit-equivalent.
  • bias_time_bin/bias_min_overlap widened from 10D/3 to 20D/1
    after confirming the stricter values were silently dropping entire
    missions for reservoirs with sparse revisit patterns.
  • merging_options config layering bug fixed: the original lookup was
    all-or-nothing (getattr(...) or DEFAULTS), so setting any override
    would silently discard every other default.
  • mission_options whitelist extended for sigma0_min, source,
    latency, short_name as each was needed.

5. Exclusion / run_config system

Files: flows.py, project.py

New feature: per-target ability to interactively review and exclude
specific passes/platforms/dates from a merge, with the result persisted
in a way that's simultaneously a log, a config, and something hand-editable.

  • {output}/{id}/run_config.yaml: one file per target, holding
    exclusions (a list of platform/orbit/date rules) and
    merging_option_overrides (per-target parameter overrides). Written
    and read by the same functions regardless of whether the change came
    from a notebook call or a hand edit — no separate "notebook state" to
    keep in sync with "the config".
  • Functions (in flows.py, exposed via Project methods so they're
    reachable without importing flows directly):
    • list_target_observations(prj, target_type, id) — summarizes
      existing observations by (platform, orbit), for reviewing what
      could be excluded.
    • exclude_from_target(prj, target_type, id, platform=, orbit=, date=, reason=)
      — adds an exclusion rule at any granularity/combination.
    • list_exclusions / remove_exclusion (by index or by matching
      criteria — e.g. platform="S3B", orbit=1517 — added after initial
      review found index-only removal asymmetric with how exclusions are
      added).
    • set_merging_option(prj, target_type, id, **kwargs) — per-target
      parameter overrides, same persistence.
    • Project._infer_target_type() — auto-resolves reservoirs vs. rivers
      when a project only configures one of them; requires an explicit
      target_type argument only when both are configured.
  • Applied in _merge_timeseries: exclusions filter the concatenated
    dataframe before any bias_correct/Kalman/svr_radial processing.
    merging_option_overrides is the highest-priority layer, above
    project-wide config, above hardcoded defaults.
  • Cache invalidation: adding/removing an exclusion, or changing a
    spatial-correction-relevant option via set_merging_option, deletes
    the cached spatial_correction_model.json (and, after section 12 was
    added, reach_slope_correction_model.json too) — forcing a fresh fit,
    since the model may have been built from data that's no longer
    included.

6. Width-based river buffer + multipolygon handling

Files: flows.py, project.py

  • _river_target_corridor now defaults to per-target buffer
    distances derived from SWORD's own width column
    (buffer = width/2 × width_buffer_factor, default factor 1.05)
    instead of one flat constant for every target.
    An explicit extraction_buffer_meters still overrides this entirely.
  • Multipolygon handling, opposite treatment for the two target types:
    • Rivers: _iter_geometry_pieces loops over each disconnected
      corridor piece and issues a separate download/query for each — both
      for ICESat-2 (exact polygon) and Sentinel-3/6 (per-piece envelope).
    • Reservoirs: _simplify_to_one_polygon collapses a multipolygon
      reservoir into one convex-hull shape, so it's still queried as a
      single target — fixed a genuine pre-existing gap in
      _download_reservoirs_icesat2, which was silently dropping every
      part of a multipolygon but the first.
    • Both helpers verified with synthetic geometries; the width-based
      math checked against real buffer-area calculations, not just
      eyeballed.

7. Plotting fixes + river corridor shading

Files: plotting.py, flows.py

  • plot_river_crossings: fixed hardcoded zoom=15zoom="auto"
    (a river AOI can span far more area than a reservoir polygon). Added
    per-target color-coding and reach midpoint markers.
  • plot_river_data: was reading directly from the raw, unfiltered
    Hydrocron CSV, disconnected from the actual clean/merge pipeline.
    Fixed to read the real merged output via a get_merged_fn callback.
  • plot_merging: References updated processing steps. Added the two newer
    optional steps, fixed subplot-count/single-file edge cases, split
    PLATFORM_COLORS (mission-level) from SATELLITE_COLORS
    (instance-level, since the real platform column values are e.g.
    "S3A" not "sentinel3").
  • Plotting now skips undersampled river targets before any of the three
    plot calls (_has_enough_observations_to_plot).
  • New: corridor_geometry parameterplot_river_crossings can
    now render the actual buffered extraction corridor (see section 6).
    generate_rivers_summaries
    computes this via the same _river_target_corridor call (same
    resolved buffer parameters) actually used for extraction.

8. Misc bug fixes

Files: project.py, rivers.yaml

  • Stale "satellites configured but have no effect" warning predated
    river support — fixed to only fire when neither reservoirs nor
    rivers is configured.
  • rivers.yaml template had SWOT's required config block commented out
    under a misleading "advanced options" header, causing a real, silent
    "SWORD downloads but nothing else does" failure mode. Updated the
    template and added an explicit warning in download_rivers() for
    this exact case.

9. Rivers: SWORD-based extraction pipeline (core)

Files: flows.py, project.py

The largest single chunk of new work — extends the
reservoir-style clean/merge pipeline to river targets (SWORD nodes/reaches).

  • _river_target_corridor (base version — see section 6 for the later
    width-based refinement), _assign_points_to_river_targets (via
    gpd.sjoin_nearest, handling point-to-line matching natively).
  • _download_rivers_icesat2/_download_rivers_sentinel: query over the
    corridor rather than per-target.
  • All three _extract_rivers_*_observations functions: split a single
    waterbody-wide download into the same per-target file structure
    reservoirs already use, so everything downstream is shared rather
    than duplicated. SWOT's version splits Hydrocron's per-waterbody CSV
    the same way, attaching a placeholder geometry purely for file-format
    consistency (nothing downstream reads it).
  • _clean_timeseries/_merge_timeseries generalized (not duplicated)
    via a target_type parameter; _target_centroid/_get_target_ids
    generalize the old reservoir-only helpers the same way.
  • DEFAULT_RIVER_MERGING_OPTIONS, kept as a separate dict from
    reservoirs' and explicitly flagged as not independently validated the
    way the reservoir defaults were (this document
    and PR_RUNBOOK.md both still describe this as unresolved).
  • Wired into project.py's create_timeseries().
  • Confirmed via direct code inspection (not assumed) that this
    architecture naturally supports multi-satellite river virtual
    stations
    _merge_timeseries loops over every configured mission
    for a target with zero target_type-specific restriction, exactly the
    same code path as reservoirs.
  • Confirmed switching feature_type between 'nodes'/'reaches' is
    a pure config change
    , not a code change — every conditional
    (SWORD file selection, target_id_col, downstream plotting/extraction)
    branches on the config value rather than hardcoding either.

10. Sentinel-6 EarthData + credential check

Files: sentinel/download.py, sentinel/__init__.py,
sentinel/preprocess.py, flows.py, project.py

  • CREODIAS only distributes Sentinel-6's Low Rate product; High Rate
    (20Hz Ku-band) is available on PO.DAAC/EarthData.
  • query_earthdata/download_earthdata mirror SWOT's existing
    earthaccess pattern rather than porting a separate client.
  • subset() extended to accept flat EarthData .nc files alongside
    CREODIAS's zipped .SEN6 directories — the actual cropping logic is
    completely unchanged and shared. Verified end-to-end against a real
    uploaded HR sample file
    (2614 rows, sane values) — the one piece of
    this feature tested against real data rather than just reviewed.
  • _sentinel6_use_earthdata()/shared _download_sentinel_for_target()
    helper, branching on mission_options["sentinel6"]["source"].
  • New: Project._require_earthdata_credentials(). Confirmed via
    earthaccess's own documentation that earthaccess.login() (called
    with no explicit strategy, same as the existing SWOT code) falls
    through environment variables → .netrcinteractive prompting
    if neither is present — hanging indefinitely in any non-interactive
    run (a scheduled job, a CLI invocation) rather than failing clearly.
    Unlike CREODIAS, which has an equivalent upfront check
    (_require_creodias_credentials), the EarthData path had none at all.
    Fixed with a matching upfront check, wired into both
    _download_reservoirs_sentinel and _download_rivers_sentinel.
    Verified the partial-credentials edge case (username set without
    password) is correctly treated as missing, not silently accepted.

11. Spatial correction tools + stale progress file fix

Files: basic_filters.py, timeseries.py, flows.py

  • apply_distance_penalty: inflates an observation's error based on
    distance from a reference point, addressing well-sampled-but-
    locationally-biased crossings that local ADM can't detect on its own.
    Off by default.
  • fit_spatial_correction_model/apply_spatial_correction: a
    genuine height correction fit from a dense source (e.g. ICESat-2),
    auto-detecting the most stable spatial axis via per-day slope
    consistency across independent day-fits (not a single pooled
    regression — confirmed empirically this matters). Off by default.
    Persisted per-target, refit only on explicit recalibrate=True.
  • Both demonstrated graphically against real numbers from an actual
    reservoir's data during review: the distance penalty recovers a
    biased Kalman estimate from 299.94m to 299.98m against a true 300.0m;
    the spatial correction's underlying slope was shown to be a real,
    reproducible pattern (5 of 6 independent day-fits agree in sign, the
    two best-supported — highest point count, highest R² — agree closely
    with each other).
  • New: stale conditional progress file fix. distance_penalty.csv/
    spatial_correction.csv are only written if that feature is active
    for the current run, but the merged_progress directory was never
    cleared beforehand (general.ifnotmakedirs only creates it if
    missing). So a run with the feature toggled off would leave a stale
    file from whenever it was last on, looking exactly like current
    output. Fixed by explicitly removing both conditional files at the
    start of merge(), before any of the current run's steps write
    anything. Verified the fix removes only these two files, leaving the
    unconditionally-rewritten step files untouched.

12. Reach slope correction (new, opt-in)

Files: flows.py

New feature addressing a real gap found during review: SWOT's own
slope/slope_u fields (from the RiverSP reach product,
hydrocron_fields.reaches) were being requested and extracted, but
silently dropped before ever reaching the merge pipeline —
PRODUCT_TIMESERIES_KEYS["swot"] never referenced them, and
Timeseries.concat() only keeps key-mapped columns.

  • _fit_reach_slope_correction: reads the raw per-target
    swot.gpkg directly (bypassing the generic Timeseries/concat()
    machinery, which is exactly where slope was being lost), takes the
    median of all available slope values for that reach (verified
    robust to an injected outlier in testing), persists it to
    reach_slope_correction_model.json following the same fit-once/
    cache/invalidate-on-exclusion-change pattern as the spatial correction
    model.
  • _apply_reach_slope_correction: for every non-SWOT row, projects
    its position onto the reach's centerline (in a local metric CRS) to
    get its along-reach distance from the geometric midpoint, then
    adjusts height by slope × distance. SWOT's own rows are left
    untouched. Verified with synthetic data: a crossing near the reach's
    start (far from the midpoint) receives a real, sensible correction
    matching the expected slope×distance magnitude; a crossing already
    near the midpoint receives essentially none; SWOT rows are completely
    unaffected either way.
  • Wired into _merge_timeseries as use_reach_slope_correction (off by
    default), gated to only ever activate when
    prj.rivers.target_id_col == "reach_id" — a no-op for node-mode
    projects, since a ~200m node doesn't have the same along-target slope
    concern a ~10km reach does.
  • Two assumptions flagged as unvalidated against real data, and
    explicitly called out in both this document and the PR description
    recommended in PR_RUNBOOK.md
    :
    1. That SWOT's own reach-level WSE is approximately
      midpoint-referenced — an evidence-based inference from the
      RiverSP processing chain (reach WSE is an aggregate of ~50 roughly
      evenly-spaced nodes, not a value evaluated at one specific point),
      not a fact directly confirmed in SWOT's product documentation.
    2. The correction's sign convention (height - slope × distance) has
      not been empirically checked against real data — confirm it
      reduces cross-mission scatter for a real reach, not increases it,
      before relying on this in production.
  • A related, corrected earlier assumption: this session's original
    reasoning for DEFAULT_RIVER_MERGING_OPTIONS's distance-penalty/
    spatial-correction defaults claimed a river target's footprint is
    "much smaller than a reservoir" — confirmed wrong for reaches
    specifically (~10km typical length), so those tools may matter just
    as much for reaches as for reservoirs. Comment corrected in section 6.

Known open items:

  • River merging_options defaults (section 9) still not
    independently validated against real river data.
  • Sentinel-6's orbit_key/stability (section 4) never re-verified
    against real S6 data the way Sentinel-3's was.
  • Reach slope correction's two core assumptions (section 12)
    explicitly unvalidated — flagged again here on purpose, since it's
    the newest and least-tested addition.

- Changed order of renaming in PLD in case users specify "lake_id" as their id as well
- Open config with utf-8 encoding to avoid bugs when editing .yaml files with Windows
- Path normalization (function in "general" + fixes to project.py) to avoid path issues when opening nc files with hdf5/netcdf4 on Windows (breaks subsetting of sentinel files)
- fixed bug in geometry with coords.pop function set as index
@KittelC KittelC closed this Jul 13, 2026
- Kalman epoch grouping + closed-form update
- svr_radial_max_iter / svr_linear_max_iter separation: linear can be capped to 5000, radial needs longer run
- Extraction skip-if-exists + Kalman grouping performance: allow overwrite = False to avoid re-reading all individual files, Kalman speed up
- Bias correction tuning (relative_orbit, bias_group_by="platform_orbit") = Group by platform and by relative orbit
Spatial correction tools (distance_penalty, spatial_correction_model) = Bias correction tools for tracks far a part - off by default
Config wiring (merging_options layering, mission_options whitelist)
Rivers:
- SWORD extraction pipeline (core)
- Sentinel-6 EarthData download (for HR data)
- Plotting fixes and Misc bug fixes (stale warning, rivers.yaml) for river processing

- Possibility to exclude tracks
- Use SWORD river width to buffer along river line when extracting other missions
- Handle multipolygon extraction
- Clean up of processing files when re-running
- Reach slope correction for river VS.
KittelC added 7 commits July 14, 2026 10:46
@KittelC

KittelC commented Jul 15, 2026

Copy link
Copy Markdown
Contributor Author

Updated to pass tests and with adequate documentation.

@KittelC
KittelC marked this pull request as ready for review July 15, 2026 14:07
@KittelC
KittelC marked this pull request as draft July 15, 2026 14:08
Updated documentation, review of flow.py decomposition

Fixed unary_union depreciation warning

Replace unary_union with union_all(), fixed depreciation warning.
@KittelC
KittelC marked this pull request as ready for review July 16, 2026 12:18
KittelC added 2 commits July 16, 2026 16:09
- Find lake using "intersects" instead of "within" when comparing to PLD (more robust if using buffered / imprecise polygon)
- Fixed potential error if not initializing first and if reservoir in CRS (reprojection order)
- Added hydroweb test
- Stop download if no lake in PLD.
USe more robust Light sqlight SWORD file to find lake_id and backfill with regional file.
Was using random sqlite from unzipped folder, which could lead to lake silently not being found because looking at wrong file.

@kzalite kzalite left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Very nice refactor and update!
I left some small comments, nothing critical.

Comment thread HydroEO/flows/_clean_engine.py Outdated
product,
{
"processing_filters": ["elevation", "MAD"],
"elevation_min_m": 0.0,

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Since these default values are repeated (min_m, max_m, mad_threshold), maybe it's worth setting them as CONSTANTS in the _constants file or at the top of the file?

prj.mission_options.get("swot", {})
.get("quality_filters", {})
.get("nodes", {})
.get("max_q", 2)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same: maybe extract as constant?

Comment thread HydroEO/flows/_river_init.py Outdated
logger.info("Kept raw SWORD zip file at %s (keep_raw_sword=True)", zip_path)


# ============================================================================

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove these commented lines

Comment thread HydroEO/flows/_river_pipeline.py Outdated
gdf.to_file(os.path.join(dst_dir, "swot.gpkg"), driver="GPKG")


# ============================================================================

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Remove these as well

Comment thread HydroEO/project.py Outdated
self.reservoirs.export_to_dfs0 = self.config["reservoirs"].get(
"export_to_dfs0", False
)
# NOTE: this was previously read via getattr(prj.reservoirs,

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This comment can be removed

KittelC added 6 commits July 17, 2026 09:59
Comparison of lake / res id failed due to str / int mismatch. Converted to int in _reser_voir_init.
SWOT fixes and delete unnecessary comments

Comparison of lake / res id failed due to str / int mismatch. Converted to int in _reser_voir_init.
CLeaned up where and how default processing parameters are defined (elev_min and max, and max_threshold) - previously defined 3 places and hardcoded in function, now called as constant and defined only in 1 place unless provided by users.
@KittelC

KittelC commented Jul 17, 2026

Copy link
Copy Markdown
Contributor Author

Reviewed comments above :

  • Removed unncessary comments, cleaned definition of constants across multiple places in the code (min elev, max elev and mad threshold defined in 3 places, now defined 1 place unless defined by user (in HydroEO/constants.py) which flows and project now import instead of redefining.
  • Fixed errors in SWOT download when reusing an already existing PLD - was reading first sqlite in folder which was sometimes regional - now reading global, light PLD and refilling reservoir ID in PLD with regional heavier file if available (prints warning message if not, telling user to redownload PLD for AOI). Fixed type error from PLD columns when comparing to user AOI (Str != int)

KittelC added 2 commits July 17, 2026 13:08
- Check overlap in area
- proximity allowance now an area overlap (in % not meters).
- currently backward compatibility, should be added.
@KittelC

KittelC commented Jul 17, 2026

Copy link
Copy Markdown
Contributor Author

Found bug in SWORD / AOI comparison:

  • now using overlap in area
  • minimum overlapping area (in %) can be set by users
  • documentation for reservoirs also updated.

@kzalite
kzalite merged commit fdbf896 into DHI:main Jul 17, 2026
5 checks passed
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants