Creating a GIS layer for "Distance from..." in R

Back to home page

Updated 14 Nov 2017 : functions readShapeSpatial and writeShapeSpatial are now deprecated, so replaced with readOGR and writeOGR .

We recently ran a workshop on Geographical Information Systems (GIS) using Quantum GIS for ecologists and wildlife researchers. For many species, distance from water, a road, forest edge, or a settlement may be an important habitat variable.

For example, we may be using automatic cameras to investigate occupancy of sites by leopards. Probability of occupancy may depend on distance from the nearest road. Given vector layers with roads and camera locations, we want to do two things:

1. produce a layer with distance-from-nearest-road as a habitat layer so that we can calculate a probability of occupancy layer, and

2. extract distance-from-nearest-road for each of our cameras in order to model our data.

We had trouble doing this in QGIS, but succeeded with R. [I've since managed to do it in QGIS: see here .] On this page we'll see

how to do a distance-from-roads layer in R ,
how to extract the distance-from-roads table in R for all the cameras, and finally
how to do a layer and table for distance from water in R, where water includes both streams (lines) and water bodies such as lakes and reservoirs (polygons).

We'll use a toy data set which you can download here . Extract all the files to a folder named "data" inside your working folder on your hard drive ( not on the desktop). 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.

I assume you have the R statistical software package installed on your computer and you are somewhat familiar with R. Open R and go to File > Change dir..., then browse to the folder containing the data folder.

Load the packages we will need, installing any you don't already have:


        library(raster)
        

        library(maptools)
        

        library(rgdal)

(You don't need to load the rgdal package, but it must be installed on your computer, and loading it is the best way to check that it is installed.)

Produce a layer with distance-from-nearest-road

Read in the roads shape file and plot it. To indicate what corresponds to a vector layer and what's a raster, I'll put .v or .r on the end of the name. The two arguments of the readOGR function are the folder containing the shape files and the name of the shape file without the .shp extension.


        
         roads.v <- readShapeSpatial("roads.shp")
         

        
        roads.v <-
        
         readOGR("data", "roads")
        
        

        plot(roads.v)

If you have several raster layers, it's important that the pixels match up exactly . So we use an existing layer as a template when creating a new raster layer. We'll use the LandUse.asc file to create our template:


        LU.r <- raster("LandUse.asc")
        

        LU.r
        

        plot(LU.r)
        

        lines(roads.v)
        

        r <- LU.r  # this will be the template
        

        r[] <- NA  # assigns all values as NA
        

        summary(r) # shows you what you have: all NA's

Now we use rasterize to create a roads raster matching the template. This takes 15 secs on my old laptop:


        roads.r <- rasterize(roads.v, r, field=1)
        

        summary(roads.r)          
				# pixels crossed by a road have "1"


        

        plot(roads.r, add=TRUE)

The yellow raster-roads should fall almost exactly on the vector-roads.

Now we have our road data in a raster file, we can use distance to create the layer we want. This can be slow if you have a raster with lots of pixels - our example has 171,600 and took 22 secs on my laptop:


        roaddist.r <- distance(roads.r)
        

        class(roaddist.r)
        

        # Check:
        

        plot(roaddist.r)
        

        lines(roads.v)

The road lines should fall neatly into the area with zero distance from road.

Now we can save to a geo-tiff file (or an ESRI ascii file) for loading into QGIS:


        writeRaster(roaddist.r, "DistFromRoad.tiff", "GTiff")

Get distance-from-nearest-road for each of our cameras

Read in the cameras shape file and plot cameras and roads together to check:


       
        
         cams.v <- readShapeSpatial("cameras.shp")
         

        
        cams.v <-
        
         readOGR("data", "cameras")
        
        

        plot(roads.v)
        

        points(cams.v)

Now extract the distance-from-road information from roaddist.r for the camera locations:


       
        distFromRoad <- extract(roaddist.r, 
				cams.v)
        

        distFromRoad

You can add this information to the attribute table for the cameras and save it in a new shape file:


       
        cams2.v <- cams.v
        

        cams2.v$roadDist <- distFromRoad
        

        head(cams2.v)
        

        writeSpatialShape(cams2.v, "cameras2")

Now load the new cameras2.shp layer back into QGIS and check the attribute table. It has a column for roadDist. It also has columns for the coordinates of the centre of the pixel containing the camera.

Getting distance-from-water

Distance from permanent water is often important, but it's trickier - at least in our example - as water includes streams on the streams layer and water bodies (lakes, reservoirs, etc) on the land use layer.

We'll start with a raster including just the water body from the land use raster, where it is coded as 5. Pixels in water.r which have a value other than 5 are replaced with NA:


       
        water.r <- LU.r
        

        water.r[LU.r != 5] <- NA
        

        summary(water.r)
        

        plot(water.r)

Now load the streams layer and plot the water body and streams together to check:


       
        
         streams.v <- readShapeSpatial("streams.shp")
         

        
        streams.v <-
        
         readOGR("data", "streams")
        
        

        plot(streams.v)  # ok.
        

        plot(water.r)
        

        lines(streams.v) # not ok!

Fixing the CRS problem

To find out why the streams didn't appear when plotted with the water body, we look first at the extents, ie. the left and right, bottom and top coordinates:


       
        extent(streams.v)
        

        extent(water.r)

They have different extents. Look carefully and you'll see a black dot in the middle of the plot (circled in blue in the diagram above): those are the streams! Check the Coordinate Reference Systems (CRS) with proj4string:


       
        proj4string(streams.v)
        

        proj4string(water.r)

streams.v has no projection information: in fact it is unprojected (ie. long-lat) in degrees with WGS84 datum. They show up fine in QGIS with On The Fly Recalculation (OTFR) enabled, but don't work for analysis. This is the hazard of using OTFR!

To fix this, we have to make the streams CRS match the other layers. Either:

1. go to the streams layer in QGIS and Save As WGS84 / Pseudo Mercator (EPSG:3857) and read into R again, or

2. convert in R as shown below:

First give streams.v a proper projection attribute for longlat WGS84, EPSG code 4326, then transform to the same CRS as the water.r layer:


       
        proj4string(streams.v) <- CRS("+init=epsg:4326")
        

        streams2.v <- spTransform(streams.v, CRS(proj4string(water.r)))

Now try the plot again: it should work this time:


       
        plot(water.r)
        

        lines(streams2.v)

Back to the main story...

Create a raster from the stream layer just as we did for the roads earlier, and plot it to check:


       
        stream2.r <- rasterize(streams2.v, r, 
				field=1)
        

        summary(stream2.r)
        

        plot(water.r)
        

        plot(stream2.r, add=TRUE)

Next we combine the water body and stream data. The water body raster has 5's for water, and we change the pixels corresponding to streams to 1's. If we need to, we can distinguish streams from other water features; for distance-from-water we only use the fact that these pixels are not NA.


       
        water2.r <- water.r
        

        water2.r[!is.na(stream2.r)] <- 1
        

        summary(water2.r)


       
        

        plot(water2.r)

The rest of the procedure is the same as for the distance-from-road information.

For the distance-from-water layer, calculate distance from each pixel to the nearest non-NA pixel in water2.r:


       
        waterdist.r <- distance(water2.r)
        

        plot(waterdist.r)
        

        lines(streams2.v)


       
        writeRaster(waterdist.r, "DistFromWater.tiff", "GTiff")

To get the distance-from-water for each camera and add it to the cameras attribute table


       
        distFromWater <- extract(waterdist.r, 
				cams.v)
        

        distFromWater
        

        

        cams3.v <- cams2.v
        

        cams3.v$waterDist <- distFromWater
        

        head(cams3.v)
        

        

        
         writeSpatialShape(cams3.v, "cameras3")
         

        
        
         writeOGR(cams3.v, "data", "cameras3", driver="ESRI Shapefile")

Load cameras3.shp into QGIS and check the attribute table.

These general methods could be extended to other distance-from-... variables, such as distance from dwellings (shown on a vector layer of points) or distance from forest edge (polygons), either outwards from the forest or inwards into the forest.

Updated 14 Nov 2017 by Mike Meredith