Поиск окрестностей (клик) по данным улиц (график)

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-полной проблемой. Обзор этой проблемы (и связанных с ней проблем) можно найти здесь . На SO есть одно описание алгоритма для его решения здесь . Насколько я могу судить, в networkx (или в python соответствующая реализация отсутствует). Я кратко попробовал этот подход, а затем отказался от него - у меня сегодня не до мозга костей для такой работы сегодня. При этом, Я буду назначать вознаграждение любому, кто может посетить эту страницу позднее, и опубликовать протестированную реализацию алгоритма, который находит SSSR в python.

I вместо этого следовали другому подходу, используя тот факт, что граф гарантированно является плоским. Вкратце, вместо того, чтобы рассматривать это как проблему с графом, мы рассматриваем это как проблему сегментации изображения. Сначала мы находим все связанные области на изображении. Затем мы определяем контур вокруг каждой области, преобразуем контуры в координатах изображения обратно в долготы и широты.

С учетом следующих определений импорта и функций:

#!/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://stackoverflow.com/a/35362787/2912349
    # https://stackoverflow.com/a/54334430/2912349
    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.")

Answer 2

Я не совсем уверен, что 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)

Answer 3

У меня нет кода, но я предполагаю, что как только я окажусь на тротуаре, если я продолжу поворачивать вправо в каждом углу, я буду циклически проходить по краям моего блока. Я не знаю библиотек, поэтому я просто поговорю здесь go.

• с вашей точки, go на север, пока не дойдете до улицы
• , поверните направо как сколько сможете и идите по улице
• по следующему углу, найдите всех звезд, выберете тот, который делает наименьший угол с вашей улицей, считая справа.
• идите по этому улица.
• поверните направо и т. д. c.

На самом деле это алгоритм, используемый для выхода из лабиринта: держите правую руку на стене и ходите. Это не работает в случае петель в лабиринте, вы просто l oop вокруг. Но это дает решение вашей проблемы.

Answer 4

Это реализация идеи Хашеми Эмада . Это работает хорошо, если исходная позиция выбрана таким образом, что существует способ сделать шаг против часовой стрелки в крутом круге. Для некоторых ребер, в частности вокруг карты, это невозможно. У меня нет идеи, как выбирать хорошие стартовые позиции или как фильтровать решения - но, может быть, у кого-то еще есть такая.

Рабочий пример (начиная с ребра (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://stackoverflow.com/a/31735642/2912349
    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')