fix: allow no-divide flowpaths (will be virtualized in nhf-builds)

src/reference_builds/pipeline/download.py

@@ -37,86 +37,6 @@ def _load_and_concat_parquet(parquet_files: list[Path]) -> gpd.GeoDataFrame:
     return pd.concat(gdfs, ignore_index=True)
 
 
-def _merge_flowpaths_without_catchments(
-    flowpaths: gpd.GeoDataFrame,
-    catchments: gpd.GeoDataFrame,
-    connectivity: pd.DataFrame,
-) -> tuple[gpd.GeoDataFrame, pd.DataFrame]:
-    """Merge flowpaths without catchments to its downstream neighbor.
-
-    FIXME: This currently breaks connectivity info as well as path length calculations among other things. Needs significant rework.
-
-    Parameters
-    ----------
-    flowpaths : gpd.GeoDataFrame
-        The flowpath geodataframe, must contain 'DnHydroSeq' and 'HydroSeq' columns
-        for upstream/downstream connectivity and 'NHDPlusID' to match with its catchment.
-    catchments : gpd.GeoDataFrame
-        The catchment geodataframe, must contain 'NHDPlusID' column
-
-    Returns
-    -------
-    tuple[gpd.GeoDataFrame, pd.DataFrame]
-        The updated flowpath geodataframe and connectivity dataframe after merging flowpaths without catchments
-    """
-    _flowpaths = flowpaths.copy()
-    _rename_mapping = {key: f"{key}_VAA" for key in connectivity.columns if key != "NHDPlusID"}
-    _connectivity = connectivity.copy().rename(columns=_rename_mapping)
-
-    # merge flowpaths with connectivity info
-    n_connect = _flowpaths["NHDPlusID"].isin(_connectivity["NHDPlusID"]).sum()
-    n_flowpaths = len(_flowpaths.index)
-    if n_connect != n_flowpaths:
-        n_missing = n_flowpaths - n_connect
-        logger.warning(
-            f"{n_missing}/{n_flowpaths} flowpaths are missing connectivity info. These flowpaths will be dropped from the reference build."
-        )
-
-    _flowpaths = _flowpaths.merge(_connectivity, on="NHDPlusID", how="inner")
-
-    # identify flowpaths without catchments
-    _has_catchment = _flowpaths["NHDPlusID"].isin(catchments["NHDPlusID"])
-    _flowpaths_without_catchments = _flowpaths[~_has_catchment]
-
-    logger.info(f"Merging {len(_flowpaths_without_catchments)} flowpaths without catchments")
-
-    for _, row in _flowpaths_without_catchments.iterrows():
-        hydroseq = row["HydroSeq_VAA"]
-        dnhydroseq = row["DnHydroSeq_VAA"]
-        downstream_flowpath = _flowpaths[_flowpaths["HydroSeq_VAA"] == dnhydroseq]
-
-        if len(downstream_flowpath) == 0 or row["TerminalFl_VAA"] == 1:
-            continue
-        elif len(downstream_flowpath) > 1:
-            logger.warning(f"Multiple downstream flowpaths found for {row['NHDPlusID']}")
-            continue
-        else:
-            # merge geometry
-            _gdf = gpd.GeoDataFrame([row, downstream_flowpath.iloc[0]], geometry="geometry")
-            merged_geom = _gdf.geometry.union_all()
-            # update geometry of downstream flowpath
-            dn_idx = downstream_flowpath.index[0]
-            _flowpaths.at[dn_idx, "geometry"] = merged_geom
-            # update connectivity
-            upstream_flowpaths = _flowpaths[_flowpaths["DnHydroSeq_VAA"] == hydroseq]
-            if len(upstream_flowpaths) >= 1:
-                _flowpaths.loc[upstream_flowpaths.index, "DnHydroSeq_VAA"] = dnhydroseq
-                _flowpaths.loc[dn_idx, "FromNode_VAA"] = upstream_flowpaths.iloc[0]["ToNode_VAA"]
-            # TODO: check if other attributes need to be updated??
-
-    # drop flowpaths without catchments after merging
-    _flowpaths = _flowpaths[_has_catchment]
-
-    # derive updated connectivity from merged flowpaths
-    connectivity = _flowpaths[list(_rename_mapping.values())].rename(
-        columns={v: k for k, v in _rename_mapping.items()}
-    )
-    connectivity["NHDPlusID"] = _flowpaths["NHDPlusID"]
-    flowpaths = _flowpaths.drop(columns=list(_rename_mapping.values()))
-
-    return flowpaths, connectivity
-
-
 def download_geoglows_data(**context: dict[str, Any]) -> dict[str, pl.DataFrame]:
     """Opens local / downloads for the reference-build process
 
@@ -196,10 +116,7 @@ def download_nhd_data(**context: dict[str, Any]) -> dict[str, pl.DataFrame]:
     _flowpaths = _validate_and_fix_geometries(data["NHDFlowline"], geom_type="flowpaths")
     catchments = _validate_and_fix_geometries(data["NHDPlusCatchment"], geom_type="divides")
 
-    _flowpaths_with_catchments = _flowpaths[_flowpaths["NHDPlusID"].isin(catchments["NHDPlusID"])]
-    flowpaths = _flowpaths_with_catchments[
-        _flowpaths_with_catchments["fcode_description"].isin(cfg.permitted_fcodes)
-    ]
+    flowpaths = _flowpaths[_flowpaths["fcode_description"].isin(cfg.permitted_fcodes)]
 
     return {
         "nhd_flowpaths": pl.from_pandas(flowpaths.to_wkb()),
@@ -253,7 +170,6 @@ def download_usgs_hf_data(**context: dict[str, Any]) -> dict[str, pl.DataFrame]:
 
     # filter/validate layers
     _flowpaths = _validate_and_fix_geometries(flowpaths, geom_type="flowpaths")
-    _flowpaths = _flowpaths[_flowpaths["comid"].isin(catchments["COMID"])]
 
     catchments = _validate_and_fix_geometries(catchments, geom_type="divides")
     catchments = _fix_divide_exclaves(catchments.to_crs("EPSG:3338")).to_crs("EPSG:4326")

src/reference_builds/utils/geometries.py

@@ -107,7 +107,7 @@ def _drop_exclaves(geom: Geometry) -> Geometry:
 
 def _find_exclaves(geom: Geometry) -> pd.Series:
     """Find and exclude non-contiguous parts of MultiPolygons, appending them to a list to be resolved later"""
-    exclaves = []
+    exclaves: list = []
     if geom.geom_type != "MultiPolygon":
         return pd.Series(data={"geometry": geom, "exclaves": exclaves}, index=["geometry", "exclaves"])