Creating a GIS layer for "Distance from..." in QGIS
|Back to home page||
In a recent post, I showed how to deal with "distance from..." data in GIS layers using the R packages for handling spatial information. The example I used there involved
We'll use the same toy data set which you can download here. Extract all the files to a folder on your hard drive. You can open the project in QGIS by double-clicking on the file "toy_example.qgs".
The example has a raster file for land use and vector layers for roads, streams and camera locations.
Let me say right away that I'm using QGIS 1.8.0 on Windows OS, either Windows XP on my laptop or Windows 7 64-bit on my desktop. I haven't tried it, but this should work fine on Linux-type OSs such as Ubuntu. It might work on Macs prior to version 10.8; there seem to be problems running QGIS on v.10.8 Macs.
The roads layer is a vector (shape file), so we first need to rasterize it with the Raster > Conversion > Rasterize plugin. Once that's done, using the Proximity plugin is easy.
Rasterizing the roads layer
If we want our new layer to cover the whole extent of the study area and we want the raster pixels to match up with those of the LandUse raster, we need to get the necessary details from the LandUse layer.
Right-click on the LandUse layer in the Layers window and select Properties, then go to the Metadata tab. Scroll down and note the Pixel Size (it's 1000, -1000) and the Layer Extent (-230000,-200000, 290000, 130000; these are xmin, xmax, ymin and ymax respectively).
Close the properties box and go to Raster > Conversion > Rasterize. For the input file, select the "roads" shape file. Leave the Attribute as "id" and enter a format and name for the output file: I used GeoTIFF and called it "roads_raster.tif". A box pops up to tell you that the file doesn't exist and you must set up the output size: ok, we'll do that, click OK.
Make sure the box next to New Size is checked. In the code panel near the bottom of the window you'll see text looking a bit like the following, but with your full file names instead of "your/path".
The two numbers after -tr are the resolution (ie. pixel size in x and y directions, both must be positive), and the four numbers after -te are the extent of the raster, both of which we got from the LandUse properties box. Click on the Help button in the Rasterize window for more details.
Now click on OK. A "Finished" box pops up and the raster layer will appear in the layers window provided you checked the "Load into canvas..." box. Close the Rasterize box - don't click OK again or you'll get another copy of the raster!
For the screen shot below, I just changed all the non-zero values to red and the 0 pixels to yellow, but any colours will do, even the blue-to-red colormap colours.
Creating the Distance from... raster
Now we are ready to do the distance-from raster. Go to Raster > Analysis > Proximity. Select roads_raster as the input file and decide on a name for the output file: I called mine "DistanceFromRoad.tif". Make sure "Dist units" are set to GEO (for GEOreferenced coordinates), check "Load into canvas..." and press OK. When finished, close the Proximity dialogue box.
In the screen shot, I've used a simple colormap for the DistanceFromRoads layer and put the roads layer on top. As you can see, they fit well.
For this you will need the Point Sampling Tool plugin. You'll probably need to fetch and install it first; see a tutorial on doing that here.
Activate the cameras and DistanceFromRoad layers. Go to Plugins > Analyses > Point sampling tool. The "Layer containing sampling points" should be "cameras". In the "Layers with fields/bands..." select all the camera fields and the DistanceFromRoad layer.
Click on the Fields tab and check the names of the fields: they are limited to 10 characters. You'll probably want to change "DistanceFr" to something more meaningful, such as "roadDist". Back in the General tab, put in a file name for the "Output point vector layer", check "Add created layer to the TOC" and click OK. A "Finished" box pops up, then close the Point Sampling Tool.
The new cameras4 layer will appear in the QGIS Layers window. Look at the attribute table and you'll see that it has the roadDist for each camera.
To check, sort by roadDist (click on the roadDist column header), select those with zero distance, and then see were they are on the map: they should be on a road.
Permanent water includes streams on the streams layer and water bodies (lakes, reservoirs, etc) on the land use layer.
Fixing the CRS of the streams layer
First we need to convert the streams vector layer to a raster layer which must match up with the LandUse raster. As we discovered when using R, the streams layer does not use the same Coordinate Reference System (CRS) as the other layers, even though it shows up nicely in QGIS when "on the fly" CRS transformation is enabled.
Go to Settings > Project Properties and the Coordinate Reference System (CRS) tab; uncheck the box next to "Enable 'on the fly' CRS transformation". The streams disappear - or at least they're reduced to a blue dot in the middle of the map!
Look at the Properties of the streams layer and on the General tab you'll see that the CRS is "EPSG:4326 - WGS 84". We need to change the CRS before we make the raster. And that's not as simple as changing the CRS description in the Properties box!
Close the Properties box, right click on the streams layer in the Layers list, and select Save As... Select a suitable name for the output file: I used "streams_projected.shp". Change the CRS from "Layer CRS" to "Selected CRS" and then Browse to "WGS 84 / Pseudo Mercator". (You could use "Project CRS", which is in fact Pseudo Mercator, but I'd prefer to see the full CRS specification up there.) Check "Add saved file..." then click OK.
The new streams_projected layer looks right without on-the-fly transformation. Remove the old streams layer.
Converting streams to a raster layer
Create a raster from the stream_projected layer just as we did for the roads earlier, using Raster > Conversion > Rasterize and putting code in the dialogue box like this:
The resulting raster should look something like this - I used a colormap style with white for 0 in the raster.
Combining the streams with water bodies from the LandUse layer
We'll use Raster > Raster Calculator to combine the two to form a new layer. Open the Raster Calculator and enter a name for the new file; I used "water_all.tif". Leave the default extents as they are.
Now we build a "Raster calculator expression" using logical operators.
On the LandUse layer, water bodies are coded as 5 (look at the Colormap tab of Properties). The pixels on the streams_raster layer with streams have values of 1, 2 or 3; those with no streams have 0.
A pixel on the new raster must show "water" if EITHER the corresponding pixel on the LandUse layer = 5 OR the corresponding pixel on the streams_raster layer > 0. Coded, it looks like this:
To avoid typos, it's best to enter this by clicking on the buttons for "(", "=", etc. or double-clicking on the name of the layer in the Raster bands box; you'll have to type the "5" and the "0". (GeoTIFF files can contain more than one raster band, and "@1" at the end of the name indicates Band 1. In this example, our raster files have only one band.)
Click on OK and the new raster should appear. It has only two values, 1 for water, 0 for no water, so I used a colormap and changed the colour for 0 to white:
Creating the Distance from... raster
Here I hit a few problems!
The procedure should be the same as for the roads. Go to Raster > Analysis > Proximity. Select water_all as the input file and enter a name for the output file, eg. "DistanceFromWater.tif". Make sure "Dist units" are set to GEO (for GEOreferenced coordinates), check "Load into canvas..." and press OK. But this is what I get:
Blue shading should be close to water, red is far from water. That's worked in the middle of the map, but the edges are shown as close to water! What has happened?
I went back to the water_all layer and used View > Identify Features to investigate the values in the raster. In the lake and exactly on the streams gave 1; away from water in the middle of the map gave 0; but a few places showed "null (no data)". Trying that with the roads_raster layer gave only zeros and small numbers (the road IDs), no nulls.
The "no data" comes from the LandUse raster, for places with no land use specified. The problem is that "no data" is coded as a number: look at the Metadata tab of Properties and it says it's -3.40282e+38, and the default in Proximity is to treat all numbers other than zero as "water" (click on Help in the Proximity dialogue box to see details). We need to be a bit cleverer!
Go to Raster > Analysis > Proximity again, select water_all as the input, and specify an output file (I used DistanceFromWater2.tif). All the pixels with water in water_all are coded "1", so check "Values" and enter 1. Click OK.
The map now looks better. I simply set the style to colormap with 20 classes, but what's happened in the lower right corner? Using the Identify tool again shows that points in this corner are up to 180 km away from water (screen shot below), but the legend only goes up to 162 km. When determining the range to use for the legend, QGIS only looks at a sample of the pixels in a raster, and here it's evidently missed the highest values!
To fix that, I went back to the Colormap tab in Properties and changed the last number in the Value column to something much bigger than 180km - and I changed the label to match.
We want to add the distance from water to the attribute table for cameras, using the file that already has the distance from road data (I called it "cameras4.shp").
Activate the cameras4 and DistanceFromWater2 layers. Go to Plugins > Analyses > Point sampling tool. The "Layer containing sampling points" should be "cameras4". In the "Layers with fields/bands..." select all the camera fields and the DistanceFromWater2 layer.
Click on the Fields tab and change "DistanceFr" to something like "waterDist". Back in the General tab, put in a file name for the "Output point vector layer" (I used "cameras5.shp"), check "Add created layer to the TOC" and click OK.
Look at the Attribute table for the new layer:
Sort by waterDist; none is zero, but select the four lowest and see where they lie on the map in relation to streams:
| Updated 17 Dec 2012 by Mike Meredith|