使用OSMnx在等时线地图上叠加设施点失败求助

问题描述 #

尝试用OSMnx绘制包含设施点的等时线地图，运行代码后仅显示变形的等时线多边形，设施点完全不显示，代码基于OSMnx官方示例编写。

原代码：

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 LineString
from shapely.geometry import Point
from shapely.geometry import Polygon

%matplotlib inline

ox.config(log_console=False, use_cache=True)

# configure the place, network type, trip times, and travel speed
place = "Menzelstr. 46, Duisburg, Germany"
network_type = "walk"
trip_times = [5, 10, 15]  # in minutes
travel_speed = 4.5  # walking speed in km/hour

# download the street network
G = ox.graph_from_address(place, network_type=network_type)

# find the centermost node and then project the graph to UTM
gdf_nodes = ox.graph_to_gdfs(G, edges=False)
x, y = gdf_nodes["geometry"].unary_union.centroid.xy
center_node = ox.distance.nearest_nodes(G, x[0], y[0])
G = ox.project_graph(G)

# add an edge attribute for time in minutes required to traverse each edge
meters_per_minute = travel_speed * 1000 / 60  # km per hour to m per minute
for _, _, _, 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(trip_times), cmap="autumn", start=0, return_hex=True)

# make the isochrone polygons
isochrone_polys = []
for trip_time in sorted(trip_times, reverse=True):
    subgraph = nx.ego_graph(G, center_node, radius=trip_time, distance="time")
    node_points = [Point((data["x"], data["y"])) for node, data in subgraph.nodes(data=True)]
    bounding_poly = gpd.GeoSeries(node_points).unary_union.convex_hull
    isochrone_polys.append(bounding_poly)

# get amenities for place
amenities = ox.geometries_from_address(place, tags={"amenity":True}, dist=1200)

# plot the network, add isochrones as colored descartes polygon patches, then add amenities
fig, ax = ox.plot_graph(
    G, show=False, close=False, edge_color="#999999", edge_alpha=0.2, node_size=0
)
for polygon, fc in zip(isochrone_polys, iso_colors):
    patch = PolygonPatch(polygon, fc=fc, ec="none", alpha=0.6, zorder=-1)
    ax.add_patch(patch)
    
amenities.plot(color="white", markersize=1, ax=ax)
plt.show()

问题分析与修复 #

1. 设施点不显示的修复 #

• 核心问题：设施点数据的坐标系和投影后的路网不匹配。ox.geometries_from_address返回的是WGS84（EPSG:4326）坐标，而ox.project_graph已将路网转为UTM坐标系，直接绘制会导致设施点偏离可视范围。
• 修复步骤：
  1. 获取设施点后，添加坐标转换代码，对齐路网坐标系：

     amenities = amenities.to_crs(G.graph['crs'])
     

  2. 原代码markersize=1过小，调整为markersize=10才能清晰显示。

2. 等时线多边形变形的修复 #

• 核心问题：用convex_hull生成的是凸包多边形，会强制把所有可达节点的最外层连起来，完全忽略路网的实际走向，导致形状严重变形。
• 修复方案（二选一）：
  • 方案一：用OSMnx内置函数生成精准等时线（推荐）
    OSMnx v1.0+提供了ox.isochrone函数，可直接生成准确的等时线多边形，替换原有的多边形生成代码：

    isochrone_polys = []
    for trip_time in sorted(trip_times, reverse=True):
        # 直接生成等时线
        isochrone_poly = ox.isochrone(G, center_node, trip_time, distance='time')
        isochrone_polys.append(isochrone_poly)
    

  • 方案二：用Alpha形状生成贴合路网的多边形
    若使用旧版OSMnx，可通过alpha_shape生成非凸多边形，更贴近实际可达范围：

    # 生成多边形时替换convex_hull为alpha_shape
    bounding_poly = ox.utils_geo.alpha_shape(gpd.GeoSeries(node_points), alpha=50)
    

    注：alpha值可根据区域大小调整，值越小多边形越贴合节点分布。

3. 完整修复后的代码 #

import geopandas as gpd
import matplotlib.pyplot as plt
import networkx as nx
import osmnx as ox
from descartes import PolygonPatch

%matplotlib inline

ox.config(log_console=False, use_cache=True)

# 配置参数
place = "Menzelstr. 46, Duisburg, Germany"
network_type = "walk"
trip_times = [5, 10, 15]  # 分钟
travel_speed = 4.5  # 步行速度 km/h

# 下载路网
G = ox.graph_from_address(place, network_type=network_type)

# 确定中心节点并投影到UTM
gdf_nodes = ox.graph_to_gdfs(G, edges=False)
x, y = gdf_nodes["geometry"].unary_union.centroid.xy
center_node = ox.distance.nearest_nodes(G, x[0], y[0])
G = ox.project_graph(G)

# 给边添加通行时间属性（分钟）
meters_per_minute = travel_speed * 1000 / 60
for _, _, _, data in G.edges(data=True, keys=True):
    data["time"] = data["length"] / meters_per_minute

# 生成等时线颜色
iso_colors = ox.plot.get_colors(n=len(trip_times), cmap="autumn", start=0, return_hex=True)

# 生成精准等时线多边形（使用OSMnx内置函数）
isochrone_polys = []
for trip_time in sorted(trip_times, reverse=True):
    isochrone_poly = ox.isochrone(G, center_node, trip_time, distance='time')
    isochrone_polys.append(isochrone_poly)

# 获取设施点并转换坐标系
amenities = ox.geometries_from_address(place, tags={"amenity":True}, dist=1200)
amenities = amenities.to_crs(G.graph['crs'])  # 对齐路网坐标系

# 绘图
fig, ax = ox.plot_graph(
    G, show=False, close=False, edge_color="#999999", edge_alpha=0.2, node_size=0
)
# 添加等时线
for polygon, fc in zip(isochrone_polys, iso_colors):
    patch = PolygonPatch(polygon.geometry.iloc[0], fc=fc, ec="none", alpha=0.6, zorder=-1)
    ax.add_patch(patch)
    
# 绘制设施点
amenities.plot(color="white", markersize=10, ax=ax)
plt.show()

内容的提问来源于stack exchange，提问作者lea_bei