Recherche de quartiers (cliques) dans les données de rue (un graphique)

Je cherche un moyen de définir automatiquement les quartiers des villes comme des polygones sur un graphique.

Ma définition d'un quartier comporte deux parties:

1. Un bloc : zone comprise entre plusieurs rues, où le nombre de rues (bords) et d'intersections (nœuds) est au minimum de trois (un triangle).
2. Un quartier : pour un bloc donné, tous les blocs directement adjacents à ce bloc et le bloc lui-même.

Voir cette illustration pour un exemple:

Par exemple B4 est un bloc défini par 7 nœuds et 6 arêtes les reliant. Comme la plupart des exemples ici, les autres blocs sont définis par 4 nœuds et 4 arêtes les reliant. De plus, le quartier de B1 comprend B2 (et vice versa) tandis que B2 comprend également B3 .

J'utilise osmnx pour obtenir des données de rue d'OSM.

1. En utilisant osmnx et networkx, comment puis-je parcourir un graphique pour trouver les nœuds et les arêtes qui définissent chaque bloc?
2. Pour chaque bloc, comment trouver les blocs adjacents?

Je travaille moi-même vers un morceau de code qui prend un graphique et une paire de coordonnées (latitude, longitude) en entrée, identifie le bloc pertinent et renvoie le polygone pour ce bloc et le quartier tel que défini ci-dessus.

Voici le code utilisé pour faire la carte:

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)

et ma tentative de trouver des cliques avec un nombre différent de nœuds et de degrés.

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)

Théorie qui pourrait être pertinente:

Énumération de tous les cycles dans un graphique non orienté

Trouver des blocs de villes à l'aide du graphique est étonnamment non trivial. Fondamentalement, cela revient à trouver le plus petit ensemble de plus petits anneaux (SSSR), qui est un problème NP-complet. Un examen de ce problème (et des problèmes connexes) peut être trouvé ici . Sur SO, il existe une description d'un algorithme pour le résoudre ici . Pour autant que je sache, il n'y a pas d'implémentation correspondante dans networkx(ou en python d'ailleurs). J'ai essayé brièvement cette approche, puis je l'ai abandonnée - mon cerveau n'est pas à la hauteur pour ce genre de travail aujourd'hui. Cela étant dit, je remettrai une prime à quiconque pourrait visiter cette page à une date ultérieure et publier une implémentation testée d'un algorithme qui trouve le SSSR en python.

Au lieu de cela, j'ai poursuivi une approche différente, en tirant parti du fait que le graphique est garanti d'être plan. En bref, au lieu de traiter cela comme un problème de graphe, nous le traitons comme un problème de segmentation d'image. Tout d'abord, nous trouvons toutes les régions connectées dans l'image. Nous déterminons ensuite le contour autour de chaque région, transformons les contours en coordonnées image en longitudes et latitudes.

Étant donné les importations et les définitions de fonction suivantes:

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

Chargez les données. Mettez en cache les importations, si vous testez cela à plusieurs reprises - sinon votre compte peut être banni. Parlant d'expérience ici.

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

Taillez les nœuds et les arêtes qui ne peuvent pas faire partie d'un cycle. Cette étape n'est pas strictement nécessaire mais donne de plus beaux contours.

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

Convertissez le tracé en image et trouvez les régions connectées:

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)

Pour chaque région étiquetée, recherchez le contour et reconvertissez les coordonnées des pixels du contour en coordonnées de données.

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

Déterminez tous les points du graphique d'origine qui se trouvent à l'intérieur (ou sur) le contour.

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)

Il est assez facile de déterminer si deux blocs sont voisins. Vérifiez simplement s'ils partagent un nœud:

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

Je ne suis pas complètement sûr que cycle_basiscela vous donnera les quartiers que vous recherchez, mais si c'est le cas, c'est une chose simple pour en obtenir le graphique de quartier:

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)

Je n'ai pas de code, mais je suppose qu'une fois que je serai sur le trottoir, si je continue à tourner à droite à chaque coin, je vais parcourir les bords de mon bloc. Je ne connais pas les bibliothèques donc je vais juste parler algo ici.

• de votre point, allez vers le nord jusqu'à atteindre une rue
• tournez à droite autant que possible et marchez dans la rue
• au coin suivant, trouvez toutes les steets, choisissez celle qui fait le plus petit angle avec votre rue en comptant à droite.
• marcher dans cette rue.
• tourner à droite, etc.

C'est en fait un algorithme à utiliser pour sortir d'un labyrinthe: gardez votre main droite sur le mur et marchez. Cela ne fonctionne pas en cas de boucles dans le labyrinthe, il suffit de boucler. Mais cela donne une solution à votre problème.

Ceci est une mise en œuvre de l'idée de Hashemi Emad . Cela fonctionne bien tant que la position de départ est choisie de telle sorte qu'il existe un moyen de faire un pas dans le sens antihoraire dans un cercle serré. Pour certains bords, en particulier autour de l'extérieur de la carte, cela n'est pas possible. Je ne sais pas comment sélectionner les bonnes positions de départ, ni comment filtrer les solutions - mais peut-être que quelqu'un d'autre en a une.

Exemple de travail (en commençant par le bord (1204573687, 4555480822)):

Exemple, où cette approche ne fonctionne pas (en commençant par le bord (1286684278, 5818325197)):

Code

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