Generating an image representing inland surface water using QGis and OSM data

Question

I'm trying to generate a raster image representing surface water for a given bounding box. I'm using OpenStreetMaps data as my input for water bodies, and I'm using QGis on Ubuntu to process this data. Eventually I would like to be able to take any bounding box as an input, and generate a raster image (2048px by 2048px) that has blue pixels where there is surface water, and black pixels where there is land (or sea; I'm not concerned with ocean water - these can be blue or black). Given that most maps/visualisation tools that use OSM data (such as openstreetmaps.org, Overpass, GeoJSON.io, etc) display inland water as solid blue polygons, I'm sure this can't be too hard to achieve!

At the moment I am trying to do this manually. Ideally this would all be accomplished with free and open source data and software. I intend to create a tool to do this automatically if and when I figure this out. Here is what I have tried:

First I grabbed some OSM data representing inland water using the Overpass API at http://overpass-turbo.eu/ I used the following query:

(
    way["natural"="water"]({bbox});
    rel["natural"="water"]({bbox});
    way["water"="lake"]({bbox});
    rel["water"="lake"]({bbox});
    way["waterway"]({bbox});
    rel["waterway"]({bbox});    
);
(._;>;);
out;

I then exported this data in GPX format and imported it as a vector layer into QGis. The lakes are all presented as outlines and the rivers as one-pixel-wide lines. I also imported a simple GeoJSON polygon that represented my bounding box, set it to have a black style, and clipped the other layer to it. I would now like to export the two layers as a 2048x2048 bitmap/png (or similar). I tried to use Raster Menu > Conversion > Rasterize on just the water layer, but this produced a blank raster. I then tried to use the Print Composer as suggested in an answer to this question but this proved to be a lot of work to get it to stretch correctly, and my polygons were not solid but just outlines so I had to fill them in manually.

I must be missing a relatively simple solution here. How can I generate a simple, 2 colour raster image representing surface water for a given bounding box?

Answer 1

I couldn't get your overpass query to work so I tried building it using the wizard. I confess I tend to use the QuickOSM plugin as I find the overpass query language a bit tricky :)

This is the query I went with:-

/*
This has been generated by the overpass-turbo wizard.
The original search was:
“natural=water or water=lake or waterway=*”
*/
[out:json][timeout:25];
// gather results
(
  // query part for: “natural=water”
  node["natural"="water"]({{bbox}});
  way["natural"="water"]({{bbox}});
  relation["natural"="water"]({{bbox}});
  // query part for: “water=lake”
  node["water"="lake"]({{bbox}});
  way["water"="lake"]({{bbox}});
  relation["water"="lake"]({{bbox}});
  // query part for: “waterway=*”
  node["waterway"]({{bbox}});
  way["waterway"]({{bbox}});
  relation["waterway"]({{bbox}});
);
// print results
out body;
>;
out skel qt;

Brought the GPX into QGIS.

However, attempting to rasterize gives the error message "not a known GDAL format". So I downloaded a GeoJSON instead - same problem. So I saved the GeoJSON to a shapefile instead, and used that.

However, I had to manually edit the GDAL command line (using the pencil icon next to the command) and include a "burn" value (here, I chose 10)

gdal_rasterize -burn 10 -ts 3000 3000 /path/to/my.shp /path/to/generated.tif

This should create a raster with 0 for the background, and 10 where there is water.

I got the same result; a black raster. The reason for this is that the default is to not stretch the histogram to the full range of values; instead it clips the bottom and top 2% of values. If fewer than 2% of pixel values are 10, it'll assume the range is 0..0, hence the black image.

If you go into Layer Properties > Rendering > Band Rendering, choose SingleBand Grey, Min/Max, Full Extent, Actual Accuracy and click Load it should show the full range of values (0 to 10).

Apply that, and you should now see your waterways.

In this case I just did one layer, the lines. You might want to use Merge Shapefiles to combine your various layers, and rasterize the combined shapefile.

If that doesn't work (e.g. mixing lines/polygons) you can either

• use Raster Calculator to mix your resulting rasters together
• use Buffer on your lines to make them polygons, then merge / rasterize

It looks a bit "gritty", with gaps. These will vanish if you zoom in. (I might have turned off anti-aliasing)

EDIT

To clarify how to fill in the lakes...

Import the GeoJSON file you exported from Overpass, and bring in the Polygon layer:-

If neccessary, save to shapefile. Rasterise this, and your lakes should be filled...

Answer 2

Have you tried using gdal_rasterize (http://www.gdal.org/gdal_rasterize.html) ?