diff --git a/.gitignore b/.gitignore index dff1f35..904cb08 100644 --- a/.gitignore +++ b/.gitignore @@ -27,4 +27,4 @@ src/snkit/_version.py *.dbf *.shx *.cpg -__init__.py +tmp.* diff --git a/src/snkit/network.py b/src/snkit/network.py index 9f27e46..b2fbbf9 100644 --- a/src/snkit/network.py +++ b/src/snkit/network.py @@ -6,7 +6,7 @@ import logging import multiprocessing import os -from typing import Any, Callable, List, Optional, Set, Tuple, Union +from typing import Any, Callable, Dict, List, Optional, Set, Tuple, Union import warnings import geopandas @@ -103,6 +103,10 @@ def __init__( edges = GeoDataFrame(geometry=[]) self.edges = edges + def copy(self) -> "Network": + """Return a Network with copies of the nodes and edges GeoDataFrames""" + return Network(self.nodes.copy(), self.edges.copy()) + def set_crs( self, crs: Optional[pyproj.CRS] = None, @@ -671,13 +675,17 @@ def geoms_to_array( return geom_arr -def concat_dedup(dfs: List[pandas.DataFrame]) -> GeoDataFrame: +def concat_dedup( + dfs: List[pandas.DataFrame], + keep: Optional[str] = "first", + subset: Optional[list[str]] = None, +) -> GeoDataFrame: """Concatenate a list of GeoDataFrames, dropping duplicate geometries - note: repeatedly drops indexes for deduplication to work """ cat = pandas.concat(dfs, axis=0, sort=False) cat.reset_index(drop=True, inplace=True) - cat_dedup = drop_duplicate_geometries(cat) + cat_dedup = drop_duplicate_geometries(cat, keep=keep, subset=subset) cat_dedup.reset_index(drop=True, inplace=True) return cat_dedup @@ -689,12 +697,14 @@ def node_connectivity_degree(node: str, network: Network) -> int: def drop_duplicate_geometries( - gdf: GeoDataFrame, keep: Optional[str] = "first" + gdf: GeoDataFrame, keep: Optional[str] = "first", subset: Optional[list[str]] = None ) -> GeoDataFrame: """Drop duplicate geometries from a dataframe""" - # as of geopandas ~0.6 this should work without explicit conversion to wkb - # discussed in https://github.com/geopandas/geopandas/issues/521 - return gdf.drop_duplicates([gdf.geometry.name]) + if subset is None: + subset = [gdf.geometry.name] + else: + subset = list(set(subset) + set([gdf.geometry.name])) + return gdf.drop_duplicates(subset=subset, keep=keep) def nearest_point_on_edges(point: Point, edges: GeoDataFrame) -> Point: @@ -1120,3 +1130,163 @@ def add_component_ids(network: Network, id_col: str = "component_id") -> Network network.nodes.loc[node_mask, id_col] = count + 1 return network + + +def merge_networks(networks: List[Network], id_col: str = "id") -> Network: + """Merge multiple networks, identifying duplicate nodes at shared locations. + + Shared nodes are matched by exact geometry. Existing node ids are preserved + where possible, but ids that would otherwise collide for distinct geometries + are renamed. Edge ids are preserved where possible, but renamed if they + would otherwise collide. If topology columns are present, edge endpoints are + remapped to the merged node ids. + + Before merging, each input network's own id_col is checked for null or + duplicated values (warning if found), null ids are filled with an empty + string or 0 depending on the column's dtype, and ids are made unique within + that network. + + Parameters + ---------- + id_col : default "id" + Column name used to identify nodes and edges + """ + + n = len(networks) + if n == 0: + warnings.warn("Merging zero networks to return empty network.") + return Network() + if n == 1: + return networks[0].copy() + + used_node_ids: Set[Any] = set() + used_edge_ids: Set[Any] = set() + + def make_unique_id(id_value: Any, used_ids: Set[Any]) -> Any: + if id_value not in used_ids: + return id_value + + suffix = 1 + candidate = f"{id_value}_{suffix}" + while candidate in used_ids: + suffix += 1 + candidate = f"{id_value}_{suffix}" + return candidate + + def prepare_ids(ids: "pandas.Series[Any]", label: str) -> "pandas.Series[Any]": + # warn about null/duplicate ids within a single network's frame, + null_mask = ids.isna() + if null_mask.any(): + warnings.warn( + f"{label}: found {int(null_mask.sum())} null {id_col}(s) before merging." + ) + duplicated_mask = ids.duplicated(keep=False) & ~null_mask + if duplicated_mask.any(): + dup_values = sorted({str(value) for value in ids[duplicated_mask]}) + warnings.warn( + f"{label}: found duplicated {id_col}s before merging: {dup_values}" + ) + + # fill nulls with a dtype-appropriate placeholder + fill_value = 0 if pandas.api.types.is_numeric_dtype(ids) else "" + filled = ids.fillna(fill_value) + + # make ids unique within that frame + seen: Set[Any] = set() + unique_values = [] + for value in filled: + unique_value = make_unique_id(value, seen) + seen.add(unique_value) + unique_values.append(unique_value) + return pandas.Series(unique_values, index=ids.index) + + # Normalise to a single geometry column name across all networks first: + # concatenating frames whose active geometry columns have different + # names would otherwise leave each network's geometry in a different + # column, so groupby/drop_duplicates would treat those rows as having + # no geometry at all rather than matching them correctly. + node_geom_col = geometry_column_name(networks[0].nodes) + any_node_has_id = any(id_col in network.nodes.columns for network in networks) + + labeled_node_frames = [] + for network_idx, network in enumerate(networks): + nodes = network.nodes.copy() + if geometry_column_name(nodes) != node_geom_col: + nodes = nodes.rename_geometry(node_geom_col) + if any_node_has_id: + if id_col not in nodes.columns: + nodes[id_col] = None + if not nodes.empty: + nodes[id_col] = prepare_ids( + nodes[id_col], f"merge_networks: nodes in network {network_idx}" + ) + nodes["_source_network"] = network_idx + nodes["_source_id"] = nodes[id_col] if any_node_has_id else None + labeled_node_frames.append(nodes) + + all_nodes = pandas.concat(labeled_node_frames, axis=0, sort=False) + all_nodes.reset_index(drop=True, inplace=True) + geom_col = node_geom_col + + # Group by geometry using the same equality/hashing that + # drop_duplicate_geometries relies on, to assign a canonical id per + # shared location and record how each network's original ids map to it. + node_id_map: Dict[Tuple[int, Any], Any] = {} + canonical_id_by_geom: Dict[Any, Any] = {} + for geom, group in all_nodes.groupby(geom_col, sort=False): + assigned_id = next( + ( + original_id + for original_id in group["_source_id"] + if original_id is not None and not pandas.isna(original_id) + ), + None, + ) + if assigned_id is not None: + assigned_id = make_unique_id(assigned_id, used_node_ids) + used_node_ids.add(assigned_id) + canonical_id_by_geom[geom] = assigned_id + for network_idx, original_id in zip( + group["_source_network"], group["_source_id"] + ): + if original_id is not None and not pandas.isna(original_id): + node_id_map[(network_idx, original_id)] = assigned_id + + nodes = drop_duplicate_geometries(all_nodes) + nodes.reset_index(drop=True, inplace=True) + if any_node_has_id: + nodes[id_col] = nodes[geom_col].map(canonical_id_by_geom) + nodes = nodes.drop(columns=["_source_network", "_source_id"]) + + edge_geom_col = geometry_column_name(networks[0].edges) + any_edge_has_id = any(id_col in network.edges.columns for network in networks) + + remapped_edge_frames: List[GeoDataFrame] = [] + for network_idx, network in enumerate(networks): + edges = network.edges.copy() + if geometry_column_name(edges) != edge_geom_col: + edges = edges.rename_geometry(edge_geom_col) + if any_edge_has_id: + if id_col not in edges.columns: + edges[id_col] = None + if not edges.empty: + edges[id_col] = prepare_ids( + edges[id_col], f"merge_networks: edges in network {network_idx}" + ) + renamed_ids = [] + for original_id in edges[id_col]: + unique_id = make_unique_id(original_id, used_edge_ids) + used_edge_ids.add(unique_id) + renamed_ids.append(unique_id) + edges[id_col] = renamed_ids + if not edges.empty and {"from_id", "to_id"}.issubset(edges.columns): + edges["from_id"] = edges["from_id"].apply( + lambda value, idx=network_idx: node_id_map.get((idx, value), value) + ) + edges["to_id"] = edges["to_id"].apply( + lambda value, idx=network_idx: node_id_map.get((idx, value), value) + ) + remapped_edge_frames.append(edges) + edges = concat_dedup(remapped_edge_frames) + + return Network(nodes, edges) diff --git a/tests/test_init.py b/tests/test_init.py index 757ec4d..2d5ebbb 100644 --- a/tests/test_init.py +++ b/tests/test_init.py @@ -4,7 +4,7 @@ import pandas as pd from pandas.testing import assert_frame_equal -from pytest import fixture, mark +from pytest import fixture, mark, warns from shapely.geometry import Point, LineString, MultiPoint, MultiLineString try: @@ -826,6 +826,197 @@ def test_add_component_ids(two_components): assert all(labelled.edges.component_id == pd.Series([2, 1, 1])) +def test_merge_zero_networks(): + merged = snkit.network.merge_networks([]) + assert merged.nodes.empty + assert merged.edges.empty + + +def test_merge_one_network(split): + merged = snkit.network.merge_networks([split]) + assert_frame_equal(merged.nodes, split.nodes) + assert_frame_equal(merged.edges, split.edges) + assert merged is not split + assert merged.nodes is not split.nodes + assert merged.edges is not split.edges + + +def test_merge_networks(split): + split_abc = snkit.Network(split.nodes[:3], split.edges[:2]) + split_cd = snkit.Network(split.nodes[2:], split.edges[2:]) + merged = snkit.network.merge_networks([split_abc, split_cd]) + assert_frame_equal(merged.nodes, split.nodes) + assert_frame_equal(merged.edges, split.edges) + + +def test_merge_networks_remaps_topology_and_renames_conflicting_ids(): + a = Point((0, 0)) + b = Point((0, 2)) + c = Point((0, 1)) + + left_nodes = GeoDataFrame(data={"id": ["node_0", "node_1"]}, geometry=[a, c]) + left_edges = GeoDataFrame( + data={"id": ["edge_0"], "from_id": ["node_0"], "to_id": ["node_1"]}, + geometry=[LineString([a, c])], + ) + + right_nodes = GeoDataFrame(data={"id": ["node_0", "node_1"]}, geometry=[c, b]) + right_edges = GeoDataFrame( + data={"id": ["edge_1"], "from_id": ["node_0"], "to_id": ["node_1"]}, + geometry=[LineString([c, b])], + ) + + merged = snkit.network.merge_networks( + [snkit.Network(left_nodes, left_edges), snkit.Network(right_nodes, right_edges)] + ) + + assert list(merged.nodes.geometry) == [a, c, b] + assert list(merged.nodes.id) == ["node_0", "node_1", "node_1_1"] + assert list(merged.edges.from_id) == ["node_0", "node_1"] + assert list(merged.edges.to_id) == ["node_1", "node_1_1"] + + +def test_merge_networks_preserves_component_ids(): + a = Point((0, 0)) + b = Point((0, 2)) + c = Point((0, 1)) + + left_nodes = GeoDataFrame( + data={"component_id": [10, 10]}, geometry=[a, c] + ) + left_edges = GeoDataFrame( + data={"component_id": [10]}, geometry=[LineString([a, c])] + ) + right_nodes = GeoDataFrame( + data={"component_id": [20, 20]}, geometry=[c, b] + ) + right_edges = GeoDataFrame( + data={"component_id": [20]}, geometry=[LineString([c, b])] + ) + + merged = snkit.network.merge_networks( + [snkit.Network(left_nodes, left_edges), snkit.Network(right_nodes, right_edges)] + ) + + assert list(merged.edges.component_id) == [10, 20] + assert list(merged.nodes.component_id) == [10, 10, 20] + + +def test_merge_networks_renames_conflicting_edge_ids(): + a = Point((0, 0)) + b = Point((0, 2)) + c = Point((0, 1)) + d = Point((0, 3)) + + left_nodes = GeoDataFrame(data={"id": ["n0", "n1"]}, geometry=[a, c]) + left_edges = GeoDataFrame( + data={"id": ["edge_0"], "from_id": ["n0"], "to_id": ["n1"]}, + geometry=[LineString([a, c])], + ) + + right_nodes = GeoDataFrame(data={"id": ["n2", "n3"]}, geometry=[b, d]) + right_edges = GeoDataFrame( + data={"id": ["edge_0"], "from_id": ["n2"], "to_id": ["n3"]}, + geometry=[LineString([b, d])], + ) + + merged = snkit.network.merge_networks( + [snkit.Network(left_nodes, left_edges), snkit.Network(right_nodes, right_edges)] + ) + + assert len(merged.edges) == 2 + assert len(set(merged.edges.id)) == 2, "edge ids should be unique after merge" + assert list(merged.edges.id) == ["edge_0", "edge_0_1"] + + +def test_merge_networks_preserves_non_identifier_column_names(): + a = Point((0, 0)) + b = Point((0, 2)) + + left_nodes = GeoDataFrame(data={"id": ["n0"], "pop 2020": [100]}, geometry=[a]) + right_nodes = GeoDataFrame(data={"id": ["n1"]}, geometry=[b]) + + merged = snkit.network.merge_networks( + [snkit.Network(left_nodes), snkit.Network(right_nodes)] + ) + + assert "pop 2020" in merged.nodes.columns + values = list(merged.nodes["pop 2020"]) + assert values[0] == 100 + assert pd.isna(values[1]) + + +def test_merge_networks_handles_mismatched_geometry_column_names(): + a = Point((0, 0)) + b = Point((0, 2)) + c = Point((0, 1)) + + left_nodes = GeoDataFrame(data={"id": ["n0", "n1"]}, geometry=[a, c]).rename_geometry( + "geom" + ) + right_nodes = GeoDataFrame(data={"id": ["n2", "n3"]}, geometry=[c, b]) + + merged = snkit.network.merge_networks( + [snkit.Network(left_nodes), snkit.Network(right_nodes)] + ) + + # 3 distinct locations (a, c, b) should survive -- c is shared, not a + # 4th node lost because left/right use differently-named geometry columns + assert len(merged.nodes) == 3 + assert set(merged.nodes.geometry) == {a, b, c} + + +def test_merge_networks_id_col_parameter(): + a = Point((0, 0)) + b = Point((0, 2)) + c = Point((0, 1)) + + left_nodes = GeoDataFrame(data={"asset_id": ["n0", "n1"]}, geometry=[a, c]) + right_nodes = GeoDataFrame(data={"asset_id": ["n0", "n1"]}, geometry=[c, b]) + + merged = snkit.network.merge_networks( + [snkit.Network(left_nodes), snkit.Network(right_nodes)], id_col="asset_id" + ) + + assert "id" not in merged.nodes.columns + assert list(merged.nodes.asset_id) == ["n0", "n1", "n1_1"] + + +def test_merge_networks_warns_and_renames_duplicate_ids_within_one_network(): + a = Point((0, 0)) + b = Point((0, 2)) + c = Point((0, 1)) + + # left has a pre-existing (invalid) duplicate id "n0" at two locations + left_nodes = GeoDataFrame(data={"id": ["n0", "n0"]}, geometry=[a, c]) + right_nodes = GeoDataFrame(data={"id": ["n2"]}, geometry=[b]) + + with warns(UserWarning, match="duplicated"): + merged = snkit.network.merge_networks( + [snkit.Network(left_nodes), snkit.Network(right_nodes)] + ) + + # ids are unique in the output -- no silent id_map corruption/overwrite + assert list(merged.nodes.id) == ["n0", "n0_1", "n2"] + + +def test_merge_networks_warns_and_fills_null_ids(): + a = Point((0, 0)) + b = Point((0, 2)) + + left_nodes = GeoDataFrame(data={"id": ["n0"]}, geometry=[a]) + right_nodes = GeoDataFrame(geometry=[b]) # no id column at all + + with warns(UserWarning, match="null"): + merged = snkit.network.merge_networks( + [snkit.Network(left_nodes), snkit.Network(right_nodes)] + ) + + # no null ids survive in the merged output + assert not merged.nodes.id.isna().any() + assert list(merged.nodes.id) == ["n0", ""] + + def test_matching_gdf_from_geoms(edge_only): expected = edge_only.edges.copy() gdf = edge_only.edges.copy()