3.6 Digital Elevation Model#
A Digital Elevation Model (DEM) is a 3D representation of a terrain’s surface, created from terrain elevation data. It represents the earth’s surface in a raster format, where each cell or pixel holds the elevation value at that location.
Used to analyze topography, such as slope, aspect, and curvature, which are essential for understanding landforms and landscape features.
This chapter covers how to create, process, and analyze DEMs using Python and shell scripts. We will walk through the process of creating a GeoTIFF file for a specific region, reprojection and resampling of DEMs, and extracting various features from DEMs
3.6.1 Characteristics#
Product/Data Type: SRTM 90m Digital Elevation Model (DEM)
Nominal Data Array Dimensions: 5° x 5° tiles
Spatial Resolution: 90 meters (at the equator)
Temporal Resolution: Single-time snapshot (data captured during the SRTM mission in 2000)
Vertical Accuracy: Less than 16 meters error
Data Format: ArcInfo ASCII and GeoTiff
Coverage: Western USA
Projection: WGS84 datum, geographic coordinate system
3.6.2 Creating a GeoTIFF Template for the Western U.S.#
Our goal is to create a GeoTIFF file that serves as a template for the western U.S. This GeoTIFF will have a specified spatial extent and resolution, and will initially contain an empty 2D array. This template can be used as a starting point for adding real elevation data later.
import os
import rasterio
from rasterio.transform import from_origin
import numpy as np
os: module is for handling file operations.
rasterio: for reading and writing raster data in GeoTIFF format.
Raster Data: Refers to any type of spatial data organized in a grid of cells or pixels. Each cell holds a value representing information about that location, such as an image or a layer of data (e.g., elevation, temperature).
GeoTIFF: Is a specific file format that stores raster data. It is an extension of the TIFF format, with added metadata for georeferencing. This metadata includes information such as the coordinate system, projection, and spatial extent, allowing the raster data to be accurately positioned on the Earth’s surface.
numpy: is for handling numerical operations and array manipulations
minx, miny, maxx, maxy = -125, 25, -100, 49
resolution = 0.036
width = int((maxx - minx) / resolution)
height = int((maxy - miny) / resolution)
data = np.zeros((height, width), dtype=np.float32)
output_filename = f"../data/dem/western_us_geotiff_template.tif"
with rasterio.open(
output_filename,
'w',
driver='GTiff',
height=height,
width=width,
count=1, # Single band
dtype=np.float32,
crs='EPSG:4326', # WGS84
transform=from_origin(minx, maxy, resolution, resolution),
) as dst:
dst.write(data, 1)
minx and maxx: Define the western and eastern boundaries.
miny and maxy: Define the southern and northern boundaries.
resolution: The distance between adjacent pixels in degrees.
Based on the spatial extent and resolution, we calculate the width and height of the DEM in pixels:
We initialize a NumPy array to hold the DEM data. This array will be empty initially but can be populated with elevation data later
Using the Rasterio library, we create the GeoTIFF file and write the empty data array to it. We specify the metadata, including the coordinate reference system (CRS) and the transformation that maps pixel coordinates to geographic coordinates
3.6.3 Reprojecting and Resampling DEMs#
Below script does
Creates a directory for storing template shapefiles.
Copies the template GeoTIFF into the working directory.
Generates a shapefile from the template TIFF to use as a clipping mask.
Reprojects and resamples the DEM using gdalwarp to match the template’s spatial extent and resolution.
Checks the details of the output file to ensure correctness.
#!/bin/bash
# This script will reproject and resample the western US DEM, clip it, to match the exact spatial extent and resolution as the template TIFF.
cd /home/chetana/gridmet_test_run
mkdir template_shp/
cp /home/chetana/western_us_geotiff_template.tif template_shp/
# Generate the template shape
gdaltindex template.shp template_shp/*.tif
# Reproject and resample the DEM
gdalwarp -s_srs EPSG:4326 -t_srs EPSG:4326 -tr 0.036 0.036 -cutline template.shp -crop_to_cutline -overwrite output_4km.tif output_4km_clipped.tif
# Check the output
gdalinfo output_4km_clipped.tif
Generated output looks like this
Processing output_4km.tif [1/1] : 0Warning 1: the source raster dataset has a SRS, but the cutline features
not. We assume that the cutline coordinates are expressed in the destination SRS.
If not, cutline results may be incorrect.
...10...20...30...40...50...60...70...80...90...100 - done.
Driver: GTiff/GeoTIFF
Files: output_4km_clipped.tif
Size is 694, 666
Coordinate System is:
GEOGCRS["WGS 84",
ENSEMBLE["World Geodetic System 1984 ensemble",
MEMBER["World Geodetic System 1984 (Transit)"],
MEMBER["World Geodetic System 1984 (G730)"],
MEMBER["World Geodetic System 1984 (G873)"],
MEMBER["World Geodetic System 1984 (G1150)"],
MEMBER["World Geodetic System 1984 (G1674)"],
MEMBER["World Geodetic System 1984 (G1762)"],
MEMBER["World Geodetic System 1984 (G2139)"],
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ENSEMBLEACCURACY[2.0]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-125.000000000000000,49.000000000000000)
Pixel Size = (0.036000000000000,-0.036000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-125.0000000, 49.0000000) (125d 0' 0.00"W, 49d 0' 0.00"N)
Lower Left (-125.0000000, 25.0240000) (125d 0' 0.00"W, 25d 1'26.40"N)
Upper Right (-100.0160000, 49.0000000) (100d 0'57.60"W, 49d 0' 0.00"N)
Lower Right (-100.0160000, 25.0240000) (100d 0'57.60"W, 25d 1'26.40"N)
Center (-112.5080000, 37.0120000) (112d30'28.80"W, 37d 0'43.20"N)
Band 1 Block=694x2 Type=Float32, ColorInterp=Gray
3.6.4 Calculating DEM Features#
Extract various terrain features from the DEM, such as slope, aspect, curvature, northness, and eastness, and save the results.
import numpy as np
import pandas as pd
from osgeo import gdal
import warnings
import rasterio
import csv
from rasterio.transform import Affine
from scipy.ndimage import sobel, gaussian_filter
def lat_lon_to_pixel(lat, lon, geotransform):
x = int((lon - geotransform[0]) / geotransform[1])
y = int((lat - geotransform[3]) / geotransform[5])
return x, y
Converts latitude and longitude coordinates to pixel coordinates using the geotransform of the raster. The geotransform provides the mapping between geographic coordinates and pixel locations.
3.6.4.1 How to calculate slope and aspect from a given dem file#
Slope: This tells us how steep the terrain is.
Aspect: This tells us the direction the slope is facing. For example, a slope might face north, south, east, or west.
def calculate_slope_aspect(dem_file):
with rasterio.open(dem_file) as dataset:
dem_data = dataset.read(1)
# transform = dataset.transform
dx, dy = np.gradient(dem_data)
slope = np.arctan(np.sqrt(dx**2 + dy**2))
slope = 90 - np.degrees(slope)
aspect = np.degrees(np.arctan2(-dy, dx))
aspect[aspect < 0] += 360
return slope, aspect
rasterio.open(dem_file): This line opens the DEM file using a library called rasterio, which is used for reading and writing geospatial data.
dataset.read(1): This reads the elevation data from the file into a 2D array called dem_data. Each element in this array represents the elevation at a specific point.
np.gradient(dem_data): This calculates the gradient of the elevation data. Think of the gradient as the rate of change of elevation. dx represents the rate of change in the horizontal direction (left to right), and dy represents the rate of change in the vertical direction (top to bottom).
np.arctan(np.sqrt(dx**2 + dy**2)): This calculates the slope in radians. The slope is found using the arctangent of the gradient’s magnitude (a combination of dx and dy).
90 - np.degrees(slope): This converts the slope from radians to degrees and adjusts it so that a flat surface has a slope of 0 degrees and a vertical surface has a slope of 90 degrees.
np.degrees(np.arctan2(-dy, dx)): This calculates the aspect in degrees. The np.arctan2 function gives the direction of the slope in radians, which we convert to degrees.
aspect[aspect < 0] += 360: This ensures that all aspect values are between 0 and 360 degrees. Sometimes, the calculated aspect can be negative, so we add 360 to these values to make them positive.
The adjustment of aspect values to the range of 0 to 360 degrees is necessary to standardize the aspect measurements. The aspect represents the compass direction that the slope faces, and compass directions are typically measured in degrees from 0 to 360:
0 degrees represents north.
90 degrees represents east.
180 degrees represents south.
270 degrees represents west. In the context of calculating aspect from the gradient:
Range of Arctan2 Output: The np.arctan2 function returns values in the range of -π to π radians, which correspond to -180 to 180 degrees. Negative values indicate directions west of north.
Positive Degree Range: To standardize these values to a positive degree range (0 to 360), we need to adjust negative aspect values.
By adding 360 degrees to any negative aspect values, we ensure that all aspect values fall within the conventional range of 0 to 360 degrees. This adjustment makes interpretation and comparison of aspect values more intuitive, aligning with standard compass directions.
Here’s a quick example to illustrate:
An aspect of -45 degrees means the slope is facing 45 degrees west of north, which is equivalent to 315 degrees in standard compass directions.
3.6.4.2 How to calculate Curvature from a dem file#
What is Curvature?
Curvature measures how much a surface deviates from being flat. In the context of terrain, curvature can tell us about features like hills, valleys, and ridges.
def calculate_curvature(dem_file, sigma=1):
with rasterio.open(dem_file) as dataset:
dem_data = dataset.read(1)
dx = sobel(dem_data, axis=1, mode='constant')
dy = sobel(dem_data, axis=0, mode='constant')
dxx = sobel(dx, axis=1, mode='constant')
dyy = sobel(dy, axis=0, mode='constant')
curvature = dxx + dyy
curvature = gaussian_filter(curvature, sigma)
return curvature
Opens the DEM file using rasterio, Reads the elevation data from the file into a 2D array called dem_data
sobel(dem_data, axis=1, mode='constant'): Applies the Sobel operator to calculate the rate of change (gradient) in the x-direction (left to right). The axis=1 means horizontal direction.
sobel(dem_data, axis=0, mode='constant'): Applies the Sobel operator to calculate the gradient in the y-direction (top to bottom). The axis=0 means vertical direction.
The Sobel operator is a method used to find the gradient of an image, highlighting changes in intensity (in this case, elevation).
dxx = sobel(dx, axis=1, mode='constant')
dyy = sobel(dy, axis=0, mode='constant')
Applies the Sobel operator again to the gradient dx to find the second derivative in the x-direction.
Applies the Sobel operator to the gradient dy to find the second derivative in the y-direction.
These second derivatives help us understand how the rate of change itself changes, which is crucial for calculating curvature.
curvature = dxx + dyy: Adds the second derivatives in both directions to get the total curvature. This gives us an idea of how the terrain bends in both the x and y directions.
gaussian_filter(curvature, sigma): Applies a Gaussian filter to smooth the curvature values. The sigma parameter controls the amount of smoothing. A higher sigma means more smoothing.
Smoothing helps reduce noise and makes the curvature data more meaningful and easier to interpret.
3.6.4.3 How to calculate gradients#
Why Northness and Eastness?
Northness and eastness describe how much a terrain slope faces north or east. They are used to understand the orientation of the terrain:
Northness: Indicates the degree to which a slope faces north.
Eastness: Indicates the degree to which a slope faces east.
Sunlight Exposure:
Temperature and Climate Effects:
Ecological and Hydrological Impacts:
def calculate_gradients(dem_file):
with rasterio.open(dem_file) as dataset:
dem_data = dataset.read(1)
dy, dx = np.gradient(dem_data, dataset.res[0], dataset.res[1])
northness = np.arctan(dy / np.sqrt(dx**2 + dy**2))
eastness = np.arctan(dx / np.sqrt(dx**2 + dy**2))
return northness, eastness
Open the DEM file and read the elevation data into a 2D array called dem_data.
dy, dx = np.gradient(dem_data, dataset.res[0], dataset.res[1])
Calculate the rate of change (gradient) of the elevation data:
dy: Gradient in the y-direction (vertical changes).
dx: Gradient in the x-direction (horizontal changes).
northness = np.arctan(dy / np.sqrt(dx**2 + dy**2))
eastness = np.arctan(dx / np.sqrt(dx**2 + dy**2))
find northness by calculating the ratio of the vertical gradient to the overall gradient and convert this ratio into an angle using arctan.
find eastness by calculating the ratio of the horizontal gradient to the overall gradient and Convert this ratio into an angle using arctan.
3.6.5 GeoTIFF to CSV Conversion#
def geotiff_to_csv(geotiff_file, csv_file, column_name):
with rasterio.open(geotiff_file) as dataset:
data = dataset.read(1)
transform = dataset.transform
height, width = data.shape
with open(csv_file, 'w', newline='') as csvfile:
csvwriter = csv.writer(csvfile)
csvwriter.writerow(['Latitude', 'Longitude', 'x', 'y', column_name])
for y in range(height):
for x in range(width):
image_value = data[y, x]
lon, lat = transform * (x, y)
csvwriter.writerow([lat, lon, x, y, image_value])
The geotiff_to_csv function reads a GeoTIFF file, extracts its pixel values and corresponding geographic coordinates, and writes this data to a CSV file. Each row in the CSV file contains the latitude, longitude, pixel x and y coordinates, and the corresponding pixel value from the GeoTIFF file. This function effectively converts raster data from a GeoTIFF into a tabular format, useful for further analysis or visualization.
3.6.6 How to save GeoTIFF with Meta Data#
def save_as_geotiff(data, output_file, src_file):
"""
Save data as a GeoTIFF file with metadata from the source file.
Args:
data (array): Data to be saved.
output_file (str): Path to the output GeoTIFF file.
src_file (str): Path to the source GeoTIFF file to inherit metadata from.
"""
with rasterio.open(src_file) as src_dataset:
profile = src_dataset.profile
transform = src_dataset.transform
# Update the data type, count, and set the transform for the new dataset
profile.update(dtype=rasterio.float32, count=1, transform=transform)
# Create the new GeoTIFF file
with rasterio.open(output_file, 'w', **profile) as dst_dataset:
# Write the data to the new GeoTIFF
dst_dataset.write(data, 1)
Opens the source GeoTIFF file using rasterio.
src_dataset.profile: Retrieves the metadata profile from the source file. This profile contains information such as the data type, coordinate reference system, width, height, etc.
src_dataset.transform: Gets the affine transformation matrix that maps pixel coordinates to geographic coordinates. This ensures that the new file will have the same geographic reference as the source file.
profile.update(dtype=rasterio.float32, count=1, transform=transform)
Updates the metadata profile for the new GeoTIFF file.
dtype=rasterio.float32: Sets the data type of the new file to 32-bit floating point. This is suitable for storing continuous data.
count=1: Sets the number of layers (bands) in the new file to 1.
transform=transform: Ensures that the new file uses the same geotransform as the source file, preserving its geographic alignment.
with rasterio.open(output_file, 'w', **profile) as dst_dataset:
dst_dataset.write(data, 1)
Creates the new GeoTIFF file with the updated profile.
Writes the provided data to the first layer (band) of the new GeoTIFF file.
3.6.7 Unleashing Terrain Insights: From DEM to CSV#
Now lets utilise all the functions we have created to convert the dem files to csv files and merge them into a single csv file consisting slope, aspect, curvature, northness, and eastness
result_dem_csv_path = "../data/dem/dem_template.csv"
result_dem_feature_csv_path = "../data/dem/dem_all.csv"
dem_file = "../data/dem/dem_file.tif"
slope_file = '../data/dem/dem_file.tif_slope.tif'
aspect_file = '../data/dem/dem_file.tif_aspect.tif'
curvature_file = '../data/dem/curvature_file.tif'
northness_file = '../data/dem/northness_file.tif'
eastness_file = '../data/dem/eastness_file.tif'
slope, aspect = calculate_slope_aspect(dem_file)
# slope = calculate_slope(dem_file)
# aspect = calculate_aspect(dem_file)
curvature = calculate_curvature(dem_file)
northness, eastness = calculate_gradients(dem_file)
# Save the slope and aspect as new GeoTIFF files
save_as_geotiff(slope, slope_file, dem_file)
save_as_geotiff(aspect, aspect_file, dem_file)
save_as_geotiff(curvature, curvature_file, dem_file)
save_as_geotiff(northness, northness_file, dem_file)
save_as_geotiff(eastness, eastness_file, dem_file)
geotiff_to_csv(dem_file, dem_file+".csv", "Elevation")
geotiff_to_csv(slope_file, slope_file+".csv", "Slope")
geotiff_to_csv(aspect_file, aspect_file+".csv", "Aspect")
geotiff_to_csv(curvature_file, curvature_file+".csv", "Curvature")
geotiff_to_csv(northness_file, northness_file+".csv", "Northness")
geotiff_to_csv(eastness_file, eastness_file+".csv", "Eastness")
# List of file paths for the CSV files
csv_files = [dem_file+".csv", slope_file+".csv", aspect_file+".csv",
curvature_file+".csv", northness_file+".csv", eastness_file+".csv"]
# Initialize an empty list to store all dataframes
dfs = []
# Read each CSV file into separate dataframes
for file in csv_files:
df = pd.read_csv(file, encoding='utf-8')
dfs.append(df)
# Merge the dataframes based on the latitude and longitude columns
merged_df = dfs[0] # Start with the first dataframe
for i in range(1, len(dfs)):
merged_df = pd.merge(merged_df, dfs[i], on=['Latitude', 'Longitude', 'x', 'y'])
# check the statistics of the columns
for column in merged_df.columns:
merged_df[column] = pd.to_numeric(merged_df[column], errors='coerce')
print(merged_df[column].describe())
# Save the merged dataframe to a new CSV file
merged_df.to_csv(result_dem_feature_csv_path, index=False)
print(f"New dem features are updated in {result_dem_feature_csv_path}")
/var/folders/tg/bv4s19f94fv3gqmxhbtgwf_80000gn/T/ipykernel_39509/2570874998.py:5: RuntimeWarning: invalid value encountered in divide
northness = np.arctan(dy / np.sqrt(dx**2 + dy**2))
/var/folders/tg/bv4s19f94fv3gqmxhbtgwf_80000gn/T/ipykernel_39509/2570874998.py:6: RuntimeWarning: invalid value encountered in divide
eastness = np.arctan(dx / np.sqrt(dx**2 + dy**2))
count 462204.000000
mean 37.030000
std 6.921275
min 25.060000
25% 31.036000
50% 37.030000
75% 43.024000
max 49.000000
Name: Latitude, dtype: float64
count 462204.00000
mean -112.52600
std 7.21226
min -125.00000
25% -118.77200
50% -112.52600
75% -106.28000
max -100.05200
Name: Longitude, dtype: float64
count 462204.000000
mean 346.500000
std 200.340552
min 0.000000
25% 173.000000
50% 346.500000
75% 520.000000
max 693.000000
Name: x, dtype: float64
count 462204.000000
mean 332.500000
std 192.257631
min 0.000000
25% 166.000000
50% 332.500000
75% 499.000000
max 665.000000
Name: y, dtype: float64
count 462204.000000
mean 1025.629083
std 808.945819
min -82.938410
25% 180.727250
50% 1018.458080
75% 1604.870600
max 4177.173000
Name: Elevation, dtype: float64
count 462204.000000
mean 20.325398
std 35.542848
min 0.044037
25% 0.521088
50% 1.558449
75% 7.970741
max 90.000000
Name: Slope, dtype: float64
count 462204.000000
mean 136.166265
std 114.602954
min -0.000000
25% 19.602830
50% 121.656285
75% 232.066403
max 359.999360
Name: Aspect, dtype: float64
count 462204.000000
mean -49.124993
std 2819.341836
min -41085.560000
25% -482.398620
50% 0.000000
75% 840.230590
max 21411.324000
Name: Curvature, dtype: float64
count 368579.000000
mean -0.012698
std 0.587101
min -0.785398
25% -0.626382
50% -0.027939
75% 0.608819
max 0.785398
Name: Northness, dtype: float64
count 368579.000000
mean -0.040279
std 0.583256
min -0.785398
25% -0.643895
50% -0.087261
75% 0.573790
max 0.785398
Name: Eastness, dtype: float64
New dem features are updated in ../data/dem/dem_all.csv