7

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?

3

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.

  • 2
    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 – mmd Dec 31 '16 at 10:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.