Pedestrian pipes

swmmanywhere/graph_utilities.py

@@ -184,8 +184,12 @@ def iterate_graphfcns(G: nx.Graph,
         G = graphfcns[function](G, addresses = addresses, **params)
         logger.info(f"graphfcn: {function} completed.")
         if verbose:
-            save_graph(graphfcns.fix_geometries(G), 
-                       addresses.model / f"{function}_graph.json")
+            save_graph(G, addresses.model / f"{function}_graph.json")
+            go.graph_to_geojson(graphfcns.fix_geometries(G),
+                                addresses.model / f"{function}_nodes.geojson",
+                                addresses.model / f"{function}_edges.geojson",
+                                G.graph['crs']
+                                )
     return G
 
 @register_graphfcn
@@ -326,7 +330,10 @@ def __call__(self,
         G = G.copy()
         lines = []
         for u, v, data in G.edges(data=True):
-            lanes = data.get('lanes',1)
+            if data.get('network_type','drive') == 'drive':
+                lanes = data.get('lanes',1)
+            else:
+                lanes = 0
             if isinstance(lanes, list):
                 lanes = sum([float(x) for x in lanes])
             else:
@@ -503,16 +510,22 @@ def __call__(self,
         """
         G = G.copy()
 
+        # Separate out streets         
+        street_edges = [(u, v, k) for u, v, k, d in G.edges(data=True, keys=True)
+                        if d.get('edge_type','street') == 'street']
+        streets = G.edge_subgraph(street_edges).copy()
+
         # Identify nodes that are within threshold of each other
-        mapping = go.merge_points([(d['x'], d['y']) for u,d in G.nodes(data=True)],
+        mapping = go.merge_points([(d['x'], d['y']) 
+                                   for u,d in streets.nodes(data=True)],
                               subcatchment_derivation.node_merge_distance)
 
         # Get indexes of node names
-        node_indices = {ix: node for ix, node in enumerate(G.nodes)}
+        node_indices = {ix: node for ix, node in enumerate(streets.nodes)}
 
         # Create a mapping of old node names to new node names
         node_names = {}
-        for ix, node in enumerate(G.nodes):
+        for ix, node in enumerate(streets.nodes):
             if ix in mapping:
                 # If the node is in the mapping, then it is mapped and 
                 # given the new coordinate (all nodes in a mapping family must
@@ -1033,28 +1046,6 @@ def __call__(self,
         # Copy graph to run shortest path on
         G_ = G.copy()
 
-        if not matched_outlets:
-            # In cases of e.g., an area with no rivers to discharge into or too
-            # small a buffer
-
-            # Identify the lowest elevation node
-            lowest_elevation_node = min(G.nodes, 
-                                    key = lambda x: G.nodes[x]['surface_elevation'])
-
-            # Create a dummy river to discharge into
-            dummy_river = {'id' : 'dummy_river',
-                           'x' : G.nodes[lowest_elevation_node]['x'] + 1,
-                           'y' : G.nodes[lowest_elevation_node]['y'] + 1}
-            G_.add_node('dummy_river')
-            nx.set_node_attributes(G_, {'dummy_river' : dummy_river})
-
-            # Update function's dicts
-            matched_outlets = {'dummy_river' : lowest_elevation_node}
-            river_points['dummy_river'] = shapely.Point(dummy_river['x'],
-                                                        dummy_river['y'])
-
-            logger.warning('No outlets found, using lowest elevation node as outlet')
-
         # Add edges between the paired river and street nodes
         for river_id, street_id in matched_outlets.items():
             # TODO instead use a weight based on the distance between the two nodes

swmmanywhere/paper/experimenter.py

@@ -137,6 +137,7 @@ def process_parameters(jobid: int,
     df = pd.DataFrame(X)
     gb = df.groupby('iter')
     n_iter = len(gb)
+    logger.info(f"{n_iter} samples created")
     flooding_results = {}
     nproc = nproc if nproc is not None else n_iter
 
@@ -167,11 +168,11 @@ def process_parameters(jobid: int,
         # Run the model
         config['model_number'] = ix
         logger.info(f"Running swmmanywhere for model {ix}")
-        address, metrics = swmmanywhere.swmmanywhere(config)
 
+        address, metrics = swmmanywhere.swmmanywhere(config)
         if metrics is None:
             raise ValueError(f"Model run {ix} failed.")
-        
+
         # Save the results
         flooding_results[ix] = {'iter': ix, 
                                 **metrics, 
@@ -206,7 +207,7 @@ def parse_arguments() -> tuple[int, int | None, Path]:
     parser = argparse.ArgumentParser(description='Process command line arguments.')
     parser.add_argument('--jobid', 
                         type=int, 
-                        default=1, 
+                        default=0, 
                         help='Job ID')
     parser.add_argument('--nproc', 
                         type=int, 
@@ -215,7 +216,7 @@ def parse_arguments() -> tuple[int, int | None, Path]:
     parser.add_argument('--config_path', 
                         type=Path, 
                         default=Path(__file__).parent.parent.parent / 'tests' /\
-                                    'test_data' / 'demo_config_sa.yml',
+                                    'test_data' / 'demo_config.yml',
                         help='Configuration file path')
 
     args = parser.parse_args()

swmmanywhere/parameters.py

@@ -85,6 +85,11 @@ class OutletDerivation(BaseModel):
 
 class TopologyDerivation(BaseModel):
     """Parameters for topology derivation."""
+    allowable_networks: list = Field(default = ['walk', 'drive'],
+                                     min_items = 1,
+                        unit = "-",
+                        description = "OSM networks to consider")
+
     weights: list = Field(default = ['chahinian_slope',
                                       'chahinian_angle',
                                       'length',
@@ -96,7 +101,8 @@ class TopologyDerivation(BaseModel):
     omit_edges: list = Field(default = ['motorway', 
                                         'motorway_link',
                                         'bridge', 
-                                        'tunnel'],
+                                        'tunnel',
+                                        'corridor'],
                         min_items = 1,
                         unit = "-",
                         description = "OSM paths pipes are not allowed under")
@@ -207,8 +213,8 @@ class HydraulicDesign(BaseModel):
 class MetricEvaluation(BaseModel):
     """Parameters for metric evaluation."""
     grid_scale: float = Field(default = 100,
-                        le = 10,
-                        ge = 5000,
+                        le = 5000,
+                        ge = 10,
                         unit = "m",
                         description = "Scale of the grid for metric evaluation")
 

swmmanywhere/prepare_data.py

@@ -89,14 +89,16 @@ def download_buildings(file_address: Path,
         logger.error(f"Error downloading data. Status code: {response.status_code}")
     return response.status_code
 
-def download_street(bbox: tuple[float, float, float, float]) -> nx.MultiDiGraph:
+def download_street(bbox: tuple[float, float, float, float],
+                    network_type = 'drive') -> nx.MultiDiGraph:
     """Get street network within a bounding box using OSMNX.
 
     [CREDIT: Taher Cheghini busn_estimator package]
 
     Args:
         bbox (tuple[float, float, float, float]): Bounding box as tuple in form 
             of (west, south, east, north) at EPSG:4326.
+        network_type (str, optional): Type of street network. Defaults to 'drive'.
 
     Returns:
         nx.MultiDiGraph: Street network with type drive and 
@@ -105,8 +107,8 @@ def download_street(bbox: tuple[float, float, float, float]) -> nx.MultiDiGraph:
     west, south, east, north = bbox
     bbox = (north, south, east, west) # not sure why osmnx uses this order
     graph = ox.graph_from_bbox(
-        bbox = bbox, network_type="drive", truncate_by_edge=True
-    )
+            bbox = bbox, network_type=network_type, truncate_by_edge=True
+        )
     return cast("nx.MultiDiGraph", graph)
 
 def download_river(bbox: tuple[float, float, float, float]) -> nx.MultiDiGraph:
@@ -122,10 +124,7 @@ def download_river(bbox: tuple[float, float, float, float]) -> nx.MultiDiGraph:
     """
     west, south, east, north = bbox
     graph = ox.graph_from_bbox(
-        north, 
-        south, 
-        east, 
-        west,
+        bbox = (north, south, east, west),
         truncate_by_edge=True, 
         custom_filter='["waterway"]')
 

swmmanywhere/preprocessing.py

@@ -221,15 +221,45 @@ def prepare_building(bbox: tuple[float, float, float, float],
 def prepare_street(bbox: tuple[float, float, float, float],
                      addresses: parameters.FilePaths,
                      target_crs: str,
-                     source_crs: str = 'EPSG:4326'):
-    """Download and reproject street graph."""
+                     source_crs: str = 'EPSG:4326',
+                     network_types = ['drive']):
+    """Download and reproject street graph.
+
+    Download the street graph within the bbox and reproject it to the UTM zone.
+    The street graph is downloaded for all network types in network_types. The
+    street graph is saved to the addresses.street directory.
+
+    Args:
+        bbox (tuple[float, float, float, float]): Bounding box coordinates in 
+            the format (minx, miny, maxx, maxy) in EPSG:4326.
+        addresses (FilePaths): Class containing the addresses of the directories.
+        target_crs (str): Target CRS to reproject the graph to.
+        source_crs (str): Source CRS of the graph.
+        network_types (list): List of network types to download. For duplicate 
+            edges, nx.compose_all selects the attributes in priority of last to
+            first. In likelihood, you want to ensure that the last network in 
+            the list is `drive`, so as to retain information about `lanes`, 
+            which is needed to calculate impervious area.
+    """
     if addresses.street.exists():
         return
-    logger.info(f'downloading street network to {addresses.street}')
-    street_network = prepare_data.download_street(bbox)
+    logger.info(f'downloading network to {addresses.street}')
+    if network_types[-1] != 'drive':
+        logger.warning("""The last network type should be `drive` to retain 
+                       `lanes` attribute, needed to calculate impervious area.""")
+    networks = []
+    for network_type in network_types:
+        network = prepare_data.download_street(bbox, network_type=network_type)
+        nx.set_edge_attributes(network, network_type, 'network_type')
+        networks.append(network)
+    street_network = nx.compose_all(networks)
+    street_network = prepare_data.download_street(bbox, network_type='drive')
+
+    # Reproject graph
     street_network = go.reproject_graph(street_network, 
                                         source_crs, 
                                         target_crs)
+
     gu.save_graph(street_network, addresses.street)
 
 def prepare_river(bbox: tuple[float, float, float, float],
@@ -248,7 +278,8 @@ def prepare_river(bbox: tuple[float, float, float, float],
 
 def run_downloads(bbox: tuple[float, float, float, float],
                   addresses: parameters.FilePaths,
-                  api_keys: dict[str, str]):
+                  api_keys: dict[str, str],
+                  network_types = ['drive']):
     """Run the data downloads.
 
     Run the precipitation, elevation, building, street and river network
@@ -260,6 +291,7 @@ def run_downloads(bbox: tuple[float, float, float, float],
             the format (minx, miny, maxx, maxy) in EPSG:4326.
         addresses (FilePaths): Class containing the addresses of the directories.
         api_keys (dict): Dictionary containing the API keys.
+        network_types (list): List of network types to download.
     """
     target_crs = go.get_utm_epsg(bbox[0], bbox[1])
 
@@ -273,7 +305,7 @@ def run_downloads(bbox: tuple[float, float, float, float],
     prepare_building(bbox, addresses, target_crs)
 
     # Download street network data
-    prepare_street(bbox, addresses, target_crs)
+    prepare_street(bbox, addresses, target_crs, network_types=network_types)
 
     # Download river network data
     prepare_river(bbox, addresses, target_crs)

swmmanywhere/swmmanywhere.py

@@ -45,13 +45,21 @@ def swmmanywhere(config: dict) -> tuple[Path, dict | None]:
     for key, val in config.get('address_overrides', {}).items():
         setattr(addresses, key, val)
 
+    # Load the parameters and perform any manual overrides
+    logger.info("Loading and setting parameters.")
+    params = parameters.get_full_parameters()
+    for category, overrides in config.get('parameter_overrides', {}).items():
+        for key, val in overrides.items():
+            setattr(params[category], key, val)
+
     # Run downloads
     logger.info("Running downloads.")
     api_keys = yaml.safe_load(config['api_keys'].open('r'))
     preprocessing.run_downloads(config['bbox'],
-                                addresses,
-                                api_keys
-                                )
+                addresses,
+                api_keys,
+                network_types = params['topology_derivation'].allowable_networks
+                )
 
     # Identify the starting graph
     logger.info("Iterating graphs.")
@@ -60,13 +68,6 @@ def swmmanywhere(config: dict) -> tuple[Path, dict | None]:
     else:
         G = preprocessing.create_starting_graph(addresses)
 
-    # Load the parameters and perform any manual overrides
-    logger.info("Loading and setting parameters.")
-    params = parameters.get_full_parameters()
-    for category, overrides in config.get('parameter_overrides', {}).items():
-        for key, val in overrides.items():
-            setattr(params[category], key, val)
-
     # Iterate the graph functions
     logger.info("Iterating graph functions.")
     G = iterate_graphfcns(G, 

tests/test_graph_utilities.py

@@ -179,26 +179,6 @@ def test_calculate_weights():
     for u, v, data in G.edges(data=True):
         assert 'weight' in data.keys()
         assert math.isfinite(data['weight'])
-
-def test_identify_outlets_no_river():
-    """Test the identify_outlets in the no river case."""
-    G, _ = load_street_network()
-    G = gu.assign_id(G)
-    G = gu.double_directed(G)
-    elev_fid = Path(__file__).parent / 'test_data' / 'elevation.tif'
-    addresses = parameters.FilePaths(base_dir = None,
-                                    project_name = None,
-                                    bbox_number = None,
-                                    model_number = None)
-    addresses.elevation = elev_fid
-    G = gu.set_elevation(G, addresses)
-    for ix, (u,v,d) in enumerate(G.edges(data=True)):
-        d['edge_type'] = 'street'
-        d['weight'] = ix
-    params = parameters.OutletDerivation()
-    G = gu.identify_outlets(G, params)
-    outlets = [(u,v,d) for u,v,d in G.edges(data=True) if d['edge_type'] == 'outlet']
-    assert len(outlets) == 1
 
 def test_identify_outlets_and_derive_topology():
     """Test the identify_outlets and derive_topology functions."""

tests/test_swmmanywhere.py

@@ -54,7 +54,7 @@ def test_swmmanywhere():
         # Set some test values
         base_dir = Path(temp_dir)
         config['base_dir'] = str(base_dir)
-        config['bbox'] = [0.05428,51.55847,0.07193,51.56726]
+        config['bbox'] = [0.05677,51.55656,0.07193,51.56726]
         config['address_overrides'] = {
             'building': str(test_data_dir / 'building.geoparquet'),
             'precipitation': str(defs_dir / 'storm.dat')