Search code examples
networkxgeopandasfoliumshapelyosmnx

How do I Plot multiple Isochrones Polygons using Python OSMnx library?

I'm trying to build isochrones for some schools, to understand the accessibility. My mode of travelling is "walk" and I wanted to create a list of at least 2 travel times, 30 minutes and 1 hour.

Here below is my code:

# import required libraries
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
import networkx as nx
import osmnx as ox
from descartes import PolygonPatch
from shapely.geometry import Point, LineString, Polygon
from IPython.display import IFrame
import folium
import itertools # will be used to perform feature joining

# %matplotlib inline
ox.__version__

# function to get isochrones
def get_isochrone(lon, lat, walk_time=[30,60], speed=4.5):
    loc = (lat, lon)
    G = ox.graph_from_point(loc, simplify=True, network_type='walk')
    # Create nodes geodataframe from Graph network (G)
    gdf_nodes = ox.graph_to_gdfs(G, edges=False)
    x, y = gdf_nodes['geometry'].unary_union.centroid.xy
    center_node = ox.get_nearest_node(G, (y[0], x[0]))
    meters_per_minute = speed * 1000 / 60 #km per hour to m per minute
    for u, v, k, data in G.edges(data=True, keys=True):
        data['time'] = data['length'] / meters_per_minute
    # get one color for each isochrone
    iso_colors = ox.plot.get_colors(n=len(walk_time), cmap="plasma", start=0, return_hex=True)
    # color the nodes according to isochrone then plot the street network
    node_colors = {}
    for walks_time, color in zip(sorted(walk_time, reverse=True), iso_colors):
        subgraph = nx.ego_graph(G, center_node, radius=walks_time, distance="time")
        for node in subgraph.nodes():
            node_colors[node] = color
        
    nc = [node_colors[node] if node in node_colors else "none" for node in G.nodes()]
    ns = [15 if node in node_colors else 0 for node in G.nodes()] 
    # make isochrone polygons
    isochrone_polys = []
    for trip_times in sorted(walk_time, reverse=True):
        subgraph = nx.ego_graph(G, center_node, radius=trip_times, distance="time")
        node_points = [Point(data['x'], data['y']) for node, data in subgraph.nodes(data=True)]
        polys = gpd.GeoSeries(node_points).unary_union.convex_hull
        isochrone_polys.append(polys)
    # adding color
    for polygon, fc in zip(isochrone_polys, iso_colors):
        patch = PolygonPatch(polygon, fc=fc, ec="none", alpha=0.6, zorder=-1)
        # isochrone_polys.add_patch(patch) 
    return isochrone_polys

When I call the get_isochrone function:

# calling the function with the coordinates
map = coordinates.apply(lambda x: get_isochrone(x.longitude, x.latitude), axis=1)

What is returned is a list of polygons with each school point having polygons equal to the number of items in the travel times list. I however noticed that the polygons returned for each point are exactly the same.

Just to verify: This script shows that they are the same:

for i in map:
    print(i[0])
    print(i[1])
    print('\n')

Here are the results: Image showing what the above script returns

As you can see, each school point returns two polygons with exactly the same coordinates, for different travel times. I have reviewed the python function several times, even changed the travel time to one item in the travel time list like 100, and I still get exactly the same coordinates. I may have missed something and I'll appreciate any help. Thanks!

In addition, I also realized that the add_patch() method in the function adds color to the plot_graph() only. How can I add the colors to the polygons in the Folium map?

Solution
    • it appears that you have used this example: https://github.com/gboeing/osmnx-examples/blob/main/notebooks/13-isolines-isochrones.ipynb
    • if all you want is polygons, then there is no need to do PolygonPatch
    • you have refactored incorrectly, for multiple walk times. You only need to generate edges once
    • pulling it all together:
      1. source some schools
      2. get_isochrone() refactored to your use case. Have changed do it returns a dict that has index and name of point / school being investigated.
      3. generate a geopandas data frame of isochrones
      4. visualise it

    data sourcing

    import osmnx as ox
import pandas as pd
import warnings
import networkx as nx
import geopandas as gpd
from shapely.geometry import Point

warnings.simplefilter(action="ignore", category=FutureWarning)
warnings.simplefilter(action="ignore", category=PendingDeprecationWarning)
ox.config(use_cache=True, log_console=False)
# get some cities
cities = ["Hereford"]  # ["Hereford", "Worcester", "Gloucester"]
cities = ox.geocode_to_gdf([{"city": c, "country": "UK"} for c in cities])

# get some schools
tags = {"amenity": "school"}
schools = pd.concat(
    [
        ox.geometries.geometries_from_polygon(r["geometry"], tags)
        for i, r in cities.iterrows()
    ]
)
schools = (
    schools.loc["way"].dropna(axis=1, thresh=len(schools) / 4).drop(columns=["nodes"])
)
# change polygon to point
schools["geometry"] = schools.to_crs(schools.estimate_utm_crs())[
    "geometry"
].centroid.to_crs(schools.crs)

    get_isochrone()

    
# function to get isochrones
def get_isochrone(
    lon, lat, walk_times=[15, 30], speed=4.5, name=None, point_index=None
):
    loc = (lat, lon)
    G = ox.graph_from_point(loc, simplify=True, network_type="walk")
    gdf_nodes = ox.graph_to_gdfs(G, edges=False)
    center_node = ox.distance.nearest_nodes(G, lon, lat)

    meters_per_minute = speed * 1000 / 60  # km per hour to m per minute
    for u, v, k, data in G.edges(data=True, keys=True):
        data["time"] = data["length"] / meters_per_minute
    polys = []
    for walk_time in walk_times:
        subgraph = nx.ego_graph(G, center_node, radius=walk_time, distance="time")
        node_points = [
            Point(data["x"], data["y"]) for node, data in subgraph.nodes(data=True)
        ]
        polys.append(gpd.GeoSeries(node_points).unary_union.convex_hull)
    info = {}
    if name:
        info["name"] = [name for t in walk_times]
    if point_index:
        info["point_index"] = [point_index for t in walk_times]
    return {**{"geometry": polys, "time": walk_times}, **info}

    integration

    WT = [5, 10, 15]
SCHOOLS = 5

# build geopandas data frame of isochrone polygons for each school
isochrones = pd.concat(
    [
        gpd.GeoDataFrame(
            get_isochrone(
                r["geometry"].x,
                r["geometry"].y,
                name=r["name"],
                point_index=i,
                walk_times=WT,
            ),
            crs=schools.crs,
        )
        for i, r in schools.head(SCHOOLS).iterrows()
    ]
)

    visualise

    warnings.filterwarnings("ignore")

gdf = isochrones.set_index(["time", "point_index"]).copy()
# remove shorter walk time from longer walk time polygon to make folium work better
for idx in range(len(WT)-1,0,-1):
    gdf.loc[WT[idx], "geometry"] = (
        gdf.loc[WT[idx]]
        .apply(
            lambda r: r["geometry"].symmetric_difference(
                gdf.loc[(WT[idx-1], r.name), "geometry"]
            ),
            axis=1,
        )
        .values
    )

m = gdf.reset_index().explore(column="time", height=300, width=500, scheme="boxplot")
schools.head(SCHOOLS).explore(m=m, marker_kwds={"radius": 3, "color": "red"})

    enter image description here

    merge overlapping polygons

    import matplotlib.cm as cm
import matplotlib.colors as colors
import folium

# merge overlapping polygons
# https://gis.stackexchange.com/questions/334459/how-to-dissolve-overlapping-polygons-using-geopandas
mergedpolys = gpd.GeoDataFrame(
    geometry=isochrones.groupby("time")["geometry"]
    .agg(lambda g: g.unary_union)
    .apply(lambda g: [g] if isinstance(g, shapely.geometry.Polygon) else g.geoms)
    .explode(),
    crs=isochrones.crs,
)

# visualize merged polygons
m = None
for i, wt in enumerate(WT[::-1]):
    m = mergedpolys.loc[[wt]].explore(
        m=m,
        color=colors.to_hex(cm.get_cmap("tab20b", len(WT))(i)),
        name=wt,
        height=300,
        width=500,
    )

m = schools.head(SCHOOLS).explore(
    m=m, marker_kwds={"radius": 3, "color": "red"}, name="schools"
)
folium.LayerControl().add_to(m)

m

    enter image description here