Finding neighbourhoods (cliques) in street data (a graph)

Question

I am looking for a way to automatically define neighbourhoods in cities as polygons on a graph.

My definition of a neighbourhood has two parts:

1. A block: An area inclosed between a number of streets, where the number of streets (edges) and intersections (nodes) is a minimum of three (a triangle).
2. A neighbourhood: For any given block, all the blocks directly adjacent to that block and the block itself.

See this illustration for an example:

E.g. B4 is block defined by 7 nodes and 6 edges connecting them. As most of the examples here, the other blocks are defined by 4 nodes and 4 edges connecting them. Also, the neighbourhood of B1 includes B2 (and vice versa) while B2 also includes B3.

I am using osmnx to get street data from OSM.

1. Using osmnx and networkx, how can I traverse a graph to find the nodes and edges that define each block?
2. For each block, how can I find the adjacent blocks?

I am working myself towards a piece of code that takes a graph and a pair of coordinates (latitude, longitude) as input, identifies the relevant block and returns the polygon for that block and the neighbourhood as defined above.

Here is the code used to make the map:

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)

and my attempt at finding cliques with different number of nodes and degrees.

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)

Theory that might be relevant:

Enumerating All Cycles in an Undirected Graph

Answer 1

Finding city blocks using the graph is surprisingly non-trivial. Basically, this amounts to finding the smallest set of smallest rings (SSSR), which is an NP-complete problem. A review of this problem (and related problems) can be found here. On SO, there is one description of an algorithm to solve it here. As far as I can tell, there is no corresponding implementation in networkx (or in python for that matter). I tried this approach briefly and then abandoned it -- my brain is not up to scratch for that kind of work today. That being said, I will award a bounty to anybody that might visit this page at a later date and post a tested implementation of an algorithm that finds the SSSR in python.

I have instead pursued a different approach, leveraging the fact that the graph is guaranteed to be planar. Briefly, instead of treating this as a graph problem, we treat this as an image segmentation problem. First, we find all connected regions in the image. We then determine the contour around each region, transform the contours in image coordinates back to longitudes and latitudes.

Given the following imports and function definitions:

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

Load the data. Do cache the imports, if testing this repeatedly -- otherwise your account can get banned. Speaking from experience here.

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 cannot be part of a cycle. This step is not strictly necessary but results in nicer contours.

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

Convert plot to image and find connected regions:

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)

For each labelled region, find the contour and convert the contour pixel coordinates back to data coordinates.

# 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 all points in the original graph that fall inside (or on) the 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)

Figuring out if two blocks are neighbors is pretty easy. Just check if they share a node:

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

Answer 2

I'm not completely sure that cycle_basis will give you the neighborhoods you seek, but if it does, it's a simple thing to get the neighborhood graph from it:

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)

Answer 3

I appreciate this question is a little bit old now but I have an alternative approach that is relatively straighforward - it does require stepping away from networkx for a moment though.

Creating the blocks

Get the projected graph:

import osmnx as ox
import geopandas as gpd
from shapely.ops import polygonize

G = ox.graph_from_address('Nørrebrogade 20, Copenhagen Municipality',
                          network_type='all', 
                          dist=500)
G_projected = ox.project_graph(G)

Convert the graph to undirected - this removes duplicate edges that would cause the subsequent polygonization to fail:

G_undirected = G_projected.to_undirected()

Extract just the edges into a GeoPandas GeoDataFrame:

G_edges_as_gdf = ox.graph_to_gdfs(G_undirected, nodes=False, edges=True)

Use polygonize from shapely.ops on the edges to create the block faces and then use these as the geometry in a new GeoDataFrame:

block_faces = list(polygonize(G_edges_as_gdf['geometry']))
blocks = gpd.GeoDataFrame(geometry=block_faces)

Plot the result:

ax = G_edges_as_gdf.plot(figsize=(10,10), color='red', zorder=0)
blocks.plot(ax=ax, facecolor='gainsboro', edgecolor='k', linewidth=2, alpha=0.5, zorder=1)

blocks created from line fragments by shapely.ops.polygonize()

Finding neighbors

PySAL does a great job of this with its spatial weights see https://pysal.org/notebooks/lib/libpysal/weights.html for further information. libpysal can be installed from conda-forge.

Here we use Rook weights to identify blocks that share an edge as in the original question. Queen weights would also include those that only share a node (i.e. meet at a street junction)

from libpysal.weights import Rook # Queen, KNN also available
w_rook = Rook.from_dataframe(blocks)

The spatial weights are just for the neighbours so we need to append the original block (index number 18 here just as an example):

block = 18
neighbors = w_rook.neighbors[block]
neighbors.append(block)
neighbors

Then we can plot using neighbors as a filter:

ax = blocks.plot(figsize=(10,10), facecolor='gainsboro', edgecolor='black')
blocks[blocks.index.isin(neighbors)].plot(ax=ax, color='red', alpha=0.5)

neighboring blocks highlighted on top of all blocks

Notes

• This doesn't resolve the sliver polygons caused by multiple road centre lines and it can be challenging to resolve ambiguities caused by pedestrianised streets and footpaths - presumably a pedestrianised street is a legitimate division between blocks but a footpath may not be.
• From memory I have in the past needed to fragment the LineStrings created from the edges further before polygonizing them - this didn't seem to be necessary with this example.
• I have left out the final step of linking the new geometries back to the Node IDs with some kind of spatial join etc.

Answer 4

I don't have a code, but I guess that once i'm on the sidewalk, if I keep turning to the right at each corner, I will cycle through the edges of my block. I don't know the libraries so I'll just talk algo here.

• from your point, go north until you reach a street
• turn right as much as you can and walk on the street
• on the next corner, find all the steets, chose the one that makes the smallest angle with your street counting from the right.
• walk on that street.
• turn right, etc.

It's actually an algorithm to use to exit a maze : keep your right hand on the wall and walk. It doesn't work in case of loops in the maze, you just loop around. But it gives a solution to your problem.

Answer 5

This is an implementation of Hashemi Emad's idea. It works well as long as the starting position is chosen such that there exists a way to step counterclockwise in a tight circle. For some edges, in particular around the outside of the map, this is not possible. I don't have an idea how to select good starting positions, or how to filter solutions -- but maybe somebody else has one.

Working example (starting with edge (1204573687, 4555480822)):

Example, where this approach does not work (starting with edge (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')