Search code examples
pythonnetworkxpolygonshapelyosmnx

How to create NetworkX graph with n nodes (OSMnx compatible) from Shapely object


Let's assume we have the following Shapely object:

from shapely import geometry

# (lat, long) coordinates
p1 = geometry.Point(-96.87,32.88)
p2 = geometry.Point(-96.77,33.00)
p3 = geometry.Point(-96.67,32.95)
p4 = geometry.Point(-96.65,32.84)
pointList = [p1, p2, p3, p4, p1]

# creat geometric object
p = geometry.Polygon([[p.x, p.y] for p in pointList])

Next we want to create OSMnx compatible NetworkX graph from this polygon. The standard way would be:

import osmnx as ox
G = ox.graph.graph_from_polygon(p, network_type='drive', simplify=False)

But instead of underlining the 'real' road network, I want to generate a 'fake' road network with some density of nodes. Each node should be connected with some number of nearest neighbors. I am not interested in buildings, lakes, etc. I just want a plain NetworkX graph inside of my polygon with some number of randomly distributed nodes and some number of edges that connect those nodes.

def f_generate_plain_polygon(p, number_of_nodes, min_number_of_connection_for_each_node):
  ...

The generated graph should have some functionality from OSMnx. For example, we could visualize such a graph:

ox.plot_graph(G, show=False, close=False, edge_color='black', bgcolor='w', edge_alpha=0.5, node_color='black')
ox.utils_graph.remove_isolated_nodes(G)

Here is an example of what function f_generate_plain_polygon should generate. It is how polygon looks like:

1

And we are adding graph with 60 nodes and at least 3 connections with nearest neighbors for each node:

2

If for some reason it is more simple to generate a graph without explicitly specifying a number of nodes but just saying that we want to generate nodes with distance e.g. 0.001 between them it is also fine.


Solution

  • Here is a possible way you can approach your problem.

    1. First, in order to use ox.plot_graph with your randomly generated graph, you'll need to know what type of graph ox.plot_graph takes as input. From the documentation here, it looks like a nx.MultiDiGraph si waht you need (doc here).

    2. Next, you can randomly generate a directed graph with H=nx.gnp_random_graph(100, 0.1, seed=None, directed=True) for instance. Note that Networkx offers a lot of graph generation functions (see here) so just pick the one that works best for you.

    3. In your next step, you can cast your graph from nx.DiGraph to nx.MultiDiGraph by using I=nx.MultiDiGraph(H).

    4. You now have a random nx.MultiDiGraph. In order to draw the graph in your polygon you need to associate to each node a x and y coordinate sampled from within the polygon. I did that using the answer in this SO post

    5. You can then link these positions to your nodes following osmnx syntax by running:

    for i,node in enumerate(I.nodes(data=True)):
      node[1]['x']=pos[i][0]
      node[1]['y']=pos[i][1]
    
    1. As of now, the only thing missing for this to work is associate a crs to your random graph. ox.plot_graph seems to be requiring it. A simple workaround is to use the one generated for the osmnx graph G in your example. I.graph={'crs':'epsg:4326'}

    2. Finally, you can plot your graph and your polygon (see here for how to plot shapely object with matplotlib):

    x,y = p.exterior.xy
    fig,ax=ox.plot_graph(I,show=False, close=False, edge_color='black', bgcolor='w', edge_alpha=0.5, node_color='black')
    ox.utils_graph.remove_isolated_nodes(I)
    ax.plot(x,y)
    ax.set_xlim([-96.9,-96.6])
    ax.set_ylim([32.82,33.02])
    
    plt.show()
    

    enter image description here

    See full code here:

    from shapely import geometry
    from shapely.geometry import Point
    import matplotlib.pyplot as plt
    import numpy as np
    import networkx as nx
    import osmnx as ox
    
    # (lat, long) coordinates
    p1 = geometry.Point(-96.87,32.88)
    p2 = geometry.Point(-96.77,33.00)
    p3 = geometry.Point(-96.67,32.95)
    p4 = geometry.Point(-96.65,32.84)
    pointList = [p1, p2, p3, p4, p1]
    
    # creat geometric object
    p = geometry.Polygon([[p.x, p.y] for p in pointList])
    
    H=nx.gnp_random_graph(100, 0.1, seed=None, directed=True) #step 2
    I=nx.MultiDiGraph(H) #step 3
    
    def random_point_in_shp(shp):  #step 4
        within = False
        while not within:
            x = np.random.uniform(shp.bounds[0], shp.bounds[2])
            y = np.random.uniform(shp.bounds[1], shp.bounds[3])
            within = shp.contains(Point(x, y))
        return [x,y]
    pos=[random_point_in_shp(p) for _ in range(100)]
    
    for i,node in enumerate(I.nodes(data=True)): #step 5 
      node[1]['x']=pos[i][0]
      node[1]['y']=pos[i][1]
    
    I.graph={'crs':'epsg:4326'} #step 6
    
    x,y = p.exterior.xy #step 7
    fig,ax=ox.plot_graph(I,show=False, close=False, edge_color='black', bgcolor='w', edge_alpha=0.5, node_color='black')
    ox.utils_graph.remove_isolated_nodes(I)
    ax.plot(x,y)
    ax.set_xlim([-96.9,-96.6])
    ax.set_ylim([32.82,33.02])
    
    plt.show()
    

    EDIT: Example with a 3-neighbors graph (60 nodes, 180 edges)

    The core of the process is the same. The difference comes from creating the node position before the graph. And then, use the lypysal weights.KNN function to connect nodes with their 3 closest neighbors as in the link you provided in comments (here). The graph can then be turned into a nx.MultiDiGraph and the rest is the same.

    from shapely import geometry
    from shapely.geometry import Point
    import matplotlib.pyplot as plt
    import numpy as np
    import networkx as nx
    import osmnx as ox
    import libpysal
    
    # (lat, long) coordinates
    p1 = geometry.Point(-96.87,32.88)
    p2 = geometry.Point(-96.77,33.00)
    p3 = geometry.Point(-96.67,32.95)
    p4 = geometry.Point(-96.65,32.84)
    pointList = [p1, p2, p3, p4, p1]
    
    # creat geometric object
    p = geometry.Polygon([[p.x, p.y] for p in pointList])
    
    def random_point_in_shp(shp):  #step 4
        within = False
        while not within:
            x = np.random.uniform(shp.bounds[0], shp.bounds[2])
            y = np.random.uniform(shp.bounds[1], shp.bounds[3])
            within = shp.contains(Point(x, y))
        return [x,y]
    
    N_nodes=60
    pos=[random_point_in_shp(p) for _ in range(60)]
    
    #create 3-neighbor graph
    kd=libpysal.cg.KDTree(np.array(pos))
    knn3=libpysal.weights.KNN(kd, k=3)
    knn_graph = knn3.to_networkx()
    K=nx.MultiDiGraph(knn_graph)
    
    
    for i,node in enumerate(K.nodes(data=True)): #step 5 
      node[1]['x']=pos[i][0]
      node[1]['y']=pos[i][1]
    
    K.graph={'crs':'epsg:4326'} #step 6
    
    x,y = p.exterior.xy #step 7
    fig,ax=ox.plot_graph(K,show=False, close=False, edge_color='black', bgcolor='w', edge_alpha=0.5, node_color='black')
    ox.utils_graph.remove_isolated_nodes(K)
    ax.plot(x,y)
    ax.set_xlim([-96.9,-96.6])
    ax.set_ylim([32.82,33.02])
    
    plt.show()
    

    enter image description here