River Catchment – Tasks
In this exercise, a river catchment area is analyzed regarding its slope and curvature. The following tasks are performed:
- Slope Visualization, Statistics and Histogram for the whole catchment area
- Analyze the slope and it’s statistics for different elevation zones within the catchment area
- Comparing slope statistics for 10m and 100m spatial resolution
- Calculating curvature without the “curvature” function in ArcGIS Pro
Study Area
As a study area, the Neubach catchment area was chosen. This little river originates in the mountains of the Dachstein massif and its slopes, before feeding into the Lammer which eventually feeds into the Salzach near Golling an der Salzach.
For the analysis, seven separated catchment areas were unified. Those areas all have the same exit point, where the Neubach drains into the Lammer. The Catchment Code (HZB_CODE) unifies those areas under the code 2-8272231-14.
The study area ranges from 850m in elevation at the drainage point until just over 2000m on the flank of the Dachstein massif.
Slope – complete Catchment Area
From the DEM, the slope per pixel is calculated in degrees. Also, a histogram is created in Python (the code can be found at the end of the post) that shows the slope of the whole study area in 5° steps as well as the mean value.
Click right on the slideshow to see the Histogram
Elevation Zone Creation
Next up, the catchment area is devided into the following four elevation zones:
- <1000m
- 1000 – 1250m
- 1250 – 1500m
- 1500 – 1750m
- > 1750m
Slope per Elevation Zone
Now, the slope average and histogram are calculated for each elevation zone.
The average slope for the different elevation zones are as follows:
- <1000m: 19.3°
- 1000 – 1250m: 25.6°
- 1250 – 1500m: 32.3°
- 1500 – 1750m: 25.5
- > 1750m: 34.9
Interestingly, the slope does not always increase the higher up the mountain you go. Zone 4 (1500 – 1750m) has a lower average slope with 25.5°, most likely doe to the flatter area in the northern zone as well as the mountain ridge in the southern area of this zone.
This observation is supported if a profile graph is drawn across all elevation zones.
Resampling to 100x100m Slope Pixels
Now, the slope .tiff is resampled to a pixel size of 100 by 100 meters to see if the slope averages change depending on the spatial resolution of the underlaying data.
Comparing the average slopes of the different elevation zones and spatial resolutions shows that the average slope does not drastically change. As a resampling algorithm to downsample the images, the bilinear method was chosen. The low difference indicates that this method yielded good results. Surely, if another resampling method such as Nearest Neighbor would have been chosen, the result would have been much less acurate.
Curvature calculation without ArcGIS curvature function
Calculating the curvature of an area is very simple with the curvature function in ArcGIS Pro, but there is another easy workaround. Since we already have the slope layer, which indicates the slope of each pixel, we can run the slope function again on the same dataset.
After the first execution, the values show the slope of the pixel. Since the curvature is the second derivative of the slope (slope-of-the-sope), the result is similar to the curvature function.
The curvature function results in pixel values from ca. -4 to 4, while the slope of the slope is still in degrees. The degrees can be scaled to the same values, which results in the same output. Unfortunately, the visualization is different for the following two curvature representations, but the areas of conformity are still visible. Clearly, the two images show the same real-life phenomenon.
Finally, as anounced, here is the code used to create all histograms in this post.
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
# Also open large images
Image.MAX_IMAGE_PIXELS = None
paths_slopes_resample=[
"slope_1000_resampled.tif",
"slope_1250_resampled.tif",
"slope_1500_resampled.tif",
"slope_1750_resampled.tif",
"slope_2050_resampled.tif",
]
paths_slopes=[
"slope_1000.tif",
"slope_1250.tif",
"slope_1500.tif",
"slope_1750.tif",
"slope_2050.tif",
]
path_slope_10 = ["slope_complete.tif"]
path_slope_100 = ["slope_complete_100.tif"]
# create clean title names from path
def extract_title(path):
name = path[:-4]
if "resample" in name:
name = name.replace("_resampled",", 100x100m")
else:
name = name+", 10x10m"
name = name.replace("_"," ")
name = name.replace("s","S")
name = name.replace("r","R")
name = name.replace("Slope","Slope in Degrees,")
name2 = name[:22]+"m Elevation Zone"+name[22:]
return(name2)
# Opening Image path, converting to numpy array, deleting nan values
def to_np(path):
im = Image.open(path)
im_array = np.array(im)
# round values to 2 decimal places
im_array = np.round(im_array, 2)
#Remove 0 (null) values
im_array[im_array == 0] = np.nan
# Remove nan values
im_array = im_array[~np.isnan(im_array)]
im_array = im_array[np.logical_not(np.isnan(im_array))]
#print((im_array))
return im_array
def graph(im_array,path,name):
# setting bns in 5 degree steps
bins=[0,5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80,85,90]
# draw histogram
n, bins, patches = plt.hist(im_array[~np.isnan(im_array)], bins=bins,color="royalblue", alpha=0.75, rwidth=0.85)
#set y grid
#plt.grid(axis='y', alpha=0.75)
# calculate mean
im_mean = im_array[~np.isnan(im_array)].mean()
# add mean line
plt.axvline(im_mean, color='red', linestyle='solid', linewidth=1)
# Set ticks, labels and title
plt.xticks(bins)
plt.xlabel("Slope (°)")
plt.ylabel('Frequency')
plt.title('Histogram\n'+name+"\nAverage Slope (°): "+str(round(im_mean,1)))
#set output name to sth computer readable and save
output_name = "histograms/"+path[:-4]+"_histogram.png"
plt.savefig(output_name,dpi=200)
print("saved: ",output_name)
# clear figure after saving
plt.clf()
plt.close("all")
for path in paths_slopes:
graph(to_np(path),path,extract_title(path))
for path in paths_slopes_resample:
graph(to_np(path),path,extract_title(path))
#Warning: Overflow in mean calculation possible
for path in path_slope_10:
graph(to_np(path),path,"Slope of full catchment area, 10x10m")
for path in path_slope_full:
graph(to_np(path),path,"Slope of full catchment area, 100x100m") 
