test and debug

swmmanywhere/geospatial_utilities.py

@@ -673,7 +673,7 @@ def derive_subbasins_streamorder(fid: Path,
         streamorder -= 1
 
     if (
-        (streamorder != streamorder_) & 
+        (streamorder != streamorder_ - 1) & 
         (os.getenv("SWMMANYWHERE_VERBOSE", "false").lower() == "true")
         ):
         logger.warning(f"""Stream order {streamorder_} resulted in no subbasins. 

swmmanywhere/graph_utilities.py

@@ -622,7 +622,8 @@ def __call__(self,
 
         # Derive road network clusters
         louv_membership = nx.community.louvain_communities(street,
-                                                           weight = 'length')
+                                                           weight = 'length',
+                                                           seed = 1)
 
         # Create gdf of street points
         street_points = gpd.GeoDataFrame(G.nodes,
@@ -686,13 +687,15 @@ def __call__(self,
         # Cut links between communities in community_omit and commuities in those
         # basins
         arcs_to_remove = []
+        street_points = street_points.reset_index().set_index('basin')
         for idx, row in community_omit.iterrows():
             community_nodes = louv_membership[int(row['community'])]
-            basin_nodes = community_basin.loc[row['basin']]['community']
+            basin_nodes = street_points.loc[[row['basin']],'id']
+            basin_nodes = set(basin_nodes).difference(community_nodes)
             arcs_to_remove.extend(
-                [(u, v) for u, v in product(community_nodes, basin_nodes)]
+                [(u, v, 0) for u, v in product(community_nodes, basin_nodes)]
                 )
-        G.remove_edges_from(arcs_to_remove)
+        G.remove_edges_from(set(G.edges).intersection(arcs_to_remove))
 
         return G
 
@@ -836,8 +839,15 @@ def __call__(self,
             keep_nodes = gpd.sjoin(nodes_gdf, 
                                    outlet_catchments[['geometry']], 
                                    predicate = 'intersects').id
-            G = G.subgraph(keep_nodes).copy()
-            G.graph = graph_
+            G_ = G.subgraph(keep_nodes).copy()
+            G_.graph = graph_
+            for node in keep_nodes:
+                path_exists = any(nx.has_path(G_, node, target) for target in outlets)
+                if path_exists:        
+                    return G_
+
+            logger.warning("""trim_to_outlets removed any paths between nodes and
+                            outlets. Returning original graph""")
             return G
 
 @register_graphfcn

tests/test_graph_utilities.py

@@ -6,6 +6,7 @@
 from __future__ import annotations
 
 import math
+import os
 import tempfile
 from pathlib import Path
 
@@ -16,6 +17,7 @@
 from swmmanywhere import parameters
 from swmmanywhere.graph_utilities import graphfcns as gu
 from swmmanywhere.graph_utilities import iterate_graphfcns, load_graph, save_graph
+from swmmanywhere.logging import logger
 
 
 def load_street_network():
@@ -365,4 +367,59 @@ def test_merge_nodes():
         node_merge_distance = 20)
     G_ = gu.merge_nodes(G, subcatchment_derivation)
     assert not set([107736,266325461,2623975694,32925453]).intersection(G_.nodes)
-    assert almost_equal(G_.nodes[25510321]['x'], 700445.0112082)
+    assert almost_equal(G_.nodes[25510321]['x'], 700445.0112082)
+
+def test_clip_to_catchments():
+    """Test the clip_to_catchments function."""
+    G, _ = load_street_network()
+    subcatchment_derivation = parameters.SubcatchmentDerivation(
+        subbasin_streamorder = 1,
+        subbasin_membership = 0.9
+    )
+    with tempfile.TemporaryDirectory() as temp_dir:
+
+        os.environ['SWMMANYWHERE_VERBOSE'] = 'true'
+        temp_path = Path(temp_dir)
+        addresses = parameters.FilePaths(base_dir = temp_path,
+                                        project_name = 'test',
+                                        bbox_number = 0,
+                                        model_number = 0)
+        addresses.nodes = addresses.base_dir / 'nodes.geojson'
+        addresses.elevation = Path(__file__).parent / 'test_data' / 'elevation.tif'
+
+        # Test clipping
+        subcatchment_derivation = parameters.SubcatchmentDerivation(
+            subbasin_streamorder = 3,
+            subbasin_membership = 0.9
+        )
+        G_ = gu.clip_to_catchments(G, 
+                                addresses=addresses,
+                                subcatchment_derivation=subcatchment_derivation)
+        assert len(G_.edges) == 30
+
+        # Test clipping with different params
+        subcatchment_derivation = parameters.SubcatchmentDerivation(
+            subbasin_streamorder = 4,
+            subbasin_membership = 0.5
+        )
+        G_ = gu.clip_to_catchments(G, 
+                                addresses=addresses,
+                                subcatchment_derivation=subcatchment_derivation)
+        assert len(G_.edges) == 32
+
+        # Check streamorder adjustment
+        os.environ['SWMMANYWHERE_VERBOSE'] = 'true'
+        logger.add(temp_path / 'test.log')
+        subcatchment_derivation = parameters.SubcatchmentDerivation(
+            subbasin_streamorder = 5,
+            subbasin_membership = 0.9
+        )
+        G_ = gu.clip_to_catchments(G, 
+                                addresses=addresses,
+                                subcatchment_derivation=subcatchment_derivation)
+        with open(temp_path / 'test.log', 'r') as f:
+            log = f.read()
+            assert 'Stream order 5 resulted in no subbasins' in log
+            assert 'Using 4 instead' in log
+        logger.remove()
+