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

10

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:

entrez la description de l'image ici

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é

tmo
la source
Problème intéressant. Vous voudrez peut-être y ajouter la balise d'algorithme. Semble que les quartiers seraient le problème le plus facile après avoir compris les blocs. En tant que quartiers, tout ce que vous recherchez, c'est un avantage partagé, n'est-ce pas? Et chaque bloc aura une liste d'arêtes ... Pour les blocs, je pense qu'il sera utile d'obtenir la direction cardinale de chaque option de rue à un nœud et de "continuer à tourner à droite" (ou à gauche) jusqu'à ce que vous ayez terminé un circuit ou atteint une impasse ou une boucle sur vous-même et revenir en arrière récursivement. Il semble qu'il y aurait cependant des cas d'angle intéressants.
Jeff H
Je pense que cette question est très similaire à votre problème non. 1. Comme vous pouvez le voir dans le lien, j'ai travaillé un peu sur le problème, et il est noueux (se révèle être NP-difficile). L'heuristique de ma réponse, cependant, pourrait encore vous donner de bons résultats.
Paul Brodersen
Étant donné que la solution que vous jugerez acceptable sera probablement également une heuristique, il peut être judicieux de définir un ensemble de données de test pour valider chaque approche. Cela signifie que, pour votre exemple de graphique, il serait bon d'avoir une annotation de tous les blocs sous une forme lisible par machine - pas seulement quelques exemples dans une image.
Paul Brodersen

Réponses:

3

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)

graphique élagué

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)

tracé des étiquettes de région

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)

tracé de contour superposé sur un graphique élagué

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)

tracé du réseau avec des nœuds appartenant au bloc en rouge

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

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)
sel-die
la source
Salut salt-die, bienvenue à SO et merci pour votre contribution. En faisant, nx.Graph(G)je perds beaucoup d'informations (directivité et type multigraph), j'ai donc du mal à vérifier votre réponse, car je ne peux pas sembler relier le nouveau graphique Ià mon graphique d'origine G.
tmo
Ce sera un peu de travail pour conserver les informations géométriques du graphique d'origine. Je vais essayer ça bientôt.
salt-die
@tmo ne fait que passer: vous devriez pouvoir utiliser la classe MultiDiGraph (qui étend Graph) dans ce cas
Théo Rubenach
1

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.

Hashemi Emad
la source
C'est une bien meilleure idée que moi. J'ajouterai une réponse avec une implémentation de votre intuition.
Paul Brodersen
0

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

entrez la description de l'image ici

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

entrez la description de l'image ici

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