From a516becb84f4ee1516ab29d5b93f2ebf8075d089 Mon Sep 17 00:00:00 2001 From: Dobson Date: Wed, 31 Jan 2024 12:54:45 +0000 Subject: [PATCH] Edit geosptail_uitilties -Load test data rather than download -test derive_rc -minor edits to fix derive subs --- swmmanywhere/geospatial_utilities.py | 21 +++---- tests/test_geospatial_utilities.py | 87 +++++++++++++++++++++++----- 2 files changed, 83 insertions(+), 25 deletions(-) diff --git a/swmmanywhere/geospatial_utilities.py b/swmmanywhere/geospatial_utilities.py index f993de47..1d0e3db6 100644 --- a/swmmanywhere/geospatial_utilities.py +++ b/swmmanywhere/geospatial_utilities.py @@ -29,7 +29,6 @@ TransformerFromCRS = lru_cache(pyproj.transformer.Transformer.from_crs) - def get_utm_epsg(x: float, y: float, crs: str | int | pyproj.CRS = 'EPSG:4326', @@ -535,7 +534,7 @@ def derive_subcatchments(G: nx.Graph, fid: Path) -> gpd.GeoDataFrame: Returns: gpd.GeoDataFrame: A GeoDataFrame containing polygons with columns: - 'geometry', 'area', 'id', and 'slope'. + 'geometry', 'area', 'id', 'width', and 'slope'. """ # Initialise pysheds grids grid = pgrid.Grid.from_raster(str(fid)) @@ -572,6 +571,9 @@ def remove_(mp): return remove_zero_area_subareas(mp, 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, @@ -580,9 +582,9 @@ def derive_rc(polys_gdf: gpd.GeoDataFrame, """Derive the RC of each subcatchment. Args: - polys_gdf (gpd.GeoDataFrame): A GeoDataFrame containing polygons with - columns: 'geometry', 'area', and 'id'. - G (nx.Graph): The input graph. + 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. building_footprints (gpd.GeoDataFrame): A GeoDataFrame containing building footprints with a 'geometry' column. @@ -592,9 +594,6 @@ def derive_rc(polys_gdf: gpd.GeoDataFrame, """ polys_gdf = polys_gdf.copy() - # TODO this should probably be in derive_subcatchments - polys_gdf['width'] = polys_gdf['area'].div(np.pi).pow(0.5) - ## Format as swmm type catchments # TODO think harder about lane widths (am I double counting here?) @@ -618,6 +617,8 @@ def derive_rc(polys_gdf: gpd.GeoDataFrame, dissolved_result['impervious_area'] = dissolved_result.geometry.area polys_gdf = pd.merge(polys_gdf, dissolved_result[['id','impervious_area']], - on = 'id') + on = 'id', + how='left').fillna(0) polys_gdf['rc'] = polys_gdf['impervious_area'] / polys_gdf['area'] * 100 - return polys_gdf \ No newline at end of file + return polys_gdf + diff --git a/tests/test_geospatial_utilities.py b/tests/test_geospatial_utilities.py index 3631abea..81bfa6df 100644 --- a/tests/test_geospatial_utilities.py +++ b/tests/test_geospatial_utilities.py @@ -7,6 +7,7 @@ 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,9 +15,15 @@ from shapely import geometry as sgeom from swmmanywhere import geospatial_utilities as go -from swmmanywhere.prepare_data import download_elevation, download_street +from swmmanywhere import graph_utilities as ge +def load_street_network(): + """Load a street network.""" + bbox = (-0.11643,51.50309,-0.11169,51.50549) + G = ge.load_graph(Path(__file__).parent / 'test_data' / 'street_graph.json') + return G, bbox + def test_interp_with_nans(): """Test the interp_interp_with_nans function.""" # Define a simple grid and values @@ -217,22 +224,15 @@ def test_burn_shape_in_raster(): def test_derive_subcatchments(): """Test the derive_subcatchments function.""" - bbox = (-0.11643,51.50309,-0.11169,51.50549) - G = download_street(bbox) - temp_fid = Path('temp.tif') + G, bbox = load_street_network() + elev_fid = Path(__file__).parent / 'test_data' / 'elevation.tif' try: - test_api_key = 'b206e65629ac0e53d599e43438560d28' - download_elevation(temp_fid, - bbox, - test_api_key) + temp_fid_r = Path('temp_reprojected.tif') crs = go.get_utm_epsg(bbox[0], bbox[1]) go.reproject_raster(crs, - temp_fid, + elev_fid, temp_fid_r) - G = go.reproject_graph(G, - 'EPSG:4326', - crs) polys = go.derive_subcatchments(G, temp_fid_r) assert 'slope' in polys.columns @@ -242,7 +242,64 @@ def test_derive_subcatchments(): assert polys.shape[0] > 0 assert polys.dropna().shape == polys.shape finally: - if temp_fid.exists(): - temp_fid.unlink() if temp_fid_r.exists(): - temp_fid_r.unlink() \ No newline at end of file + temp_fid_r.unlink() + +def test_derive_rc(): + """Test the derive_rc function.""" + G, bbox = load_street_network() + crs = go.get_utm_epsg(bbox[0], bbox[1]) + 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