gdalwarp for Raster Reprojection

Short answer

gdalwarp is GDAL's reprojection, resampling, and mosaicking engine. Unlike gdal_translate, which stays on the source pixel grid, gdalwarp builds a new pixel grid in the target CRS and resamples the source into it. You reach for it whenever the coordinate system changes, the resolution changes, several tiles must merge into one, or a raster must be clipped to a polygon. The two decisions that determine whether the output is correct are the target CRS (-t_srs) and the resampling method (-r).

A minimal reprojection from UTM to Web Mercator:

gdalwarp -t_srs EPSG:3857 -r bilinear in_utm.tif out_3857.tif

GDAL reads the source CRS from the file, so -s_srs is only needed when the source CRS is missing or wrong.

How a warp actually works

gdalwarp does not move existing pixels; it inverts. For each output pixel it computes the corresponding location in the source raster (via the PROJ coordinate operation) and then samples the source there. Because that source location almost never lands exactly on a source pixel centre, a resampling kernel interpolates a value. This is why the resampling choice is not cosmetic — it changes the numbers.

This inversion also explains the borders. When you reproject, the source's rectangular footprint becomes a rotated quadrilateral inside the new axis-aligned bounding box. The triangular corners have no source data and must be filled — with NoData if you configured it, or with zeros if you did not.

Choosing the resampling method

The default is near, which is the correct default only for categorical data. Reprojecting a DEM with the default produces a visibly stair-stepped surface and biased slope derivatives. Always pass -r explicitly.

Worked example: reproject a DEM with a controlled grid

You want a UTM 32N DEM (dem_32632.tif) in ETRS89 / UTM 32N-equivalent for a national dataset, output at exactly 10 m, aligned to a clean grid, with NoData handled:

gdalwarp \
  -t_srs EPSG:25832 \
  -tr 10 10 \
  -tap \
  -r cubic \
  -srcnodata -9999 -dstnodata -9999 \
  -co COMPRESS=DEFLATE -co PREDICTOR=3 -co TILED=YES \
  -multi -wo NUM_THREADS=ALL_CPUS \
  dem_32632.tif dem_25832.tif

What each block does:

• -tr 10 10 forces 10 m square output pixels. Without it, GDAL estimates a resolution from the warp geometry and you get an odd number like 9.873 m.
• -tap ("target aligned pixels") snaps the output extent so pixel boundaries fall on multiples of the resolution. This is essential when several rasters must align exactly for stacking or tiling — without -tap, two independently warped tiles can be offset by a fraction of a pixel and never line up.
• -srcnodata/-dstnodata ensure fill pixels are recognised on input and written as the sentinel on output, keeping the corner triangles out of your statistics.
• -multi -wo NUM_THREADS=ALL_CPUS parallelises the warp; on a large national mosaic this is a multiple-times speed-up.

Worked example: mosaic plus clip to a basin

To merge several DEM tiles and clip the result to a watershed polygon:

gdalwarp \
  -t_srs EPSG:25832 -tr 10 10 -tap -r cubic \
  -cutline basin.gpkg -cl basin_layer -crop_to_cutline \
  -dstnodata -9999 \
  -co COMPRESS=DEFLATE -co TILED=YES \
  tile_*.tif basin_dem.tif

gdalwarp accepts multiple inputs and mosaics them in one pass; later inputs overwrite earlier ones in overlap zones. -cutline supplies the clip geometry (from a GeoPackage, Shapefile, or GeoJSON), -cl names the layer, and -crop_to_cutline shrinks the output extent to the polygon's bounding box rather than keeping the full mosaic extent. Everything outside the polygon becomes -dstnodata. If the cutline is in a different CRS than -t_srs, GDAL reprojects it automatically.

Reproject vs assign: a critical distinction

gdalwarp -t_srs resamples pixels into a new CRS. If instead the raster's CRS tag is simply wrong but the pixels are already in the right system, you do not warp — you assign with gdal_translate -a_srs or gdal_edit.py -a_srs. Warping a correctly-located raster to "fix" a label resamples it needlessly and degrades it. Diagnose first with gdalinfo: if the corner lon/lat are wrong, the label is wrong (assign); if the corners are right but you need a different projection, warp.

Controlling the transformation pipeline

For high-accuracy work, the datum transformation matters as much as the projection. PROJ may offer several operations between two datums (e.g. with or without a grid shift). You can pin one with -ct (coordinate transformation pipeline string) or inspect the candidates with projinfo -s EPSG:4326 -t EPSG:25832. For most national datasets the default operation with the bundled NTv2/HTDP grids is correct, but for survey-grade reprojection between older national datums, verify which grid PROJ selected.

Common pitfalls and why they happen

• Default near resampling on a DEM. Produces blocky elevation and corrupts every downstream derivative (slope, aspect, hillshade). The kernel default exists for categorical safety, not continuous accuracy.
• Black borders in the output. The rotated footprint leaves empty corners. Missing -dstnodata stores them as zeros, which then dominate a min/max stretch and create dark wedges. Always set NoData on a warp.
• Mosaic tiles that do not align. Without -tap, independently warped tiles have sub-pixel offsets and leave seams. -tap plus a fixed -tr is the alignment guarantee.
• -s_srs used to relabel. -s_srs tells GDAL how to interpret the input before warping; it does not fix a genuinely mislabelled file's geometry. If the input CRS is wrong, fix the assignment before warping.
• Reprojecting to Web Mercator for analysis. EPSG:3857 is fine for web tiles but its area and distance distortion away from the equator makes it a poor choice for measurement. Analyse in an equal-area or appropriate UTM zone, reproject to 3857 only for display.

Validation

After a warp, run gdalinfo -stats on the output and confirm: the CRS is the requested EPSG, the pixel size matches -tr, the NoData value is present and excluded from the statistics, and the corner coordinates fall where expected. For a mosaic, spot-check seam areas at full resolution for value discontinuities. For a clip, overlay the cutline polygon in QGIS and confirm the data stops at the boundary. Keep the exact gdalwarp command with the dataset so the reprojection is reproducible.

Bathyl perspective

We treat the resampling method and -tap alignment as non-negotiable parameters on every warp, because both cause errors that are invisible in a thumbnail and only surface when someone measures a slope or stacks two layers. Reprojection is where geometric accuracy is won or lost, so it is documented, parameterised, and reproducible rather than clicked once and forgotten.

Related reading

• gdalinfo for Raster Quality Checks
• gdal_translate for Raster Conversion
• GDAL Build VRT Explained
• Shapefile vs GeoPackage vs GeoJSON
• Spatial data products

Sources

• GDAL gdalwarp
• GDAL warp resampling methods
• PROJ coordinate operations
• GDAL gdalinfo
• EPSG:3857 Web Mercator