From 1c26622a30e6380f3d8401bf38b78b2789ef9737 Mon Sep 17 00:00:00 2001 From: Quercus Hamlin Date: Wed, 15 Jul 2026 15:29:28 -0400 Subject: [PATCH] fix: allow no-divide flowpaths (will be virtualized in nhf-builds) --- src/reference_builds/pipeline/download.py | 86 +---------------------- src/reference_builds/utils/geometries.py | 2 +- 2 files changed, 2 insertions(+), 86 deletions(-) diff --git a/src/reference_builds/pipeline/download.py b/src/reference_builds/pipeline/download.py index db95018..3869251 100644 --- a/src/reference_builds/pipeline/download.py +++ b/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") diff --git a/src/reference_builds/utils/geometries.py b/src/reference_builds/utils/geometries.py index 58d9ce9..eb0dfb3 100644 --- a/src/reference_builds/utils/geometries.py +++ b/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"])