Merge pull request #98 from ImperialCollegeLondon/51-gridded-metrics

51 gridded metrics

swmmanywhere/metric_utilities.py

@@ -5,6 +5,7 @@
 """
 from collections import defaultdict
 from inspect import signature
+from itertools import product
 from typing import Callable, Optional
 
 import cytoolz.curried as tlz
@@ -14,8 +15,11 @@
 import networkx as nx
 import numpy as np
 import pandas as pd
+import shapely
 from scipy import stats
 
+from swmmanywhere.parameters import MetricEvaluation
+
 
 class MetricRegistry(dict): 
     """Registry object.""" 
@@ -30,7 +34,8 @@ def register(self, func: Callable) -> Callable:
                             "synthetic_subs": gpd.GeoDataFrame,
                             "real_subs": gpd.GeoDataFrame,
                             "synthetic_G": nx.Graph,
-                            "real_G": nx.Graph}
+                            "real_G": nx.Graph,
+                            "metric_evaluation": MetricEvaluation}
 
         sig = signature(func)
         for param, obj in sig.parameters.items():
@@ -61,7 +66,8 @@ def iterate_metrics(synthetic_results: pd.DataFrame,
                     real_results: pd.DataFrame,
                     real_subs: gpd.GeoDataFrame,
                     real_G: nx.Graph,
-                    metric_list: list[str]) -> dict[str, float]:
+                    metric_list: list[str],
+                    metric_evaluation: MetricEvaluation) -> dict[str, float]:
     """Iterate a list of metrics over a graph.
 
     Args:
@@ -72,6 +78,7 @@ def iterate_metrics(synthetic_results: pd.DataFrame,
         real_subs (gpd.GeoDataFrame): The real subcatchments.
         real_G (nx.Graph): The real graph.
         metric_list (list[str]): A list of metrics to iterate.
+        metric_evaluation (MetricEvaluation): The metric evaluation parameters.
 
     Returns:
         dict[str, float]: The results of the metrics.
@@ -87,6 +94,7 @@ def iterate_metrics(synthetic_results: pd.DataFrame,
         "real_results": real_results,
         "real_subs": real_subs,
         "real_G": real_G,
+        "metric_evaluation": metric_evaluation
     }
 
     return {m : metrics[m](**kwargs) for m in metric_list}
@@ -313,18 +321,26 @@ def edge_betweenness_centrality(G: nx.Graph,
             bt_c[n] += v
     return bt_c
 
-def align_by_subcatchment(var,
+def align_by_shape(var,
                           synthetic_results: pd.DataFrame,
                           real_results: pd.DataFrame,
-                          real_subs: gpd.GeoDataFrame,
+                          shapes: gpd.GeoDataFrame,
                           synthetic_G: nx.Graph,
                           real_G: nx.Graph) -> pd.DataFrame:
     """Align by subcatchment.
 
-    Align synthetic and real results by subcatchment and return the results.
+    Align synthetic and real results by shape and return the results.
+
+    Args:
+        var (str): The variable to align.
+        synthetic_results (pd.DataFrame): The synthetic results.
+        real_results (pd.DataFrame): The real results.
+        shapes (gpd.GeoDataFrame): The shapes to align by (e.g., grid or real_subs).
+        synthetic_G (nx.Graph): The synthetic graph.
+        real_G (nx.Graph): The real graph.
     """
-    synthetic_joined = nodes_to_subs(synthetic_G, real_subs)
-    real_joined = nodes_to_subs(real_G, real_subs)
+    synthetic_joined = nodes_to_subs(synthetic_G, shapes)
+    real_joined = nodes_to_subs(real_G, shapes)
 
     # Extract data
     real_results = extract_var(real_results, var)
@@ -347,6 +363,34 @@ def align_by_subcatchment(var,
                             )
     return results
 
+def create_grid(bbox: tuple,
+                scale: float | tuple[float,float]) -> gpd.GeoDataFrame:
+    """Create a grid of polygons.
+    
+    Create a grid of polygons based on the bounding box and scale.
+
+    Args:
+        bbox (tuple): The bounding box coordinates in the format (minx, miny, 
+            maxx, maxy).
+        scale (float | tuple): The scale of the grid. If a tuple, the scale is 
+            (dx, dy). Otherwise, the scale is dx = dy = scale.
+    
+    Returns:
+        gpd.GeoDataFrame: A geodataframe of the grid.
+    """
+    minx, miny, maxx, maxy = bbox
+
+    if isinstance(scale, tuple):
+        dx, dy = scale
+    else: 
+        dx = dy = scale
+    xmins = np.arange(minx, maxx, dx)
+    ymins = np.arange(minx, maxy, dy)
+    grid = [{'geometry' : shapely.box(x, y, x + dx, y + dy),
+             'sub_id' : i} for i, (x, y) in enumerate(product(xmins, ymins))]
+
+    return gpd.GeoDataFrame(grid)
+
 @metrics.register
 def nc_deltacon0(synthetic_G: nx.Graph,
                   real_G: nx.Graph,
@@ -647,10 +691,38 @@ def subcatchment_nse_flooding(synthetic_G: nx.Graph,
     flooding over time for each subcatchment. The metric produced is the median
     NSE across all subcatchments.
     """
-    results = align_by_subcatchment('flooding',
+    results = align_by_shape('flooding',
+                                    synthetic_results = synthetic_results,
+                                    real_results = real_results,
+                                    shapes = real_subs,
+                                    synthetic_G = synthetic_G,
+                                    real_G = real_G)
+
+    return median_nse_by_group(results, 'sub_id')
+
+@metrics.register
+def grid_nse_flooding(synthetic_G: nx.Graph,
+                            real_G: nx.Graph,
+                            synthetic_results: pd.DataFrame,
+                            real_results: pd.DataFrame,
+                            real_subs: gpd.GeoDataFrame,
+                            metric_evaluation: MetricEvaluation,
+                            **kwargs) -> float:
+    """Grid NSE flooding.
+    
+    Classify synthetic nodes to a grid and calculate the NSE of
+    flooding over time for each grid cell. The metric produced is the median
+    NSE across all grid cells.
+    """
+    scale = metric_evaluation.grid_scale
+    grid = create_grid(real_subs.total_bounds,
+                       scale)
+    grid.crs = real_subs.crs
+
+    results = align_by_shape('flooding',
                                     synthetic_results = synthetic_results,
                                     real_results = real_results,
-                                    real_subs = real_subs,
+                                    shapes = grid,
                                     synthetic_G = synthetic_G,
                                     real_G = real_G)
 

swmmanywhere/parameters.py

@@ -12,7 +12,8 @@ def get_full_parameters():
         "subcatchment_derivation": SubcatchmentDerivation(),
         "outlet_derivation": OutletDerivation(),
         "topology_derivation": TopologyDerivation(),
-        "hydraulic_design": HydraulicDesign()
+        "hydraulic_design": HydraulicDesign(),
+        "metric_evaluation": MetricEvaluation()
     }
 
 class SubcatchmentDerivation(BaseModel):
@@ -168,12 +169,18 @@ class HydraulicDesign(BaseModel):
 		description = "Depth of design storm in pipe by pipe method",
 		unit = "m")
 
-class FilePaths:
-    """Parameters for file path lookup.
+class MetricEvaluation(BaseModel):
+    """Parameters for metric evaluation."""
+    grid_scale: float = Field(default = 100,
+                        le = 10,
+                        ge = 5000,
+                        unit = "m",
+                        description = "Scale of the grid for metric evaluation")
+
 
-    TODO: this doesn't validate file paths to allow for un-initialised data
-    (e.g., subcatchments are created by a graph and so cannot be validated).
-    """
+
+class FilePaths:
+    """Parameters for file path lookup."""
 
     def __init__(self, 
                  base_dir: Path, 

swmmanywhere/swmmanywhere.py

@@ -102,7 +102,8 @@ def swmmanywhere(config: dict):
                               real_results,
                               gpd.read_file(config['real']['subcatchments']),
                               load_graph(config['real']['graph']),
-                              config['metric_list'])
+                              config['metric_list'],
+                              parameters['metric_evaluation'])
 
     return addresses.inp, metrics
 

tests/test_metric_utilities.py

@@ -8,6 +8,7 @@
 
 from swmmanywhere import metric_utilities as mu
 from swmmanywhere.graph_utilities import load_graph
+from swmmanywhere.parameters import MetricEvaluation
 
 
 def assert_close(a: float, b: float, rtol: float = 1e-3) -> None:
@@ -346,7 +347,8 @@ def test_netcomp_iterate():
                                  real_G = G,
                                  real_subs = None,
                                  real_results = None,
-                                 metric_list = netcomp_results.keys())
+                                 metric_list = netcomp_results.keys(),
+                                 metric_evaluation = MetricEvaluation())
     for metric, val in metrics.items():
         assert metric in netcomp_results
         assert np.isclose(val, 0)
@@ -358,7 +360,8 @@ def test_netcomp_iterate():
                                  real_G = G,
                                  real_subs = None,
                                  real_results = None,
-                                 metric_list = netcomp_results.keys())
+                                 metric_list = netcomp_results.keys(),
+                                 metric_evaluation = MetricEvaluation())
     for metric, val in metrics.items():
         assert metric in netcomp_results
         assert np.isclose(val, netcomp_results[metric])
@@ -446,3 +449,18 @@ def test_subcatchment_nse_flooding():
                                     real_results = results,
                                     real_subs = subs)
     assert val == 1.0
+
+    # Test gridded
+    val = mu.metrics.grid_nse_flooding(synthetic_G = G_,
+                                    synthetic_results = results_,
+                                    real_G = G,
+                                    real_results = results,
+                                    real_subs = subs,
+                                    metric_evaluation = MetricEvaluation())
+    assert val == 1.0
+
+def test_create_grid():
+    """Test the create_grid function."""
+    grid = mu.create_grid((0,0,1,1), 1/3 - 0.001)
+    assert grid.shape[0] == 16
+    assert set(grid.columns) == {'sub_id','geometry'}