Skip to content

Commit

Permalink
Merge pull request #401 from oturns/lodescodes
Browse files Browse the repository at this point in the history 
allow edge vertices in isochones
  • Loading branch information
@knaaptime
knaaptime committed Mar 5, 2024
2 parents dcf3b27 + 57259c9 commit 7a81b6c
Show file tree
Hide file tree
Showing 3 changed files with 164 additions and 71 deletions.
168 changes: 111 additions & 57 deletions geosnap/analyze/network.py
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
from warnings import warn

import geopandas as gpd
import pandas as pd
from libpysal.cg import alpha_shape_auto
Expand All @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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.
Expand All @@ -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
Expand All @@ -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
-------
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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:
Expand All @@ -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"
Expand All @@ -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)


32 changes: 21 additions & 11 deletions geosnap/io/networkio.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,21 +24,22 @@ 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"
twoway : bool, optional
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
<https://osmnx.readthedocs.io/en/stable/internals-reference.html#osmnx-speed-module>
for more information.
default_speeds : dict, optional
Expand All @@ -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
------
Expand Down Expand Up @@ -89,17 +92,21 @@ 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]],
node_x=n["x"],
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.
Expand Down Expand Up @@ -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
35 changes: 32 additions & 3 deletions geosnap/tests/test_isochrones.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"):
Expand Down Expand Up @@ -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
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

0 comments on commit 7a81b6c

Please sign in to comment.