Add streetcover to shapefile generation

swmmanywhere/geospatial_utilities.py

@@ -614,8 +614,8 @@ def remove_(mp): return remove_zero_area_subareas(mp, removed_subareas)
     return polys_gdf
 
 def derive_rc(polys_gdf: gpd.GeoDataFrame,
-              G: nx.Graph,
-              building_footprints: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
+              building_footprints: gpd.GeoDataFrame,
+              streetcover: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
     """Derive the Runoff Coefficient (RC) of each subcatchment.
 
     The runoff coefficient is the ratio of impervious area to total area. The
@@ -626,10 +626,10 @@ def derive_rc(polys_gdf: gpd.GeoDataFrame,
     Args:
         polys_gdf (gpd.GeoDataFrame): A GeoDataFrame containing polygons that
             represent subcatchments with columns: 'geometry', 'area', and 'id'. 
-        G (nx.Graph): The input graph, with node 'ids' that match polys_gdf and
-            edges with the 'id', 'width' and 'geometry' property.
         building_footprints (gpd.GeoDataFrame): A GeoDataFrame containing 
             building footprints with a 'geometry' column.
+        streetcover (gpd.GeoDataFrame): A GeoDataFrame containing street cover
+            with a 'geometry' column.
 
     Returns:
         gpd.GeoDataFrame: A GeoDataFrame containing polygons with columns:
@@ -638,23 +638,7 @@ def derive_rc(polys_gdf: gpd.GeoDataFrame,
     polys_gdf = polys_gdf.copy()
 
     ## Format as swmm type catchments 
-
-    # TODO think harder about lane widths (am I double counting here?)
-    lines = [
-        {
-            'geometry': x['geometry'].buffer(x['width'], 
-                                             cap_style=2, 
-                                             join_style=2),
-            'id': x['id']
-        }
-        for u, v, x in G.edges(data=True)
-    ]
-    lines_df = pd.DataFrame(lines)
-    lines_gdf = gpd.GeoDataFrame(lines_df, 
-                                geometry=lines_df.geometry,
-                                    crs = polys_gdf.crs)
-
-    result = gpd.overlay(lines_gdf[['geometry']], 
+    result = gpd.overlay(streetcover[['geometry']], 
                          building_footprints[['geometry']], 
                          how='union')
     result = gpd.overlay(polys_gdf, result)

swmmanywhere/graph_utilities.py

@@ -266,36 +266,55 @@ def __call__(self,
         return G
 
 @register_graphfcn
-class format_osmnx_lanes(BaseGraphFunction,
-                         required_edge_attributes = ['lanes'],
-                         adds_edge_attributes = ['width']):
-    """format_osmnx_lanes class."""
+class calculate_streetcover(BaseGraphFunction,
+                         required_edge_attributes = ['lanes']
+                         ):
+    """calculate_streetcover class."""
     # i.e., in osmnx format, i.e., empty for single lane, an int for a
     # number of lanes or a list if the edge has multiple carriageways
 
     def __call__(self,
                  G: nx.Graph, 
-                       subcatchment_derivation: parameters.SubcatchmentDerivation, 
-                       **kwargs) -> nx.Graph:
+                subcatchment_derivation: parameters.SubcatchmentDerivation, 
+                addresses: parameters.FilePaths,
+                **kwargs) -> nx.Graph:
         """Format the lanes attribute of each edge and calculates width.
 
         Args:
             G (nx.Graph): A graph
             subcatchment_derivation (parameters.SubcatchmentDerivation): A
                 SubcatchmentDerivation parameter object
+            addresses (parameters.FilePaths): A FilePaths parameter object
             **kwargs: Additional keyword arguments are ignored.
 
         Returns:
             G (nx.Graph): A graph
         """
         G = G.copy()
+        lines = []
         for u, v, data in G.edges(data=True):
             lanes = data.get('lanes',1)
             if isinstance(lanes, list):
                 lanes = sum([float(x) for x in lanes])
             else:
                 lanes = float(lanes)
-            data['width'] = lanes * subcatchment_derivation.lane_width
+            lines.append({'geometry' : data['geometry'].buffer(lanes * 
+                                subcatchment_derivation.lane_width,
+                                cap_style=2,
+                                join_style=2),
+                            'u' : u,
+                            'v' : v
+                            }
+                        )
+        lines_df = pd.DataFrame(lines)
+        lines_gdf = gpd.GeoDataFrame(lines_df, 
+                                    geometry=lines_df.geometry,
+                                        crs = G.graph['crs'])
+        if addresses.streetcover.suffix in ('.geoparquet','.parquet'):
+            lines_gdf.to_parquet(addresses.streetcover)
+        else:
+            lines_gdf.to_file(addresses.streetcover, driver='GeoJSON')
+
         return G
 
 @register_graphfcn
@@ -494,7 +513,7 @@ def __call__(self, G: nx.Graph, **kwargs) -> nx.Graph:
 
 @register_graphfcn
 class calculate_contributing_area(BaseGraphFunction,
-                                required_edge_attributes = ['id', 'geometry', 'width'],
+                                required_edge_attributes = ['id', 'geometry'],
                                 adds_edge_attributes = ['contributing_area'],
                                 adds_node_attributes = ['contributing_area']):
     """calculate_contributing_area class."""
@@ -507,7 +526,10 @@ def __call__(self, G: nx.Graph,
         
         This function calculates the contributing area for each edge. The
         contributing area is the area of the subcatchment that drains to the
-        edge. The contributing area is calculated from the elevation data.
+        edge. The contributing area is calculated from the elevation data. 
+        Runoff coefficient (RC) for each contributing area is also calculated, 
+        the RC is calculated using `addresses.buildings` and 
+        `addresses.streetcover`.        
 
         Also writes the file 'subcatchments.geojson' to addresses.subcatchments.
 
@@ -543,7 +565,12 @@ def __call__(self, G: nx.Graph,
             buildings = gpd.read_parquet(addresses.building)
         else:
             buildings = gpd.read_file(addresses.building)
-        subs_rc = go.derive_rc(subs_gdf, G, buildings)
+        if addresses.streetcover.suffix in ('.geoparquet','.parquet'):
+            streetcover = gpd.read_parquet(addresses.streetcover)
+        else:
+            streetcover = gpd.read_file(addresses.streetcover)
+
+        subs_rc = go.derive_rc(subs_gdf, buildings, streetcover)
 
         # Write subs
         # TODO - could just attach subs to nodes where each node has a list of subs

swmmanywhere/parameters.py

@@ -292,6 +292,9 @@ def _generate_elevation(self):
     def _generate_building(self):
         return self._generate_property(f'building.geo{self.extension}', 
                                         'download')
+    def _generate_streetcover(self):
+        return self._generate_property(f'streetcover.geo{self.extension}', 
+                                        'model')
     def _generate_precipitation(self):
         return self._generate_property(f'precipitation.{self.extension}', 
                                         'download')

tests/test_data/demo_config.yml

@@ -14,8 +14,8 @@ real:
 starting_graph: null
 graphfcn_list:
   - assign_id
-  - format_osmnx_lanes
   - remove_non_pipe_allowable_links
+  - calculate_streetcover
   - double_directed
   - fix_geometries
   - split_long_edges

tests/test_geospatial_utilities.py

@@ -279,20 +279,16 @@ def test_derive_rc():
                                     crs = crs)
     subs['area'] = subs.geometry.area
 
-    subs_rc = go.derive_rc(subs, G, buildings).set_index('id')
+    subs_rc = go.derive_rc(subs, buildings, buildings).set_index('id')
     assert subs_rc.loc[6277683849,'impervious_area'] == 0
     assert subs_rc.loc[107733,'impervious_area'] > 0
-    for u,v,d in G.edges(data=True):
-        d['width'] = 10
 
-    subs_rc = go.derive_rc(subs, G, buildings).set_index('id')
+    buildings.geometry = buildings.buffer(50)
+    subs_rc = go.derive_rc(subs, buildings, buildings).set_index('id')
     assert subs_rc.loc[6277683849,'impervious_area'] > 0
     assert subs_rc.loc[6277683849,'rc'] > 0
     assert subs_rc.rc.max() <= 100
 
-    for u,v,d in G.edges(data=True):
-        d['width'] = 0
-
 def test_calculate_angle():
     """Test the calculate_angle function."""
     # Test with points forming a right angle

tests/test_graph_utilities.py

@@ -52,14 +52,23 @@ def test_double_directed():
     for u, v in G.edges():
         assert (v,u) in G.edges
 
-def test_format_osmnx_lanes():
-    """Test the format_osmnx_lanes function."""
+def test_calculate_streetcover():
+    """Test the calculate_streetcover function."""
     G, _ = load_street_network()
     params = parameters.SubcatchmentDerivation()
-    G = gu.format_osmnx_lanes(G, params)
-    for u, v, data in G.edges(data=True):
-        assert 'width' in data.keys()
-        assert isinstance(data['width'], float)
+    addresses = parameters.FilePaths(base_dir = None,
+                                        project_name = None,
+                                        bbox_number = None,
+                                        model_number = None,
+                                        extension = 'json')
+    with tempfile.TemporaryDirectory() as temp_dir:
+        addresses.streetcover = Path(temp_dir) / 'streetcover.geojson'
+        _ = gu.calculate_streetcover(G, params, addresses)
+        # TODO test that G hasn't changed? or is that a waste of time?
+        assert addresses.streetcover.exists()
+        gdf = gpd.read_file(addresses.streetcover)
+        assert len(gdf) == len(G.edges)
+        assert gdf.geometry.area.sum() > 0
 
 def test_split_long_edges():
     """Test the split_long_edges function."""
@@ -82,6 +91,7 @@ def test_derive_subcatchments():
                             model_number = 1)
         addresses.elevation = Path(__file__).parent / 'test_data' / 'elevation.tif'
         addresses.building = temp_path / 'building.geojson'
+        addresses.streetcover = temp_path / 'building.geojson'
         addresses.subcatchments = temp_path / 'subcatchments.geojson'
         params = parameters.SubcatchmentDerivation()
         G, bbox = load_street_network()
@@ -275,6 +285,7 @@ def test_iterate_graphfcns():
     """Test the iterate_graphfcns function."""
     G = load_graph(Path(__file__).parent / 'test_data' / 'graph_topo_derived.json')
     params = parameters.get_full_parameters()
+    params['topology_derivation'].omit_edges = ['primary', 'bridge']
     with tempfile.TemporaryDirectory() as temp_dir:
         temp_path = Path(temp_dir)
         addresses = parameters.FilePaths(base_dir = None,
@@ -286,12 +297,13 @@ def test_iterate_graphfcns():
         addresses.model = temp_path
         G = iterate_graphfcns(G, 
                                 ['assign_id',
-                                'format_osmnx_lanes'],
+                                'remove_non_pipe_allowable_links'],
                                 params, 
                                 addresses)
         for u, v, d in G.edges(data=True):
             assert 'id' in d.keys()
-            assert 'width' in d.keys()
+        assert 'primary' not in get_edge_types(G)
+        assert len(set([d.get('bridge',None) for u,v,d in G.edges(data=True)])) == 1
 
 def test_fix_geometries():
     """Test the fix_geometries function."""