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