Calculate Cut and Fill Volume in PyQGIS

Earthworks are measured in cubic metres: how much soil must be excavated to level a building platform, how much fill a road embankment needs, how large a gravel stockpile is today compared with last month. With elevation rasters the calculation is a sum over cells — the height difference in each cell times the cell's area — and the details that matter are the reference surface, the boundary of the site, NoData handling, and the precision the inputs really support.

This recipe belongs to Terrain & Interpolation Analysis. It computes volume against a flat design level with Processing, cut and fill between two surfaces with NumPy, stockpile volumes from repeated surveys, restricts every calculation to a site boundary, and reports results with appropriate precision.

Cut and fill against a design levelA terrain profile crossed by a horizontal design level. Where existing ground is above the design level, material must be cut, shown in amber. Where it is below, material must be filled, shown in blue. Volume is the sum over cells of the height difference times the cell area, kept separately for cut and fill.Above the level: cut; below: fillcutfilldesign level — volume = Σ (Δh × cell area)

Prerequisites

  • QGIS 3.34 LTR or newer, or the QGIS 4 series.
  • Elevation rasters in a projected CRS with metre units for both horizontal and vertical values. Surveys from drones or LiDAR are typical sources; creating a DEM from a point cloud shows how to produce one.
  • A site boundary polygon.

Volume against a flat level

For a platform at a single design elevation, native:rastersurfacevolume computes the volume above or below that level directly.

import processing
from qgis.core import QgsRasterLayer, QgsVectorLayer

site = QgsVectorLayer("/data/site/platform_boundary.gpkg", "site", "ogr")
dem = processing.run("gdal:cliprasterbymasklayer", {
    "INPUT": "/data/site/survey_2026_09.tif", "MASK": site, "CROP_TO_CUTLINE": True,
    "KEEP_RESOLUTION": True, "NODATA": -9999, "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]

LEVEL = 112.40                       # design platform elevation, metres
for method, label in ((0, "above level (cut)"), (1, "below level (fill)")):
    r = processing.run("native:rastersurfacevolume", {
        "INPUT": dem, "BAND": 1, "LEVEL": LEVEL, "METHOD": method,
        "OUTPUT_HTML_FILE": "TEMPORARY_OUTPUT"})
    print(f"{label:<20} {r['VOLUME']:12,.1f} m³ over {r['AREA']:,.0f} m²")

Breakdown: Clipping to the site boundary first makes sure only cells inside the platform count — volume outside the site is irrelevant and, at the edges of a survey, often unreliable. METHOD 0 counts only cells above the level, giving the cut; 1 counts cells below, giving the fill; other methods return the net or absolute total. The algorithm reports the volume, the area of the cells that contributed and a pixel count. For a balanced earthwork, adjust the level until cut and fill are roughly equal, which a short loop over candidate levels does quickly.

Cut and fill between two surfaces

Most designs are not flat: road corridors, sloped platforms, landscaped mounds. With a design surface as a raster on the same grid, cut and fill come from the cell-by-cell difference.

Difference of two surfacesThe existing surface minus the design surface gives a difference raster. Positive cells need cutting, negative cells need filling. Multiplying by cell area and summing positive and negative parts separately gives cut and fill volumes. Both rasters must share the same grid, and NoData in either input must be excluded from both sums.existing − designexisting DEMsurveydesign surfacesame gridΔh > 0 → cutΣ Δh × areaΔh < 0 → fillΣ |Δh| × area

import numpy as np
from osgeo import gdal

def read(path):
    ds = gdal.Open(path)
    band = ds.GetRasterBand(1)
    arr = band.ReadAsArray().astype("float64")
    nd = band.GetNoDataValue()
    if nd is not None:
        arr[arr == nd] = np.nan
    return arr, ds.GetGeoTransform()

existing, gt1 = read("/data/site/existing_clip.tif")
design, gt2 = read("/data/site/design_surface.tif")
assert gt1 == gt2 and existing.shape == design.shape, "surfaces are not on the same grid"

cell_area = abs(gt1[1] * gt1[5])
diff = existing - design
valid = ~np.isnan(diff)
cut = np.nansum(np.where(diff > 0, diff, 0)) * cell_area
fill = np.nansum(np.where(diff < 0, -diff, 0)) * cell_area
print(f"cut {cut:,.0f} m³ · fill {fill:,.0f} m³ · net {cut - fill:+,.0f} m³ "
      f"over {valid.sum() * cell_area:,.0f} m²")

Breakdown: Converting NoData to NaN in both arrays and summing with nansum excludes any cell missing in either surface, so holes in the survey do not count as zero height. The assertion on geotransform and shape is essential: subtracting rasters on different grids gives nonsense without any error. Align the design surface to the survey grid first, as in resampling and aligning rasters. Cell area from the geotransform makes the calculation work at any resolution. Writing diff back to a GeoTIFF gives a cut/fill map — style it with a diverging ramp centred on zero.

Find the balanced platform level

Moving soil off site is expensive, so platforms are often set at the level where cut equals fill. With the clipped survey as an array, a quick search finds it.

site_dem, _ = read(dem)

def net_volume(level):
    d = site_dem - level
    return np.nansum(d) * cell_area           # positive: surplus cut

lo, hi = np.nanmin(site_dem), np.nanmax(site_dem)
for _ in range(40):                           # bisection
    mid = (lo + hi) / 2
    if net_volume(mid) > 0:
        lo = mid
    else:
        hi = mid
print(f"balanced level ≈ {mid:.2f} m (net {net_volume(mid):+,.0f} m³)")

Breakdown: The net volume — cut minus fill — falls steadily as the level rises, so bisection between the lowest and highest ground finds the level where it crosses zero in a few dozen steps. In practice a design level is then rounded and adjusted for drainage falls, compaction and the volume of foundations, but the balanced level is the natural starting point and shows immediately whether a proposed level implies importing or exporting material.

Measure a stockpile

A stockpile's volume is the material above the ground it sits on. The ground under the pile is not surveyed, so a base surface is constructed from the pile's edge — typically a plane or a TIN through the boundary elevations.

pile = QgsVectorLayer("/data/site/stockpile_A.gpkg", "pile", "ogr")
edge_pts = processing.run("native:pointsalonglines", {
    "INPUT": processing.run("native:polygonstolines", {"INPUT": pile, "OUTPUT": "memory:"})["OUTPUT"],
    "DISTANCE": 2, "OUTPUT": "memory:"})["OUTPUT"]
edge_z = processing.run("native:rastersampling", {
    "INPUT": edge_pts, "RASTERCOPY": "/data/site/survey_2026_09.tif",
    "COLUMN_PREFIX": "z_", "OUTPUT": "memory:"})["OUTPUT"]

base_level = float(np.median([f["z_1"] for f in edge_z.getFeatures() if f["z_1"] is not None]))
clip = processing.run("gdal:cliprasterbymasklayer", {
    "INPUT": "/data/site/survey_2026_09.tif", "MASK": pile, "CROP_TO_CUTLINE": True,
    "KEEP_RESOLUTION": True, "OUTPUT": "TEMPORARY_OUTPUT"})["OUTPUT"]
vol = processing.run("native:rastersurfacevolume", {
    "INPUT": clip, "BAND": 1, "LEVEL": base_level, "METHOD": 0,
    "OUTPUT_HTML_FILE": "TEMPORARY_OUTPUT"})["VOLUME"]
print(f"stockpile A: base {base_level:.2f} m, volume {vol:,.0f} m³")

Breakdown: Sampling the surveyed surface every two metres around the pile's outline gives the elevations where the pile meets the ground; the median is a robust flat base level that ignores a few odd edge samples. For piles on sloping ground, a flat base under- or over-counts — build a TIN from the edge points with TIN interpolation and use the two-surface method instead. Repeating the measurement on monthly surveys with the same outline and method gives a consistent inventory trend.

Compare surveys over time

Monitoring sites — quarries, landfills, construction — compare consecutive surveys to measure what moved. The two-surface method applied to two dates gives removed and added volumes.

before, gt_a = read("/data/site/survey_2026_08.tif")
after, gt_b = read("/data/site/survey_2026_09.tif")
assert gt_a == gt_b
change = after - before
noise = 0.05                                    # metres: survey vertical accuracy
removed = np.nansum(np.where(change < -noise, -change, 0)) * cell_area
added = np.nansum(np.where(change > noise, change, 0)) * cell_area
print(f"removed {removed:,.0f} m³ · added {added:,.0f} m³ (changes under {noise} m ignored)")

Breakdown: A noise threshold equal to the surveys' vertical accuracy stops random measurement error from accumulating into large fictitious volumes: over a 10-hectare site, an unbiased 3 cm error in every cell adds up to thousands of cubic metres if summed blindly in both directions. Treat the threshold as part of the method and report it with the result. Systematic offsets between surveys — a datum change, a drone flight with poor ground control — are worse; check stable areas such as roads, where the change should be zero.

Report with honest precision

Volume precision is limited by the vertical accuracy of the surfaces, not by the number of digits the computer prints.

Uncertainty from vertical accuracyA rough uncertainty for a volume is the vertical accuracy of the surface times the area. A drone survey accurate to 5 centimetres over 2 hectares gives about plus or minus 1000 cubic metres in the worst systematic case. Reporting 12,846.37 cubic metres implies false precision; 12,800 plus or minus 1,000 cubic metres is honest.Accuracy × area bounds the answervertical accuracy± 0.05 m×area2 ha=systematic bound± 1,000 m³report 12,800 ± 1,000

area_m2 = valid.sum() * cell_area
vertical_accuracy = 0.05
bound = vertical_accuracy * area_m2
print(f"cut {round(cut, -2):,.0f} m³ ± {round(bound, -2):,.0f} m³ (systematic bound)")

Breakdown: Multiplying the vertical accuracy by the area gives the volume error if the whole surface were offset by that amount — a conservative bound for systematic error. Random errors partly cancel and contribute less, but systematic ones do not. Rounding the reported volume to match — hundreds of cubic metres for a 2-hectare site surveyed to 5 cm — prevents a spreadsheet from implying precision nobody has. Always state the surfaces, dates, boundary and method alongside the number.

QGIS version compatibility

native:rastersurfacevolume has been available since QGIS 3.4 and works unchanged on 3.34 LTR, 3.40 LTR and QGIS 4, with outputs VOLUME, AREA and PIXEL_COUNT. The NumPy approach is version-independent. gdal:cliprasterbymasklayer parameters are stable across these releases.

Troubleshooting

  • Volumes are far too large. NoData was counted as a real value; check the NoData value and NaN handling.
  • Cut and fill do not match a manual check. The surfaces are not on the same grid; align them first.
  • Stockpile volume changes with the outline. The base level depends on the edge; keep the same outline across surveys.
  • Monthly changes look noisy everywhere. Apply a noise threshold and check for systematic offsets on stable ground.

Conclusion

Clip every surface to the site boundary, use native:rastersurfacevolume for flat design levels and the cell-by-cell difference on a shared grid for design surfaces, build a base from the pile edge for stockpiles, threshold changes by survey accuracy when comparing dates, and report volumes rounded to the precision the surfaces support.

Frequently Asked Questions

Can I compute volumes from contour lines? Interpolate them to a raster first, for example with TIN interpolation, then use the same methods.

Does cell size affect the result? Slightly; finer cells follow slopes more faithfully. The surface's accuracy matters more than its resolution.

How do I account for soil bulking? Multiply cut volumes by a bulking factor for the material; that is an engineering input, not a GIS one.

Can I get volumes per zone of a site? Run the calculation per polygon, or use zonal statistics on the difference raster with a sum statistic multiplied by cell area.