在街道数据（图表）中查找社区（集团）
Finding neighbourhoods (cliques) in street data (a graph)
我正在寻找一种方法来自动将城市中的社区定义为图形上的多边形。
我对社区的定义分为两部分：
- 街区：由多条街道围成的区域，其中街道（边）和交叉路口（节点）的数量至少为三个（三角形）。
- 一个社区：对于任何给定的街区，所有与该街区直接相邻的街区和街区本身。
示例请参见此图：
例如B4 是由 7 个节点和连接它们的 6 条边定义的块。正如此处的大多数示例，其他块由 4 个节点和连接它们的 4 条边定义。此外，B1 的 neighborhood 包括 B2（反之亦然），而 B2 还包括 B3.
我正在使用 osmnx 从 OSM 获取街道数据。
- 如何使用 osmnx 和 networkx 遍历图形以找到定义每个块的节点和边？
- 对于每个块，如何找到相邻的块？
我正在努力编写一段代码，将图形和一对坐标（纬度、经度）作为输入，识别相关块和 returns 该块的多边形和邻域定义如上。
这是用于制作地图的代码：
import osmnx as ox
import networkx as nx
import matplotlib.pyplot as plt
G = ox.graph_from_address('Nørrebrogade 20, Copenhagen Municipality',
network_type='all',
distance=500)
以及我尝试寻找具有不同节点数和度数的派系。
def plot_cliques(graph, number_of_nodes, degree):
ug = ox.save_load.get_undirected(graph)
cliques = nx.find_cliques(ug)
cliques_nodes = [clq for clq in cliques if len(clq) >= number_of_nodes]
print("{} cliques with more than {} nodes.".format(len(cliques_nodes), number_of_nodes))
nodes = set(n for clq in cliques_nodes for n in clq)
h = ug.subgraph(nodes)
deg = nx.degree(h)
nodes_degree = [n for n in nodes if deg[n] >= degree]
k = h.subgraph(nodes_degree)
nx.draw(k, node_size=5)
可能相关的理论：
我不完全确定 cycle_basis 会为您提供您寻找的社区，但如果可以，从中获取社区图是一件简单的事情：
import osmnx as ox
import networkx as nx
import matplotlib.pyplot as plt
G = ox.graph_from_address('Nørrebrogade 20, Copenhagen Municipality',
network_type='all',
distance=500)
H = nx.Graph(G) # make a simple undirected graph from G
cycles = nx.cycles.cycle_basis(H) # I think a cycle basis should get all the neighborhoods, except
# we'll need to filter the cycles that are too small.
cycles = [set(cycle) for cycle in cycles if len(cycle) > 2] # Turn the lists into sets for next loop.
# We can create a new graph where the nodes are neighborhoods and two neighborhoods are connected if
# they are adjacent:
I = nx.Graph()
for i, n in enumerate(cycles):
for j, m in enumerate(cycles[i + 1:], start=i + 1):
if not n.isdisjoint(m):
I.add_edge(i, j)
使用图表查找城市街区非常重要。
基本上，这相当于找到最小环的最小集合 (SSSR)，这是一个 NP 完全问题。
可以找到对这个问题（和相关问题）的评论 here。
在 SO 上，有一个算法的描述来解决它 。
据我所知，networkx（或 python 中没有相应的实现）。
我简要地尝试了这种方法，然后放弃了它——我的大脑今天还不能适应这种工作。
话虽这么说，我将奖励日后可能访问此页面的任何人，并 post 一个经过测试的算法实现，该算法在 python 中找到 SSSR。
我转而采用不同的方法，利用图形保证是平面的这一事实。
简而言之，我们不是将其视为图形问题，而是将其视为图像分割问题。
首先，我们找到图像中所有连接的区域。然后我们确定每个区域周围的轮廓，
将图像坐标中的轮廓转换回经度和纬度。
给定以下导入和函数定义：
#!/usr/bin/env python
# coding: utf-8
"""
Find house blocks in osmnx graphs.
"""
import numpy as np
import osmnx as ox
import networkx as nx
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
from skimage.measure import label, find_contours, points_in_poly
from skimage.color import label2rgb
ox.config(log_console=True, use_cache=True)
def k_core(G, k):
H = nx.Graph(G, as_view=True)
H.remove_edges_from(nx.selfloop_edges(H))
core_nodes = nx.k_core(H, k)
H = H.subgraph(core_nodes)
return G.subgraph(core_nodes)
def plot2img(fig):
# remove margins
fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
# convert to image
#
#
canvas = FigureCanvas(fig)
canvas.draw()
img_as_string, (width, height) = canvas.print_to_buffer()
as_rgba = np.fromstring(img_as_string, dtype='uint8').reshape((height, width, 4))
return as_rgba[:,:,:3]
加载数据。如果反复测试，请缓存导入——否则您的帐户可能会被禁止。
从这里的经验来说。
G = ox.graph_from_address('Nørrebrogade 20, Copenhagen Municipality',
network_type='all', distance=500)
G_projected = ox.project_graph(G)
ox.save_graphml(G_projected, filename='network.graphml')
# G = ox.load_graphml('network.graphml')
修剪不能成为循环一部分的节点和边。
这一步不是绝对必要的，但会产生更好的轮廓。
H = k_core(G, 2)
fig1, ax1 = ox.plot_graph(H, node_size=0, edge_color='k', edge_linewidth=1)
将绘图转换为图像并找到连接区域：
img = plot2img(fig1)
label_image = label(img > 128)
image_label_overlay = label2rgb(label_image[:,:,0], image=img[:,:,0])
fig, ax = plt.subplots(1,1)
ax.imshow(image_label_overlay)
对于每个标记区域，找到轮廓并将轮廓像素坐标转换回数据坐标。
# using a large region here as an example;
# however we could also loop over all unique labels, i.e.
# for ii in np.unique(labels.ravel()):
ii = np.argsort(np.bincount(label_image.ravel()))[-5]
mask = (label_image[:,:,0] == ii)
contours = find_contours(mask.astype(np.float), 0.5)
# Select the largest contiguous contour
contour = sorted(contours, key=lambda x: len(x))[-1]
# display the image and plot the contour;
# this allows us to transform the contour coordinates back to the original data cordinates
fig2, ax2 = plt.subplots()
ax2.imshow(mask, interpolation='nearest', cmap='gray')
ax2.autoscale(enable=False)
ax2.step(contour.T[1], contour.T[0], linewidth=2, c='r')
plt.close(fig2)
# first column indexes rows in images, second column indexes columns;
# therefor we need to swap contour array to get xy values
contour = np.fliplr(contour)
pixel_to_data = ax2.transData + ax2.transAxes.inverted() + ax1.transAxes + ax1.transData.inverted()
transformed_contour = pixel_to_data.transform(contour)
transformed_contour_path = Path(transformed_contour, closed=True)
patch = PathPatch(transformed_contour_path, facecolor='red')
ax1.add_patch(patch)
确定原始图中落在轮廓内（或轮廓上）的所有点。
x = G.nodes.data('x')
y = G.nodes.data('y')
xy = np.array([(x[node], y[node]) for node in G.nodes])
eps = (xy.max(axis=0) - xy.min(axis=0)).mean() / 100
is_inside = transformed_contour_path.contains_points(xy, radius=-eps)
nodes_inside_block = [node for node, flag in zip(G.nodes, is_inside) if flag]
node_size = [50 if node in nodes_inside_block else 0 for node in G.nodes]
node_color = ['r' if node in nodes_inside_block else 'k' for node in G.nodes]
fig3, ax3 = ox.plot_graph(G, node_color=node_color, node_size=node_size)
判断两个街区是否相邻非常容易。只需检查他们是否共享一个节点：
if set(nodes_inside_block_1) & set(nodes_inside_block_2): # empty set evaluates to False
print("Blocks are neighbors.")
我没有密码，但我想一旦我在人行道上，如果我在每个拐角处一直向右转，我就会骑车穿过街区的边缘。我不知道图书馆所以我就在这里谈谈算法。
- 从你所在的位置向北走，直到你到达一条街道
- 尽可能右转，在街上行走
- 在下一个拐角处，找到所有的街道，从右边数起，选择与您的街道成最小角度的街道。
- 走在那条街上。
- 右转等
它实际上是一种用于退出迷宫的算法：将右手放在墙上然后走路。它在迷宫中循环的情况下不起作用，你只是循环。但它可以解决您的问题。
这是 的实现。只要选择了起始位置，它就可以很好地工作，这样就可以在一个紧密的圆圈中逆时针前进。对于某些边缘，特别是地图外部周围，这是不可能的。我不知道如何 select 好的起始位置，或者如何过滤解决方案——但也许其他人有。
工作示例（从边缘 (1204573687, 4555480822) 开始）：
此方法不起作用的示例（从边缘 (1286684278, 5818325197) 开始）：
代码
#!/usr/bin/env python
# coding: utf-8
"""
Find house blocks in osmnx graphs.
"""
import numpy as np
import networkx as nx
import osmnx as ox
import matplotlib.pyplot as plt; plt.ion()
from matplotlib.path import Path
from matplotlib.patches import PathPatch
ox.config(log_console=True, use_cache=True)
def k_core(G, k):
H = nx.Graph(G, as_view=True)
H.remove_edges_from(nx.selfloop_edges(H))
core_nodes = nx.k_core(H, k)
H = H.subgraph(core_nodes)
return G.subgraph(core_nodes)
def get_vector(G, n1, n2):
dx = np.diff([G.nodes.data()[n]['x'] for n in (n1, n2)])
dy = np.diff([G.nodes.data()[n]['y'] for n in (n1, n2)])
return np.array([dx, dy])
def angle_between(v1, v2):
#
ang1 = np.arctan2(*v1[::-1])
ang2 = np.arctan2(*v2[::-1])
return (ang1 - ang2) % (2 * np.pi)
def step_counterclockwise(G, edge, path):
start, stop = edge
v1 = get_vector(G, stop, start)
neighbors = set(G.neighbors(stop))
candidates = list(set(neighbors) - set([start]))
if not candidates:
raise Exception("Ran into a dead end!")
else:
angles = np.zeros_like(candidates, dtype=float)
for ii, neighbor in enumerate(candidates):
v2 = get_vector(G, stop, neighbor)
angles[ii] = angle_between(v1, v2)
next_node = candidates[np.argmin(angles)]
if next_node in path:
# next_node might not be the same as the first node in path;
# therefor, we backtrack until we end back at next_node
closed_path = [next_node]
for node in path[::-1]:
closed_path.append(node)
if node == next_node:
break
return closed_path[::-1] # reverse to have counterclockwise path
else:
path.append(next_node)
return step_counterclockwise(G, (stop, next_node), path)
def get_city_block_patch(G, boundary_nodes, *args, **kwargs):
xy = []
for node in boundary_nodes:
x = G.nodes.data()[node]['x']
y = G.nodes.data()[node]['y']
xy.append((x, y))
path = Path(xy, closed=True)
return PathPatch(path, *args, **kwargs)
if __name__ == '__main__':
# --------------------------------------------------------------------------------
# load data
# # DO CACHE RESULTS -- otherwise you can get banned for repeatedly querying the same address
# G = ox.graph_from_address('Nørrebrogade 20, Copenhagen Municipality',
# network_type='all', distance=500)
# G_projected = ox.project_graph(G)
# ox.save_graphml(G_projected, filename='network.graphml')
G = ox.load_graphml('network.graphml')
# --------------------------------------------------------------------------------
# prune nodes and edges that should/can not be part of a cycle;
# this also reduces the chance of running into a dead end when stepping counterclockwise
H = k_core(G, 2)
# --------------------------------------------------------------------------------
# pick an edge and step counterclockwise until you complete a circle
# random edge
total_edges = len(H.edges)
idx = np.random.choice(total_edges)
start, stop, _ = list(H.edges)[idx]
# good edge
# start, stop = 1204573687, 4555480822
# bad edge
# start, stop = 1286684278, 5818325197
steps = step_counterclockwise(H, (start, stop), [start, stop])
# --------------------------------------------------------------------------------
# plot
patch = get_city_block_patch(G, steps, facecolor='red', edgecolor='red', zorder=-1)
node_size = [100 if node in steps else 20 for node in G.nodes]
node_color = ['crimson' if node in steps else 'black' for node in G.nodes]
fig1, ax1 = ox.plot_graph(G, node_size=node_size, node_color=node_color, edge_color='k', edge_linewidth=1)
ax1.add_patch(patch)
fig1.savefig('city_block.png')
我明白这个问题现在有点老了，但我有一个相对简单的替代方法 - 不过它确实需要暂时离开 networkx。
创建块
获取投影图：
import osmnx as ox
import geopandas as gpd
from shapely.ops import polygonize
G = ox.graph_from_address('Nørrebrogade 20, Copenhagen Municipality',
network_type='all',
dist=500)
G_projected = ox.project_graph(G)
将图形转换为无向图 - 这会删除会导致后续多边形化失败的重复边：
G_undirected = G_projected.to_undirected()
仅将边缘提取到 GeoPandas GeoDataFrame 中：
G_edges_as_gdf = ox.graph_to_gdfs(G_undirected, nodes=False, edges=True)
在边缘上使用 shapely.ops 中的 polygonize 创建块面，然后将它们用作新 GeoDataFrame 中的几何体：
block_faces = list(polygonize(G_edges_as_gdf['geometry']))
blocks = gpd.GeoDataFrame(geometry=block_faces)
绘制结果：
ax = G_edges_as_gdf.plot(figsize=(10,10), color='red', zorder=0)
blocks.plot(ax=ax, facecolor='gainsboro', edgecolor='k', linewidth=2, alpha=0.5, zorder=1)
blocks created from line fragments by shapely.ops.polygonize()
寻找邻居
PySAL 的空间权重在这方面做得很好，请参阅 https://pysal.org/notebooks/lib/libpysal/weights.html 了解更多信息。 libpysal 可以从 conda-forge 安装。
在这里，我们使用 Rook 权重来识别与原始问题中共享边的块。皇后权重还包括那些只共享一个节点的权重（即在街口相遇）
from libpysal.weights import Rook # Queen, KNN also available
w_rook = Rook.from_dataframe(blocks)
空间权重仅针对邻居，因此我们需要附加原始块（此处以索引号 18 为例）：
block = 18
neighbors = w_rook.neighbors[block]
neighbors.append(block)
neighbors
然后我们可以使用 neighbors 作为过滤器进行绘图：
ax = blocks.plot(figsize=(10,10), facecolor='gainsboro', edgecolor='black')
blocks[blocks.index.isin(neighbors)].plot(ax=ax, color='red', alpha=0.5)
neighboring blocks highlighted on top of all blocks
备注
- 这不能解决由多条道路中心线引起的狭长多边形，并且解决由步行街和人行道引起的歧义可能具有挑战性 - 大概步行街是街区之间的合法划分，但人行道可能不是.
- 根据我的记忆，我过去需要在对它们进行多边形化之前进一步分割从边缘创建的 LineStrings - 这个例子似乎没有必要。
- 我省略了通过某种空间连接等将新几何体链接回节点 ID 的最后一步。
我正在寻找一种方法来自动将城市中的社区定义为图形上的多边形。
我对社区的定义分为两部分：
- 街区：由多条街道围成的区域，其中街道（边）和交叉路口（节点）的数量至少为三个（三角形）。
- 一个社区：对于任何给定的街区，所有与该街区直接相邻的街区和街区本身。
示例请参见此图：
例如B4 是由 7 个节点和连接它们的 6 条边定义的块。正如此处的大多数示例，其他块由 4 个节点和连接它们的 4 条边定义。此外，B1 的 neighborhood 包括 B2（反之亦然），而 B2 还包括 B3.
我正在使用 osmnx 从 OSM 获取街道数据。
- 如何使用 osmnx 和 networkx 遍历图形以找到定义每个块的节点和边？
- 对于每个块，如何找到相邻的块？
我正在努力编写一段代码，将图形和一对坐标（纬度、经度）作为输入，识别相关块和 returns 该块的多边形和邻域定义如上。
这是用于制作地图的代码：
import osmnx as ox
import networkx as nx
import matplotlib.pyplot as plt
G = ox.graph_from_address('Nørrebrogade 20, Copenhagen Municipality',
network_type='all',
distance=500)
以及我尝试寻找具有不同节点数和度数的派系。
def plot_cliques(graph, number_of_nodes, degree):
ug = ox.save_load.get_undirected(graph)
cliques = nx.find_cliques(ug)
cliques_nodes = [clq for clq in cliques if len(clq) >= number_of_nodes]
print("{} cliques with more than {} nodes.".format(len(cliques_nodes), number_of_nodes))
nodes = set(n for clq in cliques_nodes for n in clq)
h = ug.subgraph(nodes)
deg = nx.degree(h)
nodes_degree = [n for n in nodes if deg[n] >= degree]
k = h.subgraph(nodes_degree)
nx.draw(k, node_size=5)
可能相关的理论：
我不完全确定 cycle_basis 会为您提供您寻找的社区，但如果可以，从中获取社区图是一件简单的事情：
import osmnx as ox
import networkx as nx
import matplotlib.pyplot as plt
G = ox.graph_from_address('Nørrebrogade 20, Copenhagen Municipality',
network_type='all',
distance=500)
H = nx.Graph(G) # make a simple undirected graph from G
cycles = nx.cycles.cycle_basis(H) # I think a cycle basis should get all the neighborhoods, except
# we'll need to filter the cycles that are too small.
cycles = [set(cycle) for cycle in cycles if len(cycle) > 2] # Turn the lists into sets for next loop.
# We can create a new graph where the nodes are neighborhoods and two neighborhoods are connected if
# they are adjacent:
I = nx.Graph()
for i, n in enumerate(cycles):
for j, m in enumerate(cycles[i + 1:], start=i + 1):
if not n.isdisjoint(m):
I.add_edge(i, j)
使用图表查找城市街区非常重要。
基本上，这相当于找到最小环的最小集合 (SSSR)，这是一个 NP 完全问题。
可以找到对这个问题（和相关问题）的评论 here。
在 SO 上，有一个算法的描述来解决它 networkx（或 python 中没有相应的实现）。
我简要地尝试了这种方法，然后放弃了它——我的大脑今天还不能适应这种工作。
话虽这么说，我将奖励日后可能访问此页面的任何人，并 post 一个经过测试的算法实现，该算法在 python 中找到 SSSR。
我转而采用不同的方法，利用图形保证是平面的这一事实。 简而言之，我们不是将其视为图形问题，而是将其视为图像分割问题。 首先，我们找到图像中所有连接的区域。然后我们确定每个区域周围的轮廓， 将图像坐标中的轮廓转换回经度和纬度。
给定以下导入和函数定义：
#!/usr/bin/env python
# coding: utf-8
"""
Find house blocks in osmnx graphs.
"""
import numpy as np
import osmnx as ox
import networkx as nx
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
from skimage.measure import label, find_contours, points_in_poly
from skimage.color import label2rgb
ox.config(log_console=True, use_cache=True)
def k_core(G, k):
H = nx.Graph(G, as_view=True)
H.remove_edges_from(nx.selfloop_edges(H))
core_nodes = nx.k_core(H, k)
H = H.subgraph(core_nodes)
return G.subgraph(core_nodes)
def plot2img(fig):
# remove margins
fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
# convert to image
#
#
canvas = FigureCanvas(fig)
canvas.draw()
img_as_string, (width, height) = canvas.print_to_buffer()
as_rgba = np.fromstring(img_as_string, dtype='uint8').reshape((height, width, 4))
return as_rgba[:,:,:3]
加载数据。如果反复测试，请缓存导入——否则您的帐户可能会被禁止。 从这里的经验来说。
G = ox.graph_from_address('Nørrebrogade 20, Copenhagen Municipality',
network_type='all', distance=500)
G_projected = ox.project_graph(G)
ox.save_graphml(G_projected, filename='network.graphml')
# G = ox.load_graphml('network.graphml')
修剪不能成为循环一部分的节点和边。 这一步不是绝对必要的，但会产生更好的轮廓。
H = k_core(G, 2)
fig1, ax1 = ox.plot_graph(H, node_size=0, edge_color='k', edge_linewidth=1)
将绘图转换为图像并找到连接区域：
img = plot2img(fig1)
label_image = label(img > 128)
image_label_overlay = label2rgb(label_image[:,:,0], image=img[:,:,0])
fig, ax = plt.subplots(1,1)
ax.imshow(image_label_overlay)
对于每个标记区域，找到轮廓并将轮廓像素坐标转换回数据坐标。
# using a large region here as an example;
# however we could also loop over all unique labels, i.e.
# for ii in np.unique(labels.ravel()):
ii = np.argsort(np.bincount(label_image.ravel()))[-5]
mask = (label_image[:,:,0] == ii)
contours = find_contours(mask.astype(np.float), 0.5)
# Select the largest contiguous contour
contour = sorted(contours, key=lambda x: len(x))[-1]
# display the image and plot the contour;
# this allows us to transform the contour coordinates back to the original data cordinates
fig2, ax2 = plt.subplots()
ax2.imshow(mask, interpolation='nearest', cmap='gray')
ax2.autoscale(enable=False)
ax2.step(contour.T[1], contour.T[0], linewidth=2, c='r')
plt.close(fig2)
# first column indexes rows in images, second column indexes columns;
# therefor we need to swap contour array to get xy values
contour = np.fliplr(contour)
pixel_to_data = ax2.transData + ax2.transAxes.inverted() + ax1.transAxes + ax1.transData.inverted()
transformed_contour = pixel_to_data.transform(contour)
transformed_contour_path = Path(transformed_contour, closed=True)
patch = PathPatch(transformed_contour_path, facecolor='red')
ax1.add_patch(patch)
确定原始图中落在轮廓内（或轮廓上）的所有点。
x = G.nodes.data('x')
y = G.nodes.data('y')
xy = np.array([(x[node], y[node]) for node in G.nodes])
eps = (xy.max(axis=0) - xy.min(axis=0)).mean() / 100
is_inside = transformed_contour_path.contains_points(xy, radius=-eps)
nodes_inside_block = [node for node, flag in zip(G.nodes, is_inside) if flag]
node_size = [50 if node in nodes_inside_block else 0 for node in G.nodes]
node_color = ['r' if node in nodes_inside_block else 'k' for node in G.nodes]
fig3, ax3 = ox.plot_graph(G, node_color=node_color, node_size=node_size)
判断两个街区是否相邻非常容易。只需检查他们是否共享一个节点：
if set(nodes_inside_block_1) & set(nodes_inside_block_2): # empty set evaluates to False
print("Blocks are neighbors.")
我没有密码，但我想一旦我在人行道上，如果我在每个拐角处一直向右转，我就会骑车穿过街区的边缘。我不知道图书馆所以我就在这里谈谈算法。
- 从你所在的位置向北走，直到你到达一条街道
- 尽可能右转，在街上行走
- 在下一个拐角处，找到所有的街道，从右边数起，选择与您的街道成最小角度的街道。
- 走在那条街上。
- 右转等
它实际上是一种用于退出迷宫的算法：将右手放在墙上然后走路。它在迷宫中循环的情况下不起作用，你只是循环。但它可以解决您的问题。
这是
工作示例（从边缘 (1204573687, 4555480822) 开始）：
此方法不起作用的示例（从边缘 (1286684278, 5818325197) 开始）：
代码
#!/usr/bin/env python
# coding: utf-8
"""
Find house blocks in osmnx graphs.
"""
import numpy as np
import networkx as nx
import osmnx as ox
import matplotlib.pyplot as plt; plt.ion()
from matplotlib.path import Path
from matplotlib.patches import PathPatch
ox.config(log_console=True, use_cache=True)
def k_core(G, k):
H = nx.Graph(G, as_view=True)
H.remove_edges_from(nx.selfloop_edges(H))
core_nodes = nx.k_core(H, k)
H = H.subgraph(core_nodes)
return G.subgraph(core_nodes)
def get_vector(G, n1, n2):
dx = np.diff([G.nodes.data()[n]['x'] for n in (n1, n2)])
dy = np.diff([G.nodes.data()[n]['y'] for n in (n1, n2)])
return np.array([dx, dy])
def angle_between(v1, v2):
#
ang1 = np.arctan2(*v1[::-1])
ang2 = np.arctan2(*v2[::-1])
return (ang1 - ang2) % (2 * np.pi)
def step_counterclockwise(G, edge, path):
start, stop = edge
v1 = get_vector(G, stop, start)
neighbors = set(G.neighbors(stop))
candidates = list(set(neighbors) - set([start]))
if not candidates:
raise Exception("Ran into a dead end!")
else:
angles = np.zeros_like(candidates, dtype=float)
for ii, neighbor in enumerate(candidates):
v2 = get_vector(G, stop, neighbor)
angles[ii] = angle_between(v1, v2)
next_node = candidates[np.argmin(angles)]
if next_node in path:
# next_node might not be the same as the first node in path;
# therefor, we backtrack until we end back at next_node
closed_path = [next_node]
for node in path[::-1]:
closed_path.append(node)
if node == next_node:
break
return closed_path[::-1] # reverse to have counterclockwise path
else:
path.append(next_node)
return step_counterclockwise(G, (stop, next_node), path)
def get_city_block_patch(G, boundary_nodes, *args, **kwargs):
xy = []
for node in boundary_nodes:
x = G.nodes.data()[node]['x']
y = G.nodes.data()[node]['y']
xy.append((x, y))
path = Path(xy, closed=True)
return PathPatch(path, *args, **kwargs)
if __name__ == '__main__':
# --------------------------------------------------------------------------------
# load data
# # DO CACHE RESULTS -- otherwise you can get banned for repeatedly querying the same address
# G = ox.graph_from_address('Nørrebrogade 20, Copenhagen Municipality',
# network_type='all', distance=500)
# G_projected = ox.project_graph(G)
# ox.save_graphml(G_projected, filename='network.graphml')
G = ox.load_graphml('network.graphml')
# --------------------------------------------------------------------------------
# prune nodes and edges that should/can not be part of a cycle;
# this also reduces the chance of running into a dead end when stepping counterclockwise
H = k_core(G, 2)
# --------------------------------------------------------------------------------
# pick an edge and step counterclockwise until you complete a circle
# random edge
total_edges = len(H.edges)
idx = np.random.choice(total_edges)
start, stop, _ = list(H.edges)[idx]
# good edge
# start, stop = 1204573687, 4555480822
# bad edge
# start, stop = 1286684278, 5818325197
steps = step_counterclockwise(H, (start, stop), [start, stop])
# --------------------------------------------------------------------------------
# plot
patch = get_city_block_patch(G, steps, facecolor='red', edgecolor='red', zorder=-1)
node_size = [100 if node in steps else 20 for node in G.nodes]
node_color = ['crimson' if node in steps else 'black' for node in G.nodes]
fig1, ax1 = ox.plot_graph(G, node_size=node_size, node_color=node_color, edge_color='k', edge_linewidth=1)
ax1.add_patch(patch)
fig1.savefig('city_block.png')
我明白这个问题现在有点老了，但我有一个相对简单的替代方法 - 不过它确实需要暂时离开 networkx。
创建块
获取投影图：
import osmnx as ox
import geopandas as gpd
from shapely.ops import polygonize
G = ox.graph_from_address('Nørrebrogade 20, Copenhagen Municipality',
network_type='all',
dist=500)
G_projected = ox.project_graph(G)
将图形转换为无向图 - 这会删除会导致后续多边形化失败的重复边：
G_undirected = G_projected.to_undirected()
仅将边缘提取到 GeoPandas GeoDataFrame 中：
G_edges_as_gdf = ox.graph_to_gdfs(G_undirected, nodes=False, edges=True)
在边缘上使用 shapely.ops 中的 polygonize 创建块面，然后将它们用作新 GeoDataFrame 中的几何体：
block_faces = list(polygonize(G_edges_as_gdf['geometry']))
blocks = gpd.GeoDataFrame(geometry=block_faces)
绘制结果：
ax = G_edges_as_gdf.plot(figsize=(10,10), color='red', zorder=0)
blocks.plot(ax=ax, facecolor='gainsboro', edgecolor='k', linewidth=2, alpha=0.5, zorder=1)
blocks created from line fragments by shapely.ops.polygonize()
寻找邻居
PySAL 的空间权重在这方面做得很好，请参阅 https://pysal.org/notebooks/lib/libpysal/weights.html 了解更多信息。 libpysal 可以从 conda-forge 安装。
在这里，我们使用 Rook 权重来识别与原始问题中共享边的块。皇后权重还包括那些只共享一个节点的权重（即在街口相遇）
from libpysal.weights import Rook # Queen, KNN also available
w_rook = Rook.from_dataframe(blocks)
空间权重仅针对邻居，因此我们需要附加原始块（此处以索引号 18 为例）：
block = 18
neighbors = w_rook.neighbors[block]
neighbors.append(block)
neighbors
然后我们可以使用 neighbors 作为过滤器进行绘图：
ax = blocks.plot(figsize=(10,10), facecolor='gainsboro', edgecolor='black')
blocks[blocks.index.isin(neighbors)].plot(ax=ax, color='red', alpha=0.5)
neighboring blocks highlighted on top of all blocks
备注
- 这不能解决由多条道路中心线引起的狭长多边形，并且解决由步行街和人行道引起的歧义可能具有挑战性 - 大概步行街是街区之间的合法划分，但人行道可能不是.
- 根据我的记忆，我过去需要在对它们进行多边形化之前进一步分割从边缘创建的 LineStrings - 这个例子似乎没有必要。
- 我省略了通过某种空间连接等将新几何体链接回节点 ID 的最后一步。