在街道数据（一个图）中寻找邻域（团伙）

Question

我正在寻找一种自动将城市中的社区定义为图形多边形的方法。

我的社区定义有两部分：

1. 一个街区：由多条街道包围的区域，其中街道数量（边缘）和路口（节点）最少为三个（三角形）。
2. 邻里：对于任何给定的街区，所有直接相邻的街区以及该街区本身。

参见以下示例说明：

例如，B4 是由7个节点和6条连接它们的边界定义的街区。就像这里大多数示例一样，其他街区由4个节点和4条连接它们的边缘定义。此外，B1 的邻里包括B2（反之亦然），而B2 也包括B3。

我使用osmnx从OSM获取街道数据。

1. 使用 osmnx 和 networkx，如何遍历图形以查找定义每个街区的节点和边缘？
2. 对于每个街区，如何找到相邻的街区？

我正在编写一段代码，该代码接受图形和一对坐标（纬度、经度）作为输入，识别相关块并返回多边形和上述邻里。

以下是用于创建地图的代码：

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)

可能相关的理论：

在无向图中枚举所有环

- tmo

---

有趣的问题。你可能想在它上面加上算法标签。看起来，当你弄清楚街区之后，邻里会是更容易的问题。对于邻里，你只需要寻找共享的边缘，对吧？每个街区都将有一个边缘列表... 对于街区，我认为获取每个节点的每个街道选项的基本方向并“保持向右转”（或向左）直到完成电路，或到达死胡同或循环回来并递归地返回会很有帮助。虽然似乎会有一些有趣的特殊情况。 - AirSquid

我认为这个问题与您的问题1非常相似。从链接中可以看出，我对此问题进行了一些工作，它很棘手（结果证明是NP-hard）。然而，我的答案中的启发式方法可能仍然可以给您提供足够好的结果。 - Paul Brodersen

由于您可能认为任何可接受的解决方案都可能是一种启发式方法，因此定义一个测试数据集以验证每种方法可能是个好主意。也就是说，对于您的示例图形，最好有一个机器可读形式的所有块的注释，而不仅仅是图像中的一些示例。 - Paul Brodersen

Answer 1

使用图形查找城市街区令人惊讶地不容易。基本上，这相当于查找最小环的最小集合（SSSR），这是一个NP完全问题。可以在此处找到有关此问题（及相关问题）的评论。在SO上，有一种算法描述来解决它此处。据我所知，networkx中没有相应的实现（或者说在Python中没有）。我尝试了这种方法，但很快就放弃了 - 我今天的大脑无法胜任那种工作。话虽如此，我会授予任何在以后访问此页面并发布在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
    # https://dev59.com/jFsW5IYBdhLWcg3wI0RH#35362787
    # https://dev59.com/jFsW5IYBdhLWcg3wI0RH#54334430
    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.")

- Paul Brodersen

Answer 2

我很感激这个问题，虽然现在有点过时了，但我有一个相对简单的替代方法——它需要暂时离开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)

通过shapely.ops.polygonize()从线段碎片创建的区块

查找邻居

PySAL通过其空间权重在这方面表现出色，有关详细信息，请参见https://pysal.org/notebooks/lib/libpysal/weights.html。可以从conda-forge安装libpysal。

在这里，我们使用Rook权重来识别共享边缘的区块，就像原始问题中一样。Queen权重也将包括仅共享节点（即在街道交叉口处相遇）的那些区块。

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)

所有区块上方突出显示相邻区块

注释

• 这并不能解决由多条道路中心线引起的小块多边形，解决由行人街道和人行道引起的歧义可能会很具有挑战性。可以假定行人街道是合法的区块分界线，但人行道可能不是。
• 据我记忆，在对边缘创建的LineStrings进行面化之前，我过去需要进一步分割它们 - 在这个示例中似乎不需要。
• 我省略了将新几何图形与某种空间连接等方式链接回节点ID的最后一步。

- Nick Bristow

---

这里有一个针对整个英国实施此策略的实现链接。 - Peter Prescott

Answer 3

我不完全确定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)

- salt-die

---

嗨，salt-die，欢迎来到SO并感谢你的贡献。当执行nx.Graph(G)时，我会丢失许多信息（定向性和多重图类型），因此我很难验证您的答案，因为我似乎无法将新图I与我的原始图G相关联。 - tmo

@tmo 只是路过：在那种情况下，你应该能够使用扩展了 Graph 的 MultiDiGraph 类。 - Théo Rubenach

Answer 4

我没有代码，但是我猜一旦我在人行道上，如果我每到一个拐角就向右转，我会穿过街区的边缘。我不知道库，所以我会谈论算法。

• 从你的位置出发，向北走直到你到达一条街道
• 尽可能向右转并沿着这条街走
• 在下一个拐角，找到所有的街道，选择与你的街道成最小角度的那个（从右数）。
• 沿着那条街走。
• 继续向右转，以此类推。

这实际上是一个迷宫中使用的算法：将右手放在墙上并走。如果迷宫中有环路，则无法解决问题，你会一直绕圈子。但它可以解决你的问题。

- Hashemi Emad

Answer 5

这是Hashemi Emad的想法的实现。只要选择的起始位置存在逆时针步进到一个紧密圆圈的方式，它就可以很好地工作。对于一些边缘，特别是在地图外部，这是不可能的。我不知道如何选择好的起始位置，或者如何过滤解决方案 - 但也许其他人有想法。

工作示例（以边缘（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):
    # https://dev59.com/41wZ5IYBdhLWcg3wLd2z#31735642
    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')

- Paul Brodersen