I have an elevation raster and I would like to get some basic terrain variables, such as slope, ruggedness, aspect or curvature.
The problem is that I don't seem to have the possibility to calculate curvature with the terrain analyst:
Can anyone tell me why this indicator isn't in the list? How can I calculate it otherwise?
Curvature is a complex terrain derivative to compute, the equation that you use depends on the resolution of your input data, as you have to ensure that the curvature results you compute can be distinguished from noise in the data.
A lot of research has been done recently on curvature calculations on high resolution LiDAR data which showed that a scaling break exists at around 2 or 3 metres resolution and above this point more different algorithms (which I am not as familiar with) need to be used. The best information about calculating topographic curvature probably comes from Hurst et al 2012 and the references therein.
The basic principle of curvature calculation, as with slope and aspect, is to pass a moving window over the elevation surface and fit the elevation values to a 6 term polynomial function, the coefficients of which will yield the slope, aspect and curvature of the centre cell of the moving window.
ArcGIS uses a 3x3 search window which will only yield good results in areas completely devoid of vegetation, which makes the tool fairly useless unless people are aware of this limitation, this may suggest why it is not present in QGIS.
The maths was derived originally (I think) in Evans (1980) and was simplified in a few pages in Principles of Geographical Information Systems (Amazon link) which I can recommend as a good guide to this kind of terrain analysis at a basic level.
One way to calculate curvature of a DEM is to convert the DEM into an ascii raster, read it into a numpy array and then perform the polynomial fitting on a moving window passing through the data. This is fairly easy to do, but is very slow to execute and needs a fair amount of optimization (these kind of operations often get ported to c++ to speed them up).
To perform the operation in QGIS you can use the GRASS plugin r.slope.aspect which is also limited by the 3x3 fixed window.
I realize this is not the simple answer you were doubtless hoping for, but I hope you understand that curvature is complex to derive in a meaningful way. All the best.
Evans, I. S. (1980), An integrated system of terrain analysis and slope mapping, Z. Geomorphol., 36, 274–295.
ESRI's version of Raster Analysis for calculating curvature might be helpful to develop a plugin for QGIS.
For each cell, a fourth-order polynomial of the form:
Z = Ax²y² + Bx²y + Cxy² + Dx² + Ey² + Fxy + Gx + Hy + I
is fit to a surface composed of a 3x3 window. The coefficients a, b, c, and so on, are calculated from this surface.
The relationships between the coefficients and the nine values of elevation for every cell numbered as shown on the diagram are as follows: Curvature values diagram Curvature values diagram
A = [(Z1 + Z3 + Z7 + Z9) / 4 - (Z2 + Z4 + Z6 + Z8) / 2 + Z5] / L4
B = [(Z1 + Z3 - Z7 - Z9) /4 - (Z2 - Z8) /2] / L3
C = [(-Z1 + Z3 - Z7 + Z9) /4 + (Z4 - Z6)] /2] / L3
D = [(Z4 + Z6) /2 - Z5] / L2
E = [(Z2 + Z8) /2 - Z5] / L2
F = (-Z1 + Z3 + Z7 - Z9) / 4L2
G = (-Z4 + Z6) / 2L
H = (Z2 - Z8) / 2L
I = Z5
The output of the Curvature tool is the second derivative of the surface—for example, the slope of the slope—such that:
Curvature = -2(D + E) * 100
Full Information and source:
http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00q90000000t000000
-2(D+E)*100 as (100/L2)*(3*Z5 - [Z2+Z4+Z6+Z8+Z5]) shows how to obtain this value as thrice the original value, Z5, minus a neighborhood sum Z2+Z4+Z6+Z8+Z5 (using a radius 1 circle), all rescaled by 100/L2. That's just three simple grid operations.
The curvature could be calculated using SAGA's module 'Terrain analysis - Morphometry ---> Slope, Aspect, Curvature'
The calculation could be done based on one of these algorithms:
LandSerf can do this .You can define windows size (3*3, 5*5, 7*7, 11*11,...) but it must be odd number. multiscale analyses. you can consider scale dependency of slope, Aspect and curvature. http://www.landserf.org/ LandSerf will fit co quadratic polynomial equation to specific predefined windows, but if you define big windows size like 50*50 it takes long time. It depends to your raster size and windows size which you define. Jo Wood wrote Landserf for his PhD thesis. It is written in Java.
Can be calculated in the different scale in LandSerf
You can also try free SAGA GIS (http://sourceforge.net/apps/trac/saga-gis/wiki) or TAS (http://www.uoguelph.ca/~hydrogeo/TAS/index.html).
Landserf gis has been written by J.Wood.you can analysis DEM, Aspect ,Slope, Plan curvature, profile curvature ,longitudinal curvature, and cross-sectional curvature,min,max, mean curvature, fractal dimension of DEM to different scale by defination of windows size .windows size should be odd.33,55, 77,99 11,11,....89*89 but it takes so long time.but result is very nice
