Merge pull request #143 from ImperialCollegeLondon/128-pipes-not-unde…

…r-roads

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
@@ -313,6 +317,10 @@ def __call__(self,
                 **kwargs) -> nx.Graph:
         """Format the lanes attribute of each edge and calculates width.
 
+        Only the `drive` network is assumed to contribute to impervious area and 
+        so others `network_types` have lanes set to 0. If no `network_type` is
+        present, the edge is assumed to be of type `drive`. 
+
         Args:
             G (nx.Graph): A graph
             subcatchment_derivation (parameters.SubcatchmentDerivation): A
@@ -326,7 +334,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:
@@ -478,15 +489,15 @@ def __call__(self,
                          )
 
 @register_graphfcn
-class merge_nodes(BaseGraphFunction):
+class merge_street_nodes(BaseGraphFunction):
     """merge_nodes class."""
     def __call__(self, 
                  G: nx.Graph, 
                  subcatchment_derivation: parameters.SubcatchmentDerivation,
                  **kwargs) -> nx.Graph:
         """Merge nodes that are close together.
 
-        This function merges nodes that are within a certain distance of each
+        Merges `street` nodes that are within a certain distance of each
         other. The distance is specified in the `node_merge_distance` attribute
         of the `subcatchment_derivation` parameter. The merged nodes are given
         the same coordinates, and the graph is relabeled with nx.relabel_nodes.
@@ -503,16 +514,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

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,47 @@ 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 'drive' in network_types and network_types[-1] != 'drive':
+        logger.warning("""The last network type should be `drive` to retain 
+                        `lanes` attribute, needed to calculate impervious area.
+                        Moving it to the last position.""")
+        network_types.pop("drive")
+        network_types.append("drive")
+    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)
+
+    # 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 +280,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 +293,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 +307,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_data/demo_config.yml

@@ -20,8 +20,8 @@ graphfcn_list:
   - remove_parallel_edges # Remove parallel edges retaining the shorter one
   - to_undirected # Convert graph to undirected to facilitate cleanup
   - split_long_edges # Set a maximum edge length
-  - merge_nodes # Merge nodes that are too close together
-  - assign_id # Remove duplicates arising from merge_nodes
+  - merge_street_nodes # Merge street nodes that are too close together
+  - assign_id # Remove duplicates arising from merge_street_nodes
   - clip_to_catchments # Clip graph to catchment subbasins
   - calculate_contributing_area # Calculate runoff coefficient
   - set_elevation # Set node elevation from DEM

tests/test_graph_utilities.py

@@ -347,12 +347,12 @@ def almost_equal(a, b, tol=1e-6):
     """Check if two numbers are almost equal."""
     return abs(a-b) < tol
 
-def test_merge_nodes():
-    """Test the merge_nodes function."""
+def test_merge_street_nodes():
+    """Test the merge_street_nodes function."""
     G, _ = load_street_network()
     subcatchment_derivation = parameters.SubcatchmentDerivation(
         node_merge_distance = 20)
-    G_ = gu.merge_nodes(G, subcatchment_derivation)
+    G_ = gu.merge_street_nodes(G, subcatchment_derivation)
     assert not set([107736,266325461,2623975694,32925453]).intersection(G_.nodes)
     assert almost_equal(G_.nodes[25510321]['x'], 700445.0112082)
 

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')