Short answer

Resampling a DEM changes the elevation surface, and because slope, aspect, hillshade, curvature, and drainage are all computed from a cell and its neighbours, resampling changes those derivatives too. To keep the result meaningful: use bilinear or cubic interpolation for continuous elevation (never nearest-neighbour), declare NoData explicitly so voids are not blended into real terrain, keep horizontal and vertical units consistent, and pick a target cell size that matches the scale of the landform you care about, not the highest number you can get. Higher resolution is not automatically better — fine grids carry noise that derivative operators amplify.

Why resampling alters terrain derivatives

A DEM is a grid of elevations. A derivative like slope is estimated from the elevation differences across a moving window (commonly the 3×3 Horn method used by GDAL and ArcGIS). The gradient therefore depends on the horizontal distance between cell centres — the cell size. Make cells larger and each slope estimate averages over more ground, smoothing peaks and filling hollows; make them smaller and the operator responds to every bump, including sensor noise.

This is why a slope figure is meaningless without its resolution. A talus slope might report a 45° maximum at 1 m LiDAR resolution and only 32° at 30 m SRTM, because the coarse grid cannot resolve the steep micro-relief. Neither number is "wrong"; they answer different questions. The rule is: report the cell size with every terrain statistic, and resample to a common grid before comparing two surfaces.

Choosing a resampling method

The method controls how new cell values are computed when the grid changes:

  • Nearest neighbour copies the closest source value. It preserves exact values and is correct for categorical rasters (land cover, geological class) — but for elevation it produces blocky stair-steps that inject false flats and cliffs, corrupting slope and hillshade. Do not use it for DEMs.
  • Bilinear interpolates from the four nearest cells. Smooth, predictable, and the safe default for continuous elevation in both up- and downsampling.
  • Cubic / cubic spline uses a 4×4 neighbourhood for a smoother surface; it can produce slight overshoot (values above/below the true range) near sharp breaks like cliff edges or building footprints in DSMs.
  • Average / mode are aggregation methods for downsampling: average is well suited to reducing a fine DEM to a coarser one because it represents the mean elevation of the cells being merged, rather than picking a single sample.

A practical guideline: upsample (finer) with bilinear or cubic; downsample (coarser) with average when you want a faithful coarse representation, or bilinear when you simply need grid alignment.

Worked example with gdalwarp

Resample a 30 m DEM to a 10 m grid in UTM 31N, with NoData handled correctly:

gdalwarp -t_srs EPSG:32631 \
  -tr 10 10 \
  -r bilinear \
  -srcnodata -9999 -dstnodata -9999 \
  -tap \
  -co COMPRESS=DEFLATE -co PREDICTOR=2 -co TILED=YES \
  dem_30m.tif dem_10m.tif

Key flags and why they matter:

  • -tr 10 10 sets the target cell size in the target CRS units (metres here).
  • -r bilinear chooses the interpolation; switch to average if you were going coarser.
  • -srcnodata/-dstnodata declare the void value so the resampler does not average -9999 into adjacent valid cells.
  • -tap ("target aligned pixels") snaps the output grid to clean coordinates, which makes later stacking and differencing exact.
  • The creation options keep the GeoTIFF compressed and tiled for downstream performance.

Then derive slope on the resampled grid:

gdaldem slope dem_10m.tif slope_deg.tif -compute_edges
gdaldem hillshade dem_10m.tif hs.tif -z 1.0 -az 315 -alt 45

-compute_edges avoids a NoData fringe one pixel wide around the dataset border.

The units trap

gdaldem slope assumes elevation (z) and horizontal (x/y) units are the same. With a DEM in a projected CRS (metres in x, y and z), slope is correct directly. With a DEM in a geographic CRS (EPSG:4326 — degrees in x/y but metres in z), the horizontal "distance" between cells is in degrees while z is in metres, so raw slope is badly wrong. Either reproject to a metric CRS first (preferred), or use the -s scale factor (gdaldem slope -s 111120) to approximate the degree-to-metre conversion. Equally, confirm whether you want slope in degrees (geometry, hazard thresholds) or percent rise (engineering, drainage) — they are not interchangeable, and gdaldem slope -p switches to percent.

Don't reach for the finest DEM by reflex

A 1 m LiDAR DSM includes tree canopy, vehicles, and survey noise; running slope on it produces a speckled, near-unusable layer for regional landform analysis. For watershed or geomorphic screening you often want a smoother landform signal — a bare-earth DTM resampled to 5–10 m, optionally low-pass filtered. Matching resolution to the question (and to the source's true accuracy, not just its pixel size) is more important than maximising pixels.

Common pitfalls and why they happen

  • Nearest-neighbour on elevation. Cause: copying a categorical-raster habit. It stair-steps the surface and ruins slope/hillshade. Use bilinear/cubic.
  • Undeclared NoData. Cause: relying on defaults. -9999 gets averaged into real cells along voids and edges, producing impossible elevations. Always pass -srcnodata/-dstnodata.
  • Slope in a geographic CRS. Cause: skipping reprojection. Degrees-vs-metres mismatch yields wrong gradients. Reproject to a metric CRS first.
  • Comparing derivatives across mismatched grids. Cause: differencing two DEMs at different cell sizes/methods. You measure the processing, not the terrain. Resample to a common grid and process identically.
  • Maximising resolution blindly. Cause: assuming finer is better. Fine grids amplify noise; match resolution to the landform scale and the data's real accuracy.

QA and validation

After resampling, inspect before trusting: check the output gdalinfo for the expected cell size, CRS, NoData value, and a plausible min/max elevation (no -9999 leaking into the stats). View the hillshade and look for striping, blocky steps, void halos, or edge artifacts. If you differenced two epochs, the difference raster should be near zero over stable ground — a non-zero bias usually means a vertical datum or unit mismatch, not real change. Where you need a defensible number, state the DEM resolution, source, vertical datum/units, and resampling method alongside it.

Bathyl perspective

Bathyl treats terrain derivatives as evidence with a stated resolution and method, not decorative relief. Each slope, aspect, or hillshade product ships with the cell size, source DEM, vertical units, and processing parameters that produced it, so a reviewer can tell what the layer supports and where interpretation becomes speculative.

Related reading

Sources