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

10

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:

insira a descrição da imagem aqui

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

tmo
fonte
Problema interessante. Você pode adicionar a tag de algoritmo a ela. Parece que os bairros seriam o problema mais fácil depois que você descobrir os blocos. Como bairros, tudo o que você procura é uma vantagem compartilhada, correto? E cada bloco terá uma lista de arestas ... Para os blocos, acho que será útil obter a direção cardinal de cada opção de rua em um nó e "continuar girando para a direita" (ou esquerda) até concluir um circuito ou alcançar um beco sem saída ou volta para si mesmo e recua de forma recursiva. Parece que haveria alguns casos interessantes, no entanto.
Jeff H
Eu acho que esta pergunta é muito semelhante ao seu problema não. 1. Como você pode ver no link, trabalhei um pouco no problema, e é um problema geral (acaba por ser difícil para o NP). A heurística na minha resposta, no entanto, ainda pode lhe dar bons resultados.
Paul Brodersen 28/02
Como qualquer solução que você considere aceitável provavelmente também será uma heurística, pode ser uma boa ideia definir um conjunto de dados de teste para validar cada abordagem. Ou seja, para o seu gráfico de exemplo, seria bom ter uma anotação de todos os blocos em formato legível por máquina - e não apenas alguns exemplos em uma imagem.
Paul Brodersen 28/02

Respostas:

3

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)

gráfico podado

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)

plotagem de rótulos de região

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)

traçado do contorno sobreposto no gráfico podado

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)

trama da rede com nós pertencentes ao bloco em vermelho

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.")
Paul Brodersen
fonte
2

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)
morrer de sal
fonte
Olá, morra, seja bem-vindo à SO e obrigado por entrar em contato. Ao fazer nx.Graph(G)isso, estou perdendo muitas informações (tipo de direcionamento e multigraph), por isso estou tendo dificuldades para verificar sua resposta, pois não consigo relacionar o novo gráfico Icom meu gráfico original G.
tmo 24/02
Será um pouco de trabalho preservar as informações geométricas do gráfico original. Vou tentar isso em breve.
salt-die
@tmo passando: você deve poder usar a classe MultiDiGraph (que estende o Graph) nesse caso
Théo Rubenach
1

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.

Hashemi Emad
fonte
Esta é uma ideia muito melhor do que eu tinha. Vou adicionar uma resposta com uma implementação da sua intuição.
Paul Brodersen
0

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

insira a descrição da imagem aqui

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

insira a descrição da imagem aqui

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')
Paul Brodersen
fonte