3

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?

  • 2
    I'm having trouble understanding where you are stuck. What I usually do, is to add My data to QGIS, render it how I need, and for a Quick export to image format, I use Project>>Save as Image. How does this not work for you Do you want to automate this? or something else? – Devdatta Tengshe Jul 6 '16 at 7:12
  • Thanks for your comment. Yes, my goal is to automate this, so I'm looking for a simpler solution than to manually edit files. If I can generate the rasters using the GUI without any manual editing, I am confident I can automate the process with the command line. Project>>Save as Image appears to save the current viewport and does not allow you to set the resolution. I need a way to directly convert (vector water data in bounding box) => (raster image of set resolution). Hope that's more clear! – Ivan Jul 10 '16 at 6:13
2
+50

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).

enter image description here

Apply that, and you should now see your waterways. enter image description here

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:-

enter image description here

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

enter image description here

  • Thanks for your detailed answer Steven! Is there a way to set the band rendering options directly through the gdal_rasterize command so that I can achieve this through the command line? Also, this does not appear to solve the problem of making the lake polygons solid, these are still just outlines? – Ivan Jul 10 '16 at 6:15
  • Make sure you choose the 'polygon' layer when you import your geojson export download. If you then save as shapefile, and rasterize as per the instructions, that should fill in your lake boundaries. – Steven Kay Jul 12 '16 at 21:56
  • @StevenKay is it possible to assign some styles to the line shapefile (like width) and then render it using gdal_rasterize? – FJ_Abbasi Apr 23 at 13:32
0

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

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.