Znajdowanie dzielnic (klik) w danych ulicznych (wykres)

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:

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

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)

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)

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)

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)

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.")

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)

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.

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

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

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