Digital Elevation Models (DEMs) store raw orthometric height values for every pixel, but the real analytical power of terrain modeling comes from computing secondary and tertiary surface derivatives. By analyzing the rate of change of elevation across a 3x3 moving kernel window, GIS algorithms derive critical geomorphological metrics: Slope (steepness), Aspect (compass direction of terrain face), Hillshade (simulated solar illumination for 3D cartographic relief), and Topographic Roughness. These terrain layers form the backbone of landslide hazard zoning, solar farm siting, hydrological runoff modeling, and forestry ecology.
📋 Prerequisites
- QGIS 3.34+ LTR installed.
- An elevation raster (e.g., 30m SRTM DEM or 12.5m ALOS PALSAR DEM).
- CRITICAL: DEM must be reprojected into a metric Projected Coordinate System (UTM).
🛠️ Technical Environment
Required Software: QGIS Raster Terrain Analysis / GDAL (Recommended: 3.34+ LTR)
Practice Dataset: ALOS PALSAR / SRTM 1 Arc-Second Global DEM
Source Portal: NASA Earthdata
CRS / Format: Projected UTM (Meters) (GeoTIFF)
Step-by-Step Workflow & Methodological Execution
Module 1: The Z-Factor and Reprojection Prerequisite
The single most common error in terrain analysis occurs when running slope or hillshade on a DEM in Geographic Coordinates (WGS 84 - EPSG:4326): • The Horizontal vs Vertical Mismatch: In EPSG:4326, the horizontal coordinates ($X, Y$) are measured in angular degrees, while elevation ($Z$) is measured in linear meters. If you calculate slope without reprojection, the algorithm divides meters by degrees, producing completely distorted, blacked-out slope maps! • The Solution: Before running terrain analysis, always reproject your DEM to a local metric Projected Coordinate System (e.g., UTM Zone 43N - EPSG:32643) via Raster -> Projections -> Warp (Reproject). If you must work in degrees, calculate the appropriate Z-Factor multiplier ($Z \approx 1 / 111,320$ at the equator).
Module 2: Calculating Slope and Aspect
1. Slope (Raster -> Analysis -> Slope): Measures the maximum rate of elevation change across adjacent cells using Horn's or Zevenbergen & Thorne's 3x3 moving kernel algorithm: • Choose units: Degrees (0° horizontal flat plain to 90° vertical cliff) or Percent Slope (rise over run $\times 100$). Percent slope is standard in civil engineering and road alignment design. 2. Aspect (Raster -> Analysis -> Aspect): Computes the compass azimuth direction that the downhill slope face points: • Expressed in degrees clockwise from 0° to 360° (0° North, 90° East, 180° South, 270° West). Flat cells are assigned -1 or 9999. • Environmental Applications: In the Northern Hemisphere, south-facing slopes receive significantly higher solar insolation, driving faster snowmelt, drier soils, and distinct vegetation communities compared to shaded north-facing slopes.
Module 3: Generating Multi-Directional Hillshades and TRI
1. Hillshade (Raster -> Analysis -> Hillshade): Simulates realistic 3D shadows by calculating hypothetical solar illumination: • Azimuth: The light source direction (standard cartographic convention is 315° - Northwest. Lighting from the south produces an optical illusion called 'relief inversion', where valleys appear as ridges!). • Altitude: The sun angle above the horizon (standard is 45°). • Multidirectional Hillshade: Merges illumination from 4 distinct directions (225°, 270°, 315°, 360°) to reveal subtle geological fault lines obscured by single-direction shadows. 2. Terrain Ruggedness Index (TRI): Quantifies topographic roughness by calculating the mean elevation difference between a focal cell and its 8 neighbors, vital for wildlife habitat modeling.
⚠️ Common Errors & Troubleshooting
❌ Slope values calculated as 89.9° across flat terrain
💡 Resolution: Z-factor error! If the DEM coordinates are in degrees (lat/long) and elevation is in meters, you must set Z-factor to ~0.000009 or reproject the DEM to UTM meters.
❌ Hillshade appears inverted (valleys look like ridges)
💡 Resolution: Default solar azimuth should be from North-West (315°); southern illumination causes the optical brain illusion of relief inversion.
💡 Expert Tips & Best Practices
- Combine a semi-transparent colored elevation ramp (40% opacity) on top of a multidirectional hillshade for stunning cartographic 3D terrain visualization.
- Use Horn's formula for smooth terrain and Zevenbergen-Thorne for steep, rugged alpine geomorphometry.
🐍 Python RichDEM / NumPy Slope & Hillshade Pipeline
import rasterio
import numpy as np
from scipy.ndimage import convolve
# Open metric UTM Elevation DEM
with rasterio.open("elevation_utm.tif") as src:
dem = src.read(1).astype('float32')
res = src.transform[0] # Pixel width in meters
profile = src.profile
# Compute gradient in X and Y using Sobel 3x3 filter kernels
kernel_x = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]]) / (8.0 * res)
kernel_y = np.array([[1, 2, 1], [0, 0, 0], [-1, -2, -1]]) / (8.0 * res)
dz_dx = convolve(dem, kernel_x)
dz_dy = convolve(dem, kernel_y)
# Calculate Slope in Radians and Degrees
slope_rad = np.arctan(np.sqrt(dz_dx**2 + dz_dy**2))
slope_deg = np.degrees(slope_rad)
# Calculate Hillshade (Sun Azimuth: 315°, Altitude: 45°)
azimuth_rad = np.radians(315.0)
altitude_rad = np.radians(45.0)
aspect_rad = np.arctan2(dz_dy, -dz_dx)
hillshade = 255.0 * (
(np.sin(altitude_rad) * np.cos(slope_rad)) +
(np.cos(altitude_rad) * np.sin(slope_rad) * np.cos(azimuth_rad - aspect_rad))
)
hillshade = np.clip(hillshade, 0, 255).astype(np.uint8)
# Export Hillshade raster
profile.update(dtype=rasterio.uint8, count=1)
with rasterio.open("calculated_hillshade.tif", "w", **profile) as dst:
dst.write(hillshade, 1)
print("Slope and Hillshade rasters generated successfully.")