Znajdowanie dzielnic (klik) w danych ulicznych (wykres)

10

Szukam sposobu automatycznego zdefiniowania dzielnic w miastach jako wielokątów na wykresie.

Moja definicja sąsiedztwa składa się z dwóch części:

  1. Blok : Obszar zawarty między wieloma ulicami, w którym liczba ulic (krawędzie) i skrzyżowań (węzły) wynosi co najmniej trzy (trójkąt).
  2. Sąsiedztwo : dla każdego bloku wszystkie bloki bezpośrednio przylegające do tego bloku i sam blok.

Zobacz tę ilustrację jako przykład:

wprowadź opis zdjęcia tutaj

Np. B4 to blok zdefiniowany przez 7 węzłów i 6 krawędzi łączących je. Jak większość przykładów tutaj, pozostałe bloki są zdefiniowane przez 4 węzły i 4 krawędzie łączące je. Ponadto sąsiedztwo z B1 obejmuje B2 (i vice versa), a B2 również B3 .

Korzystam z osmnx, aby uzyskać dane uliczne z OSM.

  1. Używając osmnx i networkx, w jaki sposób mogę przechodzić przez wykres, aby znaleźć węzły i krawędzie, które definiują każdy blok?
  2. Jak mogę znaleźć sąsiednie bloki dla każdego bloku?

Pracuję nad kawałkiem kodu, który przyjmuje na wejściu wykres i parę współrzędnych (szerokość, długość), identyfikuje odpowiedni blok i zwraca wielokąt dla tego bloku i sąsiedztwa, jak zdefiniowano powyżej.

Oto kod użyty do stworzenia mapy:

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)

i moja próba znalezienia klik z różną liczbą węzłów i stopni.

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, która może być istotna:

Zliczanie wszystkich cykli na niekierowanym wykresie

tmo
źródło
Ciekawy problem. Możesz dodać do niego znacznik algorytmu. Wygląda na to, że sąsiedztwo byłoby łatwiejszym problemem po tym, jak wymyślicie klocki. Jako dzielnice wszystko, czego szukasz, to wspólna przewaga, prawda? I każdy blok będzie miał listę krawędzi ... Myślę, że w przypadku bloków pomocne będzie uzyskanie głównego kierunku każdej opcji ulicy w węźle i „skręcanie w prawo” (lub w lewo) aż do ukończenia obwodu lub dotarcia ślepy zaułek lub zapętlić się i cofnąć rekurencyjnie. Wygląda jednak na to, że byłyby interesujące przypadki narożne.
Jeff H
Myślę, że to pytanie jest bardzo podobne do twojego problemu nie. 1. Jak widać w linku, trochę popracowałem nad problemem i jest on poważny (okazuje się, że jest trudny do NP). Heurystyka w mojej odpowiedzi może jednak dać ci wystarczająco dobre wyniki.
Paul Brodersen
Ponieważ każde rozwiązanie, które możesz uznać za akceptowalne, będzie prawdopodobnie również heurystyczne, dobrym pomysłem może być zdefiniowanie zestawu danych testowych w celu sprawdzenia poprawności każdego podejścia. Oznacza to, że na przykładowym wykresie dobrze byłoby mieć adnotację o wszystkich blokach w formie do odczytu maszynowego - a nie tylko kilka przykładów na obrazie.
Paul Brodersen

Odpowiedzi:

3

Znalezienie bloków miast za pomocą wykresu jest zaskakująco nietrywialne. Zasadniczo sprowadza się to do znalezienia najmniejszego zestawu najmniejszych pierścieni (SSSR), co stanowi problem NP-zupełny. Przegląd tego problemu (i powiązanych problemów) można znaleźć tutaj . Na tak, to jest jeden opis algorytmu rozwiązać go tutaj . O ile mi wiadomo, nie ma odpowiedniej implementacji w networkx(lub w Pythonie). Próbowałem krótko to podejście, a potem je porzuciłem - mój mózg nie jest dziś w stanie sprostać tego rodzaju pracy. Biorąc to pod uwagę, przyznam nagrodę każdemu, kto może odwiedzić tę stronę w późniejszym terminie i opublikować przetestowaną implementację algorytmu, który znajdzie SSSR w pythonie.

Zamiast tego wybrałem inne podejście, wykorzystując fakt, że wykres ma zagwarantowaną płaskość. W skrócie, zamiast traktować to jako problem z grafem, traktujemy to jako problem z segmentacją obrazu. Najpierw znajdujemy wszystkie połączone regiony na obrazie. Następnie określamy kontur wokół każdego regionu, przekształcamy kontury we współrzędnych obrazu z powrotem na długości i szerokości geograficzne.

Biorąc pod uwagę następujące definicje importu i funkcji:

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

Załaduj dane. Wykonuj buforowanie importu, jeśli testujesz to wielokrotnie - w przeciwnym razie twoje konto może zostać zbanowane. Mówiąc z doświadczenia tutaj.

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

Przycinaj węzły i krawędzie, które nie mogą być częścią cyklu. Ten krok nie jest absolutnie konieczny, ale powoduje ładniejsze kontury.

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

przycięty wykres

Konwertuj wykres na obraz i znajdź połączone regiony:

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)

wykres etykiet regionu

Dla każdego oznaczonego regionu znajdź kontur i przekonwertuj współrzędne piksela konturu z powrotem na współrzędne danych.

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

wykres konturu nałożony na przycięty wykres

Określ wszystkie punkty na oryginalnym wykresie, które mieszczą się w (lub na) konturze.

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)

wykres sieci z węzłami należącymi do bloku na czerwono

Ustalenie, czy dwa bloki są sąsiadami, jest dość łatwe. Wystarczy sprawdzić, czy współużytkują węzeł:

if set(nodes_inside_block_1) & set(nodes_inside_block_2): # empty set evaluates to False
    print("Blocks are neighbors.")
Paul Brodersen
źródło
2

Nie jestem do końca pewien, cycle_basisczy da ci to dzielnice, których szukasz, ale jeśli tak, to łatwo jest uzyskać z niego wykres sąsiedztwa:

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)
sól umrzeć
źródło
Cześć słona, witam w SO i dziękuję za włamanie się. Robiąc nx.Graph(G)tracę dużo informacji (ukierunkowanie i typ multigraficzny), więc trudno mi zweryfikować twoją odpowiedź, ponieważ nie mogę powiązać nowego wykresu Iz mój oryginalny wykres G.
tmo
Zachowanie informacji geometrycznych z oryginalnego wykresu będzie trochę trudne. Spróbuję to wkrótce.
salt-die
@tmo właśnie przechodzi: w takim przypadku powinieneś być w stanie użyć klasy MultiDiGraph (która rozszerza Graph)
Théo Rubenach
1

Nie mam kodu, ale myślę, że gdy będę na chodniku, jeśli będę nadal skręcał w prawo na każdym rogu, będę przechodził przez krawędzie mojego bloku. Nie znam bibliotek, więc porozmawiam tutaj o algo.

  • od punktu idź na północ, aż dojdziesz do ulicy
  • jak najdalej skręć w prawo i idź ulicą
  • w następnym rogu znajdź wszystkie ulice, wybierz ten, który robi najmniejszy kąt z twoją ulicą licząc od prawej.
  • chodzić po tej ulicy.
  • skręć w prawo itp.

W rzeczywistości jest to algorytm do wyjścia z labiryntu: trzymaj prawą rękę na ścianie i idź. To nie działa w przypadku pętli w labiryncie, po prostu się zapętlasz. Ale daje rozwiązanie twojego problemu.

Hashemi Emad
źródło
To znacznie lepszy pomysł niż ja. Dodam odpowiedź z implementacją twojej intuicji.
Paul Brodersen
0

To realizacja pomysłu Hashemi Emada . Działa dobrze, dopóki pozycja początkowa jest wybrana w taki sposób, że istnieje sposób, aby kroczyć w ciasnym okręgu przeciwnie do ruchu wskazówek zegara. W przypadku niektórych krawędzi, w szczególności wokół mapy, nie jest to możliwe. Nie mam pojęcia, jak wybrać dobre pozycje początkowe ani jak filtrować rozwiązania - ale może ktoś inny ma taką pozycję.

Przykład roboczy (zaczynający się od krawędzi (1204573687, 4555480822)):

wprowadź opis zdjęcia tutaj

Przykład, w którym to podejście nie działa (zaczynając od krawędzi (1286684278, 5818325197)):

wprowadź opis zdjęcia tutaj

Kod

#!/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
źródło