Localizando bairros (cliques) nos dados das ruas (um gráfico)

Estou procurando uma maneira de definir automaticamente bairros nas cidades como polígonos em um gráfico.

Minha definição de bairro tem duas partes:

1. Um quarteirão : uma área dividida entre várias ruas, em que o número de ruas (arestas) e cruzamentos (nós) é no mínimo três (um triângulo).
2. Uma vizinhança : para qualquer bloco, todos os blocos diretamente adjacentes a esse bloco e o próprio bloco.

Veja esta ilustração para um exemplo:

Por exemplo, B4 é um bloco definido por 7 nós e 6 arestas conectando-os. Como a maioria dos exemplos aqui, os outros blocos são definidos por 4 nós e 4 arestas conectando-os. Além disso, o bairro de B1 inclui B2 (e vice-versa) enquanto B2 também inclui B3 .

Estou usando o osmnx para obter dados de ruas do OSM.

1. Usando osmnx e networkx, como posso percorrer um gráfico para encontrar os nós e arestas que definem cada bloco?
2. Para cada bloco, como posso encontrar os blocos adjacentes?

Estou trabalhando em direção a um pedaço de código que usa um gráfico e um par de coordenadas (latitude, longitude) como entrada, identifica o bloco relevante e retorna o polígono para esse bloco e a vizinhança, conforme definido acima.

Aqui está o código usado para fazer o mapa:

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)

e minha tentativa de encontrar panelinhas com número diferente de nós e graus.

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)

Teoria que pode ser relevante:

Enumerando todos os ciclos em um gráfico não direcionado

Encontrar blocos urbanos usando o gráfico é surpreendentemente não trivial. Basicamente, isso equivale a encontrar o menor conjunto de anéis menores (SSSR), que é um problema completo do NP. Uma revisão deste problema (e problemas relacionados) pode ser encontrada aqui . No SO, há uma descrição de um algoritmo para resolvê-lo aqui . Até onde eu sei, não existe uma implementação correspondente em networkx(ou em python). Tentei essa abordagem brevemente e depois a abandonei - meu cérebro não está preparado para esse tipo de trabalho hoje. Dito isto, concederei uma recompensa a qualquer pessoa que possa visitar esta página posteriormente e publicarei uma implementação testada de um algoritmo que encontra o SSSR em python.

Em vez disso, segui uma abordagem diferente, aproveitando o fato de que o gráfico é garantido como plano. Resumidamente, em vez de tratar isso como um problema gráfico, tratamos isso como um problema de segmentação de imagem. Primeiro, encontramos todas as regiões conectadas na imagem. Em seguida, determinamos o contorno em torno de cada região, transformamos os contornos das coordenadas da imagem em longitudes e latitudes.

Dadas as seguintes importações e definições de função:

#!/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]

Carregue os dados. Faça o cache das importações, se estiver testando isso repetidamente - caso contrário, sua conta poderá ser banida. Falando por experiência aqui.

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')

Remova nós e arestas que não podem fazer parte de um ciclo. Esta etapa não é estritamente necessária, mas resulta em contornos mais agradáveis.

H = k_core(G, 2)
fig1, ax1 = ox.plot_graph(H, node_size=0, edge_color='k', edge_linewidth=1)

Converta plotagem em imagem e encontre regiões conectadas:

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)

Para cada região rotulada, localize o contorno e converta as coordenadas de pixel de contorno novamente em coordenadas de dados.

# 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)

Determine todos os pontos no gráfico original que caem dentro (ou sobre) o contorno.

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)

Descobrir se dois quarteirões são vizinhos é bastante fácil. Basta verificar se eles compartilham um nó:

if set(nodes_inside_block_1) & set(nodes_inside_block_2): # empty set evaluates to False
    print("Blocks are neighbors.")

Não tenho certeza absoluta de que cycle_basisirá fornecer os bairros que você procura, mas, se isso acontecer, é simples obter o gráfico do bairro:

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)

Não tenho código, mas acho que quando estiver na calçada, se continuar virando à direita em cada esquina, percorrerei as bordas do meu bloco. Eu não conheço as bibliotecas, então vou falar algo aqui.

• do seu ponto, vá para o norte até chegar a uma rua
• vire à direita o máximo que puder e ande na rua
• na próxima esquina, encontre todos os ângulos, escolha o que fizer o menor ângulo com a sua rua contando à direita.
• andar nessa rua.
• vire à direita, etc.

Na verdade, é um algoritmo a ser usado para sair de um labirinto: mantenha a mão direita na parede e ande. Não funciona no caso de loops no labirinto, você apenas dá voltas. Mas dá uma solução para o seu problema.

Esta é uma implementação da ideia de Hashemi Emad . Funciona bem desde que a posição inicial seja escolhida de forma que exista uma maneira de avançar no sentido anti-horário em um círculo apertado. Para algumas arestas, principalmente na parte externa do mapa, isso não é possível. Não faço ideia de como selecionar boas posições iniciais ou como filtrar soluções - mas talvez alguém tenha uma.

Exemplo de trabalho (começando com a borda (1204573687, 4555480822)):

Exemplo, em que essa abordagem não funciona (começando com a borda (1286684278, 5818325197)):

Código

#!/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')