How to fill water bodies with OSMnx in Python

user_1311095 0
Viewed 1398

I currently use OSMnx for a project to draw road networks in an area.

I'd now like to add water bodies so that we can clearly see which parts of an area are water and land.

So far, I've been able to identify waterbodies using the custom_filter argument to OSMnx's graph functions. I can then outline the waterbodies using the plot_graph function.

Ideally, I'd want to fill the waterbodies in (rather than only outline them). I feel like this should be possible, as in OpenStreetMap water bodies are filled, but I can't figure out how to do this with OSMnx. Does anybody have any ideas?

Here's what I have currently:

import osmnx as ox

# Get water bodies map of the New York City area
G = ox.graph_from_bbox(40.9666,40.4362,-73.6084,-74.3254, custom_filter='["natural"~"water|coastline"]', retain_all = True)

# Plot the graph in blue on a white background
ox.plot_graph(G, bgcolor='white', node_size=0, equal_aspect=True, edge_color='blue')

Which produces this image:

NYC Water Bodies Image

Do I need to somehow use a geodataframe with PlotShape? Or is plot_footprints what I need? I haven't been able to find examples of people plotting waterbodies. It seems the GDF is generally used for plotting the map of a place, and footprints is used for buildings. Though as those are both polygon-oriented plots, I feel that might be the right way.

osmnx gis openstreetmap python
2 Answers

I have been fighting with this for a while for a little project I have. What the accepted solution does is plot the fetched land on top of a water_color background, and adds the closed polygons fetched with the tag {'natural':'water'}. This works fine only if land covers the whole field of view. If you keep only one of the places in the list (say 'Manhattan, NY, USA'), then you have only this selected land in a sea of blue.

Fine if that's what you want.

What I wanted to do, and what I suspect the OP also wanted (since they fetched info from a bounding box), was to have all water-land interfaces within the bounding box. Including coastline does the job if what is needed is the contours, but coastlines are not closed polygons (and they come in separate segments) so there is no easy way to do that.

I started looking into using osmnx.geometries.linemerge and then osmnx.geo_utils.split to break a polygon along the coastlines, but I finally found that someone had done all the job already:

https://osmdata.openstreetmap.de/contact.html

This repository has all coastlines joined into polygons (either on the water side or the land side). The github repo is https://github.com/fossgis .

So I would say that the clean way of doing what the OP asks is to download and use these shapefiles. Assuming that the water polygons file has been unzipped in the working directory, this is a working example:

import osmnx as ox
import geopandas as gpd
# Bounding box
bN, bS, bE, bW = 40.9666, 40.4362, -73.6084, -74.3254
# Fetch water
G = ox.geometries.geometries_from_bbox(bN, bS, bE, bW, tags={"natural": "water"})
# Load coastline polygons
water = gpd.read_file('water-polygons-split-4326/water_polygons.shp', bbox=(bW, bN, bE, bS))
# Plot
fig, ax = ox.plot_footprints(water, bbox=(bN, bS, bE, bW),
                             color=water_color, bgcolor=land_color,
                             show=False, close=False)
ax = G.plot(ax=ax, fc=water_color, markersize=0)

Output figure

user_15608760

This isn't perfect but it gets you nearly there:

import osmnx as ox
ox.config(log_console=True, use_cache=True)

# add more items here to get all the landforms you want
places = ['Manhattan, NY, USA', 'Brooklyn, NY, USA', 'Queens, NY, USA', 'Bronx, NY, USA']
land = ox.geocode_to_gdf(places)

# get the water bodies
left, bottom, right, top = land.total_bounds
bbox = top, bottom, right, left
poly = ox.utils_geo.bbox_to_poly(*bbox)
water = ox.geometries_from_polygon(poly, tags={'natural': 'water'})

# constrain the plotting window as desired
c = land.unary_union.centroid
bbox = ox.utils_geo.bbox_from_point((c.y, c.x), dist=12000)

water_color = 'blue'
land_color = '#aaaaaa'
fig, ax = ox.plot_footprints(water, bbox=bbox,
                             color=water_color, bgcolor=water_color,
                             show=False, close=False)
ax = land.plot(ax=ax, zorder=0, fc=land_color)

osm land vs water

The key issue is that I'm currently unclear if OSM can be consistently queried for straightforward land vs water polygons (I don't normally work with land/water boundaries in my research). The places may be political boundaries, which could overlap with the water area in real life. You may want to experiment with what you query/plot as land vs water here.

user_73219420
Related