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.
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:
- 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.
- 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.
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.
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)
The lowest part of the catchment is made up of the glacial lake and the adjacent drainage channel forming the valley with its slope flanks. Due to the lake the cumulative frequency of slopes shows a steep rise close to 0° and a bimodal distribution of slopes is evident. The steep slopes of the u-shaped valley account for the higher slope values and lift the average slope of this elevation zone to approx. 25°.
This elevation zone is marked by an increase of the average slope compared to the first one. However, this can mainly be attributed to the presence of the steep slopes on the valley flanks while the parts south to the lake account for lower slope values.
Following the next zone, the present elevation zone is in terms of the covered area the second largest one. Compared to the previous one, the average slope decreases and the whole frequency distribution shifts a bit to lower values. The 5% of pixels with slope values equal to or lower than 5° is primarily located on the contiguous glacier mass in the southeast. It is surrounded by abrasively shaped steep flanks. Notwithstanding this contrast, most parts of this elevation zone show intermediate slope values between 15° and 35°.
As in the previous elevation zone, the histogram of slope values is unimodal. Furthermore, the average slope is pretty similar. However, extreme values occur less often, thus, causing the frequency distribution to be narrower. Hence, the slope values in this zone are more homogeneous. Considering the map visualisation, this does not imply that there aren’t spots with contrasting slope values next to each other.
In this elevation zone, the distribution of slope values widens again. On the one hand, there are the glacial masses located in some kind of trough structure with slope values below 10°. On the other hand, at least 20% of the area can be considered as steep slopes with values above 40°.
The highest and smallest area of the catchment is characterised by an abundance of steep slopes. More than 40% of the area exhibit slope values of more than 40° causing the average slope to reach 35°. Frost weathering and retrograde erosive expansion of the cirque form ridges with almost vertical walls in this uppermost part. Considering the continuous rise of the average slope throughout the previous elevation zones, one might also be able to detect the typical concave structure of a cirque. Nevertheless, as high slope values in lower elevation zones proof, there is always a multitude of parallel processes and corresponding shapes present in each elevation zone. This is also evident from the bimodal distribution of slope values in the current zone. Geomorphological analyses may thus be supported by stratifying calculations according to altitude zones, but still no simple conclusions can be drawn.
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.
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).
Scripts
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 data.gov.at
# Select catchment area within Nationalpark Hohe Tauern
catchment_url = "https://www.salzburg.gv.at/ogd/3ccfbe49-782d-4bb1-9b3d-f38e45b44079/Gewaesser_Einzugsgebiete.json"
wget.download(catchment_url, os.path.join(os.environ["working_dir"], "input_data"))
gs.run_command(
"v.import",
input=os.path.join(
os.environ["working_dir"], "input_data", "Gewaesser_Einzugsgebiete.json"
),
output="catchments",
)
gs.run_command(
"v.extract",
input="catchments",
output="catchment_np",
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 data.gov.at
# 5m DEM not as WMTS available, hence download in ASCII GRID format, import & cleanup
# for WMS/WMTS services r.in.wms could otherwise used
dem_url = (
"https://www.salzburg.gv.at/ogd/d585e816-1a36-4c76-b2dc-6db487d22ca3/dgm5m.zip"
)
wget.download(dem_url, os.path.join(os.environ["working_dir"], "input_data"))
zip = zipfile.ZipFile(
os.path.join(os.environ["working_dir"], "input_data", "dgm5m.zip"), "r"
)
zip.extractall(os.path.join(os.environ["working_dir"], "input_data"))
os.remove(os.path.join(os.environ["working_dir"], "input_data", "dgm5m.zip"))
gs.run_command(
"r.import",
input=os.path.join(os.environ["working_dir"], "input_data", "dgm5m", "dgm5m.asc"),
output="dem_np2",
extent="region",
resample="bicubic",
)
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")
gs.run_command(
"r.slope.aspect",
elevation="dem_np",
slope="slope_np",
format="degrees",
pcurvature="pcurve_np",
tcurvature="tcurve_np",
precision="DCELL",
overwrite=True,
)
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
gs.run_command(
"r.geomorphon",
elevation="dem_np",
forms="geomorph{}_np".format(search_rad),
search=rad_in_map_units,
)
# Aux_slope_hist function
def slope_hist(name, res=5, omit_info=False, exp_dir=export_dir):
# read slope vals
slope_raster = rasterio.open(os.path.join(exp_dir, name))
slope_vals = slope_raster.read(1)
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()
ax1.hist(
slope_vals,
bins=bins,
color="#336699",
alpha=0.9,
rwidth=0.9,
weights=(np.zeros_like(slope_vals) + 1) / slope_vals.size,
)
ax1.set_xlabel("slope (°)")
ax1.set_ylabel("relative frequency")
ax1.set_xticks(bins)
ax1.set_xlim(0, 90)
ax1.set_title(
f"average slope (°): {arr_mean}\n median slope (°): {arr_median}\n maximum slope (°): {arr_max}",
loc="right",
fontsize=10,
)
if not omit_stats:
ax1.set_title(
f"area(km²): {arr_area}\nresolution: {res}x{res}m\nnumber of pixels: {arr_size}",
loc="left",
fontsize=10,
)
# plotting cumulative frequency
steps = np.append(np.arange(0, 100, 0.001), [np.inf])
ax2.hist(
slope_vals,
bins=steps,
color="darkred",
histtype="step",
linewidth=1.5,
cumulative=True,
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)
slope_hist("slope_np.tif")
# 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)
gs.run_command(
"r.out.gdal",
input=name,
output=os.path.join(exp_dir, "{}.tif".format(name)),
format="GTiff",
overwrite=True,
)
gs.run_command("r.mask", flags="r")
export_dir = os.path.join(os.environ["working_dir"], "results")
export("dem_np")
export("slope_np")
export("pcurve_np")
export("tcurve_np")
export("hillshade_np")
export("geomorph250_np")
export("geomorph1000_np")
# 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
)
reclassify.stdin.write("\n".join(intervals).encode())
reclassify.stdin.close()
reclassify.wait()
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)
gs.run_command(
"r.out.gdal",
input="slope_np",
output=os.path.join(export_dir, "slope_np_{}.tif".format(cat)),
format="GTiff",
overwrite=True,
)
gs.run_command("r.mask", flags="r")
slope_hist("slope_np_{}.tif".format(cat))
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"]
)
reclassify.stdin.write("\n".join(intervals).encode())
reclassify.stdin.close()
reclassify.wait()
gs.mapcalc(
"{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))
gs.run_command(
"r.resamp.stats",
input="dem_np",
output="dem_{}_{}".format(res, met),
method=met,
flags="w",
)
gs.run_command("r.mask", vector="catchment_np", overwrite=True)
gs.run_command(
"r.slope.aspect",
elevation="dem_{}_{}".format(res, met),
precision="DCELL",
slope="slope_{}_{}".format(res, met),
format="degrees",
)
gs.run_command("g.region", raster="dem_np", res=5, flags="p", grow=1)
gs.run_command(
"r.resamp.interp",
input="dem_{}_{}".format(res, met),
output="dem_{}_{}".format(res, met),
method="nearest",
overwrite=True,
)
gs.run_command(
"r.resamp.interp",
input="slope_{}_{}".format(res, met),
output="slope_{}_{}".format(res, met),
method="nearest",
overwrite=True,
)
gs.mapcalc(
"{0} = {1}-{2}".format(
"dem_diff_{}_{}".format(res, met),
"dem_{}_{}".format(res, met),
"dem_np",
)
)
gs.mapcalc(
"{0} = {1}-{2}".format(
"slope_diff_{}_{}".format(res, met),
"slope_{}_{}".format(res, met),
"slope_np",
)
)
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
all_res_stats.append(res_stats)
# 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)")
ax.set_ylabel(y_label)
ax.legend(
handles=ax.get_legend_handles_labels()[0][::-1],
labels=elev_labels[::-1],
bbox_to_anchor=(1.02, 0.6),
frameon=False,
)
fig.savefig(os.path.join(exp_dir, f"res_{var}"), bbox_inches="tight", dpi=300)
else:
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)
ax1.set_ylim(val_range)
ax1.set_ylabel(y_label)
ax1.get_xaxis().set_visible(False)
ax1.text(
0.65,
0.9,
"resampling method: average",
transform=ax1.transAxes,
fontsize=8,
)
ax2.plot("res", "value", color=line_c, linewidth=1.5, data=data2)
ax2.set_ylim(val_range)
ax2.set_xlabel("resolution of resampled raster (m)")
ax2.set_ylabel(y_label)
ax2.legend(
handles=ax2.get_legend_handles_labels()[0][::-1],
labels=elev_labels[::-1],
bbox_to_anchor=(1.4, 1.45),
frameon=False,
)
ax2.text(
0.65,
0.9,
"resampling method: median",
transform=ax2.transAxes,
fontsize=8,
)
fig.savefig(
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)
