Wednesday 13 April 2011

Directional slope along changing directions: a Python plug-in for Quantum GIS

In a previous post (in Italian) I presented a stand-alone Python script whose purpose was the calculation of the directional slope along varying orientations. It can be used to calculate the slopes along directions parallel to local winds, glacier flows or other physical-environmental parameters with spatially changing orientations.

In this post, I introduce a new version, that can be used as a Quantum GIS plug-in. A similar plug-in is available for Grass, i.e. r.stream.slope. It is devoted specifically for hydrological research and considers eight possible directions for slope calculation, i.e. those corresponding to the lines joining the center of the analysed cell to the centers of the eight neighbouring cells, so that directions are discretised by angles of 45°.

In this QGIS plug-in, the orientations correspond to angles from 0° to 360°, without discretisations of their values.



Implementation

Download of the plug-in: VariableSlope.rar

Example data set: example_dataset.rar

This QGIS plug-in was created by incorporating previous Python code into the framework automatically created by the “Plugin Builder” plug-in in Quantum GIS, together with an interface created with Qt Designer.

It was developed and tested using QGIS 1.6.0 'Copiaco' on Windows Vista. It requires, among others, gdal and numpy modules, that should already be present in the default Python installation for Quantum GIS.

Unzip the plug-in in the QGIS plug-in directory and load it with the Python plug in installer. The plug-in name is “Directional slope along variable orientations”.

The plug-in window can be opened from Plugins -> Variable slope -> Variable slope.



Fig. 1. The plug-in window.


The input data sets are two rasters, one representing the elevations, the other the orientations (in degrees, from 0° to 360° or equivalently, 0° to 180° and -180° to 0°). They must share the same geometry and geographical parameters (i.e., cell sizes, number of rows and columns, coordinates of the lower-left cell).

These rasters should be already loaded in the QGIS project, in order to have them listed in the plug-in window.

Two slope calculation methods (central differences and Horn, 1981) are available. They are described in the following section.

The output format is the ARC/INFO ASCII GRID. Directional slopes are expressed by degrees (positive upwards, negative downwards). Note in Fig. 1 that the provided name for the output is without a suffix. Suffixes (simple.asc for central differences and horn.asc for Horn, 1981) will automatically be appended, depending on the chosen slope calculation methods.

If any user will encounter any problem running it, please email me at alberti.m65@gmail.com.



Theoretical bases for directional slope calculation
(Modified from http:\\www.malg.eu\directionalslope.php )

The calculation of the directional slope requires the choice of an analytical formula to evaluate it from raster data sets. This formula represents an adaptation to a discrete representation of analytical techniques for continuous and differentiable surfaces.

To determine the directional slope, we must first calculate the directional gradients along the two principal axes, x and y. From the values of dz/dx and dz/dy, one can determine the slope of the plane that approximates the surface at the local cell (i, j) using the formula (Neteler & Mitasova, 2008, eq. A.27; see also Geospatial Analysis - a comprehensive guide):

  directional slope (alpha) = arctan [(dz/dx) * sin(alpha) + (dz/dy) * cos(alpha)]

where alpha represents the orientation along which we calculate the slope. This value is the angle that the direction makes with the top of the map and varies from 0° to 360° clockwise (or equivalently, from -180° to 180°).

We can derive gradient by applying a 3x3 kernel to each cell (except those at the border). Figure 2 illustrates the adopted convention for cell indexing.



Fig. 2. Kernel for slope calculation.


Various methods are available for the estimation of the gradients along the two axes. A detailed analysis can be found in Burrough & McDonnell (1998). Two methods are implemented in this plug-in: the central differences and the Horn (1981) methods. For each gradient estimation, the former uses the values of two neighbouring cells, while the Horn (1981), suitable for rough surfaces, uses the values of six neighbouring cells.

Central differences
Gradients along the two axes are (Fig. 3):

  dz/dx = [ Elev(i, j + 1) - Elev(i, j - 1) ] / 2 cell_size

  dz/dy = [ Elev(i - 1, j) - Elev(i + 1, j) ] / 2 cell_size

In the case of cells at the grid edge, this method can be replaced by the first difference method, in which the central cell itself provides the value for the lateral, missing data.

Fig. 3. Cell weights in the central differences case.



Horn (1981) method
The gradient estimation along an axis derives from the values of six cells in a 3x3 kernel (Fig. 4). The weight applied to each cell depends on its position relative to the central cell. This method is the most suitable for rough surfaces. Formulas of the gradients are (Burrough & McDonnell, 1998, eqs. 8.9 and 8.10):

  dz/dx = [ (Elev(i-1, j + 1) + 2 Elev(i, j + 1) + Elev(i+1, j + 1)) - (Elev(i-1, j - 1) + 2 Elev(i, j - 1) + 

                 Elev(i+1, j - 1)] / 8 cell_size

  dz/dy = [ (Elev(i - 1, j-1) + 2 Elev(i - 1, j) + Elev(i - 1, j+1)) - (Elev(i + 1, j-1) + 2 Elev(i + 1, j) + 

                 Elev(i + 1, j+1)) ] / 8 cell_size


Fig. 4. Cell weights in the Horn (1981) case.