Encontrar barrios (camarillas) en datos de calles (un gráfico)

Estoy buscando una manera de definir automáticamente los barrios de las ciudades como polígonos en un gráfico.

Mi definición de vecindario tiene dos partes:

1. Un bloque : un área cerrada entre varias calles, donde la cantidad de calles (bordes) e intersecciones (nodos) es un mínimo de tres (un triángulo).
2. Un vecindario : para cualquier bloque dado, todos los bloques directamente adyacentes a ese bloque y al bloque mismo.

Vea esta ilustración para un ejemplo:

Por ejemplo, B4 es un bloque definido por 7 nodos y 6 aristas que los conectan. Como la mayoría de los ejemplos aquí, los otros bloques están definidos por 4 nodos y 4 bordes que los conectan. Además, el vecindario de B1 incluye B2 (y viceversa) mientras que B2 también incluye B3 .

Estoy usando osmnx para obtener datos de la calle de OSM.

1. Usando osmnx y networkx, ¿cómo puedo atravesar un gráfico para encontrar los nodos y bordes que definen cada bloque?
2. Para cada bloque, ¿cómo puedo encontrar los bloques adyacentes?

Estoy trabajando en un código que toma un gráfico y un par de coordenadas (latitud, longitud) como entrada, identifica el bloque relevante y devuelve el polígono para ese bloque y el vecindario como se definió anteriormente.

Aquí está el código utilizado para hacer el 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)

y mi intento de encontrar camarillas con diferente número de nodos y grados.

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)

Teoría que podría ser relevante:

Enumerar todos los ciclos en un gráfico no dirigido

Encontrar bloques de ciudades usando el gráfico es sorprendentemente no trivial. Básicamente, esto equivale a encontrar el conjunto más pequeño de anillos más pequeños (SSSR), que es un problema NP-completo. Una revisión de este problema (y problemas relacionados) se puede encontrar aquí . En SO, hay una descripción de un algoritmo para resolverlo aquí . Por lo que puedo decir, no hay una implementación correspondiente en networkx(o en Python para el caso). Probé este enfoque brevemente y luego lo abandoné: mi cerebro no está preparado para ese tipo de trabajo hoy. Dicho esto, otorgaré una recompensa a cualquiera que visite esta página más adelante y publique una implementación probada de un algoritmo que encuentre el SSSR en Python.

En cambio, he seguido un enfoque diferente, aprovechando el hecho de que se garantiza que el gráfico sea plano. Brevemente, en lugar de tratar esto como un problema gráfico, tratamos esto como un problema de segmentación de imagen. Primero, encontramos todas las regiones conectadas en la imagen. Luego determinamos el contorno alrededor de cada región, transformamos los contornos en las coordenadas de la imagen a longitudes y latitudes.

Dadas las siguientes importaciones y definiciones de funciones:

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

Cargar los datos. Guarde en caché las importaciones, si lo prueba repetidamente, de lo contrario su cuenta puede ser bloqueada. Hablando por experiencia aquí.

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

Pode nodos y bordes que no pueden ser parte de un ciclo. Este paso no es estrictamente necesario pero da como resultado contornos más bonitos.

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

Convierta la trama en imagen y encuentre regiones 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ón etiquetada, encuentre el contorno y convierta las coordenadas del píxel del contorno nuevamente en coordenadas de datos.

# 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 los puntos en el gráfico original que caen dentro (o sobre) el 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)

Averiguar si dos bloques son vecinos es bastante fácil. Solo verifique si comparten un nodo:

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

No estoy completamente seguro de que eso cycle_basiste dará los vecindarios que buscas, pero si lo hace, es algo sencillo obtener el gráfico de vecindario:

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)

No tengo un código, pero supongo que una vez que esté en la acera, si sigo girando a la derecha en cada esquina, recorreré los bordes de mi bloque. No conozco las bibliotecas, así que hablaré algo aquí.

• desde tu punto, ve hacia el norte hasta llegar a una calle
• gire a la derecha todo lo que pueda y camine por la calle
• En la siguiente esquina, encuentre todos los steets, elija el que haga el ángulo más pequeño con su calle contando desde la derecha.
• Camina por esa calle.
• gire a la derecha, etc.

En realidad, es un algoritmo para salir de un laberinto: mantén la mano derecha en la pared y camina. No funciona en caso de bucles en el laberinto, solo tienes que recorrer. Pero le da una solución a su problema.

Esta es una implementación de la idea de Hashemi Emad . Funciona bien siempre que se elija la posición inicial de tal manera que exista una forma de avanzar en sentido antihorario en un círculo cerrado. Para algunos bordes, en particular alrededor del exterior del mapa, esto no es posible. No tengo idea de cómo seleccionar buenas posiciones iniciales, o cómo filtrar soluciones, pero tal vez alguien más tenga una.

Ejemplo de trabajo (comenzando con edge (1204573687, 4555480822)):

Ejemplo, donde este enfoque no funciona (comenzando con edge (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')