Search code examples
networkxopenstreetmapgeopandasosmnx

Extract constrained polygon using OSMnx

I'm playing around with OSMnx package in order to solve the following task: - there is a point X on the map defined by latitude and longitude - we need to detect a polygon that contains that point X and is constrained by the neighboring roads - so basically the point X is inside the polygon and the neighboring roads will be borders of that polygon.

So far I've managed only to plot the visualization of the graph on the map and find the closest edge/node to point X.

In the attached image I've highlighted the area I want to extract in red.

enter image description here

Solution

  • As you are trying to find a polygon containing your point, you first need to generate polygons out of multilinestring geometry. As you did not provide your data I am downloading a sample from OSM using OSMnx.

    import osmnx as ox
import geopandas as gpd
import shapely

point = (40.742623, -73.977857)

streets_graph = ox.graph_from_point(point, distance=500, network_type='drive')
streets_graph = ox.project_graph(streets_graph)

    I have reprojected it as it is way more convenient than working with degrees, especially if you want to measure anything.

    You then have to convert OSMnx graph to geopandas GeoDataFrame.

    streets = ox.save_load.graph_to_gdfs(streets_graph, nodes=False, edges=True,
                                     node_geometry=False, fill_edge_geometry=True)

    To get some point I can work with, I'll just use the one in the centre of this geodataframe.

    point = streets.unary_union.centroid

    This is what it looks like.

    network and point overlayed

    Next you need to get polygons of your blocks defined by streets, using shapely.ops.polygonize as I suggested in the comment above and store them as GeoSeries.

    polygons = shapely.ops.polygonize(streets.geometry)
polygons = gpd.GeoSeries(polygons)

    polygonized network

    The only thing you have to do next is to find which polygon contains your point.

    target = polygons.loc[polygons.contains(point)]

    Plotting it again:

    ax = target.plot()
gpd.GeoSeries([point]).plot(ax=ax, color='r')

    final polygon

    If you want to know which streets are forming the boundary of this polygon, just intersect it with the original network. I am filtering for MultiLineString to exclude streets which intersects the polygon only in one point.

    target_streets = streets.loc[streets.intersection(target.iloc[0]).type == 'MultiLineString']

    This is how the result of that looks like.

    ax = target_streets2.plot()
gpd.GeoSeries([point]).plot(ax=ax, color='r')

    streets intersecting polygon

    Hope it helps.