Short answer

Reprojecting a raster is not the same as reprojecting a vector: every output pixel is interpolated from the input grid, so the operation resamples your data whether you intend it to or not. A clean QGIS raster reprojection means deciding three things on purpose — the target CRS, the resampling method, and the output cell size — and then validating that NoData, extent, and alignment survived the warp. The engine underneath is GDAL's gdalwarp, and you get more control by calling it directly than by accepting the defaults in the dialog.

Why raster reprojection resamples your data

A vector reprojection moves vertices: a point at one coordinate becomes the same point at a new coordinate, with no loss. A raster is a fixed grid of cells, and the source grid almost never lines up with the destination grid after the projection math is applied. GDAL has to build a new regular grid in the target CRS and then fill each new cell by sampling the old one. That sampling step is resampling, and the method you choose determines whether class codes survive, whether elevation values stay physically meaningful, and how much high-frequency detail is smoothed away.

This is why "reproject" and "resample" are effectively a single operation for rasters. You cannot reproject a 10 m grid in UTM into a geographic CRS and keep the exact original values — the cells no longer correspond one-to-one.

Choosing the resampling method

The single most consequential decision is the resampling kernel. Pick it by data type, not by habit.

  • Nearest neighbour (-r near): copies the value of the closest source cell. Use it for categorical or thematic rasters — land cover, classified lithology, soil class, any raster where the pixel value is a code. Bilinear or cubic on a class raster invents new, meaningless class numbers like "4.3" at boundaries.
  • Bilinear (-r bilinear): weighted average of the 4 nearest cells. The default for continuous data such as DEMs, temperature surfaces, and reflectance. Smooth, predictable, no overshoot.
  • Cubic (-r cubic) and cubic spline (-r cubicspline): use a 4×4 and larger neighbourhood. Sharper-looking imagery, but cubic can overshoot near steep value changes — a problem at coastlines in a DEM, where it can produce elevations below the true minimum.
  • Average / mode (-r average, -r mode): use when downsampling to a coarser cell size. average for continuous data, mode for categorical data so the most common class wins the larger cell.

If you remember one rule: categorical data uses near or mode, continuous data uses bilinear, cubic, or average.

Set the output resolution deliberately

gdalwarp will pick an output resolution for you if you omit -tr, and its guess often does not match the source. Two failure modes:

  • Too coarse: you discard real detail the source contained.
  • Too fine (oversampling): you create a bigger file whose extra pixels are pure interpolation — no new information, just disk and processing cost, plus a false impression of precision.

Compute the input ground sampling distance first (gdalinfo input.tif reports pixel size in the source CRS units), then choose an output -tr in the target CRS units that keeps the same effective resolution. For a 30 m DEM going from EPSG:4326 to EPSG:32633 (UTM, metres), 30 m is the honest output size:

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:32633 \
  -r bilinear -tr 30 30 -tap \
  -dstnodata -9999 \
  -co TILED=YES -co COMPRESS=DEFLATE \
  dem_geographic.tif dem_utm.tif

-tap (target-aligned pixels) snaps the grid to clean coordinate multiples, which keeps tiles and repeated runs consistent. Remember that in geographic CRS the cell size is in degrees (0.000277° ≈ 30 m at the equator, less elsewhere), so resolution numbers are not comparable across CRS without converting units.

The QGIS dialog vs gdalwarp

In QGIS, Raster ▸ Projections ▸ Warp (Reproject) opens the Processing wrapper for gdalwarp (gdal:warpreproject). It is convenient, but two defaults catch people out: it often leaves the output resolution blank (GDAL then estimates a cell size that is rarely round), and it does not set NoData unless you expand the advanced parameters. For production work, drive gdalwarp from the command line or paste the equivalent flags into the dialog's "Additional command-line parameters" field.

Worked example: UTM DEM to ETRS89 / a national grid

Suppose you have a 10 m DEM in WGS 84 / UTM zone 31N (EPSG:32631) and you need it in ETRS89 / LAEA Europe (EPSG:3035) for a pan-European mosaic, keeping a clean 10 m cell size:

gdalwarp \
  -s_srs EPSG:32631 \
  -t_srs EPSG:3035 \
  -tr 10 10 \
  -r bilinear \
  -tap \
  -srcnodata -9999 -dstnodata -9999 \
  -co COMPRESS=DEFLATE -co PREDICTOR=2 -co TILED=YES \
  -multi -wo NUM_THREADS=ALL_CPUS \
  dem_utm31.tif dem_laea.tif

What each flag is doing and why it matters:

  • -tr 10 10 forces an exact 10 m square cell. Without it you get an arbitrary size like 9.873 m.
  • -tap (target-aligned pixels) snaps the output grid origin to multiples of the resolution. This is the flag that makes tiles mosaic seamlessly later — every tile shares the same grid lattice.
  • -srcnodata / -dstnodata keeps your sentinel (-9999) transparent instead of being averaged into real elevations at the edges.
  • -r bilinear because a DEM is continuous; cubic risks overshoot at the coast.
  • The creation options (-co) produce a compressed, tiled GeoTIFF. PREDICTOR=2 improves DEFLATE compression on continuous integer or float data.
  • -multi -wo NUM_THREADS=ALL_CPUS parallelises the warp.

For a land cover raster the only changes are -r near (or -r mode if downsampling) and dropping PREDICTOR=2.

Why borders go black or transparent

When the source and destination CRSs are not aligned, the reprojected footprint is a rotated quadrilateral inside a rectangular raster. The corner cells outside the data have to be filled. If you did not declare NoData, GDAL fills them with 0, which renders as black for imagery or as a real elevation of 0 m in a DEM — both wrong. Always set -dstnodata so those fill cells are flagged. In QGIS you can confirm it worked under Layer Properties ▸ Transparency ▸ Additional no data value.

DEM-specific warning: reproject before you derive

If your DEM is in a geographic CRS (degrees, e.g. EPSG:4326), do not run slope, aspect, or hillshade on it directly. One degree of longitude is not the same ground distance as one degree of latitude except at the equator, so a degree-based grid has anisotropic cells and slope comes out distorted. Reproject into a projected, metric CRS with equal-area or equal-distance properties for your extent first (a UTM zone or a national grid), then derive terrain products. This is the most common cause of "my slope values look wrong" reports.

Two more DEM traps that survive a correct horizontal warp:

  • Vertical datum and units are separate from the horizontal CRS. gdalwarp does not transform vertical datums by default; if your DEM is orthometric (e.g. EGM2008) and your workflow expects ellipsoidal heights, that is a distinct correction. Keep vertical units (metres vs feet) explicit.
  • NoData and voids: always carry -dstnodata. Without it, void cells get interpolated and create phantom terrain.

Validation: prove the warp was clean

Run a short check before the layer is used downstream:

  1. gdalinfo dem_laea.tif — confirm the SRS block reads EPSG:3035, the pixel size is exactly 10×10, and the NoData value is reported.
  2. Statistics check — compare min/max/mean of input and output (gdalinfo -stats). A DEM whose minimum dropped to -9999 means NoData is leaking into the statistics; whose minimum dropped below the true terrain low means cubic overshoot.
  3. Overlay against a trusted reference — load a known-good basemap or coastline and check alignment at three or four corners, not just the centre. Misalignment that grows toward the edges usually means a missing or wrong datum transformation.
  4. Pixel count sanity — width × height should match what the new extent and resolution imply.

Common mistakes and why they happen

  • Leaving output resolution blank. GDAL estimates a non-round cell size, and every later tile uses a slightly different grid, so mosaics show seams. Fix: always set -tr and -tap.
  • Bilinear on a class raster. Produces fractional class codes that break symbology and zonal statistics. Fix: -r near.
  • No NoData declared. Black corners, polluted statistics, wrong DEM minimums. Fix: set both -srcnodata and -dstnodata.
  • Trusting on-the-fly reprojection. QGIS displays mixed CRSs aligned on screen, which hides the fact that the file on disk was never reprojected. The map looks fine; the export or the analysis is wrong.
  • Ignoring the datum step. Moving between datums (e.g. NAD27 to NAD83, or any shift involving a grid like NTv2) needs the right transformation pipeline, not just a change of EPSG code.

Bathyl perspective

We treat a reprojected raster as a derived product with a recorded recipe: source CRS, target CRS, resampling method, resolution, and the exact gdalwarp call live next to the file. That way a colleague can reproduce the grid byte-for-byte and a reviewer can see at a glance whether the resampling choice fits the data type. The discipline is cheap; silently degraded elevation data discovered three months into a project is not.

Related reading

Sources