Debugging Terracing Artefacts in Derived Terrain Rasters

Stepped, contour-like bands across gentle slopes in a hillshade almost always come from vertical quantisation in the source DEM, not from the shading parameters — so the fix is to smooth the elevation surface before deriving, and specifically not to adjust the z-factor.

Core Algorithm and Workflow

A DEM stored as 16-bit integer metres can only represent whole-metre elevations. On a slope rising by less than one metre per pixel, consecutive cells therefore hold the same value until the terrain crosses the next whole metre, at which point the value jumps. The elevation surface is a staircase, and the gradient of a staircase is zero almost everywhere with a spike at each step. Shading that gradient produces exactly what the artefact looks like: flat bands separated by hard edges that follow elevation contours.

The diagnosis follows from that mechanism:

  1. Check whether the bands follow contours. Sample elevations along a transect crossing the banding. Quantisation bands change at perfectly regular, round elevation values; real terraces do not.
  2. Read the data type and the value histogram. An integer dtype is a strong hint. A histogram of unique values showing a constant step confirms it, and gives the size of the quantum.
  3. Establish the slope at which it matters. Terracing is visible where the terrain rises by less than about one quantum per pixel. A one-metre quantum on a 30-metre grid becomes visible below roughly a two per cent slope.
Why an integer DEM produces bands on gentle ground An elevation transect across a slope rising less than one metre per pixel. The true surface is a smooth ramp with a constant gradient. The quantised surface is a staircase with flat treads and vertical risers. Below each, the derived gradient is shown: constant for the true surface, and zero with periodic spikes for the quantised one, which is what the shading renders as bands. true surface gradient constant — smooth shading quantised to whole metres gradient spikes at each riser — visible bands Both surfaces describe the same slope to within half a metre. Only one has a usable derivative.
The shading algorithm is behaving correctly in both cases. What differs is that the right-hand surface has a gradient that exists only at the risers.

Production-Ready Python Implementation

import numpy as np
import rasterio
from scipy.ndimage import gaussian_filter


def detect_quantisation(dem: np.ndarray, nodata=None) -> dict:
    """Report the vertical quantum of a DEM, if it has one."""
    values = dem[dem != nodata] if nodata is not None else dem.ravel()
    unique = np.unique(values)
    if unique.size < 3:
        return {"quantised": False, "reason": "too few distinct values to judge"}

    steps = np.diff(unique)
    step = float(np.median(steps))
    # A quantised surface has an almost constant gap between adjacent values.
    regular = bool(np.allclose(steps, step, rtol=0.02))
    return {"quantised": regular, "quantum": step if regular else None,
            "distinct_values": int(unique.size)}


def deterrace(dem_path: str, out_path: str, sigma: float = 1.0) -> dict:
    """Smooth an integer DEM just enough to remove quantisation banding."""
    with rasterio.open(dem_path) as src:
        band = src.read(1)
        profile = src.profile
        nodata = src.nodata

    report = detect_quantisation(band, nodata)
    if not report["quantised"]:
        return {"changed": False, **report}

    # Float first: smoothing an integer array re-quantises the result.
    work = band.astype("float32")
    mask = None
    if nodata is not None:
        mask = band == nodata
        work[mask] = np.nan
        # Fill voids with the local mean so the kernel does not smear nodata.
        work = np.where(np.isnan(work), np.nanmean(work), work)

    smoothed = gaussian_filter(work, sigma=sigma, mode="nearest")
    if mask is not None:
        smoothed[mask] = nodata

    profile.update(dtype="float32", compress="deflate", tiled=True)
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(smoothed.astype("float32"), 1)

    return {"changed": True, "sigma": sigma, **report}

Casting to float before smoothing is the step that is easy to omit and completely defeats the exercise: a Gaussian filter applied to an integer array rounds its output back to integers, reproducing the staircase with slightly different treads.

Filling nodata before filtering matters for the same reason -compute_edges matters when shading. A Gaussian kernel that reaches into a nodata sentinel drags that value into every cell within its radius, producing a halo of wrong elevations around every void and coastline.

Performance Tuning and Cartographic Best Practices

  • Use the smallest sigma that works. A sigma of about one pixel removes a one-unit quantum on most grids. Larger values do remove more banding and also remove genuine small landforms — gullies, levees, dune crests — which is a worse outcome than the artefact.
  • Verify against a slope histogram. Before smoothing, the slope histogram of a quantised DEM shows discrete spikes. After, it should be continuous with its tail intact. A tail that has been truncated means the smoothing was too strong and real steep terrain was flattened.
  • Smooth the DEM, not the shading. Smoothing the derived hillshade blurs every real edge along with the artefact, softening ridgelines. Where the source cannot be modified, a bilateral filter on the shading preserves edges much better than a Gaussian.
  • Prefer a float source. Where a float DEM is available, use it. All of this is remediation for a data decision made upstream, and the remediation is never as good as not needing it.
  • Tile with an overlap. The filter reads neighbours, so a tiled deterracing pass needs an overlap of at least three sigma to avoid seams — the same requirement as the shading pass itself.
Checking that smoothing removed the artefact and not the terrain Two slope histograms. Before smoothing, the distribution consists of discrete spikes at regular slope values with gaps between them, and a tail extending to 40 degrees. After smoothing with a sigma of one pixel, the distribution is continuous and the tail still reaches 40 degrees. A third panel shows an over-smoothed case where sigma is four and the tail has been truncated at 22 degrees. before discrete spikes, tail to 40° sigma 1 — correct continuous, tail intact sigma 4 — too much tail truncated at 22° The tail is the diagnostic. Removing the spikes is easy; removing them while keeping the steep terrain is what distinguishes a correct sigma from a convenient one. Assert the 99th percentile slope is preserved to within a couple of degrees.
Any sigma removes the banding. Only a small one removes it without also removing the cliffs, and the slope histogram is what tells the two cases apart.
Three causes of banding, and the test that separates them A decision path from an observed band pattern. If band boundaries land on regular round elevation values the cause is source quantisation, fixed by smoothing the DEM. If the bands align with a grid rather than with contours the cause is resampling, fixed by re-deriving from the native grid. If boundaries are irregular and follow the landscape the bands are real terrain and must not be smoothed away. bands observed sample a transect boundaries at regular round elevations source quantisation — smooth the DEM before deriving bands align with a grid, not with contours resampling artefact — re-derive from the native grid boundaries irregular, following the landscape genuine terraces — do not smooth; this is the subject The third case is why an automatic deterracing pass should report what it changed, not run silently.
Only the first two are defects. Applying the fix for them to the third erases the landform the map was made to show, which no downstream check will notice.

Integration and Next Steps

Deterracing belongs in the DEM preparation step, ahead of every derivative — hillshade, slope, aspect and contours all inherit the artefact otherwise. Run it once per source raster and cache the result, so the per-sheet or per-tile pipeline consumes a clean surface. The assertion that the 99th percentile slope survives makes a good contract test of the kind described in the fundamentals overview: it is cheap, it runs on every rebuild, and it fails loudly when a new source arrives with a different quantum.

Frequently Asked Questions

How do I tell terracing from real terraced terrain?

Real terraces — agricultural benches, river strandlines, lava flows — follow the landscape and occur at irregular elevations that reflect its history. Quantisation bands occur at perfectly regular intervals matching the DEM’s vertical quantum, usually every whole metre or foot, and their boundaries land on round numbers. Extract elevations along a transect crossing the bands and look at where the values change: a constant step at round values is quantisation, and anything else is terrain.

Why does raising the z-factor make terracing worse?

Because the artefact is a discontinuity in the gradient and the z-factor multiplies the gradient. A one-metre riser spread across a fifty-metre run is a two per cent slope change, which is barely visible; exaggerate ten-fold and it becomes a twenty per cent change, which shades as a hard edge. Exaggeration amplifies whatever is in the derivative, and on a quantised surface a large part of what is in the derivative is the quantisation.

Does resampling the DEM remove terracing?

Partly. Bilinear or cubic resampling to a finer grid interpolates between quantised values and reduces the visible banding, but it also invents elevations that were never measured and can introduce artefacts aligned to the resampling grid. Where the DEM has to be resampled anyway it is a reasonable side benefit; as a deliberate fix, an explicit smoothing kernel is preferable because its strength is a single documented parameter that can be tuned and asserted against.

Should I smooth the DEM or the derived hillshade?

The DEM. Smoothing the shading blurs the artefact and every genuine edge together, so ridgelines soften at the same rate as the bands disappear. Smoothing the elevation surface removes the staircase before the gradient is taken, leaving the derivative clean and the ridges sharp. When the source raster genuinely cannot be modified — a read-only published product, say — a bilateral filter on the shading is the better fallback, because it smooths within regions while preserving the strong edges that carry the terrain’s structure.

Back to Terrain and Hillshade Automation