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:
- Get the bounds of GeoTIFF.
- Query the Overpass API for all features with "natural=water".
- Initialise a numpy array with same dimensions as GeoTIFF.
- Translate Lat/Lon from OSM features to indexes in numpy array.
- 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?