From caf5b77b58b69ca9a749127efdc88bfc03ab48ac Mon Sep 17 00:00:00 2001 From: eli knaap Date: Mon, 4 Mar 2024 16:20:45 -0800 Subject: [PATCH 1/6] allow edge vertices in isochones --- geosnap/analyze/network.py | 121 +++++++++++++++++++++++-------------- geosnap/io/networkio.py | 15 +++-- 2 files changed, 86 insertions(+), 50 deletions(-) diff --git a/geosnap/analyze/network.py b/geosnap/analyze/network.py index 723c8d56..3aebd3d4 100644 --- a/geosnap/analyze/network.py +++ b/geosnap/analyze/network.py @@ -1,3 +1,5 @@ +from warnings import warn + import geopandas as gpd import pandas as pd from libpysal.cg import alpha_shape_auto @@ -6,11 +8,20 @@ def _geom_to_hull(geom, ratio, allow_holes): - return concave_hull(MultiPoint(geom.tolist()), ratio=ratio, allow_holes=allow_holes) + if isinstance(geom, list): + return concave_hull(MultiPoint(geom), ratio=ratio, allow_holes=allow_holes) + elif isinstance(geom, gpd.GeoDataFrame): + return concave_hull( + MultiPoint(geom.geometry.get_coordinates().values), + ratio=ratio, + allow_holes=allow_holes, + ) + def _geom_to_alpha(geom): return alpha_shape_auto(geom.get_coordinates()[["x", "y"]].values) + def _points_to_poly(df, column, hull="shapely", ratio=0.2, allow_holes=False): if hull == "libpysal": output = df.groupby(column)["geometry"].apply(_geom_to_alpha) @@ -134,9 +145,9 @@ def isochrones_from_id( df = matrix[matrix.cost <= distance] nodes = node_df[node_df.index.isin(df.destination.tolist())] if hull == "libpysal": - alpha = _geom_to_alpha(nodes.geometry) + alpha = _geom_to_alpha(nodes.geometry.tolist()) elif hull == "shapely": - alpha = _geom_to_hull(nodes.geometry, ratio=ratio, allow_holes=allow_holes) + alpha = _geom_to_hull(nodes.geometry.tolist(), ratio=ratio, allow_holes=allow_holes) else: raise ValueError( f"`algorithm must be either 'alpha' or 'hull' but {hull} was passed" @@ -161,45 +172,51 @@ def isochrones_from_gdf( hull="shapely", ratio=0.2, allow_holes=False, + use_edges=True, ): """Create travel isochrones for several origins simultaneously - Parameters - ---------- - origins : geopandas.GeoDataFrame - a geodataframe containing the locations of origin point features - threshold: float - maximum travel distance to define the isochrone, measured in the same - units as edges_df in the pandana.Network object. If the network was - created with pandana this is usually meters; if it was created with - urbanaccess this is usually travel time in minutes. - network : pandana.Network - pandana Network instance for calculating the shortest path isochrone - for each origin feature - network_crs : str, int, pyproj.CRS (optional) - the coordinate system used to store x and y coordinates in the passed - pandana network. If None, the network is assumed to be stored in the - same CRS as the origins geodataframe - reindex : bool - if True, use the dataframe index as the origin and destination IDs - (rather than the node_ids of the pandana.Network). Default is True - hull : str, {'libpysal', 'shapely'} - Which method to generate container polygons (concave hulls) for destination - points. If 'libpysal', use `libpysal.cg.alpha_shape_auto` to create the - concave hull, else if 'shapely', use `shapely.concave_hull`. - Default is libpysal - ratio : float - ratio keyword passed to `shapely.concave_hull`. Only used if - `hull='shapely'`. Default is 0.3 - allow_holes : bool - keyword passed to `shapely.concave_hull` governing whether holes are - allowed in the resulting polygon. Only used if `hull='shapely'`. - Default is False. - - Returns - ------- - GeoPandas.DataFrame - polygon geometries with the isochrones for each origin point feature + Parameters + ---------- + origins : geopandas.GeoDataFrame + a geodataframe containing the locations of origin point features + threshold: float + maximum travel distance to define the isochrone, measured in the same + units as edges_df in the pandana.Network object. If the network was + created with pandana this is usually meters; if it was created with + urbanaccess this is usually travel time in minutes. + network : pandana.Network + pandana Network instance for calculating the shortest path isochrone + for each origin feature + network_crs : str, int, pyproj.CRS (optional) + the coordinate system used to store x and y coordinates in the passed + pandana network. If None, the network is assumed to be stored in the + same CRS as the origins geodataframe + reindex : bool + if True, use the dataframe index as the origin and destination IDs + (rather than the node_ids of the pandana.Network). Default is True + hull : str, {'libpysal', 'shapely'} + Which method to generate container polygons (concave hulls) for destination + points. If 'libpysal', use `libpysal.cg.alpha_shape_auto` to create the + concave hull, else if 'shapely', use `shapely.concave_hull`. + Default is libpysal + ratio : float + ratio keyword passed to `shapely.concave_hull`. Only used if + `hull='shapely'`. Default is 0.3 + allow_holes : bool + keyword passed to `shapely.concave_hull` governing whether holes are + allowed in the resulting polygon. Only used if `hull='shapely'`. + Default is False. + use_edges: bool + If true, use vertices from the Network.edge_df to make the polygon more + accurate by adhering to roadways. Requires that the 'geometry' column be + available on the Network.edges_df, most commonly by using + `geosnap.io.get_network_from_gdf` + + Returns + ------- + GeoPandas.DataFrame + polygon geometries with the isochrones for each origin point feature """ if network_crs is None: @@ -221,16 +238,30 @@ def isochrones_from_gdf( reindex=False, drop_nonorigins=False, ) + if (use_edges) and ("geometry" not in network.edges_df.columns): + warn( + "use_edges is True, but edge geometries are not available " + "in the Network object. Try recreating with " + "geosnap.io.get_network_from_gdf" + ) alphas = [] for origin in matrix.origin.unique(): do = matrix[matrix.origin == origin] - dest_pts = destinations.loc[do["destination"]] + dest_pts = gpd.GeoDataFrame(destinations.loc[do["destination"]]) + if use_edges: + dest_pts = dest_pts.geometry.tolist() + else: + edges = network.edges_df.copy() + roads = edges[ + (edges["to"].isin(dest_pts.index.values)) + & (edges["from"].isin(dest_pts.index.values)) + ] + dest_pts = roads + if hull == "libpysal": - alpha = _geom_to_alpha(dest_pts.geometry) + alpha = _geom_to_alpha(dest_pts) elif hull == "shapely": - alpha = _geom_to_hull( - dest_pts.geometry, ratio=ratio, allow_holes=allow_holes - ) + alpha = _geom_to_hull(dest_pts, ratio=ratio, allow_holes=allow_holes) else: raise ValueError( f"`algorithm must be either 'alpha' or 'hull' but {hull} was passed" @@ -245,5 +276,3 @@ def isochrones_from_gdf( if reindex: df = df.rename(index=mapper) return gpd.GeoDataFrame(df, crs=network_crs) - - diff --git a/geosnap/io/networkio.py b/geosnap/io/networkio.py index bc61130c..eaa1d13a 100644 --- a/geosnap/io/networkio.py +++ b/geosnap/io/networkio.py @@ -89,9 +89,10 @@ def get_network_from_gdf( n, e = ox.utils_graph.graph_to_gdfs(graph) if output_crs is not None: n = _reproject_osm_nodes(n, input_crs=4326, output_crs=output_crs) + e = e.to_crs(output_crs) e = e.reset_index() - return pdna.Network( + net= pdna.Network( edge_from=e["u"], edge_to=e["v"], edge_weights=e[[impedance]], @@ -99,7 +100,10 @@ def get_network_from_gdf( node_y=n["y"], twoway=twoway, ) + # keep the geometries on hand, since we have them already + net.edges_df = gpd.GeoDataFrame(net.edges_df, geometry=e.geometry, crs=output_crs) + return net def project_network(network, output_crs=None, input_crs=4326): """Reproject a pandana.Network object into another coordinate system. @@ -128,14 +132,17 @@ def project_network(network, output_crs=None, input_crs=4326): # take original x,y coordinates and convert into geopandas.Series, then reproject nodes = _reproject_osm_nodes(network.nodes_df, input_crs, output_crs) + edges = network.edges_df.copy() + if 'geometry' in edges: + edges = edges.to_crs(output_crs) # reinstantiate the network (needs to rebuild the tree) net = pdna.Network( node_x=nodes["x"], node_y=nodes["y"], - edge_from=network.edges_df["from"], - edge_to=network.edges_df["to"], - edge_weights=network.edges_df[network.impedance_names], + edge_from=edges["from"], + edge_to=network.edges["to"], + edge_weights=network.edges[network.impedance_names], twoway=network._twoway, ) return net From 90bcd72b334b6a2b930b2789bee9adcde0de20b0 Mon Sep 17 00:00:00 2001 From: eli knaap Date: Mon, 4 Mar 2024 16:27:51 -0800 Subject: [PATCH 2/6] typo in conditional logic --- geosnap/analyze/network.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geosnap/analyze/network.py b/geosnap/analyze/network.py index 3aebd3d4..cc788616 100644 --- a/geosnap/analyze/network.py +++ b/geosnap/analyze/network.py @@ -248,7 +248,7 @@ def isochrones_from_gdf( for origin in matrix.origin.unique(): do = matrix[matrix.origin == origin] dest_pts = gpd.GeoDataFrame(destinations.loc[do["destination"]]) - if use_edges: + if use_edges is False: dest_pts = dest_pts.geometry.tolist() else: edges = network.edges_df.copy() From c47bb469f2338327df97ccaa2b39d66437699b95 Mon Sep 17 00:00:00 2001 From: eli knaap Date: Mon, 4 Mar 2024 18:35:07 -0800 Subject: [PATCH 3/6] naming fix --- geosnap/analyze/network.py | 32 ++++++++++++++++++-------------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/geosnap/analyze/network.py b/geosnap/analyze/network.py index cc788616..f7a4e20b 100644 --- a/geosnap/analyze/network.py +++ b/geosnap/analyze/network.py @@ -19,6 +19,9 @@ def _geom_to_hull(geom, ratio, allow_holes): def _geom_to_alpha(geom): + if isinstance(geom, list): + return alpha_shape_auto(gpd.GeoSeries(geom).get_coordinates()[["x", "y"]].values) + return alpha_shape_auto(geom.get_coordinates()[["x", "y"]].values) @@ -67,7 +70,7 @@ def pdna_to_adj(origins, network, threshold, reindex=True, drop_nonorigins=True) # map node ids in the network to index in the gdf mapper = dict(zip(node_ids, origins.index.values)) - namer = {"source": "origin", "distance": "cost"} + namer = {"source": "origin", network.impedance_names[0]: "cost"} adj = network.nodes_in_range(node_ids, threshold) adj = adj.rename(columns=namer) @@ -248,15 +251,16 @@ def isochrones_from_gdf( for origin in matrix.origin.unique(): do = matrix[matrix.origin == origin] dest_pts = gpd.GeoDataFrame(destinations.loc[do["destination"]]) - if use_edges is False: - dest_pts = dest_pts.geometry.tolist() - else: - edges = network.edges_df.copy() - roads = edges[ - (edges["to"].isin(dest_pts.index.values)) - & (edges["from"].isin(dest_pts.index.values)) - ] - dest_pts = roads + if use_edges is True: + if "geometry" not in network.edges_df.columns: + pass + else: + edges = network.edges_df.copy() + roads = edges[ + (edges["to"].isin(dest_pts.index.values)) + & (edges["from"].isin(dest_pts.index.values)) + ] + dest_pts = roads if hull == "libpysal": alpha = _geom_to_alpha(dest_pts) @@ -271,8 +275,8 @@ def isochrones_from_gdf( alpha["distance"] = threshold alpha["origin"] = origin alphas.append(alpha) - df = pd.concat(alphas, ignore_index=True) - df = df.set_index("origin") - if reindex: - df = df.rename(index=mapper) + df = pd.concat(alphas, ignore_index=True) + df = df.set_index("origin") + if reindex: + df = df.rename(index=mapper, errors='raise') return gpd.GeoDataFrame(df, crs=network_crs) From 6d80a19d20464dfa7a3b2cf816e6016f64620818 Mon Sep 17 00:00:00 2001 From: eli knaap Date: Mon, 4 Mar 2024 18:49:01 -0800 Subject: [PATCH 4/6] edges in from_id --- geosnap/analyze/network.py | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/geosnap/analyze/network.py b/geosnap/analyze/network.py index f7a4e20b..f7474862 100644 --- a/geosnap/analyze/network.py +++ b/geosnap/analyze/network.py @@ -85,7 +85,7 @@ def pdna_to_adj(origins, network, threshold, reindex=True, drop_nonorigins=True) def isochrones_from_id( - origin, network, threshold, hull="shapely", ratio=0.2, allow_holes=False + origin, network, threshold, hull="shapely", ratio=0.2, allow_holes=False, use_edges=True ): """Create travel isochrone(s) from a single origin using a pandana network. @@ -147,10 +147,20 @@ def isochrones_from_id( # select the nodes within each threshold distance and take their alpha shape df = matrix[matrix.cost <= distance] nodes = node_df[node_df.index.isin(df.destination.tolist())] + if use_edges is True: + if "geometry" not in network.edges_df.columns: + pass + else: + edges = network.edges_df.copy() + roads = edges[ + (edges["to"].isin(nodes.index.values)) + & (edges["from"].isin(nodes.index.values)) + ] + nodes = roads if hull == "libpysal": - alpha = _geom_to_alpha(nodes.geometry.tolist()) + alpha = _geom_to_alpha(nodes) elif hull == "shapely": - alpha = _geom_to_hull(nodes.geometry.tolist(), ratio=ratio, allow_holes=allow_holes) + alpha = _geom_to_hull(nodes, ratio=ratio, allow_holes=allow_holes) else: raise ValueError( f"`algorithm must be either 'alpha' or 'hull' but {hull} was passed" From 5cd95315a3b9b4baf92cb2edf1a5be49f68034bb Mon Sep 17 00:00:00 2001 From: eli knaap Date: Tue, 5 Mar 2024 12:17:33 -0800 Subject: [PATCH 5/6] add tests --- geosnap/analyze/network.py | 35 +++++++++++++++++++++----------- geosnap/io/networkio.py | 25 +++++++++++++---------- geosnap/tests/test_isochrones.py | 31 +++++++++++++++++++++++++--- 3 files changed, 65 insertions(+), 26 deletions(-) diff --git a/geosnap/analyze/network.py b/geosnap/analyze/network.py index f7474862..46322256 100644 --- a/geosnap/analyze/network.py +++ b/geosnap/analyze/network.py @@ -20,7 +20,9 @@ def _geom_to_hull(geom, ratio, allow_holes): def _geom_to_alpha(geom): if isinstance(geom, list): - return alpha_shape_auto(gpd.GeoSeries(geom).get_coordinates()[["x", "y"]].values) + return alpha_shape_auto( + gpd.GeoSeries(geom).get_coordinates()[["x", "y"]].values + ) return alpha_shape_auto(geom.get_coordinates()[["x", "y"]].values) @@ -85,7 +87,14 @@ def pdna_to_adj(origins, network, threshold, reindex=True, drop_nonorigins=True) def isochrones_from_id( - origin, network, threshold, hull="shapely", ratio=0.2, allow_holes=False, use_edges=True + origin, + network, + threshold, + hull="shapely", + ratio=0.2, + allow_holes=False, + use_edges=True, + network_crs=4326 ): """Create travel isochrone(s) from a single origin using a pandana network. @@ -95,11 +104,12 @@ def isochrones_from_id( A single or list of node id(s) from a `pandana.Network.nodes_df` to serve as isochrone origins network : pandana.Network - A pandana network object + A pandana network object (preferably created with + geosnap.io.get_network_from_gdf) threshold : int or list A single or list of threshold distances for which isochrones will be - computed. These are in the - same units as edges from the pandana.Network.edge_df + computed. These are in the same impedance units as edges from the + pandana.Network.edge_df hull : str, {'libpysal', 'shapely'} Which method to generate container polygons (concave hulls) for destination points. If 'libpysal', use `libpysal.cg.alpha_shape_auto` to create the @@ -112,6 +122,9 @@ def isochrones_from_id( keyword passed to `shapely.concave_hull` governing whether holes are allowed in the resulting polygon. Only used if `algorithm='hull'`. Default is False. + network_crs : int or pyproj.CRS + which coordinate system the pandana.Network node coordinates are stored in. + Default is 4326 Returns ------- @@ -126,7 +139,7 @@ def isochrones_from_id( node_df = gpd.GeoDataFrame( network.nodes_df, geometry=gpd.points_from_xy(network.nodes_df.x, network.nodes_df.y), - crs=4326, + crs=network_crs, ) maxdist = max(threshold) if isinstance(threshold, list) else threshold @@ -166,7 +179,7 @@ def isochrones_from_id( f"`algorithm must be either 'alpha' or 'hull' but {hull} was passed" ) - alpha = gpd.GeoDataFrame(geometry=pd.Series(alpha), crs=4326) + alpha = gpd.GeoDataFrame(geometry=pd.Series(alpha), crs=network_crs) alpha["distance"] = distance dfs.append(alpha) @@ -194,10 +207,8 @@ def isochrones_from_gdf( origins : geopandas.GeoDataFrame a geodataframe containing the locations of origin point features threshold: float - maximum travel distance to define the isochrone, measured in the same - units as edges_df in the pandana.Network object. If the network was - created with pandana this is usually meters; if it was created with - urbanaccess this is usually travel time in minutes. + maximum travel cost to define the isochrone, measured in the same + impedance units as edges_df in the pandana.Network object. network : pandana.Network pandana Network instance for calculating the shortest path isochrone for each origin feature @@ -288,5 +299,5 @@ def isochrones_from_gdf( df = pd.concat(alphas, ignore_index=True) df = df.set_index("origin") if reindex: - df = df.rename(index=mapper, errors='raise') + df = df.rename(index=mapper, errors="raise") return gpd.GeoDataFrame(df, crs=network_crs) diff --git a/geosnap/io/networkio.py b/geosnap/io/networkio.py index eaa1d13a..537f62f7 100644 --- a/geosnap/io/networkio.py +++ b/geosnap/io/networkio.py @@ -24,9 +24,9 @@ def get_network_from_gdf( Parameters ---------- gdf : geopandas.GeoDataFrame - dataframe covering the study area of interest; this should be stored in lat/long - (epsg 4326). Note the first step is to take the unary union of this dataframe, - which is expensive, so large dataframes may be time-consuming + dataframe covering the study area of interest; Note the first step is to take + the unary union of this dataframe, which is expensive, so large dataframes may + be time-consuming. The network will inherit the CRS from this dataframe network_type : str, {"all_private", "all", "bike", "drive", "drive_service", "walk"} the type of network to collect from OSM (passed to `osmnx.graph_from_polygon`) by default "walk" @@ -34,11 +34,12 @@ def get_network_from_gdf( Whether to treat the pandana.Network as directed or undirected. For a directed network, use `twoway=False` (which is the default). For an undirected network (e.g. a walk network) where travel can flow in both directions, the network is more - efficient when twoway=True. This has implications for drive networks or multimodal - networks where impedance is different depending on travel direction. + efficient when twoway=True but forces the impedance to be equal in both + directions. This has implications for auto or multimodal + networks where impedance is generally different depending on travel direction. add_travel_times : bool, default=False whether to use posted travel times from OSM as the impedance measure (rather - than network-distance). Speeds are based on max drive speeds, see + than network-distance). Speeds are based on max posted drive speeds, see for more information. default_speeds : dict, optional @@ -50,7 +51,9 @@ def get_network_from_gdf( ------- pandana.Network a pandana.Network object with node coordinates stored in the same system as the - input geodataframe + input geodataframe. If add_travel_times is True, the network impedance + is travel time measured in seconds (assuming automobile travel speeds); else + the impedance is travel distance measured in meters Raises ------ @@ -92,7 +95,7 @@ def get_network_from_gdf( e = e.to_crs(output_crs) e = e.reset_index() - net= pdna.Network( + net = pdna.Network( edge_from=e["u"], edge_to=e["v"], edge_weights=e[[impedance]], @@ -133,7 +136,7 @@ def project_network(network, output_crs=None, input_crs=4326): # take original x,y coordinates and convert into geopandas.Series, then reproject nodes = _reproject_osm_nodes(network.nodes_df, input_crs, output_crs) edges = network.edges_df.copy() - if 'geometry' in edges: + if "geometry" in edges.columns: edges = edges.to_crs(output_crs) # reinstantiate the network (needs to rebuild the tree) @@ -141,8 +144,8 @@ def project_network(network, output_crs=None, input_crs=4326): node_x=nodes["x"], node_y=nodes["y"], edge_from=edges["from"], - edge_to=network.edges["to"], - edge_weights=network.edges[network.impedance_names], + edge_to=edges["to"], + edge_weights=edges[[network.impedance_names[0]]], twoway=network._twoway, ) return net diff --git a/geosnap/tests/test_isochrones.py b/geosnap/tests/test_isochrones.py index 4ff2db66..0d87f875 100644 --- a/geosnap/tests/test_isochrones.py +++ b/geosnap/tests/test_isochrones.py @@ -3,14 +3,14 @@ isochrones_from_gdf, ) from geosnap import DataStore -from geosnap.io import get_acs, get_network_from_gdf +from geosnap.io import get_acs, get_network_from_gdf, project_network import pandana as pdna import geopandas as gpd import os import pytest import sys from numpy.testing import assert_almost_equal - +import osmnx as ox def get_data(): if not os.path.exists("./41740.h5"): @@ -83,4 +83,29 @@ def test_network_constructor(): tracts = get_acs(DataStore(), county_fips='48301', level='tract', years=2015) walk_net = get_network_from_gdf(tracts) # this will grow depending on the size of the OSM network when tested... - assert walk_net.edges_df.shape[0] > 6000 \ No newline at end of file + assert walk_net.edges_df.shape[0] > 6000 + +def test_isos_with_edges(): + tracts = get_acs(DataStore(), county_fips='48301', level='tract', years=2015) + walk_net = get_network_from_gdf(tracts) + type(walk_net) + facilities = ox.features.features_from_polygon( + tracts.unary_union, {"amenity": "fuel"} +) + #facilities = facilities[facilities.geometry.type == "Point"] + alpha = isochrones_from_gdf( + facilities, network=walk_net, threshold=2000, use_edges=True +) + print(alpha.area.round(8)) + # this will grow depending on the size of the OSM network when tested... + assert alpha.area.round(8).iloc[0] == 0.00026001 + +def test_project_network(): + tracts = get_acs(DataStore(), county_fips='48301', level='tract', years=2015) + walk_net = get_network_from_gdf(tracts) + # this will grow depending on the size of the OSM network when tested... + tracts = tracts.to_crs(tracts.estimate_utm_crs()) + walk_net = project_network(walk_net, output_crs=tracts.crs) + nodes = walk_net.get_node_ids(tracts.centroid.x, tracts.centroid.y) + print(nodes) + assert nodes[0] == 7876436325 \ No newline at end of file From 57259c9357bce6d7e87ad2022c9a1b99f1c7e401 Mon Sep 17 00:00:00 2001 From: eli knaap Date: Tue, 5 Mar 2024 12:45:06 -0800 Subject: [PATCH 6/6] skip test on win --- geosnap/tests/test_isochrones.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/geosnap/tests/test_isochrones.py b/geosnap/tests/test_isochrones.py index 0d87f875..679be8f2 100644 --- a/geosnap/tests/test_isochrones.py +++ b/geosnap/tests/test_isochrones.py @@ -85,6 +85,10 @@ def test_network_constructor(): # this will grow depending on the size of the OSM network when tested... assert walk_net.edges_df.shape[0] > 6000 +@pytest.mark.skipif( + sys.platform.startswith("win"), + reason="skipping test on windows because of dtype issue", +) def test_isos_with_edges(): tracts = get_acs(DataStore(), county_fips='48301', level='tract', years=2015) walk_net = get_network_from_gdf(tracts) @@ -99,7 +103,7 @@ def test_isos_with_edges(): print(alpha.area.round(8)) # this will grow depending on the size of the OSM network when tested... assert alpha.area.round(8).iloc[0] == 0.00026001 - + def test_project_network(): tracts = get_acs(DataStore(), county_fips='48301', level='tract', years=2015) walk_net = get_network_from_gdf(tracts)