Masters programme | E-portfolio
Semester I & II
Spatial Analyses

Terrain analyses

Qualitative and quantitative assessment of terrain in a mountain catchment area. Analysis of the influence of different DEM resolutions on the calculation of elevations and slopes.

Table of Contents

Introduction & Study area

The following analyses focus on a water catchment area at the southern border of the state of Salzburg, Austria. The catchment leads up to the Großvenediger mountain and comprises an area of approx. 22 km2.

AoI-Delineation based on official data on catchments, orthophoto as baseimage derived from

As is evident from the map above, slight inaccuracies are inherent to the catchment delineation. However, for the purpose of the present analyses this generalised demarcation is sufficient. A Digital Elevation Model (DEM) with 5m spatial resolution is available for the area and builds the foundation for the subsequent analyses:

  1. Describing the terrain of the given area in a detailed manner using standard geomorphometric measures such as elevation, slope and curvature but also more advanced methods to capture geomorphological forms.
  2. Analysing the influence of different resampling resolutions on the quality of elevation data and derived slope information.

To get a bit more familiar with the Area of Interest (AoI), the scene view and the hillshade representation below provide a visual impression of the glacially formed catchment. Apparently, the main part of the catchment constitutes a glacial cirque with a small lake at its low point. The trough valley that follows from there impressively shows the erosive effect of former glacier masses. The recent glacier tongues are distinctive in the hillshade due to their surface homogeneity.

Hillshade with sun evelation angle of 30° and azimuth of 270° (east)

A. General geomorphometric analyses

A.1 Elevation

The lowest point of the catchment is located in the through valley in the northwest. Contrary, the highest point at 3655m a.s.l. is right below the summit of the Großvenediger (south east). The maximum difference in elevation is approx. 1700m.

Elevations based on the DEM5

A.2 Slope

The average slope in the AoI amounts to 26° and some parts are nearly vertical (maximum slope 83°). Throughout the whole catchment, the histogram of slopes shows a unimodal, right skewed distribution. In general, low slope values are calculated for the valley, the glacial lake and the glacier masses themselves. Next to these areas, steep slopes are visible and often a product of some sort of erosional activity.  Without further investigation, it can be conjectured that steep slopes tend to occur more often at higher altitudes. However, it is also evident that steep slopes are not limited to these areas. In order to investigate this a bit more in depth, slopes differentiated according to elevation zones may be analysed. Therefore a reclassification of elevations with four equally spaced elevation zones of 250m and two extended ones at the bottom and top of the range was carried out. For each elevation zone at least some spots with slope angles close to 80° exist. It is notable that the observation of right skewed frequency distributions holds equally true for all elevation zones, even though the multimodality of some of the distributions obscures this. The assumption of increasing slopes at higher altitudes can be confirmed partially. More details are provided below.

Slopes overlaid on the orthophoto (left) and compared to elevations (right)

Histogram of slopes including summarising statistics

A.3 Curvature

The curvature represents the second derivative of a surface, thus, describing the change in slope. Geometrically spoken, determining the curvature corresponds to fitting a circle to the shape of the surface at a given point. The smaller the radius of the circle that best approximates the shape of the surface, the higher the curvature. Speaking of a circle in a 1-dimensional sense stems from the fact that curvature is usually calculated along a specific direction, i.e. a certain plane cutting through the surface. This plane may be oriented in line with the steepest slope, then referred to as profile curvature. Alternatively, tangential curvature uses the plane perpendicular to the slope as a reference line. Curvature may also be expressed as the mean curvature along both perpendicular directions corresponding to the idea of fitting a sphere to the surface. However, calculating curvature along each direction individually may be more useful since profile and tangential curvature are of different significance depending on the purpose of the analysis. To characterise the acceleration and deceleration of flow down the surface by force of gravity, profile curvature may be applied. Tangential curvature is suited to characterise convergence and divergence of flow across the surface (for extended explanations and visualisations go here).

Below both profile and tangential curvature are presented with their magnitude classified into 5 intervals. Deep blue areas are concavely shaped with radii of curvature of less than 50m. Contrary, rad areas represent convex forms. Using the swipe tool to compare both measures of curvature, one immediately sees that tangential curvature dominates. All the narrow channel structures running down the slope contribute to tangential curvature but not to profile curvature. Areas of strong (positive or negative) tangential curvature correspond to areas of steep slopes. This is not equally true for profile curvature. Clearly visible in the eastern part of the catchment, profile curvature allows to delineate the convexly shaped mountain ridges. Furthermore, the shore of the central lake represents the case of concavely shaped footslopes, which are equally well delineated. The glaciers appear as low curvature areas along both directions. Their radius of curvature do not exceed 100m. Similar to the hillshade representation, they are identifiable as larger uniform areas in the maps below.

Profile curvature (left) and tangential curvature (right)

A.4 Automated mapping of forms

As indicated in the previous descriptions, geomorphometric measures may be utilised to identify geomorphological forms. Peaks are unique in the way that they are surrounded by areas of lower elevation, ridges are concavely shaped, footslopes are convex forms and so on. To automate the delineation of these and other terrain forms the relatively new concept of geomorphons may be applied. This concept describes the visibility neighbourhood of a pixel in form of a 8-tuple pattern of relative elevations. The tuple {+,-,-,-,-,-,-,0}, for example, says that six of the neighbouring pixels are at lower elevations, one at an equal one and one at a higher one. This tuple may then be automatically classified as representing the terrestrial form of a ridge. The finite number of 498 possible different tuples (geomorphons) constitutes a comprehensive set of morphological terrain types including the 10 most popular ones as shown below.

Forms represented by geomorphons, Jasiewicz & Stepinski 2013

The identification of forms is, of course, a scale-dependent issue. The concept of geomorphons accounts for this by using self-adapted neighbourhoods instead of relying on fixed-sized neighbourhoods. The neighbouring pixels forming the 8-tuple are determined within a prespecified search distance of the focus pixel based on the line-of-sight principle. Therefore, choosing an adequate search radius forms one critical point during the parametrisation of the algorithm. Higher search radii enable the identification of larger landforms, smaller ones allow local scale delineations of large-scale patterns. There is a variety of other parameters – including at least the flatness threshold – which may be modified to ensure a proper delineation of forms. However, exploring the broad range of parametrisation possibilities aiming at reaching the best-suited results is not the focus of the present analyses. 

The figures in the following are only intended to provide a first insight into the possibilities for automated mapping of landforms based on DEM data. The results of applying an implementation of the geomorphons algorithm under two different search radii are presented. The maximum search radius of 1000m obviously shows a more uniform delineation of landforms over larger areas. The valley gains in width compared to the classification based on the 250m radius, as with the larger search radius for more focus pixels all visible neighbours appear to be higher. Conversely, the line of sight principle means that the delineation of ridges, for example, changes little. Areas around a ridge are not identified as part of the ridge even with a larger search radius, as the higher points of the nearby ridge are still the closest visible neighbours. While parts of the glacial lake or its shore are identified as footslopes with the smaller search radius, higher points are increasingly found within the larger search radius leading to classification as a valley area or even pit.

Form identification based on search radii of 250m (left) and 1000m (right) 

B. Influence of resampling resolution on the derivation of elevation & slope information

Testing the accuracy of elevation and slope data as a function of resampling resolution is a necessary prerequisite for choosing appropriate DEM input data for all kinds of analyses. The tolerable error in the DEM input data, of course, varies with regard to the scale of interest that applies for the analysis to be carried out. But as soon as one is able to specify this acceptable error level, it might be advisable to go with the coarsest resolution DEM that still meets the requirements. This not only saves computing resources and time but also financial resources if money would otherwise be spent unnecessarily on the acquisition of higher-resolution input data. In the present case with open data availability for the relatively small catchment area, these considerations do not play a role, but for analyses in other larger areas such considerations are important.

In order to balance the aim of obtaining reliable data on the accuracy of geomorphometric measures for a range of different resampling resolutions one the one hand and the practical need of reducing the computational effort on the other hand, the following analysis design was applied: The original DEM5 was resampled more than 20 times to get DEMs with resolutions between 5m and 500m. Resampling was done using the average method. A selection of four of the resampled DEMs is shown below. For each of the resampled DEMs the same subsequent calculations were performed. First, slope values were calculated and, second, the mean absolute error (MAE) for both elevation and slope was determined by performing a comparison to the original DEM5 and the derived slope dataset, respectively. The resulting quasi-continous dataset on the resampling resolution dependent accuracy of elevations and slopes can than be visualised as line plots.

Unsurprisingly, the MAE of elevation increases linearly with decreasing resolutions of the resampled raster up to the maximum MAE of approx. 55m for the DEM500. When differentiating the MAE according to the elevation zones, it becomes apparent that the error is highest for the elevation zone 3250-3750m, followed by the zones of altitudes between 1750-2500m. The strong average slopes in these zones explain this pattern. Fluctuations in the course of the curves from resolution of about 300m upwards are particularly noticeable for the smaller zones. At these resolutions only a few pixels cover the smaller elevation zones. At a resolution of 500m, the entire AoI is covered by only 81 pixels. This also causes the mean elevation per zone to deviate from their true value.

Regarding slope measures, inaccuracies also rise with decreasing resampling resolutions. However, different to the evolution of elevation errors, the slope MAE does not show a linear pattern. The curves flatten at lower resampling resolutions since MAE are already well above 10°. The degradation of slope information is worse when moving from a DEM5 to a DEM10 compared to the difference between a DEM400 and DEM500. It is interesting to see, that resampling the DEM5 to a DEM10 already introduces a MAE of about 2.5°. Of course, depending on the application, the single slope pixel values and the corresponding MAE may not be as important. If instead the average slope for a whole area is of interest, the introduced error is a bit lower as the figure representing the average slope per elevation zone indicates. In both cases, however, the error is particular large within high curvature areas. This is particularly evident when representing the information not in a spatially aggregated form of a line plot but as a map. Underneath the following plots, this is done for a selection of three resampling resolutions. For slope biases, an almost perfect positive correlation with profile curvature is obvious (compare to section A.3).

Finally, in order to further analyse the influence of different resampling methods on the derivation of geomorphometric measures, the previous analyses were repeated using the median instead of the average as a resampling method. Note that choosing common resampling methods such as „nearest“ , „bilinear“ or „bicubic“ does not make a lot of sense in the present case of resampling to rasters with resolutions that are coarser than the source raster by orders of magnitude. Applying the median method for resampling, the results turned out out to be very similar to the ones already presented. Therefore, further visualisations are omitted.


All analyses presented above were carried out in a GRASS GIS environment. Histogram visualisations were python-based and maps were produced with the help of ArcGIS. The main scripts for following along are shown below.

Setup & Data import

As mentioned above, projection problems exist when importing the DEM. These are probably due to an unreadable projection file as the attempt to create a mapset on the basis of the file using g.proj fails. Also when testing the data import with QGIS, the automatic recognition of the coordinate reference system (CRS) fails. Using manual specification of the CRS (here EPSG:31258), the data is georeferenced correctly. While this manual georeferencing works in QGIS, the corresponding standard approach in GRASS GIS (create new mapset with source CRS -> import data into mapset -> r.proj to import into mapset with target CRS) does not work for reasons that cannot be explained further.

					## Project setup - shell commands
# grass78 -c epsg:25833 "grassdata/terrain_analyses" -e
# grass78 "grassdata/terrain_analyses/PERMANENT"

# working_dir='/home/felix/Desktop/assignment_terrain'
# export working_dir

## open python shell within grass project
## import packages
# python3
from tqdm import tqdm
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import zipfile as zipfile
import wget
import rasterio
import grass.script as gs
import os

# Import data
# A. Water catchment areas from
# Select catchment area within Nationalpark Hohe Tauern
catchment_url = "", os.path.join(os.environ["working_dir"], "input_data"))
        os.environ["working_dir"], "input_data", "Gewaesser_Einzugsgebiete.json"
    where="HZB_CODE = '2  8272 19  1 0 0 0 0'",
gs.run_command("g.region", vector="catchment_np", res=5, flags="ap", grow=200)

# B. DEM (5m res) from
# 5m DEM not as WMTS available, hence download in ASCII GRID format, import & cleanup
# for WMS/WMTS services could otherwise used
dem_url = (
), os.path.join(os.environ["working_dir"], "input_data"))
zip = zipfile.ZipFile(
    os.path.join(os.environ["working_dir"], "input_data", ""), "r"
zip.extractall(os.path.join(os.environ["working_dir"], "input_data"))
os.remove(os.path.join(os.environ["working_dir"], "input_data", ""))
    input=os.path.join(os.environ["working_dir"], "input_data", "dgm5m", "dgm5m.asc"),
del zip

Analyses - part A

					# Analyses I: Characterise AOI terrain using high-res DEM
# Visualising aoi using geomorphometric measures (slope, curvature, etc)

# Calculate geomorphometric measures
gs.run_command("r.relief", input="dem_np", output="hillshade_np")

for search_rad in [250, 1000]:
    idx_res = [x for x in gs.parse_command("g.region", flags="g").keys()].index("nsres")
    res = int([x for x in gs.parse_command("g.region", flags="g").values()][idx_res])
    rad_in_map_units = search_rad / res

# Aux_slope_hist function
def slope_hist(name, res=5, omit_info=False, exp_dir=export_dir):
    # read slope vals
    slope_raster =, name))
    slope_vals =
    slope_vals = slope_vals[~np.isnan(slope_vals)]
    # calc descriptive statistics
    arr_size = slope_vals.size
    arr_area = np.round(arr_size * res * res / 1000000, 2)
    arr_mean = np.round(np.mean(slope_vals), 2)
    arr_median = np.round(np.median(slope_vals), 2)
    arr_max = np.round(np.max(slope_vals), 2)
    # set fig layout
    fig, ax1 = plt.subplots()
    ax2 = ax1.twinx()
    # plotting relative frequency
    bins = np.arange(0, 95, 5).tolist()
        weights=(np.zeros_like(slope_vals) + 1) / slope_vals.size,
    ax1.set_xlabel("slope (°)")
    ax1.set_ylabel("relative frequency")
    ax1.set_xlim(0, 90)
        f"average slope (°): {arr_mean}\n median slope (°): {arr_median}\n maximum slope (°): {arr_max}",
    if not omit_stats:
            f"area(km²): {arr_area}\nresolution: {res}x{res}m\nnumber of pixels: {arr_size}",
    # plotting cumulative frequency
    steps = np.append(np.arange(0, 100, 0.001), [np.inf])
        weights=(np.zeros_like(slope_vals) + 1) / slope_vals.size,
    ax2.set_ylabel("cumulative frequency")
    fig.savefig(os.path.join(exp_dir, f'hist_{name.split(".")[0]}'), dpi=300)


# Export data for further vis
# Aux_export function
def export(name, exp_dir=export_dir):
    gs.run_command("r.mask", vector="catchment_np", overwrite=True)
        output=os.path.join(exp_dir, "{}.tif".format(name)),
    gs.run_command("r.mask", flags="r")

export_dir = os.path.join(os.environ["working_dir"], "results")

# Calculate slope per elevation zone
# Reclassify dem - elevation zones of 250m intervals
intervals = [
    "1750 thru 2250 = 17502250",
    "2250 thru 2500 = 22502500",
    "2500 thru 2750 = 25002750",
    "2750 thru 3000 = 27503000",
    "3000 thru 3250 = 30003250",
    "3250 thru 3750 = 32503750",

gs.run_command("r.mask", vector="catchment_np", overwrite=True)
reclassify = gs.feed_command(
    "r.reclass", input="dem_np", rules="-", output="dem_reclass_5", overwrite=True
gs.mapcalc("dem_reclass_5 = int(dem_reclass_5)", overwrite=True)
gs.run_command("r.mask", flags="r")

for cat in [x.rsplit(" ", 1)[1] for x in intervals]:
    gs.run_command("r.mask", raster="dem_reclass_5", maskcats=int(cat), overwrite=True)
        output=os.path.join(export_dir, "slope_np_{}.tif".format(cat)),
    gs.run_command("r.mask", flags="r")

Analyses - part B

					# Analyses II: Evaluating influence of resampling with different res
# Define tested resampling resolutions
resolutions = np.arange(10, 50, 5).tolist()
resolutions.extend(np.arange(50, 100, 10))
resolutions.extend(np.arange(100, 505, 20))

# Define tested resampling methods
methods = ["average", "median"]

# Perform resampling, slope calculation and diff calculation to DEM5/Slope5
res_specific_layers = {}
for res in tqdm(resolutions):
    single_res_layers = {}
    single_res_layers["lyr_reclass"] = "dem_reclass_{}".format(res)
    gs.run_command("g.region", raster="dem_np", res=res, flags="p", grow=1)
    reclassify = gs.feed_command(
        "r.reclass", input="dem_np", rules="-", output=single_res_layers["lyr_reclass"]
        "{0} = int({0})".format(single_res_layers["lyr_reclass"]), overwrite=True
    single_res_layers["lyr_dem"], single_res_layers["lyr_slope"] = [], []
    single_res_layers["lyr_dem_diff"], single_res_layers["lyr_slope_diff"] = [], []
    for met in methods:
        gs.run_command("g.region", raster="dem_np", res=res, flags="p", grow=1)
        single_res_layers["lyr_dem"].append("dem_{}_{}".format(res, met))
        single_res_layers["lyr_slope"].append("slope_{}_{}".format(res, met))
        single_res_layers["lyr_dem_diff"].append("dem_diff_{}_{}".format(res, met))
        single_res_layers["lyr_slope_diff"].append("slope_diff_{}_{}".format(res, met))
            output="dem_{}_{}".format(res, met),
        gs.run_command("r.mask", vector="catchment_np", overwrite=True)
            elevation="dem_{}_{}".format(res, met),
            slope="slope_{}_{}".format(res, met),
        gs.run_command("g.region", raster="dem_np", res=5, flags="p", grow=1)
            input="dem_{}_{}".format(res, met),
            output="dem_{}_{}".format(res, met),
            input="slope_{}_{}".format(res, met),
            output="slope_{}_{}".format(res, met),
            "{0} = {1}-{2}".format(
                "dem_diff_{}_{}".format(res, met),
                "dem_{}_{}".format(res, met),
            "{0} = {1}-{2}".format(
                "slope_diff_{}_{}".format(res, met),
                "slope_{}_{}".format(res, met),
        gs.run_command("r.mask", flags="r")
    res_specific_layers[res] = single_res_layers

# Calculate stats per zone
all_res_stats = []
gs.run_command("r.mask", raster="dem_np", overwrite=True)

for res, specific_layers in tqdm(res_specific_layers.items()):
    # get layers per resolution & export them
    gs.run_command("g.region", raster="dem_np", res=5, flags="p", grow=1)
    reclass_layer = [layer for layer in specific_layers.values()][0]
    terrain_layers = [layer for layer in specific_layers.values()][1:]
    terrain_layers = [layer for sublist in terrain_layers for layer in sublist]
    list(map(export, terrain_layers))
    gs.run_command("g.region", raster="dem_np", res=5, flags="p", grow=1)
    for lyr in terrain_layers:
        res_stats = {}
        res_stats["res"] = res
        res_stats["metric"] = lyr.rsplit("_", 2)[0]
        res_stats["resample"] = lyr.rsplit("_", 2)[-1]
        # calculate mean abs error (MAE) per elevation zone using reclass layer
        lyr_zonal_stats = gs.parse_command(
            "r.univar", map=lyr, zones=reclass_layer, flags="t"
        idx_mae = [i.split("|") for i in lyr_zonal_stats.keys()][0].index("mean_of_abs")
        lyr_zonal_means = [i.split("|")[idx_mae] for i in lyr_zonal_stats.keys()][1:]
        lyr_zonal_means = [float(val) for val in lyr_zonal_means]
        zones = gs.parse_command("r.describe", map=reclass_layer, flags="1nd").keys()
        zones = ["elevzone_{}".format(zone) for zone in zones]
        for zone_name, zone_mean in zip(zones, lyr_zonal_means):
            res_stats[zone_name] = zone_mean

# Merge stat data & export as csv
res_stats = pd.DataFrame(all_res_stats)
res_stats = pd.melt(
    res_stats, id_vars=["res", "metric", "resample"], var_name="elev_zone"
res_stats["elevation"] = [int(x.rsplit("_")[1]) for x in res_stats["elev_zone"]]
res_stats = res_stats.sort_values(["res", "elevation"])
res_stats.to_csv(os.path.join(export_dir, "res_stats.csv"))

# Define aux_function for plotting resolution dependent accuracies
def resolution_plots(var, y_label, compare_resample_method=False, exp_dir=export_dir):
    terrain_colors = np.array(
            [38, 115, 0],
            [153, 230, 0],
            [255, 230, 0],
            [230, 189, 115],
            [168, 112, 0],
            [212, 141, 107],
    terrain_colors = terrain_colors / 255
    elev_labels = [
        str(x)[0:4] + "m - " + str(x)[4:] + "m"
        for x in np.unique(res_stats["elevation"])
    if not compare_resample_method:
        fig, ax = plt.subplots()
        for i, zone in enumerate(np.unique(res_stats["elev_zone"])):
            data = res_stats[res_stats["metric"] == var]
            data = data[(data["elev_zone"] == zone) & (data["resample"] == "average")]
            line_c = terrain_colors[i]
            ax.plot("res", "value", color=line_c, linewidth=1.5, data=data)
            ax.set_xlabel("resolution of resampled raster (m)")
                bbox_to_anchor=(1.02, 0.6),
        fig.savefig(os.path.join(exp_dir, f"res_{var}"), bbox_inches="tight", dpi=300)
        val_range = [
            min(res_stats[res_stats["metric"] == var]["value"]),
            max(res_stats[res_stats["metric"] == var]["value"]),
        fig, [ax1, ax2] = plt.subplots(nrows=2)
        for i, zone in enumerate(np.unique(res_stats["elev_zone"])):
            data = res_stats[res_stats["metric"] == var]
            data1 = data[(data["elev_zone"] == zone) & (data["resample"] == "average")]
            data2 = data[(data["elev_zone"] == zone) & (data["resample"] == "median")]
            line_c = terrain_colors[i]
            ax1.plot("res", "value", color=line_c, linewidth=1.5, data=data1)
                "resampling method: average",
            ax2.plot("res", "value", color=line_c, linewidth=1.5, data=data2)
            ax2.set_xlabel("resolution of resampled raster (m)")
                bbox_to_anchor=(1.4, 1.45),
                "resampling method: median",
            os.path.join(exp_dir, f"res_II_{var}"), bbox_inches="tight", dpi=300

resolution_plots("dem", "mean elevation (m)")
resolution_plots("slope", "mean slope (°)")
resolution_plots("slope", "mean slope (°)", compare_resample_method=True)
resolution_plots("dem_diff", "elevation - mean absolute error (m)")
resolution_plots("slope_diff", "slope - mean absolute error (°)")

# Plot slope histograms for a selection of resampled rasters
slope_hist("slope_10_average.tif", omit_info=True)
slope_hist("slope_100_average.tif", omit_info=True)
slope_hist("slope_500_average.tif", omit_info=True)