How to Compare DSMs in Python: Analyze Surface Elevation Changes and Generate Heat Maps

Introduction

Digital Surface Models (DSMs) represent the Earth’s surface, including all features such as buildings, trees, and other objects above the ground. They are widely used in geospatial analysis for applications like urban planning, flood modeling, and vegetation monitoring. By comparing two DSMs from different time periods, we can detect changes in surface elevation caused by construction, deforestation, or natural processes.

The data and source code for this project can be accessed via the following link:
https://github.com/geospatialearning/compare_dsm_with_heatmap .
The datasets used in this analysis were downloaded from the UK Environment Agency’s open survey data portal:
https://environment.data.gov.uk/survey .

In this tutorial, you’ll learn how to:

Load and process DSMs using Python.
Compute elevation differences between two DSMs.
Visualize changes with a heat map.
Calculate volumetric changes (fill and cut).
Save the results as a new GeoTIFF file and an image.

Why Compare DSMs?

Comparing DSMs allows us to:

Monitor Urban Development: Track the addition or removal of structures, such as buildings or infrastructure.
Assess Environmental Changes: Study deforestation, reforestation, or changes in land use.
Analyze Natural Processes: Detect erosion, sediment deposition, or subsidence.
Support Decision-Making: Provide actionable insights for urban planners, environmental scientists, and engineers.

Step-by-Step Workflow

Step 1: Import Required Libraries

We’ll use Python libraries like rasterio, numpy, and matplotlib to process and visualize the DSMs:

import rasterio
import numpy as np
import matplotlib.pyplot as plt

rasterio: Reads and processes raster data.
numpy: Performs numerical operations on arrays.
matplotlib: Creates visualizations like heat maps.

Step 2: Define Paths to Input DSMs

Specify the file paths for the two DSMs:

dsm_2018_path = "national_lidar_programme_dsm-2018-1-TQ07nw/DSM_TQ0075_P_10752_20180118_20180118.tif"
dsm_2023_path = "national_lidar_programme_dsm-2023-1-TQ07nw/DSM_TQ0075_P_12757_20230109_20230315.tif"

Replace these paths with the actual locations of your DSM files.

Step 3: Load and Read the DSMs

Use rasterio to open the DSMs and read their elevation data as NumPy arrays:

with rasterio.open(dsm_2018_path) as src_2018, rasterio.open(dsm_2023_path) as src_2023:
    dsm_2018 = src_2018.read(1).astype(np.float32)  # Read as float32 to avoid overflow
    dsm_2023 = src_2023.read(1).astype(np.float32)
    profile = src_2018.profile  # Metadata for saving output

src_2018.read(1): Reads the first band of the raster as a NumPy array.
.astype(np.float32): Ensures the data type is suitable for numerical calculations.
profile: Stores metadata (e.g., resolution, CRS, transform) for later use.

Step 4: Calculate the Elevation Difference

Compute the difference between the 2023 and 2018 DSMs:

difference = dsm_2023 - dsm_2018

Positive values indicate areas where the surface has risen (e.g., new buildings or vegetation growth).
Negative values indicate areas where the surface has lowered (e.g., demolition or deforestation).

Step 5: Handle NoData and Invalid Values

Mask out invalid values (NoData, Inf) to prevent errors during calculations:

nodata_value = profile['nodata']
if nodata_value is not None:
    difference[(dsm_2018 == nodata_value) | (dsm_2023 == nodata_value)] = np.nan
difference[np.isinf(difference)] = np.nan

np.nan: Represents missing or invalid data.
np.isinf(): Detects infinite values.

Step 6: Visualize the Heat Map

Create a heat map to visualize surface elevation changes:

plt.figure(figsize=(10, 8))
plt.imshow(difference, cmap='coolwarm', vmin=-2, vmax=2)
plt.colorbar(label="Elevation Change (m)")
plt.title("Heat Map of Surface Elevation Changes (2018 vs 2023)")
plt.savefig("Surface_Elevation_Change_Heatmap.png", dpi=300, bbox_inches='tight')  # Save the figure
plt.show()

cmap='coolwarm': Uses a diverging color scheme (blue for cuts, red for fills).
vmin and vmax: Set the range of values displayed in the heat map.
plt.savefig(): Saves the heat map as an image file.

Step 7: Calculate Volumetric Changes

Compute the total volume of surface rise (fill) and surface lowering (cut):

transform = profile['transform']
pixel_area = transform[0] * -transform[4]  # Pixel width × height

volume_raise = np.nansum(np.where(difference > 0, difference, 0)) * pixel_area
volume_cut = np.nansum(np.where(difference < 0, difference, 0)) * pixel_area

print(f"Total Volume of Surface Rise (Fill): {volume_raise:.2f} m³")
print(f"Total Volume of Surface Lowering (Cut): {volume_cut:.2f} m³")

np.nansum(): Ignores NaN values during summation.
pixel_area: Multiplies the elevation difference by the area of each pixel to calculate volume.

Step 8: Save the Difference Raster

Save the computed difference as a new GeoTIFF file:

profile.update(dtype=rasterio.float32, nodata=np.nan)
output_path = "difference_2018_2023.tif"
with rasterio.open(output_path, 'w', **profile) as dst:
    dst.write(difference.astype(rasterio.float32), 1)

profile.update(): Updates metadata for the new raster.
dst.write(): Writes the difference array to the file.

Expected Outputs

Heat Map:

A visual representation of surface elevation changes saved as Surface_Elevation_Change_Heatmap.png.

Volumetric Analysis:

Printed volumes of surface rise (fill) and surface lowering (cut) in cubic meters.

Difference Raster:

A new GeoTIFF file (difference_2018_2023.tif) containing the elevation differences.

Applications of DSM Comparison

Urban Planning: Identify new constructions or demolished structures.
Environmental Monitoring: Track changes in vegetation or land use.
Disaster Management: Assess damage caused by natural disasters like landslides or floods.
Infrastructure Development: Monitor progress at construction sites.

Troubleshooting Tips

Mismatched Resolutions or Extents:

Ensure both DSMs have the same resolution and extent. Use tools like gdalwarp to align them if needed.

Invalid Values:

Check for NaN or Inf values in the input DSMs and handle them appropriately.

Overflow Errors:

Use float32 or float64 data types to avoid numerical overflow.

Conclusion

This tutorial demonstrates how to compare two DSMs to analyze surface elevation changes over time. By following these steps, you can:

Visualize changes using a heat map.
Quantify volumetric changes (surface rise and lowering).
Save the results for further analysis or reporting.

Feel free to adapt this workflow to other datasets or applications, such as monitoring urban development, analyzing vegetation changes, or assessing flood risks. Let me know if you need further assistance!

I hope this tutorial will create a good foundation for you. If you want tutorials on another GIS topic or you have any queries, please send an mail at contact@spatial-dev.guru.