Updating workflows for reservoirs and rivers

AGENTS.md

@@ -16,7 +16,9 @@ This document provides codebase knowledge for AI coding agents (Claude, Copilot,
 
 2. **Rivers** (`rivers` branch)
    - SWOT Hydrocron API download (public, no credentials needed)
+   - Multi-satellite data download (SWOT Lake SP, ICESat-2 ATL13, Sentinel-3, Sentinel-6)
    - SWORD v17b database integration for river node/reach matching
+   - Full timeseries extraction, cleaning, filtering, and multi-mission merge aligned to SWORD reaches or nodes
    - Output: per-river Hydrocron timeseries (CSV)
 
 3. **SWOT Raster Tiles** (`swot_raster` branch)

CITATION.cff

@@ -4,22 +4,32 @@ message: >-
   HydroEO's core timeseries pipeline additionally implements a
   methodology described in Schwatke et al. (2015) -- see the
   "references" section below and HydroEO/utils/filters/basic_filters.py
-  -- please cite that publication as well when using the Kalman filter /
-  outlier rejection pipeline it implements.
+  -- please cite that publication as well when using the Kalman filter,
+  SVR, or windowed ADM-based outlier/uncertainty weighting it
+  implements.
 title: "HydroEO"
 version: "0.1.0"
 license: MIT
 repository-code: "https://github.com/DHI/HydroEO"
 authors:
-  - family-names: "Chewning"
-    given-names: "Connor"
-    affiliation: "DHI"
   - family-names: "Kittel"
     given-names: "Cecile"
     affiliation: "DHI"
+  - family-names: "Chewning"
+    given-names: "Connor"
+    affiliation: "DHI"
   - family-names: "Zalite"
     given-names: "Karlis"
     affiliation: "DHI"
+  - family-names: "Monica"
+    given-names: "Coppo Frias"
+    affiliation: "DHI"
+  - family-names: "Sarah"
+    given-names: "Franze"
+    affiliation: "DHI"
+  - family-names: "Paul"
+    given-names: "Senty"
+    affiliation: "DHI"
 references:
   - type: article
     authors:
@@ -41,5 +51,7 @@ references:
     notes: >-
       HydroEO's core timeseries pipeline (see
       HydroEO/utils/filters/basic_filters.py) follows the windowed
-      outlier rejection, multi-mission bias correction, and Kalman
-      filter approach described in this publication.
+      ADM-based outlier/uncertainty weighting, SVR, and Kalman filter
+      approach described in this publication. The multi-mission bias
+      correction and mission-combination logic are this codebase's own
+      contributions, not sourced from this paper.

HydroEO/constants.py

@@ -88,6 +88,13 @@
 
 SUPPORTED_CLEAN_FILTERS = ["elevation", "MAD", "daily_mean", "hampel", "rolling_median"]
 
+# Project.__sat_init and HydroEO.flows._clean_engine both
+# reference these constants
+DEFAULT_PROCESSING_FILTERS = ["elevation", "MAD"]
+DEFAULT_ELEVATION_MIN_M = 0.0
+DEFAULT_ELEVATION_MAX_M = 8000.0
+DEFAULT_MAD_THRESHOLD = 5.0
+
 # ============================================================================
 # Mission Defaults
 # ============================================================================
@@ -96,46 +103,49 @@
     "swot": {
         "download": False,
         "process": False,
-        "pld_match_max_distance_m": 100.0,
+        # Minimum PLD-lake overlap as a percentage (0-100) of the
+        # reservoir's own area. Previously used pld_match_max_distance_m.
+        # NOTE: Note backward compatible.
+        "pld_match_min_overlap_pct": 10.0,
         "exclude_obs_id_values": ["no_data"],
         "hydrocron_fields": SWOT_DEFAULT_HYDROCRON_FIELDS,
         "quality_filters": SWOT_DEFAULT_QUALITY_FILTERS,
-        "processing_filters": ["elevation", "MAD"],
-        "elevation_min_m": 0.0,
-        "elevation_max_m": 8000.0,
-        "mad_threshold": 5.0,
+        "processing_filters": DEFAULT_PROCESSING_FILTERS,
+        "elevation_min_m": DEFAULT_ELEVATION_MIN_M,
+        "elevation_max_m": DEFAULT_ELEVATION_MAX_M,
+        "mad_threshold": DEFAULT_MAD_THRESHOLD,
     },
     "icesat2": {
         "download": False,
         "process": False,
         "atl13_fields": ICESAT2_DEFAULT_FIELDS,
         "atl13": {"pass_invalid": False, "beams": [], "spots": []},
         "track_keys": ICESAT2_SUPPORTED_TRACK_KEYS,
-        "processing_filters": ["elevation", "MAD"],
-        "elevation_min_m": 0.0,
-        "elevation_max_m": 8000.0,
-        "mad_threshold": 5.0,
+        "processing_filters": DEFAULT_PROCESSING_FILTERS,
+        "elevation_min_m": DEFAULT_ELEVATION_MIN_M,
+        "elevation_max_m": DEFAULT_ELEVATION_MAX_M,
+        "mad_threshold": DEFAULT_MAD_THRESHOLD,
     },
     "sentinel3": {
         "download": False,
         "process": False,
         "subset_file_id": "enhanced_measurement.nc",
         "sigma0_max": 1e5,
         "download_threads": 1,
-        "processing_filters": ["elevation", "MAD"],
-        "elevation_min_m": 0.0,
-        "elevation_max_m": 8000.0,
-        "mad_threshold": 5.0,
+        "processing_filters": DEFAULT_PROCESSING_FILTERS,
+        "elevation_min_m": DEFAULT_ELEVATION_MIN_M,
+        "elevation_max_m": DEFAULT_ELEVATION_MAX_M,
+        "mad_threshold": DEFAULT_MAD_THRESHOLD,
     },
     "sentinel6": {
         "download": False,
         "process": False,
         "subset_file_id": "enhanced_measurement.nc",
         "sigma0_max": 1e5,
         "download_threads": 1,
-        "processing_filters": ["elevation", "MAD"],
-        "elevation_min_m": 0.0,
-        "elevation_max_m": 8000.0,
-        "mad_threshold": 5.0,
+        "processing_filters": DEFAULT_PROCESSING_FILTERS,
+        "elevation_min_m": DEFAULT_ELEVATION_MIN_M,
+        "elevation_max_m": DEFAULT_ELEVATION_MAX_M,
+        "mad_threshold": DEFAULT_MAD_THRESHOLD,
     },
-}
+}

HydroEO/downloaders/hydroweb.py

@@ -30,8 +30,12 @@
 
 
 def download_PLD(
-    download_dir: str, bounds: list, raw_pld_path: str = None, keep_raw: bool = True
-):
+    download_dir: str,
+    bounds: list,
+    raw_pld_path: str = None,
+    keep_raw: bool = True,
+    continent_codes: list = None,
+    ):
     """Download and subset SWOT Prior Lake Database.
 
     Parameters
@@ -44,7 +48,13 @@ def download_PLD(
         Path to existing PLD zip file or extracted folder. If provided, skips download.
     keep_raw : bool, default True
         If False, delete raw zip and temp extraction folder after subset is created.
-    """
+    continent_codes : list[str], optional
+        Continent-tile suffixes (e.g. ["AS", "AU"]) identifying which
+        full-schema PLD files (containing res_id, not just lake_id +
+        geometry) to use for backfilling res_id only. 
+        Inspect your downloaded PLD folder to see which ones are actually present. 
+        If omitted, all full-schema tiles found are used.
+     """
     # create download directory if needed
     general.ifnotmakedirs(download_dir)
 
@@ -59,10 +69,12 @@ def download_PLD(
     # Track whether the zip was downloaded by HydroEO (vs user-provided)
     hydroweb_managed_zip = True
 
-    # Check if extracted files already exist
-    extracted_files_exist = os.path.isdir(extracted_dir) and any(
-        f.endswith(".sqlite") for f in os.listdir(extracted_dir)
-    )
+    def _dir_has_pld_files(d):
+        return any(f.endswith((".sqlite", ".gpkg")) for f in os.listdir(d))
+
+     # Check if extracted files already exist
+    extracted_files_exist = os.path.isdir(extracted_dir) and _dir_has_pld_files(
+        extracted_dir)
 
     # Handle user-provided raw_pld_path
     if raw_pld_path is not None and os.path.exists(raw_pld_path):
@@ -74,24 +86,23 @@ def download_PLD(
         elif os.path.isdir(raw_pld_path):
             # User provided a directory
             logger.info("Using provided PLD directory: %s", raw_pld_path)
-            # Check if it contains .sqlite files directly
-            if any(f.endswith(".sqlite") for f in os.listdir(raw_pld_path)):
+            # Check if it contains PLD files (.sqlite or .gpkg) directly
+            if _dir_has_pld_files(raw_pld_path):
                 extracted_dir = raw_pld_path
                 extracted_files_exist = True
             else:
                 # Check for nested SWOT_PRIOR_LAKE_DATABASE structure
                 nested_path = os.path.join(
                     raw_pld_path, "SWOT_PRIOR_LAKE_DATABASE", "SWOT_PRIOR_LAKE_DATABASE"
                 )
-                if os.path.isdir(nested_path) and any(
-                    f.endswith(".sqlite") for f in os.listdir(nested_path)
-                ):
+                if os.path.isdir(nested_path) and _dir_has_pld_files(nested_path):
                     extracted_dir = nested_path
                     extracted_files_exist = True
                     unzipped_dir = raw_pld_path  # track parent for deletion logic
                 else:
                     logger.warning(
-                        "Provided directory does not contain .sqlite files: %s",
+                        "Provided directory does not contain .sqlite or "
+                        ".gpkg PLD files: %s",
                         raw_pld_path,
                     )
                     return
@@ -149,28 +160,130 @@ def download_PLD(
     # Get list of downloaded files
     downloaded_files = os.listdir(extracted_dir)
 
-    # now clean up and merge lake datafiles
+    def _read_pld_file(filepath, is_sqlite):
+        if is_sqlite:
+            # SWOT PLD .sqlite tiles store lake_id as the feature ID
+            g = gpd.read_file(filepath, layer="lake", fid_as_index=True)
+            g.columns = [c.lower() for c in g.columns]
+            g.index.name = "lake_id"
+            g = g.reset_index()
+        else:
+            # .gpkg PLD files already carry lake_id as a normal column.
+            g = gpd.read_file(filepath)
+            g.columns = [c.lower() for c in g.columns]
+        return g
+
     logger.info("Merging products and removing temporary files")
-    gdf_list = list()
-
-    for file in downloaded_files:
-        if file.endswith(".sqlite"):
-            filepath = os.path.join(extracted_dir, file)
-            logger.info("found %s", filepath)
-
-            # Read layer directly as GeoDataFrame and normalize column names
-            gdf = gpd.read_file(filepath, layer="lake", fid_as_index=True)
-            gdf.columns = [c.lower() for c in gdf.columns]
-            gdf.index.name = "lake_id"
-            gdf = gdf.reset_index()[["lake_id", "res_id", "geometry"]]
-            gdf = gdf[["lake_id", "res_id", "geometry"]]
-
-            # Filter to bounds and append
-            gdf = gdf.loc[gdf.within(shapely.Polygon.from_bounds(*bounds))]
-            gdf_list.append(gdf)
-
-    # once we have processed all files, concatenate them
-    gdf = pd.concat(gdf_list).reset_index(drop=True)
+    light_files = [
+        f
+        for f in downloaded_files
+        if "_light" in f.lower() and f.lower().endswith((".gpkg", ".sqlite"))
+    ]
+    all_tile_files = [
+        f
+        for f in downloaded_files
+        if f not in light_files and (f.endswith(".sqlite") or f.endswith(".gpkg"))
+    ]
+    bbox_poly = shapely.Polygon.from_bounds(*bounds)
+
+    if light_files:
+        if len(light_files) > 1:
+            logger.warning(
+                "Multiple '_light' PLD files found (%s) -- using all of "
+                "them, but normally there should be exactly one.",
+                ", ".join(light_files),
+            )
+        coverage_list = [
+            _read_pld_file(os.path.join(extracted_dir, f), f.lower().endswith(".sqlite"))
+            for f in light_files
+        ]
+        coverage_gdf = pd.concat(coverage_list, ignore_index=True)
+        if "res_id" not in coverage_gdf.columns:
+            coverage_gdf["res_id"] = None
+        coverage_gdf = coverage_gdf[["lake_id", "res_id", "geometry"]]
+
+        # Find lake using intersection with PLD
+        coverage_gdf = coverage_gdf.loc[
+            coverage_gdf.intersects(bbox_poly)
+        ].reset_index(drop=True)
+
+        backfill_tile_files = all_tile_files
+        if continent_codes:
+            backfill_tile_files = [
+                f
+                for f in all_tile_files
+                if any(code.lower() in f.lower() for code in continent_codes)
+            ]
+            unmatched_codes = [
+                code
+                for code in continent_codes
+                if not any(code.lower() in f.lower() for f in all_tile_files)
+            ]
+            if unmatched_codes:
+                logger.warning(
+                    "hydroweb.continent_codes %s requested for PLD res_id "
+                    "backfill, but no matching tile file was found among "
+                    "the downloaded PLD data (%s). Re-download the PLD "
+                    "if you expect a tile for these regions.",
+                    unmatched_codes,
+                    ", ".join(downloaded_files),
+                )
+
+        if backfill_tile_files and len(coverage_gdf) > 0:
+            res_id_map = {}
+            for file in backfill_tile_files:
+                filepath = os.path.join(extracted_dir, file)
+                logger.info("found %s (for res_id backfill)", filepath)
+                tile_gdf = _read_pld_file(filepath, file.endswith(".sqlite"))
+                if "res_id" in tile_gdf.columns:
+                    res_id_map.update(dict(zip(tile_gdf["lake_id"], tile_gdf["res_id"])))
+            if res_id_map:
+                coverage_gdf["res_id"] = coverage_gdf["lake_id"].map(
+                    res_id_map
+                ).combine_first(coverage_gdf["res_id"])
+
+        still_missing = int(coverage_gdf["res_id"].isna().sum())
+        if still_missing > 0:
+            logger.warning(
+                "%d of %d PLD lake(s) within this project's bounds could "
+                "not have res_id backfilled from the available "
+                "tile file(s) (%s). Set hydroweb.continent_codes in "
+                "your config to the correct region code(s) "
+                "and/or re-download the PLD to include the missing tile for backfilling res_id.",
+                still_missing,
+                len(coverage_gdf),
+                ", ".join(all_tile_files) if all_tile_files else "none present",
+            )
+    else:
+        logger.warning(
+            "No global '_light' PLD file found among the downloaded PLD "
+            "data (%s) -- falling back to per-continent/region "
+            "tile files. If this PLD download "
+            "Check missing_res_id and redownload PLD if data is unexpectedly missing. " \
+            "If you expect a global '_light' file, re-download the PLD to include it.",
+            ", ".join(downloaded_files),
+        )
+        coverage_list = [
+            _read_pld_file(os.path.join(extracted_dir, f), f.endswith(".sqlite"))
+            for f in all_tile_files
+        ]
+        if coverage_list:
+            coverage_gdf = pd.concat(coverage_list, ignore_index=True)
+        else:
+            coverage_gdf = gpd.GeoDataFrame(
+                {"lake_id": [], "res_id": [], "geometry": []}, crs="EPSG:4326"
+            )
+        if "res_id" not in coverage_gdf.columns:
+            coverage_gdf["res_id"] = None
+        coverage_gdf = coverage_gdf[["lake_id", "res_id", "geometry"]]
+        coverage_gdf = coverage_gdf.loc[
+            coverage_gdf.intersects(bbox_poly)
+        ].reset_index(drop=True)
+
+    gdf = coverage_gdf
+
+    gdf["res_id"] = gdf["res_id"].astype("float64")
+
 
     # save the concatenated dataframe as GPKG
     export_path = os.path.join(download_dir, "PLD_subset.gpkg")