Note

Interactive online version: Binder badge

Simplification of roundabouts#

Analysis using street network can be sensitive to the actual geometry representation. For example, one roundabout can be typically represented as 4 or more nodes, while it is topologically only one. Hence measurements attempting to capture topology (e.g. centrality) or node density will be skewed.

The following code exemplifies a simple process for road network simplification of roundabouts using momepy.roundabout_simplification().

For this example we will fetch some data from OpenStreetMap (using osmnx). However, most other data sources should also work as long as their topology has been corrected and can be stored in a GeoDataFrame.

[1]:
import geopandas as gpd
import matplotlib.pyplot as plt
import momepy as mm
import osmnx as ox

Load data#

Using osmnx download a neighborhood and reproject to its local CRS.

Two things are required to achieve good results: 1. Reproject the network to a projected CRS (in meters) 1. Transform the graph to an undirected graph: - This helps to remove overlapping LineStrings once moving to GeoDataFrame.

[2]:
place = "Chamberi, Madrid"

G = ox.graph_from_place(
    place, network_type="drive", simplify=True, buffer_dist=200
)
G_projected = ox.project_graph(G)

edges = ox.graph_to_gdfs(
    ox.get_undirected(G_projected),  # prevents some (semi)duplicate geometries
    nodes=False,
    edges=True,
    node_geometry=False,
    fill_edge_geometry=True,
)

edges.head(3)
/var/folders/2f/fhks6w_d0k556plcv3rfmshw0000gn/T/ipykernel_93551/3487061367.py:3: FutureWarning: The buffer_dist argument has been deprecated and will be removed in the v2.0.0 release. Buffer your query area directly, if desired. See the OSMnx v2 migration guide: https://github.com/gboeing/osmnx/issues/1123
  G = ox.graph_from_place(
/Users/martin/miniforge3/envs/momepy/lib/python3.11/site-packages/osmnx/graph.py:381: FutureWarning: The buffer_dist argument has been deprecated and will be removed in the v2.0.0 release. Buffer your results directly, if desired. See the OSMnx v2 migration guide: https://github.com/gboeing/osmnx/issues/1123
  gdf_place = geocoder.geocode_to_gdf(
/Users/martin/miniforge3/envs/momepy/lib/python3.11/site-packages/osmnx/graph.py:392: FutureWarning: The 'unary_union' attribute is deprecated, use the 'union_all()' method instead.
  polygon = gdf_place["geometry"].unary_union
/var/folders/2f/fhks6w_d0k556plcv3rfmshw0000gn/T/ipykernel_93551/3487061367.py:9: FutureWarning: The `get_undirected` function is deprecated and will be removed in the v2.0.0 release. Replace it with `convert.to_undirected` instead. See the OSMnx v2 migration guide: https://github.com/gboeing/osmnx/issues/1123
  ox.get_undirected(G_projected),  # prevents some (semi)duplicate geometries
[2]:
osmid oneway name highway reversed length from to geometry maxspeed lanes bridge width access tunnel junction
u v key
21990706 21990729 0 390868364 True Calle de Mejía Lequerica residential False 60.664 21990706 21990729 LINESTRING (440818.925 4475309.73, 440776.64 4... NaN NaN NaN NaN NaN NaN NaN
21990729 25906109 0 44884571 True Calle de Mejía Lequerica residential False 35.206 21990729 25906109 LINESTRING (440776.64 4475353.246, 440750.521 ... NaN NaN NaN NaN NaN NaN NaN
25906107 0 163743121 True Calle de la Beneficencia residential False 65.736 21990729 25906107 LINESTRING (440776.64 4475353.246, 440825.335 ... NaN NaN NaN NaN NaN NaN NaN
[3]:
ax = edges.plot(figsize=(8, 12), color="red")
ax.set_axis_off()
../../_images/user_guide_preprocessing_roundabout_simplification_4_0.png

Default behavior#

[4]:
edges_output = mm.roundabout_simplification(edges)
[5]:
ax = edges.plot(figsize=(8, 12), color="red")
edges_output.plot(ax=ax, figsize=(8, 12), color="black")

ax.set_axis_off()
../../_images/user_guide_preprocessing_roundabout_simplification_7_0.png

Using the output column 'simplification_group' we can count the number of roundabouts that were simplified:

[6]:
edges_output.simplification_group.nunique()
[6]:
7

In total 7 roundabouts were simplified, and most of their incoming edges were also considered during the process.

Let’s look at them individually.

[7]:
simplification_groups = edges_output.simplification_group.unique()[1:]

fig, axs = plt.subplots(1, len(simplification_groups), figsize=(25, 7))
axs = axs.flatten()

for i in simplification_groups:
    mask = edges_output.simplification_group == i
    minx, miny, maxx, maxy = edges_output[mask].geometry.total_bounds

    edges.plot(ax=axs[i], color="red")
    edges_output.plot(ax=axs[i], color="black")

    axs[i].set_xlim(minx, maxx)
    axs[i].set_ylim(miny, maxy)
    axs[i].set_axis_off()

plt.show()
../../_images/user_guide_preprocessing_roundabout_simplification_11_0.png

It’s worth investigating why some roundabouts are still missing from simplification. Since the selection of roundabouts is based on the resulting polygons after polygonizing the road network let’s use the circular compactness (measurable as momepy.CircularCompactness) parameter of the resulting polygons as the main selection attribute.

[8]:
from shapely.ops import polygonize

polys = gpd.GeoDataFrame(
    geometry=list(polygonize(edges.geometry)), crs=edges.crs
)
circom_serie = mm.CircularCompactness(polys).series
polys.loc[:, "circom"] = circom_serie

The circom_threshold parameter (default = 0.7) establishes the limit at which roundabouts should be selected.

[9]:
circom_threshold = 0.7
mask = circom_serie > circom_threshold
colors = ["green" if m else "#dbdada" for m in mask]

fig, axs = plt.subplots(1, 2, figsize=(20, 8))

polys.plot(ax=axs[0], column="circom", cmap="PiYG")
axs[0].set_axis_off()
axs[0].title.set_text("Circular Compactness values")

polys.plot(ax=axs[1], color=colors)
axs[1].set_axis_off()
axs[1].title.set_text("Polygons above 'circom_threshold' ")
../../_images/user_guide_preprocessing_roundabout_simplification_15_0.png

Let’s see the following two false negatives:

[10]:
false_negatives = [[312, 314, 304], 162]


fig, axs = plt.subplots(1, 2, figsize=(20, 6))


edges.plot(ax=axs[0], color="black")
polys.plot(ax=axs[0], column="circom", cmap="PiYG")
minx, miny, maxx, maxy = polys.loc[false_negatives[0]].geometry.total_bounds
axs[0].set_xlim(minx - 75, maxx + 75)
axs[0].set_ylim(miny - 75, maxy + 75)
axs[0].set_axis_off()
axs[0].title.set_text("Roundabout indexes:\n" + str(false_negatives[0]))

edges.plot(ax=axs[1], color="black")
polys.plot(ax=axs[1], column="circom", cmap="PiYG")
minx, miny, maxx, maxy = polys.loc[false_negatives[1]].geometry.bounds
label_str = str(round(polys.loc[false_negatives[1]].circom, 3))

axs[1].set_xlim(minx - 75, maxx + 75)
axs[1].set_ylim(miny - 75, maxy + 75)
axs[1].set_axis_off()
axs[1].title.set_text("Roundabout index: 162\n 'circom': " + label_str)
../../_images/user_guide_preprocessing_roundabout_simplification_17_0.png

The above helps to exemplify two main reasons why some roundabouts are not being selected: 1. They are crossed by other roads; which breaks them into multiple polygons - eg. index: [314, 312, 304] 1. They don’t meet the circom_threshold - eg. index: [162] –> (0.676)

Adapting parameters#

Finally, let’s also explore what the outcome would look like if we lower the circom_threshold and include_adjacent is set to False

[11]:
edges_output_65 = mm.roundabout_simplification(
    edges, circom_threshold=0.65, include_adjacent=False
)

m = edges.explore(color="black", tiles="CartoDB positron")
edges_output.explore(m=m, color="red")
edges_output_65.explore(m=m, color="blue")
m
[11]:
Make this Notebook Trusted to load map: File -> Trust Notebook