在街道数据（图表）中查找社区（集团）

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)

可能相关的理论：

枚举无向图中的所有循环

Answer 1

使用该图查找城市街区出奇地非平凡。基本上，这相当于找到最小环的最小集合 (SSSR)，这是一个 NP 完全问题。可以在此处找到对此问题（和相关问题）的评论。就这么，有一个算法的一个描述来解决它在这里。据我所知，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
    # /sf/answers/2475395121/
    # /sf/answers/3803410131/
    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.")