Volumetric analysis is a powerful technique used to estimate the volume of natural and man-made features from elevation data. Whether measuring stockpiles at a mine, calculating the capacity of an excavation, or monitoring changes in terrain over time, accurate volume calculations provide valuable insights for planning, resource management, and quality control. This tutorial walks through the steps involved in the workflow using the CATALYST Professional Python API, with a Digital Surface Model (DSM), Digital Terrain Model (DTM), and polygon boundaries that define the area of interest (AOI).
The process involves generating a height difference raster from the DSM and DTM, identifying the elevation differences within each stockpile boundary, and calculating the total volume based on the height difference and pixel area.
To accompany this tutorial, a Jupyter Notebook that implements the workflow using the CATALYST Python API alongside open-source Python libraries is provided. This notebook serves as a practical example that can be followed, modified, and adapted to different datasets.
CATALYST Volumetric Analysis Jupyter Notebook - Stockpile_Volumetric_Analysis.ipynb
Before performing volumetric analysis, both DSM and DTM elevation data are required. The process of creating these datasets is covered in detail in various tutorials, which are linked below. If these files are already available, you can proceed to the next section to calculate the height difference.
Once the DSM and DTM rasters have been generated, the next step is to create a raster containing the elevation difference between the surface and underlying terrain. This is done by subtracting the DTM from the DSM on a pixel-by-pixel basis. The resulting image values correspond to the height of the objects above the ground, which will be the primary input for the volume calculations.
The Change Detection tool is well-suited for this step, as it computes the difference between a working raster and reference raster. In this scenario, the DSM is used as the working raster and the DTM as the reference raster. When automating the workflow in Python, the CHDETOP algorithm provides the equivalent functionality, allowing the difference raster to be generated programmatically.
To generate the height difference raster in Focus:
The stockpile segments are extracted through image segmentation, which groups neighbouring pixels with similar characteristics into distinct regions to represent individual objects. This way, each stockpile is isolated from the surrounding terrain and neighbouring stockpiles, allowing for easy volume calculation. The EXPOLRAS algorithm is used here to extract and refine the vectors.
To extract the stockpile polygons:
With the stockpile boundaries extracted and the height difference raster generated, the final step is to calculate the volume of each stockpile. This is accomplished by combining the height information from the raster with the area represented by the image resolution. For each pixel contained within a segment, the product of the height difference and the pixel area is calculated. The sum of these values across all pixels within a segment represents the estimated stockpile volume.
This can be done in Python with the help of the numPy, rasterio, and geopandas open-source libraries along with the CATALYST Python API.
import numpy as np
import rasterio
from rasterio.mask import mask
import geopandas as gpd
import pci
from pci.api import datasource as ds
# Calculate Volume of stockpile segments
stockpile_gpd = gpd.read_file(stockpile_poly)
height_diff = rasterio.open(height_diff)
volumes = []
with ds.open_dataset(dsm, ds.eAM_READ) as dataset:
resolution = dataset.resolution_in_meters
pixel_area = resolution[0] * resolution[1]
for idx, row in stockpile_gpd.iterrows():
# Clip difference raster to stockpile polygon
clipped, _ = mask(height_diff, [row.geometry], crop=True, filled=False)
# Remove band dimension
delta_height = clipped[0]
# Sum height values
height_sum = np.sum(delta_height)
# Calculate volume
volume = height_sum * pixel_area
volumes.append(volume)
# Save to file
stockpile_gpd["volume_m3"] = volumes
stockpile_gpd.to_file(os.path.join(stockpile_vol_dir, "stockpile_volume.shp"))
height_diff.close()In this sample script, stockpile_poly represents the stockpile segments extracted using the EXPOLRAS algorithm, height_diff represents the height difference raster generated with the change detection tool.
The image resolution is obtained from the CATALYST Python API dataset object's resolution_in_meters tuple property. The x and y resolution values are multiplied to determine the area represented by a single pixel. Next, the height_sum parameter accumulates the height difference values for all pixels within a stockpile segment. Multiplying this accumulated height by pixel_area produces the estimated volume for that stockpile.
Once the volume has been calculated for each segment, the resulting values are added as a new attribute to the stockpile geopandas object. The updated dataset is then exported as a shapefile, allowing the calculated volumes to be viewed, analyzed, or used in other GIS applications.
The estimated volume of the stockpile, as shown in Figure 1, is 416,799.78 cubic metres (m3).