Merge branch 'main' into 128-pipes-not-under-roads

swmmanywhere/geospatial_utilities.py

@@ -699,6 +699,7 @@ def remove_(mp): return remove_zero_area_subareas(mp, removed_subareas)
     polys_gdf['width'] = polys_gdf['area'].div(np.pi).pow(0.5)
     return polys_gdf
 
+
 def derive_rc(subcatchments: gpd.GeoDataFrame,
               building_footprints: gpd.GeoDataFrame,
               streetcover: gpd.GeoDataFrame) -> gpd.GeoDataFrame:

swmmanywhere/misc/debug_derive_rc.py

@@ -0,0 +1,80 @@
+"""Debugging the derive_rc_alt function.
+
+This is an alternative function for derive_rc in geospatial_utilities/derive_rc
+it is used to double check that they are performing correctly. It may also be
+more computationally efficient.
+"""
+from __future__ import annotations
+
+import geopandas as gpd
+import pandas as pd
+import shapely
+
+
+def derive_rc_alt(subcatchments: gpd.GeoDataFrame,
+                  building_footprints: gpd.GeoDataFrame,
+                  streetcover: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
+    """Derive the Runoff Coefficient (RC) of each subcatchment (alt)."""
+    building_footprints = building_footprints.copy()
+    streetcover = streetcover.copy()
+    subcatchments = subcatchments.copy()
+
+    # Derive impervious without overlaps
+    intersecting_geom_sidx, intersecting_geom_bidx = \
+        building_footprints.sindex.query(streetcover.geometry, predicate='intersects')
+    intersecting_building_geom = \
+        building_footprints.iloc[intersecting_geom_bidx].geometry
+    intersecting_street_geom = streetcover.iloc[intersecting_geom_sidx].geometry
+
+    if intersecting_geom_sidx.shape[0] > 0:
+        unified_geom = \
+            intersecting_building_geom.unary_union.union(
+                intersecting_street_geom.unary_union)
+        unified_geom = gpd.GeoDataFrame(geometry=[unified_geom], 
+                                        crs=building_footprints.crs)
+    else:
+        unified_geom = gpd.GeoDataFrame(geometry=[], crs=building_footprints.crs)
+    building_footprints = building_footprints.drop(
+        building_footprints.index[intersecting_geom_bidx],
+                                                   axis = 0)
+    streetcover = streetcover.drop(streetcover.index[intersecting_geom_sidx],
+                                   axis = 0)
+
+    new_geoms = [g for g in [unified_geom, building_footprints, streetcover] 
+                 if g.shape[0] > 0]
+    # Create the "unified" impervious geometries
+    impervious = gpd.GeoDataFrame(pd.concat(new_geoms), 
+                                    crs=building_footprints.crs)
+    subcat_tree = subcatchments.sindex
+    bf_pidx, sb_pidx = subcat_tree.query(impervious.geometry,
+                                         predicate='intersects')
+    sb_idx = subcatchments.iloc[sb_pidx].index
+
+    # Calculate impervious area and runoff coefficient (rc)
+    subcatchments["impervious_area"] = 0.0
+
+    # Calculate all intersection-impervious geometries
+    intersection_area = shapely.intersection(
+                subcatchments.iloc[sb_pidx].geometry.to_numpy(), 
+                impervious.iloc[bf_pidx].geometry.to_numpy())
+
+    # Indicate which catchment each intersection is part of 
+    intersections = pd.DataFrame([{'sb_idx': ix,
+                                  'impervious_geometry': ia}
+                                  for ix, ia in zip(sb_idx, intersection_area)]
+                                  )
+
+    # Aggregate by catchment
+    areas = (
+        intersections
+        .groupby('sb_idx')
+        .apply(shapely.ops.unary_union)
+        .apply(shapely.area)
+    )
+
+    # Store as impervious area in subcatchments
+    subcatchments["impervious_area"] = 0
+    subcatchments.loc[areas.index, "impervious_area"] = areas
+    subcatchments["rc"] = subcatchments["impervious_area"] / \
+        subcatchments.geometry.area * 100
+    return subcatchments

tests/test_geospatial_utilities.py

@@ -18,13 +18,20 @@
 
 from swmmanywhere import geospatial_utilities as go
 from swmmanywhere import graph_utilities as ge
+from swmmanywhere.misc.debug_derive_rc import derive_rc_alt
 
 
 def load_street_network():
     """Load a street network."""
     G = ge.load_graph(Path(__file__).parent / 'test_data' / 'street_graph.json')
     return G
 
+def almost_equal(a, b, tol=1e-6):
+    """Check if two numbers are almost equal."""
+    if hasattr(a, 'shape'):
+        return ((a - b) < tol).all().all()
+    return abs(a-b) < tol
+
 def test_interp_with_nans():
     """Test the interp_interp_with_nans function."""
     # Define a simple grid and values
@@ -127,11 +134,6 @@ def test_reproject_raster():
         fid.unlink(missing_ok=True)
         new_fid.unlink(missing_ok=True)
 
-
-def almost_equal(a, b, tol=1e-6):
-    """Check if two numbers are almost equal."""
-    return abs(a-b) < tol
-
 def test_get_transformer():
     """Test the get_transformer function."""
     # Test a northern hemisphere point
@@ -278,6 +280,9 @@ def test_derive_rc():
                             (700329, 5709883),
                             (700351, 5709883)])]
 
+    streetcover = [d['geometry'].buffer(5) for u,v,d in G.edges(data=True)]
+    streetcover = gpd.GeoDataFrame(geometry = streetcover, crs = crs)
+
     subs = gpd.GeoDataFrame(data = {'id' : [107733,
                                             1696030874,
                                             6277683849]
@@ -286,16 +291,49 @@ def test_derive_rc():
                                     crs = crs)
     subs['area'] = subs.geometry.area
 
+    # Test no RC
     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
 
-    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
+    # Test alt method
+    subs_rc_alt = derive_rc_alt(subs, buildings, buildings).set_index('id')
+    assert almost_equal(
+        subs_rc[['area','impervious_area','rc']],
+        subs_rc_alt[['area','impervious_area','rc']]
+        )
+
+    # Test some RC
+    subs_rc = go.derive_rc(subs, buildings, streetcover).set_index('id')
+    assert almost_equal(subs_rc.loc[6277683849,'impervious_area'],
+                        1092.452579)
+    assert almost_equal(subs_rc.loc[6277683849,'rc'],
+                        21.770677)
     assert subs_rc.rc.max() <= 100
 
+    # Test alt method
+    subs_rc_alt = derive_rc_alt(subs, buildings, streetcover).set_index('id')
+    assert almost_equal(
+        subs_rc[['area','impervious_area','rc']],
+        subs_rc_alt[['area','impervious_area','rc']]
+        )
+
+    # Test intersecting buildings and streets
+    buildings = buildings.overlay(streetcover.dissolve(), how = 'union')
+    subs_rc = go.derive_rc(subs, buildings, streetcover).set_index('id')
+    assert almost_equal(subs_rc.loc[6277683849,'impervious_area'],
+                        1092.452579)
+    assert almost_equal(subs_rc.loc[6277683849,'rc'],
+                        21.770677)
+    assert subs_rc.rc.max() <= 100
+
+    # Test alt method
+    subs_rc_alt = derive_rc_alt(subs, buildings, streetcover).set_index('id')
+    assert almost_equal(
+        subs_rc[['area','impervious_area','rc']],
+        subs_rc_alt[['area','impervious_area','rc']]
+        )
+
 def test_calculate_angle():
     """Test the calculate_angle function."""
     # Test with points forming a right angle