Tag Archives: Numba

Create a slope layer from a DEM

Leave a reply

To calculate slope from a DEM the gdaldem utility can be used:

gdaldem slope -of KEA in_dem.kea out_slope.kea

There is an option to pass in a scale to convert the horizontal spacing from degrees to metres. However, if the DEM covers a range of latitudes it is desirable to convert the horizontal spacing on a per-pixel basis. To accomplish this I wrote a Python script, with RIOS used to read and write out data in blocks, making it efficient to process large datasets.

Implementation
As the slope calculation requires looping through a block of data, a task which is slow in pure Python, two strategies are used to improve the speed:

  1. Numba is used to compile the Python code – which is a lot faster than pure Python
  2. Fortran code is used to perform the actual slope calculation, compiled as a Python module using f2py

To provide a comparison of the time required for each implementation slope was calculated from a 1800 x 1800 pixel DEM using each method:

Implementation Time (s)
Pure Python 136.7
Numba 2.6
Fortran 2.1

Tests using Python 3.4 with Numba 0.13.4 running under OS X with 1.8 GHz Intel i5 processor. The gfortran Fortran compiler was used. Average runtime over 3 runs.

So the Fortran version is slightly faster than the Numba version but only slightly and required a lot more effort to implement. Both are significantly faster than the pure Python version. The code provides a nice example of applying more complicated algorithms to images with RIOS used to handle the I/O.

To account for noise in the DEM there is also a version of the slope calculation which will use least squares fitting to fit a plane over a window of pixels and calculate the slope from this. As the Python version requires calls to the NumPy linear fitting code there is no improvement using Numba. The Fortran version uses Accelerate under OS X or ATLAS under Linux for the plane fitting.

The code is available to download from https://github.com/MiXIL/calcSlopeDegrees

Usage
To calculate slope from a DEM where the horizontal sloping is in degrees:

calcSlopeDegrees.py --spacing_degrees inDEM.kea outSlope.kea

To calculate slope by fitting a plane over 9 x 9 window:

calcSlopeDegrees.py --plane_ls --window_size 9 inDEM.kea outSlope.kea