diff --git a/geosnap/analyze/network.py b/geosnap/analyze/network.py index 723c8d56..46322256 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,25 @@ 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): + 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) + 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) @@ -56,7 +72,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) @@ -71,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 + 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. @@ -81,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 @@ -98,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 ------- @@ -112,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 @@ -133,16 +160,26 @@ 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) + alpha = _geom_to_alpha(nodes) elif hull == "shapely": - alpha = _geom_to_hull(nodes.geometry, 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" ) - 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) @@ -161,45 +198,49 @@ 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 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 + 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 +262,31 @@ 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 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.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" @@ -240,10 +296,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) - - diff --git a/geosnap/io/networkio.py b/geosnap/io/networkio.py index bc61130c..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 ------ @@ -89,9 +92,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 +103,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 +135,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.columns: + 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=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..679be8f2 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,33 @@ 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 + +@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) + 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