Projecting Rasters: Why the Resampling Method Matters

GIS Concepts · background · by Jeff Jenness

Vector data project cleanly. A point, a polyline or a polygon is a set of precisely defined vertices, and each vertex can be transformed into a new coordinate system and back again without changing. Rasters are different. A raster is a set of rows and columns, and when it is projected it has to become a different set of rows and columns, with (almost always) square cells, laid out horizontally and vertically in the new coordinate system. The new cells do not line up with the old ones, so every new cell has to be handed a value by some rule, and the choice of rule — the resampling method — decides whether your projected raster is faithful or quietly corrupted. Project a raster out and back, and you will almost never get the original raster.

This page condenses a topic Jeff teaches at more length in the Projections and Datums module of his GIS training course, where the lab exercises let you watch a raster change as it is projected out and back.

The short version Use nearest neighbor (or majority) for categorical rasters such as land cover, where inventing a value that never existed would be nonsense. Use bilinear interpolation for continuous surfaces such as elevation, and do it before you derive slope, aspect, curvature, topographic position or roughness from them. A DEM projected with nearest neighbor looks fine; its derivatives do not.

Two grids that do not line up

The green grid below is the original raster; the blue grid is the same ground after projection into a coordinate system whose axes point a little differently. Each blue cell overlaps several green cells, and no blue cell overlaps exactly one. The question every projection has to answer is: what value goes in each blue cell?

A green square grid overlaid by a blue grid of the same cell size rotated by a few degrees, so every blue cell overlaps parts of several green cells
The original raster (green) and the projected raster's cells (blue). Every new cell straddles several old ones.

Four ways to fill the new cells

It is easier to think about with each cell drawn as a dot at its center. Green dots are the original cell centers; blue dots are the centers of the new cells. ArcGIS offers four rules for assigning a value to each blue dot.

Green and blue dot grids; one blue dot is highlighted together with the single nearest green dot
Nearest neighbor. The blue cell takes the value of the one nearest green cell. The new raster can only contain values that existed in the original, which is what you want for categorical data.
Green and blue dot grids; one blue dot is highlighted together with the four nearest green dots around it
Bilinear interpolation. The blue cell takes the distance-weighted average of the four nearest green cells. Best for continuous surfaces such as a DEM.

The other two are variations on these. Cubic convolution averages the sixteen nearest original cells instead of four, which suits continuous surfaces but smooths them noticeably. Majority looks at a 4 × 4 block of original cells and takes the most common value, which, like nearest neighbor, guarantees the new raster holds only values the original had, and so is a second good choice for categorical rasters.

The nearest-neighbor problem: repeats and skips

Nearest neighbor is the default in many tools, and for a categorical raster it is the right default. The trouble comes when it is applied to a continuous surface. Look again at the two dot grids. Because the blue grid is slightly rotated or slightly differently spaced, there are places where two blue cells share the same nearest green cell, and places where a green cell is nearest to no blue cell at all. The original values get repeated in some new cells and skipped in others, and because the two grids drift in and out of phase, the repeats and skips fall on a regular cycle across the raster.

Green and blue dot grids with four red ellipses circling places where two blue dots share one nearest green dot, or a green dot has no blue dot nearby
Circled: places where one original value will be copied into two new cells, and where an original value will be dropped. The pattern recurs at a regular spacing across the whole raster.

For a categorical raster this is harmless; a land cover class repeated in one extra cell is still that class. For an elevation surface it means the projected DEM contains a faint but perfectly regular lattice of little terraces: pairs of identical cells where the ground should have been rising, and small jumps where a value was skipped. The elevations themselves are barely affected. A nearest-neighbor pull can move a value by at most half a cell diagonal, about 21 m on a 30-m grid, which in steep ground is a few meters of elevation and in gentle ground is nothing. If all you need is elevation, you will never notice. But slope is the difference between neighboring cells, and a difference of zero followed by a double-sized jump is exactly what the repeats and skips produce.

Why the DEM looks fine and the slope does not

Here is the effect on real data. The 30-m DEM from the Corridor Design Tutorial was projected from UTM Zone 12N to the North America Albers Equal Area Conic system twice, once with nearest neighbor and once with bilinear interpolation, and slope was then computed from each projected DEM.

Four panels: the two projected DEMs on top, indistinguishable from each other; the slope from the nearest-neighbor DEM below left, covered in a fine cross-hatched lattice of stripes; the slope from the bilinear DEM below right, clean
Top: the DEM projected with nearest neighbor (left) and bilinear interpolation (right). You cannot tell them apart. Bottom: slope derived from each. The nearest-neighbor version is overlaid with a lattice of stripes that has nothing to do with the terrain; the bilinear version is clean. The stripes run diagonally here because the Albers grid is rotated a few degrees from the UTM grid. Between a projected system and a geographic one they typically run horizontally and vertically.

Everything derived from the surface inherits the lattice: aspect, curvature, the topographic position index, hillshades, and every roughness measure, all of which are built from differences between neighboring cells. Some of them, such as curvature and the ruggedness indices, amplify it, because they are differences of differences. A slope-position or landform classification made from such a DEM will carry stripes of “steep slope” and “ridge” marching across flat valley floors, and a habitat model built on that factor will carry them too.

Hydrology suffers in a different way. The lattice of little terraces and jumps acts like a network of shallow channels laid across the surface, and flow-direction and flow-accumulation tools, which send each cell's water to its lowest neighbor, will follow those channels. Streams modeled from a nearest-neighbor DEM run in straight lines along the artifact grid, turn at right angles, and pool in places where no water would ever pool, and the topographic wetness index and every watershed delineation built on top of them go wrong with them.

What to do

In these tools The Clip Data to Analysis Area tool projects rasters with nearest neighbor, always, because it handles categorical and continuous rasters in the same run and nearest neighbor is the one method that is safe for both. So if a DEM arrives in a different coordinate system from your analysis, project it yourself with bilinear interpolation first and hand the projected copy to the clip tool, which then clips it with no further resampling. Every tool in the Topographic Analysis group — slope position, landform, TPI, the roughness indices — assumes the DEM it is given was produced this way, and the Topographic Wetness Index tool most of all, since it computes flow direction and flow accumulation from the DEM and every false channel becomes a false stream.

Credits

Adapted from the Projections and Datums module of Jeff Jenness's GIS training course, where the dot-grid diagrams above are used to teach resampling and where a lab exercise has students project a land cover raster out to geographic coordinates and back again to watch how much of it changes. The slope demonstration was computed for this page from the corridor tutorial's DEM.