diff --git a/dev-requirements.txt b/dev-requirements.txt index ee78d0e0..88ec4101 100644 --- a/dev-requirements.txt +++ b/dev-requirements.txt @@ -12,6 +12,8 @@ affine==2.4.0 # via # pysheds # rasterio +annotated-types==0.6.0 + # via pydantic attrs==23.2.0 # via # fiona @@ -197,6 +199,10 @@ pre-commit==3.6.0 # via swmmanywhere (pyproject.toml) pyarrow==14.0.2 # via swmmanywhere (pyproject.toml) +pydantic==2.5.3 + # via swmmanywhere (pyproject.toml) +pydantic-core==2.14.6 + # via pydantic pyparsing==3.1.1 # via # matplotlib @@ -282,7 +288,10 @@ tqdm==4.66.1 # cdsapi # swmmanywhere (pyproject.toml) typing-extensions==4.9.0 - # via mypy + # via + # mypy + # pydantic + # pydantic-core tzdata==2023.4 # via pandas urllib3==2.1.0 diff --git a/pyproject.toml b/pyproject.toml index 3024c537..dedfb435 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,6 +27,7 @@ dependencies = [ # TODO definitely don't need all of these "osmnx", "pandas", "pyarrow", + "pydantic", "pysheds", "pyswmm", "PyYAML", diff --git a/requirements.txt b/requirements.txt index 86670e7f..0770bb21 100644 --- a/requirements.txt +++ b/requirements.txt @@ -12,6 +12,8 @@ affine==2.4.0 # via # pysheds # rasterio +annotated-types==0.6.0 + # via pydantic attrs==23.2.0 # via # fiona @@ -157,6 +159,10 @@ pillow==10.2.0 # scikit-image pyarrow==14.0.2 # via swmmanywhere (pyproject.toml) +pydantic==2.5.3 + # via swmmanywhere (pyproject.toml) +pydantic-core==2.14.6 + # via pydantic pyparsing==3.1.1 # via # matplotlib @@ -223,6 +229,10 @@ tqdm==4.66.1 # via # cdsapi # swmmanywhere (pyproject.toml) +typing-extensions==4.9.0 + # via + # pydantic + # pydantic-core tzdata==2023.4 # via pandas urllib3==2.1.0 diff --git a/swmmanywhere/geospatial_utilities.py b/swmmanywhere/geospatial_utilities.py index 6749f4c7..c63cecbe 100644 --- a/swmmanywhere/geospatial_utilities.py +++ b/swmmanywhere/geospatial_utilities.py @@ -6,28 +6,34 @@ @author: Barnaby Dobson """ +import math +from copy import deepcopy from functools import lru_cache -from typing import Optional +from pathlib import Path +from typing import List, Optional import geopandas as gpd import networkx as nx import numpy as np import pandas as pd import pyproj +import pysheds import rasterio as rst import rioxarray +from pysheds import grid as pgrid from rasterio import features from scipy.interpolate import RegularGridInterpolator from shapely import geometry as sgeom +from shapely import ops as sops from shapely.strtree import STRtree +from tqdm import tqdm TransformerFromCRS = lru_cache(pyproj.transformer.Transformer.from_crs) - def get_utm_epsg(x: float, y: float, crs: str | int | pyproj.CRS = 'EPSG:4326', - datum_name: str = "WGS 84"): + datum_name: str = "WGS 84") -> str: """Get the UTM CRS code for a given coordinate. Note, this function is taken from GeoPandas and modified to use @@ -112,13 +118,13 @@ def interp_with_nans(xy: tuple[float,float], def interpolate_points_on_raster(x: list[float], y: list[float], - elevation_fid: str) -> list[float ]: + elevation_fid: Path) -> list[float ]: """Interpolate points on a raster. Args: x (list): X coordinates. y (list): Y coordinates. - elevation_fid (str): Filepath to elevation raster. + elevation_fid (Path): Filepath to elevation raster. Returns: elevation (float): Elevation at point. @@ -129,16 +135,16 @@ def interpolate_points_on_raster(x: list[float], data[data == src.nodata] = None # Get the raster's coordinates - x = np.linspace(src.bounds.left, src.bounds.right, src.width) - y = np.linspace(src.bounds.bottom, src.bounds.top, src.height) + x_grid = np.linspace(src.bounds.left, src.bounds.right, src.width) + y_grid = np.linspace(src.bounds.bottom, src.bounds.top, src.height) # Define grid - xx, yy = np.meshgrid(x, y) + xx, yy = np.meshgrid(x_grid, y_grid) grid = np.vstack([xx.ravel(), yy.ravel()]).T values = data.ravel() # Define interpolator - interp = RegularGridInterpolator((y,x), + interp = RegularGridInterpolator((y_grid,x_grid), np.flipud(data), method='linear', bounds_error=False, @@ -147,14 +153,14 @@ def interpolate_points_on_raster(x: list[float], return [interp_with_nans((y_, x_), interp, grid, values) for x_, y_ in zip(x,y)] def reproject_raster(target_crs: str, - fid: str, - new_fid: Optional[str] = None): + fid: Path, + new_fid: Optional[Path] = None) -> None: """Reproject a raster to a new CRS. Args: target_crs (str): Target CRS in EPSG format (e.g., EPSG:32630). - fid (str): Filepath to the raster to reproject. - new_fid (str, optional): Filepath to save the reprojected raster. + fid (Path): Filepath to the raster to reproject. + new_fid (Path, optional): Filepath to save the reprojected raster. Defaults to None, which will just use fid with '_reprojected'. """ # Open the raster @@ -165,7 +171,7 @@ def reproject_raster(target_crs: str, # Define the output filepath if new_fid is None: - new_fid = fid.replace('.tif','_reprojected.tif') + new_fid = Path(str(fid.with_suffix('')) + '_reprojected.tif') # Save the reprojected raster reprojected.rio.to_raster(new_fid) @@ -295,15 +301,15 @@ def nearest_node_buffer(points1: dict[str, sgeom.Point], def burn_shape_in_raster(geoms: list[sgeom.LineString], depth: float, - raster_fid: str, - new_raster_fid: str): + raster_fid: Path, + new_raster_fid: Path) -> None: """Burn a depth into a raster along a list of shapely geometries. Args: geoms (list): List of Shapely geometries. depth (float): Depth to carve. - raster_fid (str): Filepath to input raster. - new_raster_fid (str): Filepath to save the carved raster. + raster_fid (Path): Filepath to input raster. + new_raster_fid (Path): Filepath to save the carved raster. """ with rst.open(raster_fid) as src: # read data @@ -332,4 +338,332 @@ def burn_shape_in_raster(geoms: list[sgeom.LineString], crs=src.crs, transform=src.transform, nodata = src.nodata) as dest: - dest.write(data, 1) \ No newline at end of file + dest.write(data, 1) + +def condition_dem(grid: pysheds.sgrid.sGrid, + dem: pysheds.sview.Raster) -> pysheds.sview.Raster: + """Condition a DEM with pysheds. + + Args: + grid (pysheds.sgrid.sGrid): The grid object. + dem (pysheds.sview.Raster): The input DEM. + + Returns: + pysheds.sview.Raster: The conditioned DEM. + """ + # Fill pits, depressions, and resolve flats in the DEM + pit_filled_dem = grid.fill_pits(dem) + flooded_dem = grid.fill_depressions(pit_filled_dem) + inflated_dem = grid.resolve_flats(flooded_dem) + + return inflated_dem + +def compute_flow_directions(grid: pysheds.sgrid.sGrid, + inflated_dem: pysheds.sview.Raster) \ + -> tuple[pysheds.sview.Raster, tuple]: + """Compute flow directions. + + Args: + grid (pysheds.sgrid.sGrid): The grid object. + inflated_dem (pysheds.sview.Raster): The input DEM. + + Returns: + pysheds.sview.Raster: Flow directions. + tuple: Direction mapping. + """ + dirmap = (64, 128, 1, 2, 4, 8, 16, 32) + flow_dir = grid.flowdir(inflated_dem, dirmap=dirmap) + return flow_dir, dirmap + +def calculate_flow_accumulation(grid: pysheds.sgrid.sGrid, + flow_dir: pysheds.sview.Raster, + dirmap: tuple) -> pysheds.sview.Raster: + """Calculate flow accumulation. + + Args: + grid (pysheds.sgrid.sGrid): The grid object. + flow_dir (pysheds.sview.Raster): Flow directions. + dirmap (tuple): Direction mapping. + + Returns: + pysheds.sview.Raster: Flow accumulations. + """ + flow_acc = grid.accumulation(flow_dir, dirmap=dirmap) + return flow_acc + +def delineate_catchment(grid: pysheds.sgrid.sGrid, + flow_acc: pysheds.sview.Raster, + flow_dir: pysheds.sview.Raster, + dirmap: tuple, + G: nx.Graph) -> gpd.GeoDataFrame: + """Delineate catchments. + + Args: + grid (pysheds.sgrid.Grid): The grid object. + flow_acc (pysheds.sview.Raster): Flow accumulations. + flow_dir (pysheds.sview.Raster): Flow directions. + dirmap (tuple): Direction mapping. + G (nx.Graph): The input graph with nodes containing 'x' and 'y'. + + Returns: + gpd.GeoDataFrame: A GeoDataFrame containing polygons with columns: + 'geometry', 'area', and 'id'. Sorted by area in descending order. + """ + polys = [] + # Iterate over the nodes in the graph + for id, data in tqdm(G.nodes(data=True), total=len(G.nodes)): + # Snap the node to the nearest grid cell + x, y = data['x'], data['y'] + grid_ = deepcopy(grid) + x_snap, y_snap = grid_.snap_to_mask(flow_acc > 5, (x, y)) + + # Delineate the catchment + catch = grid_.catchment(x=x_snap, + y=y_snap, + fdir=flow_dir, + dirmap=dirmap, + xytype='coordinate', + algorithm = 'recursive' + ) + # n.b. recursive algorithm is not recommended, but crashes with a seg + # fault occasionally otherwise. + + grid_.clip_to(catch) + + # Polygonize the catchment + shapes = grid_.polygonize() + catchment_polygon = sops.unary_union([sgeom.shape(shape) for shape, + value in shapes]) + + # Add the catchment to the list + polys.append({'id': id, + 'geometry': catchment_polygon, + 'area': catchment_polygon.area}) + polys = sorted(polys, key=lambda d: d['area'], reverse=True) + return gpd.GeoDataFrame(polys, crs = grid.crs) + +def remove_intersections(polys: gpd.GeoDataFrame) -> gpd.GeoDataFrame: + """Remove intersections from polygons. + + Args: + polys (gpd.GeoDataFrame): A geodataframe with columns id and geometry. + + Returns: + gpd.GeoDataFrame: A geodataframe with columns id and geometry. + """ + result_polygons = polys.copy() + for i in tqdm(range(len(polys))): + for j in range(i + 1, len(polys)): + if polys.iloc[i]['geometry'].intersects(polys.iloc[j]['geometry']): + polyi = result_polygons.iloc[i]['geometry'] + polyj = polys.iloc[j]['geometry'] + result_polygons.at[polys.index[i], + 'geometry'] = polyi.difference(polyj) + return result_polygons + +def remove_zero_area_subareas(mp: sgeom.MultiPolygon, + removed_subareas: List[sgeom.Polygon]) \ + -> sgeom.MultiPolygon: + """Remove subareas with zero area from a multipolygon. + + Args: + mp (sgeom.MultiPolygon): A multipolygon. + removed_subareas (List[sgeom.Polygon]): A list of removed subareas. + + Returns: + sgeom.MultiPolygon: A multipolygon with zero area subareas removed. + """ + if not hasattr(mp, 'geoms'): + return mp + + largest = max(mp.geoms, key=lambda x: x.area) + removed = [subarea for subarea in mp.geoms if subarea != largest] + removed_subareas.extend(removed) + return largest + +def attach_unconnected_subareas(polys_gdf: gpd.GeoDataFrame, + unconnected_subareas: List[sgeom.Polygon]) \ + -> gpd.GeoDataFrame: + """Attach unconnected subareas to the nearest polygon. + + Args: + polys_gdf (gpd.GeoDataFrame): A GeoDataFrame containing polygons with + columns: 'geometry', 'area', and 'id'. + unconnected_subareas (List[sgeom.Polygon]): A list of subareas that are + not attached to others. + + Returns: + gpd.GeoDataFrame: A GeoDataFrame containing polygons. + """ + tree = STRtree(polys_gdf.geometry) + for subarea in unconnected_subareas: + nearest_poly = tree.nearest(subarea) + ind = polys_gdf.index[nearest_poly] + new_poly = polys_gdf.loc[ind, 'geometry'].union(subarea) + if hasattr(new_poly, 'geoms'): + new_poly = max(new_poly.geoms, key=lambda x: x.area) + polys_gdf.at[ind, 'geometry'] = new_poly + return polys_gdf + +def calculate_slope(polys_gdf: gpd.GeoDataFrame, + grid: pysheds.sgrid.sGrid, + cell_slopes: np.ndarray) -> gpd.GeoDataFrame: + """Calculate the average slope of each polygon. + + Args: + polys_gdf (gpd.GeoDataFrame): A GeoDataFrame containing polygons with + columns: 'geometry', 'area', and 'id'. + grid (pysheds.sgrid.sGrid): The grid object. + cell_slopes (np.ndarray): The slopes of each cell in the grid. + + Returns: + gpd.GeoDataFrame: A GeoDataFrame containing polygons with an added + 'slope' column. + """ + polys_gdf['slope'] = None + for idx, row in polys_gdf.iterrows(): + mask = features.geometry_mask([row.geometry], + grid.shape, + grid.affine, + invert=True) + average_slope = cell_slopes[mask].mean().mean() + polys_gdf.loc[idx, 'slope'] = max(float(average_slope), 0) + return polys_gdf + +def derive_subcatchments(G: nx.Graph, fid: Path) -> gpd.GeoDataFrame: + """Derive subcatchments from the nodes on a graph and a DEM. + + Args: + G (nx.Graph): The input graph with nodes containing 'x' and 'y'. + fid (Path): Filepath to the DEM. + + Returns: + gpd.GeoDataFrame: A GeoDataFrame containing polygons with columns: + 'geometry', 'area', 'id', 'width', and 'slope'. + """ + # Initialise pysheds grids + grid = pgrid.Grid.from_raster(str(fid)) + dem = grid.read_raster(str(fid)) + + # Condition the DEM + inflated_dem = condition_dem(grid, dem) + + # Compute flow directions + flow_dir, dirmap = compute_flow_directions(grid, inflated_dem) + + # Calculate slopes + cell_slopes = grid.cell_slopes(dem, flow_dir) + + # Calculate flow accumulations + flow_acc = calculate_flow_accumulation(grid, flow_dir, dirmap) + + # Delineate catchments + polys = delineate_catchment(grid, flow_acc, flow_dir, dirmap, G) + + # Remove intersections + result_polygons = remove_intersections(polys) + + # Convert to GeoDataFrame + polys_gdf = result_polygons.dropna(subset=['geometry']) + polys_gdf = polys_gdf[~polys_gdf['geometry'].is_empty] + + # Remove zero area subareas and attach to nearest polygon + removed_subareas: List[sgeom.Polygon] = list() # needed for mypy + def remove_(mp): return remove_zero_area_subareas(mp, removed_subareas) + polys_gdf['geometry'] = polys_gdf['geometry'].apply(remove_) + polys_gdf = attach_unconnected_subareas(polys_gdf, removed_subareas) + + # Calculate area and slope + polys_gdf['area'] = polys_gdf.geometry.area + polys_gdf = calculate_slope(polys_gdf, grid, cell_slopes) + + # Calculate width + polys_gdf['width'] = polys_gdf['area'].div(np.pi).pow(0.5) + return polys_gdf + +def derive_rc(polys_gdf: gpd.GeoDataFrame, + G: nx.Graph, + building_footprints: gpd.GeoDataFrame) -> gpd.GeoDataFrame: + """Derive the RC of each subcatchment. + + 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. + + Returns: + gpd.GeoDataFrame: A GeoDataFrame containing polygons with columns: + 'geometry', 'area', 'id', 'impervious_area', and 'rc'. + """ + polys_gdf = polys_gdf.copy() + + ## Format as swmm type catchments + + # TODO think harder about lane widths (am I double counting here?) + lines = [] + for u, v, x in G.edges(data=True): + lines.append({'geometry' : x['geometry'].buffer(x['width'], + cap_style = 2, + join_style=2), + 'id' : x['id']}) + 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']], + building_footprints[['geometry']], + how='union') + result = gpd.overlay(polys_gdf, result) + + dissolved_result = result.dissolve(by='id').reset_index() + dissolved_result['impervious_area'] = dissolved_result.geometry.area + polys_gdf = pd.merge(polys_gdf, + dissolved_result[['id','impervious_area']], + on = 'id', + how='left').fillna(0) + polys_gdf['rc'] = polys_gdf['impervious_area'] / polys_gdf['area'] * 100 + return polys_gdf + +def calculate_angle(point1: tuple[float,float], + point2: tuple[float,float], + point3: tuple[float,float]) -> float: + """Calculate the angle between three points. + + Calculate the angle between the vectors formed by (point1, + point2) and (point2, point3) + + Args: + point1 (tuple): The first point (x,y). + point2 (tuple): The second point (x,y). + point3 (tuple): The third point (x,y). + + Returns: + float: The angle between the three points in degrees. + """ + vector1 = (point1[0] - point2[0], point1[1] - point2[1]) + vector2 = (point3[0] - point2[0], point3[1] - point2[1]) + + dot_product = vector1[0] * vector2[0] + vector1[1] * vector2[1] + magnitude1 = math.sqrt(vector1[0]**2 + vector1[1]**2) + magnitude2 = math.sqrt(vector2[0]**2 + vector2[1]**2) + + if magnitude1 * magnitude2 == 0: + # Avoid division by zero + return float('inf') + + cosine_angle = dot_product / (magnitude1 * magnitude2) + + # Ensure the cosine value is within the valid range [-1, 1] + cosine_angle = min(max(cosine_angle, -1), 1) + + # Calculate the angle in radians + angle_radians = math.acos(cosine_angle) + + # Convert angle to degrees + angle_degrees = math.degrees(angle_radians) + + return angle_degrees \ No newline at end of file diff --git a/swmmanywhere/graph_utilities.py b/swmmanywhere/graph_utilities.py new file mode 100644 index 00000000..ad81ce75 --- /dev/null +++ b/swmmanywhere/graph_utilities.py @@ -0,0 +1,994 @@ +# -*- coding: utf-8 -*- +"""Created on 2024-01-26. + +@author: Barney +""" +import json +import tempfile +from abc import ABC, abstractmethod +from collections import defaultdict +from heapq import heappop, heappush +from itertools import product +from pathlib import Path +from typing import Any, Callable, Dict, Hashable, List, Optional + +import geopandas as gpd +import networkx as nx +import numpy as np +import osmnx as ox +import pandas as pd +import shapely +from tqdm import tqdm + +from swmmanywhere import geospatial_utilities as go +from swmmanywhere import parameters + + +def load_graph(fid: Path) -> nx.Graph: + """Load a graph from a file saved with save_graph. + + Args: + fid (Path): The path to the file + + Returns: + G (nx.Graph): A graph + """ + json_data = json.loads(fid.read_text()) + + G = nx.node_link_graph(json_data,directed=True) + for u, v, data in G.edges(data=True): + if 'geometry' in data: + geometry_coords = data['geometry'] + line_string = shapely.LineString(shapely.wkt.loads(geometry_coords)) + data['geometry'] = line_string + return G +def _serialize_line_string(obj): + if isinstance(obj, shapely.LineString): + return obj.wkt + else: + return obj +def save_graph(G: nx.Graph, + fid: Path) -> None: + """Save a graph to a file. + + Args: + G (nx.Graph): A graph + fid (Path): The path to the file + """ + json_data = nx.node_link_data(G) + + with open(fid, 'w') as file: + json.dump(json_data, + file, + default = _serialize_line_string) + + +class BaseGraphFunction(ABC): + """Base class for graph functions. + + On a SWMManywhere project the intention is to iterate over a number of + graph functions. Each graph function may require certain attributes to + be present in the graph. Each graph function may add attributes to the + graph. This class provides a framework for graph functions to check + their requirements and additions a-priori when the list is provided. + """ + + required_edge_attributes: List[str] = list() + adds_edge_attributes: List[str] = list() + required_node_attributes: List[str] = list() + adds_node_attributes: List[str] = list() + def __init_subclass__(cls, + required_edge_attributes: Optional[List[str]] = None, + adds_edge_attributes: Optional[List[str]] = None, + required_node_attributes : Optional[List[str]] = None, + adds_node_attributes : Optional[List[str]] = None + ): + """Set the required and added attributes for the subclass.""" + cls.required_edge_attributes = required_edge_attributes if \ + required_edge_attributes else [] + cls.adds_edge_attributes = adds_edge_attributes if \ + adds_edge_attributes else [] + cls.required_node_attributes = required_node_attributes if \ + required_node_attributes else [] + cls.adds_node_attributes = adds_node_attributes if \ + adds_node_attributes else [] + + @abstractmethod + def __call__(self, + G: nx.Graph, + *args, + **kwargs) -> nx.Graph: + """Run the graph function.""" + return G + + def validate_requirements(self, + edge_attributes: set, + node_attributes: set) -> None: + """Validate that the graph has the required attributes.""" + for attribute in self.required_edge_attributes: + assert attribute in edge_attributes, f"{attribute} not in edge attributes" + + for attribute in self.required_node_attributes: + assert attribute in node_attributes, f"{attribute} not in node attributes" + + + def add_graphfcn(self, + edge_attributes: set, + node_attributes: set) -> tuple[set, set]: + """Add the attributes that the graph function adds.""" + self.validate_requirements(edge_attributes, node_attributes) + edge_attributes = edge_attributes.union(self.adds_edge_attributes) + node_attributes = node_attributes.union(self.adds_node_attributes) + return edge_attributes, node_attributes + +class GraphFunctionRegistry(dict): + """Registry object.""" + + def register(self, cls): + """Register a graph function.""" + if cls.__name__ in self: + raise ValueError(f"{cls.__name__} already in the graph functions registry!") + + self[cls.__name__] = cls() + return cls + + def __getattr__(self, name): + """Get a graph function from the graphfcn dict.""" + try: + return self[name] + except KeyError: + raise AttributeError(f"{name} NOT in the graph functions registry!") + + +graphfcns = GraphFunctionRegistry() + +def register_graphfcn(cls) -> Callable: + """Register a graph function. + + Args: + cls (Callable): A class that inherits from BaseGraphFunction + + Returns: + cls (Callable): The same class + """ + graphfcns.register(cls) + return cls + +def get_osmid_id(data: dict) -> Hashable: + """Get the ID of an edge.""" + id_ = data.get('osmid', data.get('id')) + if isinstance(id_, list): + id_ = id_[0] + return id_ + +@register_graphfcn +class assign_id(BaseGraphFunction, + required_edge_attributes = ['osmid'], + adds_edge_attributes = ['id'] + ): + """assign_id class.""" + + def __call__(self, + G: nx.Graph, + **kwargs) -> nx.Graph: + """Assign an ID to each edge. + + This function takes a graph and assigns an ID to each edge. The ID is + assigned to the 'id' attribute of each edge. Needed because some edges + have 'osmid', some have a list of 'osmid', others have 'id'. + + Args: + G (nx.Graph): A graph + **kwargs: Additional keyword arguments are ignored. + + Returns: + G (nx.Graph): The same graph with an ID assigned to each edge + """ + for u, v, data in G.edges(data=True): + data['id'] = get_osmid_id(data) + return G + +@register_graphfcn +class format_osmnx_lanes(BaseGraphFunction, + required_edge_attributes = ['lanes'], + adds_edge_attributes = ['width']): + """format_osmnx_lanes 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: + """Format the lanes attribute of each edge and calculates width. + + Args: + G (nx.Graph): A graph + subcatchment_derivation (parameters.SubcatchmentDerivation): A + SubcatchmentDerivation parameter object + **kwargs: Additional keyword arguments are ignored. + + Returns: + G (nx.Graph): A graph + """ + G = G.copy() + 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 + return G + +@register_graphfcn +class double_directed(BaseGraphFunction, + required_edge_attributes = ['id']): + """double_directed class.""" + + def __call__(self, G: nx.Graph, **kwargs) -> nx.Graph: + """Create a 'double directed graph'. + + This function duplicates a graph and adds reverse edges to the new graph. + These new edges share the same data as the 'forward' edges but have a new + 'id'. An undirected graph is not suitable because it removes data of one of + the edges if there are edges in both directions between two nodes + (necessary to preserve, e.g., consistent 'width'). + + Args: + G (nx.Graph): A graph + **kwargs: Additional keyword arguments are ignored. + + Returns: + G (nx.Graph): A graph + """ + #TODO the geometry is left as is currently - should be reveresed, however + # in original osmnx geometry there are some incorrectly directed ones + # someone with more patience might check start and end Points to check + # which direction the line should be going in... + G_new = G.copy() + for u, v, data in G.edges(data=True): + if (v, u) not in G.edges: + reverse_data = data.copy() + reverse_data['id'] = f"{data['id']}.reversed" + G_new.add_edge(v, u, **reverse_data) + return G_new + +@register_graphfcn +class split_long_edges(BaseGraphFunction, + required_edge_attributes = ['id', 'geometry', 'length']): + """split_long_edges class.""" + + def __call__(self, + G: nx.Graph, + subcatchment_derivation: parameters.SubcatchmentDerivation, + **kwargs) -> nx.Graph: + """Split long edges into shorter edges. + + This function splits long edges into shorter edges. The edges are split + into segments of length 'max_street_length'. The first and last segment + are connected to the original nodes. Intermediate segments are connected + to newly created nodes. The 'geometry' of the original edge must be + a LineString. + + Args: + G (nx.Graph): A graph + subcatchment_derivation (parameters.SubcatchmentDerivation): A + SubcatchmentDerivation parameter object + **kwargs: Additional keyword arguments are ignored. + + Returns: + graph (nx.Graph): A graph + """ + #TODO refactor obviously + max_length = subcatchment_derivation.max_street_length + graph = G.copy() + edges_to_remove = [] + edges_to_add = [] + nodes_to_add = [] + maxlabel = max(graph.nodes) + 1 + ll = 0 + + def create_new_edge_data(line, data, id_): + new_line = shapely.LineString(line) + new_data = data.copy() + new_data['id'] = id_ + new_data['length'] = new_line.length + new_data['geometry'] = shapely.LineString([(x[0], x[1]) + for x in new_line.coords]) + return new_data + + for u, v, data in graph.edges(data=True): + line = data['geometry'] + length = data['length'] + if ((u, v) not in edges_to_remove) & ((v, u) not in edges_to_remove): + if length > max_length: + new_points = [shapely.Point(x) + for x in ox.utils_geo.interpolate_points(line, + max_length)] + if len(new_points) > 2: + for ix, (start, end) in enumerate(zip(new_points[:-1], + new_points[1:])): + new_data = create_new_edge_data([start, + end], + data, + f"{data['id']}.{ix}") + if (v,u) in graph.edges: + # Create reversed data + data_r = graph.get_edge_data(v, u).copy()[0] + id_ = f"{data_r['id']}.{ix}" + new_data_r = create_new_edge_data([end, start], + data_r.copy(), + id_) + if ix == 0: + # Create start to first intermediate + edges_to_add.append((u, maxlabel + ll, new_data.copy())) + nodes_to_add.append((maxlabel + ll, + {'x': + new_data['geometry'].coords[-1][0], + 'y': + new_data['geometry'].coords[-1][1]})) + + if (v, u) in graph.edges: + # Create first intermediate to start + edges_to_add.append((maxlabel + ll, + u, + new_data_r.copy())) + + ll += 1 + elif ix == len(new_points) - 2: + # Create last intermediate to end + edges_to_add.append((maxlabel + ll - 1, + v, + new_data.copy())) + if (v, u) in graph.edges: + # Create end to last intermediate + edges_to_add.append((v, + maxlabel + ll - 1, + new_data_r.copy())) + else: + nodes_to_add.append((maxlabel + ll, + {'x': + new_data['geometry'].coords[-1][0], + 'y': + new_data['geometry'].coords[-1][1]})) + # Create N-1 intermediate to N intermediate + edges_to_add.append((maxlabel + ll - 1, + maxlabel + ll, + new_data.copy())) + if (v, u) in graph.edges: + # Create N intermediate to N-1 intermediate + edges_to_add.append((maxlabel + ll, + maxlabel + ll - 1, + new_data_r.copy())) + ll += 1 + edges_to_remove.append((u, v)) + if (v, u) in graph.edges: + edges_to_remove.append((v, u)) + + for u, v in edges_to_remove: + if (u, v) in graph.edges: + graph.remove_edge(u, v) + + for node in nodes_to_add: + graph.add_node(node[0], **node[1]) + + for edge in edges_to_add: + graph.add_edge(edge[0], edge[1], **edge[2]) + + return graph + +@register_graphfcn +class calculate_contributing_area(BaseGraphFunction, + required_edge_attributes = ['id', 'geometry', 'width'], + adds_edge_attributes = ['contributing_area'], + adds_node_attributes = ['contributing_area']): + """calculate_contributing_area class.""" + + def __call__(self, G: nx.Graph, + subcatchment_derivation: parameters.SubcatchmentDerivation, + addresses: parameters.FilePaths, + **kwargs) -> nx.Graph: + """Calculate the contributing area for each edge. + + 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. + + Also writes the file 'subcatchments.geojson' to addresses.subcatchments. + + Args: + G (nx.Graph): A graph + subcatchment_derivation (parameters.SubcatchmentDerivation): A + SubcatchmentDerivation parameter object + addresses (parameters.FilePaths): An FilePaths parameter object + **kwargs: Additional keyword arguments are ignored. + + Returns: + G (nx.Graph): A graph + """ + G = G.copy() + + # Carve + # TODO I guess we don't need to keep this 'carved' file.. + # maybe could add verbose/debug option to keep it + with tempfile.TemporaryDirectory() as temp_dir: + temp_fid = Path(temp_dir) / "carved.tif" + go.burn_shape_in_raster([d['geometry'] for u,v,d in G.edges(data=True)], + subcatchment_derivation.carve_depth, + addresses.elevation, + temp_fid) + + # Derive + subs_gdf = go.derive_subcatchments(G,temp_fid) + + # RC + buildings = gpd.read_file(addresses.building) + subs_rc = go.derive_rc(subs_gdf, G, buildings) + + # Write subs + # TODO - could just attach subs to nodes where each node has a list of subs + subs_rc.to_file(addresses.subcatchments, driver='GeoJSON') + + # Assign contributing area + imperv_lookup = subs_rc.set_index('id').impervious_area.to_dict() + + # Set node attributes + nx.set_node_attributes(G, 0.0, 'contributing_area') + nx.set_node_attributes(G, imperv_lookup, 'contributing_area') + + # Prepare edge attributes + edge_attributes = {edge: G.nodes[edge[0]]['contributing_area'] + for edge in G.edges} + + # Set edge attributes + nx.set_edge_attributes(G, edge_attributes, 'contributing_area') + return G + +@register_graphfcn +class set_elevation(BaseGraphFunction, + required_node_attributes = ['x', 'y'], + adds_node_attributes = ['elevation']): + """set_elevation class.""" + + def __call__(self, G: nx.Graph, + addresses: parameters.FilePaths, + **kwargs) -> nx.Graph: + """Set the elevation for each node. + + This function sets the elevation for each node. The elevation is + calculated from the elevation data. + + Args: + G (nx.Graph): A graph + addresses (parameters.FilePaths): An FilePaths parameter object + **kwargs: Additional keyword arguments are ignored. + + Returns: + G (nx.Graph): A graph + """ + G = G.copy() + x = [d['x'] for x, d in G.nodes(data=True)] + y = [d['y'] for x, d in G.nodes(data=True)] + elevations = go.interpolate_points_on_raster(x, + y, + addresses.elevation) + elevations_dict = {id_: elev for id_, elev in zip(G.nodes, elevations)} + nx.set_node_attributes(G, elevations_dict, 'elevation') + return G + +@register_graphfcn +class set_surface_slope(BaseGraphFunction, + required_node_attributes = ['elevation'], + adds_edge_attributes = ['surface_slope']): + """set_surface_slope class.""" + + def __call__(self, G: nx.Graph, + **kwargs) -> nx.Graph: + """Set the surface slope for each edge. + + This function sets the surface slope for each edge. The surface slope is + calculated from the elevation data. + + Args: + G (nx.Graph): A graph + **kwargs: Additional keyword arguments are ignored. + + Returns: + G (nx.Graph): A graph + """ + G = G.copy() + # Compute the slope for each edge + slope_dict = {(u, v, k): (G.nodes[u]['elevation'] - G.nodes[v]['elevation']) + / d['length'] for u, v, k, d in G.edges(data=True, + keys=True)} + + # Set the 'surface_slope' attribute for all edges + nx.set_edge_attributes(G, slope_dict, 'surface_slope') + return G + +@register_graphfcn +class set_chahinan_angle(BaseGraphFunction, + required_node_attributes = ['x','y'], + adds_edge_attributes = ['chahinan_angle']): + """set_chahinan_angle class.""" + + def __call__(self, G: nx.Graph, + **kwargs) -> nx.Graph: + """Set the Chahinan angle for each edge. + + This function sets the Chahinan angle for each edge. The Chahinan angle is + calculated from the geometry of the edge and weighted according to the + angle (based on: https://doi.org/10.1016/j.compenvurbsys.2019.101370) + + Args: + G (nx.Graph): A graph + **kwargs: Additional keyword arguments are ignored. + + Returns: + G (nx.Graph): A graph + """ + # TODO - in a double directed graph, not sure how meaningful this is + # TODO could probably refactor + G = G.copy() + for u,v,d in G.edges(data=True): + min_weight = float('inf') + for node in G.successors(v): + if node == u: + continue + + p1 = (G.nodes[u]['x'], G.nodes[u]['y']) + p2 = (G.nodes[v]['x'], G.nodes[v]['y']) + p3 = (G.nodes[node]['x'], G.nodes[node]['y']) + angle = go.calculate_angle(p1,p2,p3) + chahinan_weight = np.interp(angle, + [0, 90, 135, 180, 225, 270, 360], + [1, 0.2, 0.7, 0, 0.7, 0.2, 1] + ) + min_weight = min(chahinan_weight, min_weight) + if min_weight == float('inf'): + min_weight = 0 + d['chahinan_angle'] = min_weight + return G + +@register_graphfcn +class calculate_weights(BaseGraphFunction, + required_edge_attributes = + parameters.TopologyDerivation().weights, + adds_edge_attributes = ['weight']): + """calculate_weights class.""" + # TODO.. I guess if someone defines their own weights, this will need + # to change, will want an automatic way to do that... + + def __call__(self, G: nx.Graph, + topo_derivation: parameters.TopologyDerivation, + **kwargs) -> nx.Graph: + """Calculate the weights for each edge. + + This function calculates the weights for each edge. The weights are + calculated from the edge attributes. + + Args: + G (nx.Graph): A graph + topo_derivation (parameters.TopologyDerivation): A TopologyDerivation + parameter object + **kwargs: Additional keyword arguments are ignored. + + Returns: + G (nx.Graph): A graph + """ + # Calculate bounds to normalise between + bounds: Dict[Any, List[float]] = defaultdict(lambda: [np.Inf, -np.Inf]) + + for (u, v, d), w in product(G.edges(data=True), + topo_derivation.weights): + bounds[w][0] = min(bounds[w][0], d.get(w, np.Inf)) + bounds[w][1] = max(bounds[w][1], d.get(w, -np.Inf)) + + G = G.copy() + for u, v, d in G.edges(data=True): + total_weight = 0 + for attr, bds in bounds.items(): + # Normalise + weight = (d[attr] - bds[0]) / (bds[1] - bds[0]) + # Exponent + weight = weight ** getattr(topo_derivation,f'{attr}_exponent') + # Scaling + weight = weight * getattr(topo_derivation,f'{attr}_scaling') + # Sum + total_weight += weight + # Set + d['weight'] = total_weight + return G + +@register_graphfcn +class identify_outlets(BaseGraphFunction, + required_edge_attributes = ['length', 'edge_type'], + required_node_attributes = ['x', 'y']): + """identify_outlets class.""" + + def __call__(self, G: nx.Graph, + outlet_derivation: parameters.OutletDerivation, + **kwargs) -> nx.Graph: + """Identify outlets in a combined river-street graph. + + This function identifies outlets in a combined river-street graph. An + outlet is a node that is connected to a river and a street. + + Args: + G (nx.Graph): A graph + outlet_derivation (parameters.OutletDerivation): An OutletDerivation + parameter object + **kwargs: Additional keyword arguments are ignored. + + Returns: + G (nx.Graph): A graph + """ + G = G.copy() + river_points = {} + street_points = {} + + # Get the points for each river and street node + for u, v, d in G.edges(data=True): + upoint = shapely.Point(G.nodes[u]['x'], G.nodes[u]['y']) + vpoint = shapely.Point(G.nodes[v]['x'], G.nodes[v]['y']) + if d['edge_type'] == 'river': + river_points[u] = upoint + river_points[v] = vpoint + else: + street_points[u] = upoint + street_points[v] = vpoint + + # Pair up the river and street nodes + matched_outlets = go.nearest_node_buffer(river_points, + street_points, + outlet_derivation.river_buffer_distance) + + # Copy graph to run shortest path on + G_ = G.copy() + + # 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 + G_.add_edge(street_id, river_id, + **{'length' : outlet_derivation.outlet_length, + 'edge_type' : 'outlet', + 'id' : f'{street_id}-{river_id}-outlet'}) + + # Add edges from the river nodes to a waste node + for river_node in river_points.keys(): + if G.out_degree(river_node) == 0: + G_.add_edge(river_node, + 'waste', + **{'length' : 0, + 'edge_type' : 'outlet', + 'id' : f'{river_node}-waste-outlet'}) + + # Set the length of the river edges to 0 - from a design perspective + # once water is in the river we don't care about the length - since it + # costs nothing + for u,v,d in G_.edges(data=True): + if d['edge_type'] == 'river': + d['length'] = 0 + + # Find shortest path to identify only 'efficient' outlets. The MST + # makes sense here over shortest path as each node is only allowed to + # be visited once - thus encouraging fewer outlets. In shortest path + # nodes near rivers will always just pick their nearest river node. + G_ = nx.minimum_spanning_tree(G_.to_undirected(), + weight = 'length') + + # Retain the shortest path outlets in the original graph + for u,v,d in G_.edges(data=True): + if (d['edge_type'] == 'outlet') & (v != 'waste'): + G.add_edge(u,v,**d) + + return G + +@register_graphfcn +class derive_topology(BaseGraphFunction, + required_edge_attributes = ['edge_type', # 'rivers' and 'streets' + 'weight'], + adds_node_attributes = ['outlet', 'shortest_path']): + """derive_topology class.""" + + + def __call__(self, G: nx.Graph, + **kwargs) -> nx.Graph: + """Derive the topology of a graph. + + Runs a djiikstra-based algorithm to identify the shortest path from each + node to its nearest outlet (weighted by the 'weight' edge value). The + returned graph is one that only contains the edges that feature on the + shortest paths. + + Args: + G (nx.Graph): A graph + **kwargs: Additional keyword arguments are ignored. + + Returns: + G (nx.Graph): A graph + """ + G = G.copy() + + # Identify outlets + outlets = [u for u,v,d in G.edges(data=True) if d['edge_type'] == 'outlet'] + + # Remove non-street edges/nodes and unconnected nodes + nodes_to_remove = [] + for u, v, d in G.edges(data=True): + if d['edge_type'] != 'street': + if d['edge_type'] == 'outlet': + nodes_to_remove.append(v) + else: + nodes_to_remove.append(u) + nodes_to_remove.append(v) + + isolated_nodes = list(nx.isolates(G)) + + for u in set(nodes_to_remove).union(isolated_nodes): + G.remove_node(u) + + # Check for negative cycles + if nx.negative_edge_cycle(G, weight = 'weight'): + raise ValueError('Graph contains negative cycle') + + # Initialize the dictionary with infinity for all nodes + shortest_paths = {node: float('inf') for node in G.nodes} + + # Initialize the dictionary to store the paths + paths: dict[Hashable,list] = {node: [] for node in G.nodes} + + # Set the shortest path length to 0 for outlets + for outlet in outlets: + shortest_paths[outlet] = 0 + paths[outlet] = [outlet] + + # Initialize a min-heap with (distance, node) tuples + heap = [(0, outlet) for outlet in outlets] + while heap: + # Pop the node with the smallest distance + dist, node = heappop(heap) + + # For each neighbor of the current node + for neighbor, _, edge_data in G.in_edges(node, data=True): + # Calculate the distance through the current node + alt_dist = dist + edge_data['weight'] + # If the alternative distance is shorter + + if alt_dist >= shortest_paths[neighbor]: + continue + + # Update the shortest path length + shortest_paths[neighbor] = alt_dist + # Update the path + paths[neighbor] = paths[node] + [neighbor] + # Push the neighbor to the heap + heappush(heap, (alt_dist, neighbor)) + + edges_to_keep = set() + for path in paths.values(): + # Assign outlet + outlet = path[0] + for node in path: + G.nodes[node]['outlet'] = outlet + G.nodes[node]['shortest_path'] = shortest_paths[node] + + # Store path + for i in range(len(path) - 1): + edges_to_keep.add((path[i+1], path[i])) + + # Remvoe edges not on paths + new_graph = G.copy() + for u,v in G.edges(): + if (u,v) not in edges_to_keep: + new_graph.remove_edge(u,v) + return new_graph + +def design_pipe(ds_elevation: float, + chamber_floor: float, + edge_length: float, + pipe_design: parameters.HydraulicDesign, + Q: float + ) -> nx.Graph: + """Design a pipe. + + This function designs a pipe by iterating over a range of diameters and + depths. It returns the diameter and depth of the pipe that minimises the + cost function, while also maintaining or minimising feasibility parameters + associated with: surcharging, velocity and filling ratio. + + Args: + ds_elevation (float): The downstream elevationq + chamber_floor (float): The elevation of the chamber floor + edge_length (float): The length of the edge + pipe_design (parameters.HydraulicDesign): A HydraulicDesign parameter + object + Q (float): The flow rate + + Returns: + diam (float): The diameter of the pipe + depth (float): The depth of the pipe + """ + designs = product(pipe_design.diameters, + np.linspace(pipe_design.min_depth, + pipe_design.max_depth, + 10) # TODO should 10 be a param? + ) + pipes = [] + for diam, depth in designs: + A = (np.pi * diam ** 2 / 4) + n = 0.012 # mannings n + R = A / (np.pi * diam) # hydraulic radius + # TODO... presumably need to check depth > (diam + min_depth) + + elev_diff = chamber_floor - (ds_elevation - depth) + slope = elev_diff / edge_length + # Always pick a pipe that is feasible without surcharging + # if available + surcharge_feasibility = 0.0 + # Use surcharged elevation + while slope <= 0: + surcharge_feasibility += 0.05 + slope = (chamber_floor + surcharge_feasibility - + (ds_elevation - depth)) / edge_length + # TODO could make the feasibility penalisation increase + # when you get above surface_elevation[node]... but + # then you'd need a feasibility tracker and an offset + # tracker + v = (slope ** 0.5) * (R ** (2/3)) / n + filling_ratio = Q / (v * A) + # buffers from: https://www.polypipe.com/sites/default/files/Specification_Clauses_Underground_Drainage.pdf + average_depth = (depth + chamber_floor) / 2 + V = edge_length * (diam + 0.3) * (average_depth + 0.1) + cost = 1.32 / 2000 * (9579.31 * diam ** 0.5737 + 1153.77 * V**1.31) + v_feasibility = max(pipe_design.min_v - v, 0) + \ + max(v - pipe_design.max_v, 0) + fr_feasibility = max(filling_ratio - pipe_design.max_fr, 0) + """ + TODO shear stress... got confused here + density = 1000 + dyn_visc = 0.001 + hydraulic_diameter = 4 * (A * filling_ratio**2) / \ + (np.pi * diam * filling_ratio) + Re = density * v * 2 * (diam / 4) * (filling_ratio ** 2) / dyn_visc + fd = 64 / Re + shear_stress = fd * density * v**2 / fd + shear_feasibility = max(min_shear - shear_stress, 0) + """ + slope = (chamber_floor - (ds_elevation - depth)) / edge_length + pipes.append({'diam' : diam, + 'depth' : depth, + 'slope' : slope, + 'v' : v, + 'fr' : filling_ratio, + # 'tau' : shear_stress, + 'cost' : cost, + 'v_feasibility' : v_feasibility, + 'fr_feasibility' : fr_feasibility, + 'surcharge_feasibility' : surcharge_feasibility, + # 'shear_feasibility' : shear_feasibility + }) + + pipes_df = pd.DataFrame(pipes).dropna() + if pipes_df.shape[0] > 0: + ideal_pipe = pipes_df.sort_values(by=['surcharge_feasibility', + 'v_feasibility', + 'fr_feasibility', + # 'shear_feasibility', + 'cost'], + ascending = True).iloc[0] + return ideal_pipe.diam, ideal_pipe.depth + else: + raise Exception('something odd - no non nan pipes') + +def process_successors(G: nx.Graph, + node: Hashable, + surface_elevations: dict[Hashable, float], + chamber_floor: dict[Hashable, float], + edge_diams: dict[tuple[Hashable,Hashable,int], float], + pipe_design: parameters.HydraulicDesign + ) -> None: + """Process the successors of a node. + + This function processes the successors of a node. It designs a pipe to the + downstream node and sets the diameter and downstream invert level of the + pipe. It also sets the downstream invert level of the downstream node. It + returns None but modifies the edge_diams and chamber_floor dictionaries. + + Args: + G (nx.Graph): A graph + node (Hashable): A node + surface_elevations (dict): A dictionary of surface elevations keyed by + node + chamber_floor (dict): A dictionary of chamber floor elevations keyed by + node + edge_diams (dict): A dictionary of pipe diameters keyed by edge + pipe_design (parameters.HydraulicDesign): A HydraulicDesign parameter + object + """ + for ix, ds_node in enumerate(G.successors(node)): + edge = G.get_edge_data(node,ds_node,0) + # Find contributing area with ancestors + # TODO - could do timearea here if i hated myself enough + anc = nx.ancestors(G,node).union([node]) + tot = sum([G.nodes[anc_node]['contributing_area'] for anc_node in anc]) + + M3_PER_HR_TO_M3_PER_S = 1 / 60 / 60 + Q = tot * pipe_design.precipitation * M3_PER_HR_TO_M3_PER_S + + # Design the pipe to find the diameter and invert depth + diam, depth = design_pipe(surface_elevations[ds_node], + chamber_floor[node], + edge['length'], + pipe_design, + Q + ) + edge_diams[(node,ds_node,0)] = diam + chamber_floor[ds_node] = surface_elevations[ds_node] - depth + if ix > 0: + print('''a node has multiple successors, + not sure how that can happen if using shortest path + to derive topology''') + +@register_graphfcn +class pipe_by_pipe(BaseGraphFunction, + required_edge_attributes = ['length', 'elevation'], + required_node_attributes = ['contributing_area', 'elevation'], + adds_edge_attributes = ['diameter'], + adds_node_attributes = ['chamber_floor_elevation']): + """pipe_by_pipe class.""" + # If doing required_graph_attributes - it would be something like 'dendritic' + + def __call__(self, + G: nx.Graph, + pipe_design: parameters.HydraulicDesign, + **kwargs + )->nx.Graph: + """Pipe by pipe hydraulic design. + + Starting from the most upstream node, design a pipe to the downstream node + specifying a diameter and downstream invert level. A range of diameters and + invert levels are tested (ranging between conditions defined in + pipe_design). From the tested diameters/inverts, a selection is made based + on each pipe's satisfying feasibility constraints on: surcharge velocity, + filling ratio, (and shear stress - not currently implemented). Prioritising + feasibility in this order it identifies pipes with the preferred feasibility + level. If multiple pipes are feasible, it picks the lowest cost pipe. Once + the feasible pipe is identified, the diameter and downstream invert are set + and then the next downstream pipe can be designed. + + This approach is based on the pipe-by-pipe design proposed in: + https://doi.org/10.1016/j.watres.2021.117903 + + Args: + G (nx.Graph): A graph + pipe_design (parameters.HydraulicDesign): A HydraulicDesign parameter + object + **kwargs: Additional keyword arguments are ignored. + + Returns: + G (nx.Graph): A graph + """ + G = G.copy() + surface_elevations = {n : d['elevation'] for n, d in G.nodes(data=True)} + topological_order = list(nx.topological_sort(G)) + chamber_floor = {} + edge_diams: dict[tuple[Hashable,Hashable,int],float] = {} + # Iterate over nodes in topological order + for node in tqdm(topological_order): + # Check if there's any nodes upstream, if not set the depth to min_depth + if len(nx.ancestors(G,node)) == 0: + chamber_floor[node] = surface_elevations[node] - pipe_design.min_depth + + process_successors(G, + node, + surface_elevations, + chamber_floor, + edge_diams, + pipe_design + ) + + nx.function.set_edge_attributes(G, edge_diams, "diameter") + nx.function.set_node_attributes(G, chamber_floor, "chamber_floor_elevation") + return G \ No newline at end of file diff --git a/swmmanywhere/parameters.py b/swmmanywhere/parameters.py new file mode 100644 index 00000000..2f66dcf0 --- /dev/null +++ b/swmmanywhere/parameters.py @@ -0,0 +1,267 @@ +# -*- coding: utf-8 -*- +"""Created on 2024-01-26. + +@author: Barney +""" + +from pathlib import Path + +import numpy as np +from pydantic import BaseModel, Field, model_validator + + +class SubcatchmentDerivation(BaseModel): + """Parameters for subcatchment derivation.""" + lane_width: float = Field(default = 3.5, + ge = 2.0, + le = 5.0, + unit = "m", + description = "Width of a road lane.") + + carve_depth: float = Field(default = 2.0, + ge = 1.0, + le = 3.0, + unit = "m", + description = "Depth of road/river carve for flow accumulation.") + + max_street_length: float = Field(default = 60.0, + ge = 20.0, + le = 100.0, + unit = "m", + description = "Distance to split streets into segments.") + +class OutletDerivation(BaseModel): + """Parameters for outlet derivation.""" + max_river_length: float = Field(default = 30.0, + ge = 5.0, + le = 100.0, + unit = "m", + description = "Distance to split rivers into segments.") + + river_buffer_distance: float = Field(default = 150.0, + ge = 50.0, + le = 300.0, + unit = "m", + description = "Buffer distance to link rivers to streets.") + + outlet_length: float = Field(default = 40.0, + ge = 10.0, + le = 600.0, + unit = "m", + description = "Length to discourage street drainage into river buffers.") + +class TopologyDerivation(BaseModel): + """Parameters for topology derivation.""" + weights: list = Field(default = ['surface_slope', + 'chahinan_angle', + 'length', + 'contributing_area'], + min_items = 1, + unit = "-", + description = "Weights for topo derivation") + + surface_slope_scaling: float = Field(default = 1, + le = 1, + ge = 0, + unit = "-", + description = "Constant to apply to surface slope in topo derivation") + + chahinan_angle_scaling: float = Field(default = 1, + le = 1, + ge = 0, + unit = "-", + description = "Constant to apply to chahinan angle in topo derivation") + + length_scaling: float = Field(default = 1, + le = 1, + ge = 0, + unit = "-", + description = "Constant to apply to length in topo derivation") + + contributing_area_scaling: float = Field(default = 1, + le = 1, + ge = 0, + unit = "-", + description = "Constant to apply to contributing area in topo derivation") + + surface_slope_exponent: float = Field(default = 1, + le = 2, + ge = -2, + unit = "-", + description = "Exponent to apply to surface slope in topo derivation") + + chahinan_angle_exponent: float = Field(default = 1, + le = 2, + ge = -2, + unit = "-", + description = "Exponent to apply to chahinan angle in topo derivation") + + length_exponent: float = Field(default = 1, + le = 2, + ge = -2, + unit = "-", + description = "Exponent to apply to length in topo derivation") + + contributing_area_exponent: float = Field(default = 1, + le = 2, + ge = -2, + unit = "-", + description = "Exponent to apply to contributing area in topo derivation") + + + + @model_validator(mode='after') + def check_weights(cls, values): + """Check that weights have associated scaling and exponents.""" + for weight in values.weights: + if not hasattr(values, f'{weight}_scaling'): + raise ValueError(f"Missing {weight}_scaling") + if not hasattr(values, f'{weight}_exponent'): + raise ValueError(f"Missing {weight}_exponent") + return values + +# TODO move this to tests and run it if we're happy with this way of doing things +class NewTopo(TopologyDerivation): + """Demo for changing weights that should break the validator.""" + weights: list = Field(default = ['surface_slope', + 'chahinan_angle', + 'length', + 'contributing_area', + 'test'], + min_items = 1, + unit = "-", + description = "Weights for topo derivation") + +class HydraulicDesign(BaseModel): + """Parameters for hydraulic design.""" + diameters: list = Field(default = np.linspace(0.15,3,int((3-0.15)/0.075) + 1), + min_items = 1, + unit = "m", + description = """Diameters to consider in + pipe by pipe method""") + max_fr: float = Field(default = 0.8, + upper_limit = 1, + lower_limit = 0, + unit = "-", + description = "Maximum filling ratio in pipe by pipe method") + min_shear: float = Field(default = 2, + upper_limit = 3, + lower_limit = 0, + unit = "Pa", + description = "Minimum wall shear stress in pipe by pipe method") + min_v: float = Field(default = 0.75, + upper_limit = 2, + lower_limit = 0, + unit = "m/s", + description = "Minimum velocity in pipe by pipe method") + max_v: float = Field(default = 5, + upper_limit = 10, + lower_limit = 3, + unit = "m/s", + description = "Maximum velocity in pipe by pipe method") + min_depth: float = Field(default = 0.5, + upper_limit = 1, + lower_limit = 0, + unit = "m", + description = "Minimum excavation depth in pipe by pipe method") + max_depth: float = Field(default = 5, + upper_limit = 10, + lower_limit = 2, + unit = "m", + description = "Maximum excavation depth in pipe by pipe method") + precipitation: float = Field(default = 0.006, + upper_limit = 0.010, + lower_limit = 0.001, + description = "Depth of design storm in pipe by pipe method", + unit = "m") + +class FilePaths: + """Parameters for file path lookup. + + 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). + """ + + def __init__(self, + base_dir: Path, + project_name: str, + bbox_number: int, + model_number: str, + extension: str='json'): + """Initialise the class.""" + self.base_dir = base_dir + self.project_name = project_name + self.bbox_number = bbox_number + self.model_number = model_number + self.extension = extension + + def __getattr__(self, name): + """Fetch the address.""" + return self._fetch_address(name) + + def _generate_path(self, *subdirs): + """Generate a path.""" + return self.base_dir.joinpath(*subdirs) + + def _generate_property(self, + property_name: str, + location: str): + """Generate a property. + + Check if the property exists in the class, otherwise generate it. + + Args: + property_name (str): Name of the folder/file. + location (str): Name of the folder that the property_name exists + in. + + Returns: + Path: Path to the property. + """ + if property_name in self.__dict__.keys(): + return self.__dict__[property_name] + else: + return self._generate_path(self.project_name, + getattr(self, location), + property_name) + + def _fetch_address(self, name): + """Fetch the address. + + Generate a path to the folder/file described by name. If the + folder/file has already been set, then it will be returned. Otherwise + it will be generated according to the default structure defined below. + + Args: + name (str): Name of the folder/file. + + Returns: + Path: Path to the folder/file. + """ + if name == 'project': + return self._generate_path(self.project_name) + elif name == 'national': + return self._generate_property('national', 'project') + elif name == 'bbox': + return self._generate_property(f'bbox_{self.bbox_number}', + 'project') + elif name == 'model': + return self._generate_property(f'model_{self.model_number}', + 'bbox') + elif name == 'subcatchments': + return self._generate_property(f'subcatchments.{self.extension}', + 'model') + elif name == 'download': + return self._generate_property('download', + 'bbox') + elif name == 'elevation': + return self._generate_property('elevation.tif', 'download') + elif name == 'building': + return self._generate_property(f'building.{self.extension}', + 'download') + elif name == 'precipitation': + return self._generate_property(f'precipitation.{self.extension}', + 'download') + else: + raise AttributeError(f"Attribute {name} not found") + \ No newline at end of file diff --git a/swmmanywhere/prepare_data.py b/swmmanywhere/prepare_data.py index 951e5f38..bdbbdafb 100644 --- a/swmmanywhere/prepare_data.py +++ b/swmmanywhere/prepare_data.py @@ -52,7 +52,7 @@ def get_country(x: float, data = yaml.safe_load(file) # Get country ISO code from coordinates - location = geolocator.reverse("{0}, {1}".format(y, x), exactly_one=True) + location = geolocator.reverse(f"{y}, {x}", exactly_one=True) iso_country_code = location.raw.get("address", {}).get("country_code") iso_country_code = iso_country_code.upper() @@ -219,7 +219,7 @@ def download_precipitation(bbox: tuple[float, float, float, float], } c = cdsapi.Client(url='https://cds.climate.copernicus.eu/api/v2', - key='{0}:{1}'.format(username, api_key)) + key=f'{username}:{api_key}') # Get data c.retrieve('reanalysis-era5-single-levels', request, diff --git a/tests/test_data/elevation.tif b/tests/test_data/elevation.tif new file mode 100644 index 00000000..c1ec0382 Binary files /dev/null and b/tests/test_data/elevation.tif differ diff --git a/tests/test_data/graph_topo_derived.json b/tests/test_data/graph_topo_derived.json new file mode 100644 index 00000000..f129673b --- /dev/null +++ b/tests/test_data/graph_topo_derived.json @@ -0,0 +1 @@ +{"directed": true, "multigraph": true, "graph": {"created_date": "2024-01-31 12:27:58", "created_with": "OSMnx 1.8.1", "crs": "EPSG:32630", "simplified": true}, "nodes": [{"y": 5709902.110360867, "x": 700277.5569639901, "street_count": 3, "outlet": 12354833, "shortest_path": 51, "elevation": 0, "contributing_area": 0, "id": 107733}, {"y": 5709939.408081475, "x": 700277.6837391303, "street_count": 3, "outlet": 12354833, "shortest_path": 51, "elevation": 1, "contributing_area": 1, "id": 107734}, {"y": 5709950.664565533, "x": 700311.2287216682, "street_count": 3, "outlet": 12354833, "shortest_path": 55, "elevation": 2, "contributing_area": 2, "id": 107735}, {"y": 5709945.54203236, "x": 700328.4358706471, "street_count": 3, "outlet": 12354833, "shortest_path": 61, "elevation": 3, "contributing_area": 3, "id": 107736}, {"y": 5709897.32580564, "x": 700345.9401207004, "street_count": 3, "outlet": 12354833, "shortest_path": 111, "elevation": 4, "contributing_area": 4, "id": 107737}, {"y": 5709878.507201574, "x": 700314.3472000923, "street_count": 3, "outlet": 12354833, "shortest_path": 52, "elevation": 5, "contributing_area": 5, "id": 107738}, {"y": 5710021.9827727135, "x": 700465.4044789741, "street_count": 4, "outlet": 25472347, "shortest_path": 60, "elevation": 6, "contributing_area": 6, "id": 109753}, {"y": 5709694.8741118815, "x": 700044.4345900133, "highway": "turning_circle", "street_count": 3, "outlet": 12354833, "shortest_path": 0, "elevation": 7, "contributing_area": 7, "id": 12354833}, {"y": 5709859.5404103035, "x": 700260.646795129, "street_count": 3, "outlet": 12354833, "shortest_path": 52, "elevation": 8, "contributing_area": 8, "id": 21392086}, {"y": 5709943.283586768, "x": 700130.0316394513, "street_count": 3, "outlet": 25472347, "shortest_path": 0, "elevation": 9, "contributing_area": 9, "id": 25472347}, {"y": 5709918.64998135, "x": 700182.2264369107, "street_count": 1, "outlet": 25472347, "shortest_path": 25, "elevation": 10, "contributing_area": 10, "id": 25472348}, {"y": 5709665.057833289, "x": 700141.1675964806, "street_count": 4, "outlet": 12354833, "shortest_path": 21, "elevation": 11, "contributing_area": 11, "id": 25472373}, {"y": 5709861.753318606, "x": 700486.534735846, "street_count": 3, "outlet": 12354833, "shortest_path": 201, "elevation": 12, "contributing_area": 12, "id": 25472468}, {"y": 5710133.23541422, "x": 700406.4117507454, "street_count": 3, "outlet": 25472347, "shortest_path": 27, "elevation": 13, "contributing_area": 13, "id": 25472854}, {"y": 5709836.784993341, "x": 700419.6382007168, "street_count": 3, "outlet": 12354833, "shortest_path": 166, "elevation": 14, "contributing_area": 14, "id": 25472893}, {"y": 5709805.458949697, "x": 700446.664728181, "street_count": 3, "outlet": 12354833, "shortest_path": 74, "elevation": 15, "contributing_area": 15, "id": 25510321}, {"y": 5709878.072050431, "x": 700330.9734851549, "street_count": 3, "outlet": 12354833, "shortest_path": 69, "elevation": 16, "contributing_area": 16, "id": 32925453}, {"y": 5710484.769616377, "x": 699976.6805692167, "street_count": 1, "outlet": 12354833, "shortest_path": 63, "elevation": 17, "contributing_area": 17, "id": 248122264}, {"y": 5709800.829107233, "x": 700450.6767434096, "street_count": 3, "outlet": 12354833, "shortest_path": 113, "elevation": 18, "contributing_area": 18, "id": 266325461}, {"y": 5709931.922096681, "x": 700343.3585345388, "street_count": 3, "outlet": 12354833, "shortest_path": 71, "elevation": 19, "contributing_area": 19, "id": 770549936}, {"y": 5709959.908314285, "x": 700369.9469864559, "street_count": 3, "outlet": 12354833, "shortest_path": 70, "elevation": 20, "contributing_area": 20, "id": 1696030874}, {"y": 5710395.710567596, "x": 700021.9512225334, "street_count": 3, "outlet": 12354833, "shortest_path": 54, "elevation": 21, "contributing_area": 21, "id": 1881001588}, {"y": 5709815.86155232, "x": 700437.6921531621, "street_count": 2, "outlet": 12354833, "shortest_path": 112, "elevation": 22, "contributing_area": 22, "id": 2623975694}, {"y": 5709873.623166366, "x": 700362.1435140749, "street_count": 3, "outlet": 12354833, "shortest_path": 112, "elevation": 23, "contributing_area": 23, "id": 6277683849}], "links": [{"osmid": [157379106, 5035977, 724068972, 1068268847, 1068268848, 323063155, 371146076], "oneway": true, "lanes": "2", "ref": "A3200", "name": "York Road", "highway": "primary", "maxspeed": "20 mph", "reversed": false, "length": 275.046, "geometry": "LINESTRING (700141.1675964806 5709665.057833289, 700143.5157680553 5709677.065346105, 700144.8474277792 5709679.601053864, 700155.976942681 5709700.673855638, 700160.7768450535 5709709.19243132, 700172.0393874501 5709729.179248003, 700185.4803273565 5709752.960108612, 700191.8141976679 5709764.167181919, 700209.5905451218 5709797.66218869, 700216.9778449582 5709810.525491558, 700224.2875415371 5709823.24098586, 700228.962895601 5709831.387259596, 700233.8496671842 5709839.98729696, 700236.7225068416 5709844.499128756, 700241.8766513304 5709850.548563866, 700261.0016506243 5709870.979407257, 700265.1289820187 5709876.632025824, 700267.7934926444 5709880.4341245955, 700270.6096913561 5709884.4426481845, 700272.9768698463 5709889.791986443, 700274.1428695226 5709894.938038215, 700277.5569639901 5709902.110360867)", "width": 0, "id": "157379106.reversed", "edge_type": "street", "weight": 2, "source": 107733, "target": 25472373, "key": 0}, {"osmid": 4040995, "oneway": true, "lanes": "3", "ref": "A301", "name": "Tenison Way", "highway": "primary", "maxspeed": "20 mph", "junction": "circular", "reversed": false, "length": 37.45, "geometry": "LINESTRING (700277.5569639901 5709902.110360867, 700276.5294887233 5709907.1921423655, 700276.0456861216 5709918.397616499, 700277.6837391303 5709939.408081475)", "width": 0, "id": "4040995.reversed", "edge_type": "street", "weight": 5, "source": 107734, "target": 107733, "key": 0}, {"osmid": [253278314, 24935919], "oneway": true, "lanes": ["2", "1"], "ref": "A301", "highway": "primary", "maxspeed": "20 mph", "junction": "circular", "reversed": false, "length": 38.699000000000005, "geometry": "LINESTRING (700277.6837391303 5709939.408081475, 700281.0620739642 5709943.082457772, 700283.3080724685 5709945.865862207, 700286.8370136549 5709949.245527317, 700291.5898057348 5709951.548802148, 700295.19372216 5709952.147561552, 700298.571576143 5709953.015790332, 700302.8354936898 5709952.738624432, 700306.8351147252 5709951.938800403, 700311.2287216682 5709950.664565533)", "width": 0, "id": "253278314.reversed", "edge_type": "street", "weight": 7, "source": 107735, "target": 107734, "key": 0}, {"osmid": 24935923, "oneway": true, "lanes": "3", "ref": "A301", "highway": "primary", "maxspeed": "20 mph", "junction": "circular", "reversed": false, "length": 17.9, "geometry": "LINESTRING (700311.2287216682 5709950.664565533, 700328.4358706471 5709945.54203236)", "width": 0, "id": "24935923.reversed", "edge_type": "street", "weight": 11, "source": 107736, "target": 107735, "key": 0}, {"osmid": [228563659, 24935916], "oneway": true, "lanes": "3", "ref": "A301", "highway": "primary", "maxspeed": "20 mph", "junction": "circular", "reversed": false, "length": 24.528999999999996, "geometry": "LINESTRING (700345.9401207004 5709897.32580564, 700343.9882393219 5709893.429296876, 700342.7545134749 5709891.409623581, 700341.5918455586 5709889.5263809515, 700338.673571965 5709885.992611689, 700337.8972975462 5709885.060000548, 700330.9734851549 5709878.072050431)", "width": 0, "id": 228563659, "edge_type": "street", "weight": 12, "source": 107737, "target": 32925453, "key": 0}, {"osmid": 24935915, "oneway": true, "lanes": "2", "ref": "A301", "highway": "primary", "maxspeed": "20 mph", "junction": "circular", "reversed": false, "length": 43.962999999999994, "geometry": "LINESTRING (700314.3472000923 5709878.507201574, 700306.2843434939 5709882.42051824, 700299.5883693597 5709885.508079189, 700294.6641010015 5709888.253542988, 700289.9725552127 5709892.322180368, 700284.8256592196 5709897.363912976, 700277.5569639901 5709902.110360867)", "width": 0, "id": 24935915, "edge_type": "street", "weight": 16, "source": 107738, "target": 107733, "key": 0}, {"osmid": 3115950, "name": "Cornwall Road", "highway": "unclassified", "maxspeed": "20 mph", "oneway": false, "reversed": true, "length": 125.803, "geometry": "LINESTRING (700465.4044789741 5710021.9827727135, 700461.4439979517 5710029.710301515, 700458.6490779112 5710035.022915902, 700456.6658586246 5710038.7863448765, 700444.9310978653 5710061.072771147, 700430.5406586338 5710086.205261883, 700416.6051415108 5710113.371273859, 700407.1722169924 5710131.751011451, 700406.4117507454 5710133.23541422)", "width": 0, "id": 3115950, "edge_type": "street", "weight": 18, "source": 109753, "target": 25472854, "key": 0}, {"osmid": [725226531, 307619438, 1068268846, 1068268849, 323063156, 42507450], "oneway": true, "lanes": "2", "ref": "A3200", "name": "York Road", "highway": "primary", "maxspeed": "20 mph", "reversed": false, "length": 228.89699999999996, "geometry": "LINESTRING (700260.646795129 5709859.5404103035, 700253.727705374 5709851.372400684, 700250.164089411 5709846.755359551, 700246.2016118084 5709841.855335842, 700242.4905154916 5709837.277025176, 700233.1152255475 5709823.722810672, 700219.3327136479 5709798.569795305, 700216.9147153466 5709794.153864184, 700212.9726675999 5709786.971906629, 700204.4524102686 5709770.411498058, 700185.1938145162 5709735.37705126, 700169.0912214798 5709706.658454267, 700156.226779984 5709684.114146786, 700153.8828750349 5709679.578684956, 700151.8880911011 5709676.059180303, 700150.8420187996 5709674.035825833, 700141.1675964806 5709665.057833289)", "width": 0, "id": 725226531, "edge_type": "street", "weight": 23, "source": 21392086, "target": 25472373, "key": 0}, {"osmid": 4253584, "oneway": true, "name": "Concert Hall Approach", "highway": "living_street", "maxspeed": "20 mph", "reversed": false, "length": 57.56, "geometry": "LINESTRING (700130.0316394513 5709943.283586768, 700182.2264369107 5709918.64998135)", "width": 0, "id": "4253584.reversed", "edge_type": "street", "weight": 28, "source": 25472348, "target": 25472347, "key": 0}, {"osmid": 4253560, "name": "Chicheley Street", "highway": "residential", "maxspeed": "20 mph", "width": 0, "oneway": false, "reversed": false, "length": 101.119, "geometry": "LINESTRING (700141.1675964806 5709665.057833289, 700131.2293895038 5709668.763878157, 700127.7724262028 5709669.897037534, 700123.7282148367 5709671.118403503, 700082.5852502533 5709684.718811786, 700059.2118402369 5709691.993324356, 700056.0186594417 5709692.958758778, 700050.1884291518 5709694.532943944, 700044.4345900133 5709694.8741118815)", "id": 4253560, "edge_type": "street", "weight": 29, "source": 25472373, "target": 12354833, "key": 0}, {"osmid": 4253650, "name": "Exton Street", "highway": "unclassified", "maxspeed": "20 mph", "oneway": false, "reversed": true, "length": 71.213, "geometry": "LINESTRING (700486.534735846 5709861.753318606, 700454.4513036673 5709848.849677348, 700429.4004307707 5709839.998849332, 700426.7727709545 5709839.137886363, 700419.6382007168 5709836.784993341)", "width": 0, "id": 4253650, "edge_type": "street", "weight": 32, "source": 25472468, "target": 25472893, "key": 0}, {"osmid": [243937860, 223569640, 2423885, 1222144887, 376614608, 545538484, 545538485, 1222144886, 753458071, 545538486, 425225428, 1222144890, 753458070], "name": ["Upper Ground", "Belvedere Road"], "highway": ["residential", "unclassified"], "maxspeed": "20 mph", "oneway": false, "reversed": [false, true], "length": 342.27299999999997, "geometry": "LINESTRING (700406.4117507454 5710133.23541422, 700387.069359366 5710122.928634459, 700372.9230384012 5710116.011803224, 700369.1779954545 5710114.2381824525, 700363.0648608486 5710111.36889115, 700349.9293818473 5710105.249243745, 700327.3277679619 5710097.052294423, 700311.837243554 5710091.441121343, 700285.8403761584 5710080.526944091, 700264.3797259327 5710070.88315461, 700260.658661233 5710068.676280703, 700255.3352328222 5710065.515354177, 700249.9942922055 5710061.741291486, 700245.8203461394 5710058.859569429, 700242.7033414934 5710056.665404704, 700224.1362021959 5710044.307535186, 700221.4609648484 5710042.364655568, 700218.4947517125 5710040.399164544, 700199.8158378075 5710027.5247511575, 700180.0977083285 5710009.810042807, 700167.8430164627 5709996.4987158, 700161.3400735646 5709988.881739712, 700154.4322770827 5709980.7811199045, 700142.0556038325 5709963.333798693, 700132.9058145428 5709947.762005055, 700130.0316394513 5709943.283586768)", "width": 0, "id": 243937860, "edge_type": "street", "weight": 34, "source": 25472854, "target": 25472347, "key": 0}, {"osmid": [1207319057, 256768302], "oneway": true, "lanes": "1", "ref": "A301", "name": "Waterloo Road", "highway": "primary", "maxspeed": "20 mph", "reversed": false, "length": 29.652, "geometry": "LINESTRING (700419.6382007168 5709836.784993341, 700427.245395453 5709833.856050073, 700430.6112952617 5709829.924490856, 700432.30694177 5709828.1652200045, 700435.93212929 5709824.355256594, 700437.6921531621 5709815.86155232)", "width": 0, "id": 1207319057, "edge_type": "street", "weight": 37, "source": 25472893, "target": 2623975694, "key": 0}, {"osmid": [437211562, 4041071, 409351991, 359842682, 323062908], "oneway": true, "name": "Mepham Street", "highway": "unclassified", "maxspeed": "20 mph", "reversed": false, "length": 202.61500000000004, "geometry": "LINESTRING (700260.646795129 5709859.5404103035, 700264.606979534 5709848.115749547, 700270.2421185977 5709844.262482439, 700276.2468099121 5709841.258961956, 700280.3130668337 5709840.005186374, 700290.7579613369 5709836.486469852, 700317.3910231055 5709825.499927625, 700334.9227590943 5709819.911372949, 700341.4646550727 5709818.265384199, 700350.8427009505 5709815.89617834, 700383.4195720793 5709807.683451332, 700407.4894184885 5709803.36655895, 700438.1292499496 5709800.5675001955, 700442.1522457196 5709802.875500894, 700446.664728181 5709805.458949697)", "width": 0, "id": "437211562.reversed", "edge_type": "street", "weight": 40, "source": 25510321, "target": 21392086, "key": 0}, {"osmid": 207115006, "oneway": true, "lanes": "3", "ref": "A301", "highway": "primary", "maxspeed": "20 mph", "junction": "circular", "reversed": false, "length": 16.852, "geometry": "LINESTRING (700330.9734851549 5709878.072050431, 700325.9959171354 5709876.884555874, 700321.5250918429 5709877.120134854, 700314.3472000923 5709878.507201574)", "width": 0, "id": 207115006, "edge_type": "street", "weight": 41, "source": 32925453, "target": 107738, "key": 0}, {"osmid": [225025827, 200598186, 200598189, 228564878, 928025711, 35900943, 944217725], "bridge": "yes", "oneway": true, "lanes": ["2", "1"], "ref": "A301", "name": "Waterloo Bridge", "highway": "primary", "maxspeed": "20 mph", "reversed": false, "length": 630.8910000000001, "geometry": "LINESTRING (699976.6805692167 5710484.769616377, 699995.0223578266 5710453.266369763, 700014.2757750297 5710422.968325912, 700023.9057410788 5710408.003189852, 700042.265373482 5710374.4741320945, 700154.1424509434 5710191.497236617, 700176.0592502693 5710155.5588213885, 700180.562426344 5710148.153163577, 700195.976710407 5710122.748642229, 700201.8983423435 5710115.499165202, 700239.9790628456 5710051.914439671, 700276.6209293818 5709991.179535505, 700279.6394946901 5709986.710817551, 700285.3315518312 5709977.537060874, 700292.9918154946 5709966.91541046, 700298.4768209211 5709960.172166666, 700301.9810694954 5709956.780495522, 700305.95837635 5709953.552255062, 700311.2287216682 5709950.664565533)", "width": 0, "id": 225025827, "edge_type": "street", "weight": 44, "source": 248122264, "target": 107735, "key": 0}, {"osmid": 699097039, "lanes": "4", "ref": "A301", "name": "Waterloo Road", "highway": "primary", "maxspeed": "20 mph", "oneway": false, "reversed": true, "length": 6.116, "geometry": "LINESTRING (700450.6767434096 5709800.829107233, 700446.664728181 5709805.458949697)", "width": 0, "id": 699097039, "edge_type": "street", "weight": 45, "source": 266325461, "target": 25510321, "key": 0}, {"osmid": 228562216, "oneway": true, "lanes": "2", "ref": "A301", "highway": "primary", "maxspeed": "20 mph", "junction": "circular", "reversed": false, "length": 20.893, "geometry": "LINESTRING (700328.4358706471 5709945.54203236, 700333.3858597063 5709943.3766637165, 700337.0812959784 5709941.128379782, 700341.0533185115 5709936.274174798, 700341.9863743042 5709934.818844887, 700343.3585345388 5709931.922096681)", "width": 0, "id": "228562216.reversed", "edge_type": "street", "weight": 47, "source": 770549936, "target": 107736, "key": 0}, {"osmid": [83801193, 251642745], "oneway": true, "ref": "A3200", "name": "Stamford Street", "highway": "primary", "maxspeed": "20 mph", "reversed": false, "length": 44.416, "geometry": "LINESTRING (700328.4358706471 5709945.54203236, 700334.077684354 5709945.920587545, 700339.4473318071 5709946.856317236, 700344.5164533333 5709948.715572506, 700348.4576651177 5709950.452366831, 700352.1663099085 5709952.447236735, 700356.7277008864 5709955.021407143, 700360.3472058966 5709957.513862129, 700369.9469864559 5709959.908314285)", "width": 0, "id": "83801193.reversed", "edge_type": "street", "weight": 51, "source": 1696030874, "target": 107736, "key": 0}, {"osmid": [200596577, 224637155, 948685543, 948685544, 670348489, 23013971, 4253331], "oneway": true, "lanes": ["2", "1"], "ref": "A301", "name": "Waterloo Bridge", "highway": "primary", "maxspeed": "20 mph", "reversed": false, "length": 526.335, "bridge": "yes", "geometry": "LINESTRING (700277.6837391303 5709939.408081475, 700275.7128986834 5709948.672989107, 700275.6714573554 5709953.069868358, 700275.1963905097 5709961.235701989, 700274.8376413824 5709970.152202144, 700274.4762318035 5709974.380560487, 700273.746804977 5709978.249198908, 700272.3189850443 5709982.557972666, 700270.3512829997 5709986.812039583, 700243.8802948985 5710030.354161469, 700218.2029132167 5710070.876579723, 700212.394118875 5710080.725057999, 700179.172642919 5710137.051966281, 700177.2608309576 5710142.7113350835, 700172.1774497676 5710151.085170908, 700040.4501411039 5710365.171293509, 700029.6364467517 5710382.840209466, 700021.9512225334 5710395.710567596)", "width": 0, "id": "200596577.reversed", "edge_type": "street", "weight": 52, "source": 1881001588, "target": 107734, "key": 0}, {"osmid": 256768303, "lanes": "4", "ref": "A301", "name": "Waterloo Road", "highway": "primary", "maxspeed": "20 mph", "oneway": false, "reversed": false, "length": 13.715, "geometry": "LINESTRING (700437.6921531621 5709815.86155232, 700446.664728181 5709805.458949697)", "width": 0, "id": 256768303, "edge_type": "street", "weight": 53, "source": 2623975694, "target": 25510321, "key": 0}, {"osmid": 228568372, "oneway": true, "lanes": "2", "ref": "A301", "name": "Waterloo Road", "highway": "primary", "maxspeed": "20 mph", "reversed": false, "length": 32.080999999999996, "geometry": "LINESTRING (700362.1435140749 5709873.623166366, 700356.7236303644 5709872.37364107, 700356.0138620024 5709872.04496849, 700351.4285008686 5709872.36508479, 700346.2485579546 5709873.6639322005, 700338.4782655792 5709875.974101273, 700330.9734851549 5709878.072050431)", "width": 0, "id": 228568372, "edge_type": "street", "weight": 56, "source": 6277683849, "target": 32925453, "key": 0}]} \ No newline at end of file diff --git a/tests/test_data/street_graph.json b/tests/test_data/street_graph.json new file mode 100644 index 00000000..b1f9a2e0 --- /dev/null +++ b/tests/test_data/street_graph.json @@ -0,0 +1 @@ +{"directed": true, "multigraph": true, "graph": {"created_date": "2024-01-31 12:27:58", "created_with": "OSMnx 1.8.1", "crs": "EPSG:32630", "simplified": true}, "nodes": [{"y": 5709902.110360867, "x": 700277.5569639901, "street_count": 3, "id": 107733}, {"y": 5709939.408081475, "x": 700277.6837391303, "street_count": 3, "id": 107734}, {"y": 5709950.664565533, "x": 700311.2287216682, "street_count": 3, "id": 107735}, {"y": 5709945.54203236, "x": 700328.4358706471, "street_count": 3, "id": 107736}, {"y": 5709897.32580564, "x": 700345.9401207004, "street_count": 3, "id": 107737}, {"y": 5709878.507201574, "x": 700314.3472000923, "street_count": 3, "id": 107738}, {"y": 5710021.9827727135, "x": 700465.4044789741, "street_count": 4, "id": 109753}, {"y": 5709694.8741118815, "x": 700044.4345900133, "highway": "turning_circle", "street_count": 3, "id": 12354833}, {"y": 5709859.5404103035, "x": 700260.646795129, "street_count": 3, "id": 21392086}, {"y": 5709943.283586768, "x": 700130.0316394513, "street_count": 3, "id": 25472347}, {"y": 5709918.64998135, "x": 700182.2264369107, "street_count": 1, "id": 25472348}, {"y": 5709665.057833289, "x": 700141.1675964806, "street_count": 4, "id": 25472373}, {"y": 5709861.753318606, "x": 700486.534735846, "street_count": 3, "id": 25472468}, {"y": 5710133.23541422, "x": 700406.4117507454, "street_count": 3, "id": 25472854}, {"y": 5709836.784993341, "x": 700419.6382007168, "street_count": 3, "id": 25472893}, {"y": 5709805.458949697, "x": 700446.664728181, "street_count": 3, "id": 25510321}, {"y": 5709878.072050431, "x": 700330.9734851549, "street_count": 3, "id": 32925453}, {"y": 5710484.769616377, "x": 699976.6805692167, "street_count": 1, "id": 248122264}, {"y": 5709800.829107233, "x": 700450.6767434096, "street_count": 3, "id": 266325461}, {"y": 5709931.922096681, "x": 700343.3585345388, "street_count": 3, "id": 770549936}, {"y": 5709959.908314285, "x": 700369.9469864559, "street_count": 3, "id": 1696030874}, {"y": 5710395.710567596, "x": 700021.9512225334, "street_count": 3, "id": 1881001588}, {"y": 5709815.86155232, "x": 700437.6921531621, "street_count": 2, "id": 2623975694}, {"y": 5709873.623166366, "x": 700362.1435140749, "street_count": 3, "id": 6277683849}], "links": [{"osmid": 4040995, "oneway": true, "lanes": "3", "ref": "A301", "name": "Tenison Way", "highway": "primary", "maxspeed": "20 mph", "junction": "circular", "reversed": false, "length": 37.45, "geometry": "LINESTRING (700277.5569639901 5709902.110360867, 700276.5294887233 5709907.1921423655, 700276.0456861216 5709918.397616499, 700277.6837391303 5709939.408081475)", "width": 0, "id": 4040995, "source": 107733, "target": 107734, "key": 0}, {"osmid": [200596577, 224637155, 948685543, 948685544, 670348489, 23013971, 4253331], "oneway": true, "lanes": ["2", "1"], "ref": "A301", "name": "Waterloo Bridge", "highway": "primary", "maxspeed": "20 mph", "reversed": false, "length": 526.335, "bridge": "yes", "geometry": "LINESTRING (700277.6837391303 5709939.408081475, 700275.7128986834 5709948.672989107, 700275.6714573554 5709953.069868358, 700275.1963905097 5709961.235701989, 700274.8376413824 5709970.152202144, 700274.4762318035 5709974.380560487, 700273.746804977 5709978.249198908, 700272.3189850443 5709982.557972666, 700270.3512829997 5709986.812039583, 700243.8802948985 5710030.354161469, 700218.2029132167 5710070.876579723, 700212.394118875 5710080.725057999, 700179.172642919 5710137.051966281, 700177.2608309576 5710142.7113350835, 700172.1774497676 5710151.085170908, 700040.4501411039 5710365.171293509, 700029.6364467517 5710382.840209466, 700021.9512225334 5710395.710567596)", "width": 0, "id": [200596577, 224637155, 948685543, 948685544, 670348489, 23013971, 4253331], "source": 107734, "target": 1881001588, "key": 0}, {"osmid": [253278314, 24935919], "oneway": true, "lanes": ["2", "1"], "ref": "A301", "highway": "primary", "maxspeed": "20 mph", "junction": "circular", "reversed": false, "length": 38.699000000000005, "geometry": "LINESTRING (700277.6837391303 5709939.408081475, 700281.0620739642 5709943.082457772, 700283.3080724685 5709945.865862207, 700286.8370136549 5709949.245527317, 700291.5898057348 5709951.548802148, 700295.19372216 5709952.147561552, 700298.571576143 5709953.015790332, 700302.8354936898 5709952.738624432, 700306.8351147252 5709951.938800403, 700311.2287216682 5709950.664565533)", "width": 0, "id": [253278314, 24935919], "source": 107734, "target": 107735, "key": 0}, {"osmid": 24935923, "oneway": true, "lanes": "3", "ref": "A301", "highway": "primary", "maxspeed": "20 mph", "junction": "circular", "reversed": false, "length": 17.9, "geometry": "LINESTRING (700311.2287216682 5709950.664565533, 700328.4358706471 5709945.54203236)", "width": 0, "id": 24935923, "source": 107735, "target": 107736, "key": 0}, {"osmid": [83801193, 251642745], "oneway": true, "ref": "A3200", "name": "Stamford Street", "highway": "primary", "maxspeed": "20 mph", "reversed": false, "length": 44.416, "geometry": "LINESTRING (700328.4358706471 5709945.54203236, 700334.077684354 5709945.920587545, 700339.4473318071 5709946.856317236, 700344.5164533333 5709948.715572506, 700348.4576651177 5709950.452366831, 700352.1663099085 5709952.447236735, 700356.7277008864 5709955.021407143, 700360.3472058966 5709957.513862129, 700369.9469864559 5709959.908314285)", "width": 0, "id": [83801193, 251642745], "source": 107736, "target": 1696030874, "key": 0}, {"osmid": 228562216, "oneway": true, "lanes": "2", "ref": "A301", "highway": "primary", "maxspeed": "20 mph", "junction": "circular", "reversed": false, "length": 20.893, "geometry": "LINESTRING (700328.4358706471 5709945.54203236, 700333.3858597063 5709943.3766637165, 700337.0812959784 5709941.128379782, 700341.0533185115 5709936.274174798, 700341.9863743042 5709934.818844887, 700343.3585345388 5709931.922096681)", "width": 0, "id": 228562216, "source": 107736, "target": 770549936, "key": 0}, {"osmid": [228563659, 24935916], "oneway": true, "lanes": "3", "ref": "A301", "highway": "primary", "maxspeed": "20 mph", "junction": "circular", "reversed": false, "length": 24.528999999999996, "geometry": "LINESTRING (700345.9401207004 5709897.32580564, 700343.9882393219 5709893.429296876, 700342.7545134749 5709891.409623581, 700341.5918455586 5709889.5263809515, 700338.673571965 5709885.992611689, 700337.8972975462 5709885.060000548, 700330.9734851549 5709878.072050431)", "width": 0, "id": [228563659, 24935916], "source": 107737, "target": 32925453, "key": 0}, {"osmid": [4041072, 670353288, 1163014149, 1163014150], "oneway": true, "lanes": "1", "ref": "A301", "name": "Waterloo Road", "highway": "primary", "maxspeed": "20 mph", "reversed": false, "length": 29.020999999999997, "geometry": "LINESTRING (700345.9401207004 5709897.32580564, 700349.2380955268 5709892.823608096, 700351.2189495779 5709887.534484578, 700353.5553747164 5709883.272725233, 700355.2608794149 5709880.7343362635, 700357.6691066612 5709877.82280663, 700362.1435140749 5709873.623166366)", "width": 0, "id": [4041072, 670353288, 1163014149, 1163014150], "source": 107737, "target": 6277683849, "key": 0}, {"osmid": [228568056, 5035978, 670631404], "oneway": true, "lanes": "2", "ref": "A3200", "name": "York Road", "highway": "primary", "maxspeed": "20 mph", "reversed": false, "length": 60.64600000000001, "geometry": "LINESTRING (700314.3472000923 5709878.507201574, 700310.6387952155 5709877.7372764675, 700305.7560426467 5709877.845259409, 700299.4779484478 5709878.443827132, 700292.0008729227 5709878.783519834, 700286.8243099415 5709878.412237795, 700282.0695676366 5709877.389475359, 700277.0818624147 5709875.578042671, 700272.9440150785 5709872.129823205, 700262.6592841079 5709861.9025773555, 700260.646795129 5709859.5404103035)", "width": 0, "id": [228568056, 5035978, 670631404], "source": 107738, "target": 21392086, "key": 0}, {"osmid": 24935915, "oneway": true, "lanes": "2", "ref": "A301", "highway": "primary", "maxspeed": "20 mph", "junction": "circular", "reversed": false, "length": 43.962999999999994, "geometry": "LINESTRING (700314.3472000923 5709878.507201574, 700306.2843434939 5709882.42051824, 700299.5883693597 5709885.508079189, 700294.6641010015 5709888.253542988, 700289.9725552127 5709892.322180368, 700284.8256592196 5709897.363912976, 700277.5569639901 5709902.110360867)", "width": 0, "id": 24935915, "source": 107738, "target": 107733, "key": 0}, {"osmid": 3115950, "name": "Cornwall Road", "highway": "unclassified", "maxspeed": "20 mph", "oneway": false, "reversed": true, "length": 125.803, "geometry": "LINESTRING (700465.4044789741 5710021.9827727135, 700461.4439979517 5710029.710301515, 700458.6490779112 5710035.022915902, 700456.6658586246 5710038.7863448765, 700444.9310978653 5710061.072771147, 700430.5406586338 5710086.205261883, 700416.6051415108 5710113.371273859, 700407.1722169924 5710131.751011451, 700406.4117507454 5710133.23541422)", "width": 0, "id": 3115950, "source": 109753, "target": 25472854, "key": 0}, {"osmid": [1015370840, 1068268844, 1068268845], "ref": "A3200", "name": "Stamford Street", "highway": "primary", "maxspeed": "20 mph", "oneway": false, "reversed": true, "length": 113.57300000000001, "geometry": "LINESTRING (700465.4044789741 5710021.9827727135, 700453.098476616 5710014.035996378, 700450.5354608949 5710012.420353288, 700445.3293230028 5710008.929791383, 700410.4907315228 5709986.274230827, 700382.0277860094 5709967.756897977, 700369.9469864559 5709959.908314285)", "width": 0, "id": [1015370840, 1068268844, 1068268845], "source": 109753, "target": 1696030874, "key": 0}, {"osmid": [240228802, 2731434, 371146645], "name": "Belvedere Road", "highway": "residential", "maxspeed": "20 mph", "oneway": false, "reversed": false, "length": 263.92800000000005, "geometry": "LINESTRING (700044.4345900133 5709694.8741118815, 700045.1858113637 5709708.7896099, 700045.753666956 5709712.7205270715, 700048.2667245618 5709723.899323567, 700054.3328123755 5709746.14201081, 700060.8472276325 5709768.469191039, 700070.128409068 5709800.470777224, 700081.6009791871 5709835.231271199, 700085.2635685451 5709846.332917027, 700087.2385952899 5709852.46845025, 700091.0908391678 5709862.107698213, 700095.4972763445 5709872.848934459, 700097.9855905913 5709878.303173154, 700103.8378423693 5709891.150335212, 700105.4044098433 5709894.953602542, 700115.1380090672 5709913.8778502615, 700130.0316394513 5709943.283586768)", "width": 0, "id": [240228802, 2731434, 371146645], "source": 12354833, "target": 25472347, "key": 0}, {"osmid": 4253560, "name": "Chicheley Street", "highway": "residential", "maxspeed": "20 mph", "width": 0, "oneway": false, "reversed": true, "length": 101.119, "geometry": "LINESTRING (700044.4345900133 5709694.8741118815, 700050.1884291518 5709694.532943944, 700056.0186594417 5709692.958758778, 700059.2118402369 5709691.993324356, 700082.5852502533 5709684.718811786, 700123.7282148367 5709671.118403503, 700127.7724262028 5709669.897037534, 700131.2293895038 5709668.763878157, 700141.1675964806 5709665.057833289)", "id": 4253560, "source": 12354833, "target": 25472373, "key": 0}, {"osmid": [437211562, 4041071, 409351991, 359842682, 323062908], "oneway": true, "name": "Mepham Street", "highway": "unclassified", "maxspeed": "20 mph", "reversed": false, "length": 202.61500000000004, "geometry": "LINESTRING (700260.646795129 5709859.5404103035, 700264.606979534 5709848.115749547, 700270.2421185977 5709844.262482439, 700276.2468099121 5709841.258961956, 700280.3130668337 5709840.005186374, 700290.7579613369 5709836.486469852, 700317.3910231055 5709825.499927625, 700334.9227590943 5709819.911372949, 700341.4646550727 5709818.265384199, 700350.8427009505 5709815.89617834, 700383.4195720793 5709807.683451332, 700407.4894184885 5709803.36655895, 700438.1292499496 5709800.5675001955, 700442.1522457196 5709802.875500894, 700446.664728181 5709805.458949697)", "width": 0, "id": [437211562, 4041071, 409351991, 359842682, 323062908], "source": 21392086, "target": 25510321, "key": 0}, {"osmid": [725226531, 307619438, 1068268846, 1068268849, 323063156, 42507450], "oneway": true, "lanes": "2", "ref": "A3200", "name": "York Road", "highway": "primary", "maxspeed": "20 mph", "reversed": false, "length": 228.89699999999996, "geometry": "LINESTRING (700260.646795129 5709859.5404103035, 700253.727705374 5709851.372400684, 700250.164089411 5709846.755359551, 700246.2016118084 5709841.855335842, 700242.4905154916 5709837.277025176, 700233.1152255475 5709823.722810672, 700219.3327136479 5709798.569795305, 700216.9147153466 5709794.153864184, 700212.9726675999 5709786.971906629, 700204.4524102686 5709770.411498058, 700185.1938145162 5709735.37705126, 700169.0912214798 5709706.658454267, 700156.226779984 5709684.114146786, 700153.8828750349 5709679.578684956, 700151.8880911011 5709676.059180303, 700150.8420187996 5709674.035825833, 700141.1675964806 5709665.057833289)", "width": 0, "id": [725226531, 307619438, 1068268846, 1068268849, 323063156, 42507450], "source": 21392086, "target": 25472373, "key": 0}, {"osmid": 4253584, "oneway": true, "name": "Concert Hall Approach", "highway": "living_street", "maxspeed": "20 mph", "reversed": false, "length": 57.56, "geometry": "LINESTRING (700130.0316394513 5709943.283586768, 700182.2264369107 5709918.64998135)", "width": 0, "id": 4253584, "source": 25472347, "target": 25472348, "key": 0}, {"osmid": [2731434, 240228802, 371146645], "name": "Belvedere Road", "highway": "residential", "maxspeed": "20 mph", "oneway": false, "reversed": true, "length": 263.928, "geometry": "LINESTRING (700130.0316394513 5709943.283586768, 700115.1380090672 5709913.8778502615, 700105.4044098433 5709894.953602542, 700103.8378423693 5709891.150335212, 700097.9855905913 5709878.303173154, 700095.4972763445 5709872.848934459, 700091.0908391678 5709862.107698213, 700087.2385952899 5709852.46845025, 700085.2635685451 5709846.332917027, 700081.6009791871 5709835.231271199, 700070.128409068 5709800.470777224, 700060.8472276325 5709768.469191039, 700054.3328123755 5709746.14201081, 700048.2667245618 5709723.899323567, 700045.753666956 5709712.7205270715, 700045.1858113637 5709708.7896099, 700044.4345900133 5709694.8741118815)", "width": 0, "id": [2731434, 240228802, 371146645], "source": 25472347, "target": 12354833, "key": 0}, {"osmid": [243937860, 223569640, 2423885, 1222144887, 376614608, 753458071, 425225428, 545538485, 545538486, 753458070, 1222144886, 545538484, 1222144890], "name": ["Upper Ground", "Belvedere Road"], "highway": ["residential", "unclassified"], "maxspeed": "20 mph", "oneway": false, "reversed": [false, true], "length": 342.27299999999997, "geometry": "LINESTRING (700130.0316394513 5709943.283586768, 700132.9058145428 5709947.762005055, 700142.0556038325 5709963.333798693, 700154.4322770827 5709980.7811199045, 700161.3400735646 5709988.881739712, 700167.8430164627 5709996.4987158, 700180.0977083285 5710009.810042807, 700199.8158378075 5710027.5247511575, 700218.4947517125 5710040.399164544, 700221.4609648484 5710042.364655568, 700224.1362021959 5710044.307535186, 700242.7033414934 5710056.665404704, 700245.8203461394 5710058.859569429, 700249.9942922055 5710061.741291486, 700255.3352328222 5710065.515354177, 700260.658661233 5710068.676280703, 700264.3797259327 5710070.88315461, 700285.8403761584 5710080.526944091, 700311.837243554 5710091.441121343, 700327.3277679619 5710097.052294423, 700349.9293818473 5710105.249243745, 700363.0648608486 5710111.36889115, 700369.1779954545 5710114.2381824525, 700372.9230384012 5710116.011803224, 700387.069359366 5710122.928634459, 700406.4117507454 5710133.23541422)", "width": 0, "id": [243937860, 223569640, 2423885, 1222144887, 376614608, 753458071, 425225428, 545538485, 545538486, 753458070, 1222144886, 545538484, 1222144890], "source": 25472347, "target": 25472854, "key": 0}, {"osmid": 4253560, "name": "Chicheley Street", "highway": "residential", "maxspeed": "20 mph", "width": 0, "oneway": false, "reversed": false, "length": 101.119, "geometry": "LINESTRING (700141.1675964806 5709665.057833289, 700131.2293895038 5709668.763878157, 700127.7724262028 5709669.897037534, 700123.7282148367 5709671.118403503, 700082.5852502533 5709684.718811786, 700059.2118402369 5709691.993324356, 700056.0186594417 5709692.958758778, 700050.1884291518 5709694.532943944, 700044.4345900133 5709694.8741118815)", "id": 4253560, "source": 25472373, "target": 12354833, "key": 0}, {"osmid": [157379106, 5035977, 724068972, 1068268847, 1068268848, 323063155, 371146076], "oneway": true, "lanes": "2", "ref": "A3200", "name": "York Road", "highway": "primary", "maxspeed": "20 mph", "reversed": false, "length": 275.046, "geometry": "LINESTRING (700141.1675964806 5709665.057833289, 700143.5157680553 5709677.065346105, 700144.8474277792 5709679.601053864, 700155.976942681 5709700.673855638, 700160.7768450535 5709709.19243132, 700172.0393874501 5709729.179248003, 700185.4803273565 5709752.960108612, 700191.8141976679 5709764.167181919, 700209.5905451218 5709797.66218869, 700216.9778449582 5709810.525491558, 700224.2875415371 5709823.24098586, 700228.962895601 5709831.387259596, 700233.8496671842 5709839.98729696, 700236.7225068416 5709844.499128756, 700241.8766513304 5709850.548563866, 700261.0016506243 5709870.979407257, 700265.1289820187 5709876.632025824, 700267.7934926444 5709880.4341245955, 700270.6096913561 5709884.4426481845, 700272.9768698463 5709889.791986443, 700274.1428695226 5709894.938038215, 700277.5569639901 5709902.110360867)", "width": 0, "id": [157379106, 5035977, 724068972, 1068268847, 1068268848, 323063155, 371146076], "source": 25472373, "target": 107733, "key": 0}, {"osmid": 4253650, "name": "Exton Street", "highway": "unclassified", "maxspeed": "20 mph", "oneway": false, "reversed": true, "length": 71.213, "geometry": "LINESTRING (700486.534735846 5709861.753318606, 700454.4513036673 5709848.849677348, 700429.4004307707 5709839.998849332, 700426.7727709545 5709839.137886363, 700419.6382007168 5709836.784993341)", "width": 0, "id": 4253650, "source": 25472468, "target": 25472893, "key": 0}, {"osmid": 3115950, "name": "Cornwall Road", "highway": "unclassified", "maxspeed": "20 mph", "oneway": false, "reversed": false, "length": 125.803, "geometry": "LINESTRING (700406.4117507454 5710133.23541422, 700407.1722169924 5710131.751011451, 700416.6051415108 5710113.371273859, 700430.5406586338 5710086.205261883, 700444.9310978653 5710061.072771147, 700456.6658586246 5710038.7863448765, 700458.6490779112 5710035.022915902, 700461.4439979517 5710029.710301515, 700465.4044789741 5710021.9827727135)", "width": 0, "id": 3115950, "source": 25472854, "target": 109753, "key": 0}, {"osmid": [243937860, 223569640, 2423885, 1222144887, 376614608, 545538484, 545538485, 1222144886, 753458071, 545538486, 425225428, 1222144890, 753458070], "name": ["Upper Ground", "Belvedere Road"], "highway": ["residential", "unclassified"], "maxspeed": "20 mph", "oneway": false, "reversed": [false, true], "length": 342.27299999999997, "geometry": "LINESTRING (700406.4117507454 5710133.23541422, 700387.069359366 5710122.928634459, 700372.9230384012 5710116.011803224, 700369.1779954545 5710114.2381824525, 700363.0648608486 5710111.36889115, 700349.9293818473 5710105.249243745, 700327.3277679619 5710097.052294423, 700311.837243554 5710091.441121343, 700285.8403761584 5710080.526944091, 700264.3797259327 5710070.88315461, 700260.658661233 5710068.676280703, 700255.3352328222 5710065.515354177, 700249.9942922055 5710061.741291486, 700245.8203461394 5710058.859569429, 700242.7033414934 5710056.665404704, 700224.1362021959 5710044.307535186, 700221.4609648484 5710042.364655568, 700218.4947517125 5710040.399164544, 700199.8158378075 5710027.5247511575, 700180.0977083285 5710009.810042807, 700167.8430164627 5709996.4987158, 700161.3400735646 5709988.881739712, 700154.4322770827 5709980.7811199045, 700142.0556038325 5709963.333798693, 700132.9058145428 5709947.762005055, 700130.0316394513 5709943.283586768)", "width": 0, "id": [243937860, 223569640, 2423885, 1222144887, 376614608, 545538484, 545538485, 1222144886, 753458071, 545538486, 425225428, 1222144890, 753458070], "source": 25472854, "target": 25472347, "key": 0}, {"osmid": 4253650, "name": "Exton Street", "highway": "unclassified", "maxspeed": "20 mph", "oneway": false, "reversed": false, "length": 71.213, "geometry": "LINESTRING (700419.6382007168 5709836.784993341, 700426.7727709545 5709839.137886363, 700429.4004307707 5709839.998849332, 700454.4513036673 5709848.849677348, 700486.534735846 5709861.753318606)", "width": 0, "id": 4253650, "source": 25472893, "target": 25472468, "key": 0}, {"osmid": [725226529, 670353289, 725226540, 107677719], "lanes": "4", "ref": "A301", "name": "Waterloo Road", "highway": "primary", "maxspeed": "20 mph", "oneway": false, "reversed": true, "length": 68.741, "geometry": "LINESTRING (700419.6382007168 5709836.784993341, 700415.9275751942 5709840.647271395, 700409.9926658429 5709846.281364492, 700399.9372265527 5709854.191464778, 700387.1601167364 5709862.0386844715, 700378.9869482361 5709866.10343341, 700370.5343034118 5709869.678336915, 700362.1435140749 5709873.623166366)", "width": 0, "id": [725226529, 670353289, 725226540, 107677719], "source": 25472893, "target": 6277683849, "key": 0}, {"osmid": [1207319057, 256768302], "oneway": true, "lanes": "1", "ref": "A301", "name": "Waterloo Road", "highway": "primary", "maxspeed": "20 mph", "reversed": false, "length": 29.652, "geometry": "LINESTRING (700419.6382007168 5709836.784993341, 700427.245395453 5709833.856050073, 700430.6112952617 5709829.924490856, 700432.30694177 5709828.1652200045, 700435.93212929 5709824.355256594, 700437.6921531621 5709815.86155232)", "width": 0, "id": [1207319057, 256768302], "source": 25472893, "target": 2623975694, "key": 0}, {"osmid": 256768303, "lanes": "4", "ref": "A301", "name": "Waterloo Road", "highway": "primary", "maxspeed": "20 mph", "oneway": false, "reversed": true, "length": 13.715, "geometry": "LINESTRING (700446.664728181 5709805.458949697, 700437.6921531621 5709815.86155232)", "width": 0, "id": 256768303, "source": 25510321, "target": 2623975694, "key": 0}, {"osmid": 699097039, "lanes": "4", "ref": "A301", "name": "Waterloo Road", "highway": "primary", "maxspeed": "20 mph", "oneway": false, "reversed": false, "length": 6.116, "geometry": "LINESTRING (700446.664728181 5709805.458949697, 700450.6767434096 5709800.829107233)", "width": 0, "id": 699097039, "source": 25510321, "target": 266325461, "key": 0}, {"osmid": 207115006, "oneway": true, "lanes": "3", "ref": "A301", "highway": "primary", "maxspeed": "20 mph", "junction": "circular", "reversed": false, "length": 16.852, "geometry": "LINESTRING (700330.9734851549 5709878.072050431, 700325.9959171354 5709876.884555874, 700321.5250918429 5709877.120134854, 700314.3472000923 5709878.507201574)", "width": 0, "id": 207115006, "source": 32925453, "target": 107738, "key": 0}, {"osmid": [225025827, 200598186, 200598189, 228564878, 928025711, 35900943, 944217725], "bridge": "yes", "oneway": true, "lanes": ["2", "1"], "ref": "A301", "name": "Waterloo Bridge", "highway": "primary", "maxspeed": "20 mph", "reversed": false, "length": 630.8910000000001, "geometry": "LINESTRING (699976.6805692167 5710484.769616377, 699995.0223578266 5710453.266369763, 700014.2757750297 5710422.968325912, 700023.9057410788 5710408.003189852, 700042.265373482 5710374.4741320945, 700154.1424509434 5710191.497236617, 700176.0592502693 5710155.5588213885, 700180.562426344 5710148.153163577, 700195.976710407 5710122.748642229, 700201.8983423435 5710115.499165202, 700239.9790628456 5710051.914439671, 700276.6209293818 5709991.179535505, 700279.6394946901 5709986.710817551, 700285.3315518312 5709977.537060874, 700292.9918154946 5709966.91541046, 700298.4768209211 5709960.172166666, 700301.9810694954 5709956.780495522, 700305.95837635 5709953.552255062, 700311.2287216682 5709950.664565533)", "width": 0, "id": [225025827, 200598186, 200598189, 228564878, 928025711, 35900943, 944217725], "source": 248122264, "target": 107735, "key": 0}, {"osmid": 699097039, "lanes": "4", "ref": "A301", "name": "Waterloo Road", "highway": "primary", "maxspeed": "20 mph", "oneway": false, "reversed": true, "length": 6.116, "geometry": "LINESTRING (700450.6767434096 5709800.829107233, 700446.664728181 5709805.458949697)", "width": 0, "id": 699097039, "source": 266325461, "target": 25510321, "key": 0}, {"osmid": 228562217, "oneway": true, "lanes": "4", "ref": "A301", "highway": "primary", "maxspeed": "20 mph", "junction": "circular", "reversed": false, "length": 34.67, "geometry": "LINESTRING (700343.3585345388 5709931.922096681, 700345.9401207004 5709897.32580564)", "width": 0, "id": 228562217, "source": 770549936, "target": 107737, "key": 0}, {"osmid": [251642744, 157380564], "oneway": true, "ref": "A3200", "name": "Stamford Street", "highway": "primary", "maxspeed": "20 mph", "reversed": false, "length": 38.856, "lanes": "1", "geometry": "LINESTRING (700369.9469864559 5709959.908314285, 700363.816512409 5709952.895947482, 700359.7617173692 5709950.163598803, 700355.2675727794 5709946.5898859985, 700352.9588338891 5709944.516642271, 700348.8771816592 5709939.823396935, 700347.6618187202 5709937.514923466, 700343.3585345388 5709931.922096681)", "width": 0, "id": [251642744, 157380564], "source": 1696030874, "target": 770549936, "key": 0}, {"osmid": [1015370840, 1068268844, 1068268845], "ref": "A3200", "name": "Stamford Street", "highway": "primary", "maxspeed": "20 mph", "oneway": false, "reversed": false, "length": 113.57300000000001, "geometry": "LINESTRING (700369.9469864559 5709959.908314285, 700382.0277860094 5709967.756897977, 700410.4907315228 5709986.274230827, 700445.3293230028 5710008.929791383, 700450.5354608949 5710012.420353288, 700453.098476616 5710014.035996378, 700465.4044789741 5710021.9827727135)", "width": 0, "id": [1015370840, 1068268844, 1068268845], "source": 1696030874, "target": 109753, "key": 0}, {"osmid": 256768303, "lanes": "4", "ref": "A301", "name": "Waterloo Road", "highway": "primary", "maxspeed": "20 mph", "oneway": false, "reversed": false, "length": 13.715, "geometry": "LINESTRING (700437.6921531621 5709815.86155232, 700446.664728181 5709805.458949697)", "width": 0, "id": 256768303, "source": 2623975694, "target": 25510321, "key": 0}, {"osmid": [1163014824, 1163014825, 256768300], "oneway": true, "lanes": "2", "ref": "A301", "name": "Waterloo Road", "highway": "primary", "maxspeed": "20 mph", "reversed": false, "length": 29.411, "geometry": "LINESTRING (700437.6921531621 5709815.86155232, 700429.4368226664 5709819.188040226, 700426.1114026883 5709822.798269364, 700423.1837722174 5709826.190360595, 700419.6382007168 5709836.784993341)", "width": 0, "id": [1163014824, 1163014825, 256768300], "source": 2623975694, "target": 25472893, "key": 0}, {"osmid": [670353289, 725226529, 725226540, 107677719], "lanes": "4", "ref": "A301", "name": "Waterloo Road", "highway": "primary", "maxspeed": "20 mph", "oneway": false, "reversed": false, "length": 68.741, "geometry": "LINESTRING (700362.1435140749 5709873.623166366, 700370.5343034118 5709869.678336915, 700378.9869482361 5709866.10343341, 700387.1601167364 5709862.0386844715, 700399.9372265527 5709854.191464778, 700409.9926658429 5709846.281364492, 700415.9275751942 5709840.647271395, 700419.6382007168 5709836.784993341)", "width": 0, "id": [670353289, 725226529, 725226540, 107677719], "source": 6277683849, "target": 25472893, "key": 0}, {"osmid": 228568372, "oneway": true, "lanes": "2", "ref": "A301", "name": "Waterloo Road", "highway": "primary", "maxspeed": "20 mph", "reversed": false, "length": 32.080999999999996, "geometry": "LINESTRING (700362.1435140749 5709873.623166366, 700356.7236303644 5709872.37364107, 700356.0138620024 5709872.04496849, 700351.4285008686 5709872.36508479, 700346.2485579546 5709873.6639322005, 700338.4782655792 5709875.974101273, 700330.9734851549 5709878.072050431)", "width": 0, "id": 228568372, "source": 6277683849, "target": 32925453, "key": 0}]} \ No newline at end of file diff --git a/tests/test_geospatial_utilities.py b/tests/test_geospatial_utilities.py index fcca0cdc..0a0a899c 100644 --- a/tests/test_geospatial_utilities.py +++ b/tests/test_geospatial_utilities.py @@ -4,9 +4,10 @@ @author: Barney """ -import os +from pathlib import Path from unittest.mock import MagicMock, patch +import geopandas as gpd import networkx as nx import numpy as np import rasterio as rst @@ -14,8 +15,14 @@ from shapely import geometry as sgeom from swmmanywhere import geospatial_utilities as go +from swmmanywhere import graph_utilities as ge +def load_street_network(): + """Load a street network.""" + G = ge.load_graph(Path(__file__).parent / 'test_data' / 'street_graph.json') + return G + def test_interp_with_nans(): """Test the interp_interp_with_nans function.""" # Define a simple grid and values @@ -62,10 +69,12 @@ def test_interpolate_points_on_raster(mock_rst_open): y = [0.25, 0.75] # Call the function - result = go.interpolate_points_on_raster(x, y, 'fake_path') + result = go.interpolate_points_on_raster(x, + y, + Path('fake_path')) - # [3,2] feels unintuitive but it's because rasters measure from the top - assert result == [3.0, 2.0] + # [2.75, 2.25] feels unintuitive but it's because rasters measure from the top + assert result == [2.75, 2.25] def test_get_utm(): """Test the get_utm_epsg function.""" @@ -94,29 +103,27 @@ def create_raster(fid): def test_reproject_raster(): """Test the reproject_raster function.""" # Create a mock raster file - fid = 'test.tif' + fid = Path('test.tif') try: create_raster(fid) # Define the input parameters target_crs = 'EPSG:32630' - new_fid = 'test_reprojected.tif' + new_fid = Path('test_reprojected.tif') # Call the function go.reproject_raster(target_crs, fid) # Check if the reprojected file exists - assert os.path.exists(new_fid) + assert new_fid.exists() # Check if the reprojected file has the correct CRS with rst.open(new_fid) as src: assert src.crs.to_string() == target_crs finally: # Regardless of test outcome, delete the temp file - if os.path.exists(fid): - os.remove(fid) - if os.path.exists(new_fid): - os.remove(new_fid) + fid.unlink(missing_ok=True) + new_fid.unlink(missing_ok=True) def almost_equal(a, b, tol=1e-6): @@ -190,8 +197,8 @@ def test_burn_shape_in_raster(): # Define the input parameters depth = 1.0 - raster_fid = 'input.tif' - new_raster_fid = 'output.tif' + raster_fid = Path('input.tif') + new_raster_fid = Path('output.tif') try: create_raster(raster_fid) @@ -207,7 +214,105 @@ def test_burn_shape_in_raster(): assert (data != data_).any() finally: # Regardless of test outcome, delete the temp file - if os.path.exists(raster_fid): - os.remove(raster_fid) - if os.path.exists(new_raster_fid): - os.remove(new_raster_fid) \ No newline at end of file + raster_fid.unlink(missing_ok=True) + new_raster_fid.unlink(missing_ok=True) + +def test_derive_subcatchments(): + """Test the derive_subcatchments function.""" + G = load_street_network() + elev_fid = Path(__file__).parent / 'test_data' / 'elevation.tif' + + polys = go.derive_subcatchments(G, elev_fid) + assert 'slope' in polys.columns + assert 'area' in polys.columns + assert 'geometry' in polys.columns + assert 'id' in polys.columns + assert polys.shape[0] > 0 + assert polys.dropna().shape == polys.shape + assert polys.crs == G.graph['crs'] + + +def test_derive_rc(): + """Test the derive_rc function.""" + G = load_street_network() + crs = G.graph['crs'] + eg_bldg = sgeom.Polygon([(700291,5709928), + (700331,5709927), + (700321,5709896), + (700293,5709900), + (700291,5709928)]) + buildings = gpd.GeoDataFrame(geometry = [eg_bldg], + crs = crs) + subs = [sgeom.Polygon([(700262, 5709928), + (700262, 5709883), + (700351, 5709883), + (700351, 5709906), + (700306, 5709906), + (700306, 5709928), + (700262, 5709928)]), + sgeom.Polygon([(700306, 5709928), + (700284, 5709928), + (700284, 5709950), + (700374, 5709950), + (700374, 5709906), + (700351, 5709906), + (700306, 5709906), + (700306, 5709928)]), + sgeom.Polygon([(700351, 5709883), + (700351, 5709906), + (700374, 5709906), + (700374, 5709883), + (700396, 5709883), + (700396, 5709816), + (700329, 5709816), + (700329, 5709838), + (700329, 5709883), + (700351, 5709883)])] + + subs = gpd.GeoDataFrame(data = {'id' : [107733, + 1696030874, + 6277683849] + }, + geometry = subs, + crs = crs) + subs['area'] = subs.geometry.area + + subs_rc = go.derive_rc(subs, G, 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') + 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 + point1 = (0, 0) + point2 = (1, 0) + point3 = (1, 1) + assert go.calculate_angle(point1, point2, point3) == 90 + + # Test with points forming a straight line + point1 = (0, 0) + point2 = (1, 0) + point3 = (2, 0) + assert go.calculate_angle(point1, point2, point3) == 180 + + # Test with points forming an angle of 45 degrees + point1 = (0, 0) + point2 = (1, 0) + point3 = (0, 1) + assert almost_equal(go.calculate_angle(point1, point2, point3), 45) + + # Test with points forming an angle of 0 degrees + point1 = (0, 0) + point2 = (1, 0) + point3 = (0, 0) + assert go.calculate_angle(point1, point2, point3) == 0 \ No newline at end of file diff --git a/tests/test_graph_utilities.py b/tests/test_graph_utilities.py new file mode 100644 index 00000000..87cb0f55 --- /dev/null +++ b/tests/test_graph_utilities.py @@ -0,0 +1,224 @@ +# -*- coding: utf-8 -*- +"""Created on 2024-01-26. + +@author: Barney +""" +import math +import tempfile +from pathlib import Path + +import geopandas as gpd +import networkx as nx +from shapely import geometry as sgeom + +from swmmanywhere import parameters +from swmmanywhere.graph_utilities import graphfcns as gu +from swmmanywhere.graph_utilities import load_graph, save_graph + + +def load_street_network(): + """Load a street network.""" + bbox = (-0.11643,51.50309,-0.11169,51.50549) + G = load_graph(Path(__file__).parent / 'test_data' / 'street_graph.json') + return G, bbox + +def test_save_load(): + """Test the save_graph and load_graph functions.""" + # Load a street network + G,_ = load_street_network() + with tempfile.TemporaryDirectory() as temp_dir: + # Save the graph + save_graph(G, Path(temp_dir) / 'test_graph.json') + # Load the graph + G_new = load_graph(Path(temp_dir) / 'test_graph.json') + # Check if the loaded graph is the same as the original graph + assert nx.is_isomorphic(G, G_new) + +def test_assign_id(): + """Test the assign_id function.""" + G, _ = load_street_network() + G = gu.assign_id(G) + for u, v, data in G.edges(data=True): + assert 'id' in data.keys() + assert isinstance(data['id'], int) + +def test_double_directed(): + """Test the double_directed function.""" + G, _ = load_street_network() + G = gu.assign_id(G) + G = gu.double_directed(G) + for u, v in G.edges(): + assert (v,u) in G.edges + +def test_format_osmnx_lanes(): + """Test the format_osmnx_lanes 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) + +def test_split_long_edges(): + """Test the split_long_edges function.""" + G, _ = load_street_network() + G = gu.assign_id(G) + max_length = 20 + params = parameters.SubcatchmentDerivation(max_street_length = max_length) + G = gu.split_long_edges(G, params) + for u, v, data in G.edges(data=True): + assert data['length'] <= (max_length * 2) + +def test_derive_subcatchments(): + """Test the derive_subcatchments function.""" + with tempfile.TemporaryDirectory() as temp_dir: + temp_path = Path(temp_dir) + addresses = parameters.FilePaths(base_dir = temp_path, + project_name = 'test', + bbox_number = 1, + extension = 'json', + model_number = 1) + addresses.elevation = Path(__file__).parent / 'test_data' / 'elevation.tif' + addresses.building = temp_path / 'building.geojson' + addresses.subcatchments = temp_path / 'subcatchments.geojson' + params = parameters.SubcatchmentDerivation() + G, bbox = load_street_network() + + # mock up buildings + eg_bldg = sgeom.Polygon([(700291.346,5709928.922), + (700331.206,5709927.815), + (700321.610,5709896.444), + (700293.192,5709900.503), + (700291.346,5709928.922)]) + gdf = gpd.GeoDataFrame(geometry = [eg_bldg], + crs = G.graph['crs']) + gdf.to_file(addresses.building, driver='GeoJSON') + + G = gu.calculate_contributing_area(G, params, addresses) + for u, v, data in G.edges(data=True): + assert 'contributing_area' in data.keys() + assert isinstance(data['contributing_area'], float) + + for u, data in G.nodes(data=True): + assert 'contributing_area' in data.keys() + assert isinstance(data['contributing_area'], float) + +def test_set_elevation_and_slope(): + """Test the set_elevation and set_surface_slope function.""" + G, _ = load_street_network() + with tempfile.TemporaryDirectory() as temp_dir: + temp_path = Path(temp_dir) + addresses = parameters.FilePaths(base_dir = temp_path, + project_name = 'test', + bbox_number = 1, + extension = 'json', + model_number = 1) + addresses.elevation = Path(__file__).parent / 'test_data' / 'elevation.tif' + G = gu.set_elevation(G, addresses) + for id_, data in G.nodes(data=True): + assert 'elevation' in data.keys() + assert math.isfinite(data['elevation']) + assert data['elevation'] > 0 + + G = gu.set_surface_slope(G) + for u, v, data in G.edges(data=True): + assert 'surface_slope' in data.keys() + assert math.isfinite(data['surface_slope']) + +def test_chahinan_angle(): + """Test the chahinan_angle function.""" + G, _ = load_street_network() + G = gu.set_chahinan_angle(G) + for u, v, data in G.edges(data=True): + assert 'chahinan_angle' in data.keys() + assert math.isfinite(data['chahinan_angle']) + +def test_calculate_weights(): + """Test the calculate_weights function.""" + G, _ = load_street_network() + params = parameters.TopologyDerivation() + for weight in params.weights: + for ix, (u,v,data) in enumerate(G.edges(data=True)): + data[weight] = ix + + G = gu.calculate_weights(G, params) + for u, v, data in G.edges(data=True): + assert 'weight' in data.keys() + assert math.isfinite(data['weight']) + +def test_identify_outlets_and_derive_topology(): + """Test the identify_outlets and derive_topology functions.""" + G, _ = load_street_network() + G = gu.assign_id(G) + G = gu.double_directed(G) + for ix, (u,v,d) in enumerate(G.edges(data=True)): + d['edge_type'] = 'street' + d['weight'] = ix + + params = parameters.OutletDerivation(river_buffer_distance = 300) + dummy_river1 = sgeom.LineString([(699913.878,5709769.851), + (699932.546,5709882.575)]) + dummy_river2 = sgeom.LineString([(699932.546,5709882.575), + (700011.524,5710060.636)]) + dummy_river3 = sgeom.LineString([(700011.524,5710060.636), + (700103.427,5710169.052)]) + + G.add_edge('river1', 'river2', **{'length' : 10, + 'edge_type' : 'river', + 'id' : 'river1-to-river2', + 'geometry' : dummy_river1}) + G.add_edge('river2', 'river3', **{'length' : 10, + 'edge_type' : 'river', + 'id' : 'river2-to-river3', + 'geometry' : dummy_river2}) + + G.add_edge('river3', 'river4', **{'length' : 10, + 'edge_type' : 'river', + 'id' : 'river3-to-river4', + 'geometry' : dummy_river3}) + + G.nodes['river1']['x'] = 699913.878 + G.nodes['river1']['y'] = 5709769.851 + G.nodes['river2']['x'] = 699932.546 + G.nodes['river2']['y'] = 5709882.575 + G.nodes['river3']['x'] = 700011.524 + G.nodes['river3']['y'] = 5710060.636 + G.nodes['river4']['x'] = 700103.427 + G.nodes['river4']['y'] = 5710169.052 + + # Test outlet derivation + G_ = G.copy() + 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) == 2 + + # Test topo derivation + G_ = gu.derive_topology(G_) + assert len(G_.edges) == 22 + + # Test outlet derivation parameters + G_ = G.copy() + params.outlet_length = 600 + 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_pipe_by_pipe(): + """Test the pipe_by_pipe function.""" + G = load_graph(Path(__file__).parent / 'test_data' / 'graph_topo_derived.json') + for ix, (u,d) in enumerate(G.nodes(data=True)): + d['elevation'] = ix + d['contributing_area'] = ix + + params = parameters.HydraulicDesign() + + G = gu.pipe_by_pipe(G, params) + for u, v, d in G.edges(data=True): + assert 'diameter' in d.keys() + assert d['diameter'] in params.diameters + + for u, d in G.nodes(data=True): + assert 'chamber_floor_elevation' in d.keys() + assert math.isfinite(d['chamber_floor_elevation']) + \ No newline at end of file