Create bitmap from OSM data and GeoTIFF in Python

Question

I want to create a bitmap (as a 2D numpy array) for a given satellite image in GeoTIFF format. The bitmap should have the same dimensions as the GeoTIFF and should get the value 1 at every pixel where the GeoTIFF has water.

My current approach would be:

1. Get the bounds of GeoTIFF.
2. Query the Overpass API for all features with "natural=water".
3. Initialise a numpy array with same dimensions as GeoTIFF.
4. Translate Lat/Lon from OSM features to indexes in numpy array.
5. Use indexes from step 4 to set entries in numpy array to 1.

My problem lies at step 5. Water areas are often closed ways which is simply a list of nodes. This means that in my numpy array I would have to identify all indices which lie in the polygon defined by the given nodes. This sounds like a problem suitable for Shapely but I looked at the documentation for Polygons and couldn't find anything.

I could loop through every possible index in my numpy array and check if that point lies in my polygon but that seems to be a bit cumbersome for me.

So my question is: Is there an easy solution to the problem in step 5?

Does this whole process seem reasonable or is there a more straightforward way to create a bitmap for the occurrence of water in a GeoTIFF?

Answer 1

This process is often called "rasterization" and your approach is right on. Shapely doesn't do this, but its sibling package, Rasterio, does your steps 1, 3, 4, 5: https://mapbox.github.io/rasterio/topics/features.html#burning-shapes-into-a-raster. If you can get waterbody polygons from OSM in GeoJSON format (using http://wiki.openstreetmap.org/wiki/Overpass_turbo/GeoJSON, maybe), you can burn the polygons into a numpy array using Rasterio's features.rasterize method: https://mapbox.github.io/rasterio/api/rasterio.features.html#rasterio.features.rasterize.

• mmd: Checking for natural=water will not suffice if you're in the middle of the ocean (see overpass-turbo.eu/s/kWS as an example), as it would not return any data at all. Besides, including waterway=* may also be needed. As a more general advise, you may want to take a look at tools like OSMCoastLine, or maybe the OpenStreetMap-Carto rendering style on osm.org for some inspiration. You can even download ready-made water polygons here: openstreetmapdata.com/data/water-polygons