10.2 Raster Resampling Techniques: Nearest Neighbor, Bilinear Interpolation & Cubic Convolution

10.2 Raster Resampling Techniques: Nearest Neighbor, Bilinear Interpolation & Cubic Convolution

Core Principle: Resampling is an unavoidable mathematical necessity whenever a raster undergoes geometric manipulation—including coordinate reprojection, rotation, georeferencing rectification, or cell size coarsening/refinement. Because the cell centers of the transformed output grid rarely coincide with input cell centers, an interpolation algorithm must determine the output cell value. Matching the resampling technique to the underlying data measurement scale (categorical/nominal vs. continuous/ratio) is critical: using interpolation on discrete land cover classes produces nonsensical data corruption, while using Nearest Neighbor on elevation surfaces introduces severe stair-stepped aliasing.

---

1. The Computational Necessity of Raster Resampling

A raster dataset is a rigid, regular tessellation of square or rectangular grid cells arranged in orthogonal rows and columns. When a raster is subjected to any geometric transformation—such as:

• Reprojecting from a geographic CRS (e.g., WGS84 latitude/longitude in angular degrees) to a projected CRS (e.g., UTM Zone 18N in linear meters),
• Rectifying an unreferenced image using georeferencing transformation equations,
• Rotating a grid to match a flight track or street orientation, or
• Rescaling cell resolution (e.g., upsampling 30-meter Landsat data to 10-meter pixels, or downsampling 1-meter LiDAR elevation to 5-meter cells),

the original grid geometry is altered. The new output raster grid possesses its own regular, orthogonal cell spacing oriented to the target coordinate system. Consequently, the cell centers of the target output grid almost never match the cell centers of the source input grid.

   INPUT RASTER GRID (Source Space)               OUTPUT RASTER GRID (Target Space)
   Rotated / Transformed                           Orthogonal to Target CRS
       +-------+-------+                              +---------+---------+
      /       /       /                               |         |         |
     / (x1)  / (x2)  /                                |    ?    |    ?    |
    +-------+-------+         Transformation          | [Cell]  | [Cell]  |
   /       /       /      ---------------------->     | (Out 1) | (Out 2) |
  / (x3)  / (x4)  /                                   +---------+---------+
 +-------+-------+                                    |         |         |
                                                      |    ?    |    ?    |
                                                      +---------+---------+

Forward Mapping vs. Backward (Inverse) Mapping

Two fundamental computational paradigms exist for raster transformation:

1. Forward Mapping (Source-to-Target):
   • The algorithm iterates through every cell in the input raster, calculates its new spatial location in target space using the transformation equations, and deposits the pixel value into the destination grid.
   • The Fatal Flaws: Because the geometric transformation involves scaling, rotation, and shearing, mapped input cell centers do not land uniformly on target cell centers. This produces gaps (unassigned null holes) where no input pixel landed, and collisions (overlapping contention) where multiple input pixels compete for the same output cell. Forward mapping requires costly post-processing hole-filling algorithms and is rarely used for standard raster transformation.
2. Backward Mapping (Target-to-Source / Inverse Mapping):
   • The algorithm iterates systematically through every cell in the output raster grid from column 1 to $N$ and row 1 to $M$.
   • For each output cell center coordinate $(X_{out}, Y_{out})$, the inverse transformation equations compute the exact corresponding floating-point coordinate $(u_{in}, v_{in})$ back within the source raster coordinate space.
   • Because $(u_{in}, v_{in})$ is almost always a fractional decimal (e.g., column 142.38, row 89.71), it falls between the discrete integer pixel centers of the source image.
   • A resampling algorithm then inspects the surrounding source pixels and interpolates or assigns a definitive cell value to the output pixel.
   • The Operational Advantage: Backward mapping guarantees that every single cell in the output raster receives exactly one calculated value, completely preventing holes, data voids, and multi-pixel overlap collisions.

                         BACKWARD MAPPING PIPELINE
   
   1. For each Output Cell (Row, Col) in Target Raster:
      [ Output Cell Center: (X_out, Y_out) ]
                    |
                    v  (Apply Inverse Transformation Equations)
   2. Map backward into Source Raster Coordinate Space:
      [ Fractional Source Location: (u_in = 142.38, v_in = 89.71) ]
                    |
                    v  (Select Resampling Kernel)
   3. Sample Neighboring Source Pixels:
      +-------------------------------------------------------------+
      |  Nearest Neighbor:       Bilinear Interpolation:   Cubic:   |
      |  Takes single closest    Weighted average of       Weighted |
      |  cell center (1x1)       4 nearest centers (2x2)   16 cells |
      +-------------------------------------------------------------+
                    |
                    v
   4. Assign interpolated value to Output Cell (Row, Col)

---

2. Nearest Neighbor (NN) Resampling

Nearest Neighbor (NN) resampling is the simplest and computationally fastest resampling algorithm. It operates on a $1 \times 1$ sample space.

Algorithmic Mechanics

Given the fractional coordinate $(u, v)$ in the source raster, the algorithm computes the Euclidean distance to the four surrounding integer pixel centers and selects the single cell whose center is physically closest:

$V_{out} = V(\text{round}(u), \text{round}(v))$

No mathematical averaging, weighted interpolation, or value modification occurs. The output pixel is assigned the unaltered, raw digital number (DN) of that single closest neighbor.

                     NEAREST NEIGHBOR KERNEL (1x1)
   
        (col, row)               (col+1, row)
            +-------------------------+
            |                         |
            |       (u, v)            |
            |         *               |
            |     (Target)            |
            |    Distance d1 < d2     |
            |                         |
            +-------------------------+
        (col, row+1)             (col+1, row+1)
        
   Output value is assigned directly from the top-left pixel (col, row)
   because distance d1 to its center is the shortest.

Core Properties and Applications

• Value Preservation: Because output values are copied directly from source cells without mathematical manipulation, no new pixel values are ever created. The original data range, unique value list, frequency histogram, and discrete classifications remain completely intact.
• Mandatory for Categorical and Discrete Data: Nearest Neighbor is the only standard resampling method appropriate for categorical, nominal, or ordinal datasets.
  • Land Use / Land Cover (LULC) rasters (e.g., National Land Cover Database - NLCD).
  • Soil taxonomic classes (e.g., Soil Survey Geographic Database - SSURGO).
  • Zoning boundaries, parcel identification numbers (PINs), administrative district codes.
  • Binary suitability masks ($0 = \text{unsuitable}, 1 = \text{suitable}$).
  • Raster attribute tables (RAT) with text strings or categorical lookup keys.
• Positional Accuracy & Spatial Artifacts: Because Nearest Neighbor snaps to the nearest cell center, it introduces a horizontal positional displacement of up to $\pm 0.5$ of a pixel width (half-pixel shift). When applied to continuous imagery or rotated grids, it produces noticeable visual artifacts: stair-stepping (aliasing), jagged edges along linear features, and blocky pixel duplication.
• Computational Speed: Extremely fast ($O(1)$ lookup per output cell center, zero floating-point arithmetic operations).

---

3. Bilinear Interpolation (BI)

Bilinear Interpolation (BI) calculates output cell values using a distance-weighted average of the 4 nearest cell centers ($2 \times 2$ pixel neighborhood).

Algorithmic Mechanics

The algorithm performs linear interpolation sequentially along two orthogonal dimensions (first along the horizontal row axis, then along the vertical column axis, or vice versa).

                   BILINEAR INTERPOLATION KERNEL (2x2)
   
        P00 (u0, v0)                 P10 (u1, v0)
            *---------------------------*
            |      |                    |
            |      | Δv                 |
            |------* (u, v) ------------|  <-- Interpolate between P00 and P10
            |      |                    |      to find intermediate R1
            |      | 1 - Δv             |  <-- Interpolate between P01 and P11
            |      |                    |      to find intermediate R2
            *---------------------------*
        P01 (u0, v1)                 P11 (u1, v1)
        
   Final value at (u, v) is interpolated vertically between R1 and R2.

Let the fractional coordinates within the bounding unit square be $\Delta u = u - u_0$ and $\Delta v = v - v_0$, where $0 \le \Delta u, \Delta v < 1$, and let the four neighboring cell values be $P_{00}, P_{10}, P_{01}, P_{11}$:

1. Interpolate horizontally across row 0: $R_1 = (1 - \Delta u) P_{00} + \Delta u \cdot P_{10}$
2. Interpolate horizontally across row 1: $R_2 = (1 - \Delta u) P_{01} + \Delta u \cdot P_{11}$
3. Interpolate vertically between $R_1$ and $R_2$: $V_{out} = (1 - \Delta v) R_1 + \Delta v \cdot R_2$

Expanding into the full bilinear polynomial equation:

$V_{out} = (1 - \Delta u)(1 - \Delta v) P_{00} + \Delta u (1 - \Delta v) P_{10} + (1 - \Delta u) \Delta v \cdot P_{01} + \Delta u \Delta v \cdot P_{11}$

Core Properties and Applications

• Surface Continuity: Produces a continuous, smooth output surface free from the jagged stair-stepping of Nearest Neighbor.
• Ideal for Continuous Surfaces: Highly recommended for continuous ratio or interval datasets where values vary smoothly across space:
  • Digital Elevation Models (DEMs), Digital Surface Models (DSMs), and bathymetric grids.
  • Slope, aspect, and topographic curvature surfaces.
  • Atmospheric data: ambient temperature, barometric pressure, precipitation grids.
  • Continuous vegetation indices: Normalized Difference Vegetation Index (NDVI), Leaf Area Index (LAI).
  • Radiometrically calibrated geophysical datasets (gravity, magnetic anomalies).
• The Low-Pass Filtering (Smoothing) Effect: Because Bilinear Interpolation averages four neighboring values, it acts as a spatial low-pass filter. It smooths out high-frequency noise but also dampens localized extremes: sharp mountain peaks are lowered slightly, narrow ravines are raised, and crisp feature boundaries are softened.
• Range Boundedness: The output value is mathematically guaranteed to fall strictly within the minimum and maximum range of the four contributing input cells: $\min(P_{ij}) \le V_{out} \le \max(P_{ij})$. It will never generate values outside this local range.

[!CAUTION] Exam Trap: Applying Bilinear Interpolation to Categorical Land Cover. Never apply Bilinear Interpolation to discrete or categorical data! If an output cell falls between a cell with code $11$ (Water) and a cell with code $41$ (Deciduous Forest), bilinear interpolation will calculate a weighted average of approximately $26$. In standard classification schemes, code $26$ may represent Low-Intensity Residential or an unassigned null class. The output raster becomes corrupted with meaningless intermediate values and the raster attribute table is destroyed.

---

4. Cubic Convolution (CC)

Cubic Convolution (CC) is an advanced resampling technique that evaluates a $4 \times 4$ kernel of the 16 nearest cell centers. It fits a piecewise cubic polynomial surface that closely approximates the theoretically ideal sinc reconstruction filter $(\sin(\pi x) / (\pi x))$.

                     CUBIC CONVOLUTION KERNEL (4x4)
   
            *       *       *       *
            (0,3)   (1,3)   (2,3)   (3,3)
            
            *       *       *       *
            (0,2)   (1,2)   (2,2)   (3,2)
                        (u, v)
            *       *     x *       *
            (0,1)   (1,1)   (2,1)   (3,1)
            
            *       *       *       *
            (0,0)   (1,0)   (2,0)   (3,0)
            
   Evaluates 16 neighbor cells using a cubic spline weighting function.

Algorithmic Mechanics

The continuous cubic weighting function $W(d)$ evaluates the distance $d$ from the fractional coordinate $(u, v)$ to each of the 16 cell centers (typically evaluated for $|d| < 2$ with parameter $a = -0.5$):

$W(d) = \begin{cases} (a + 2)|d|^3 - (a + 3)|d|^2 + 1 & \text{for } |d| \le 1 \\ a|d|^3 - 5a|d|^2 + 8a|d| - 4a & \text{for } 1 < |d| < 2 \\ 0 & \text{for } |d| \ge 2 \end{cases}$

Notice that for distances between $1$ and $2$ pixel units, the weighting function $W(d)$ dips into negative values. These negative sidelobes act as a sharpening mechanism, enhancing edges and high-frequency spatial gradients.

Core Properties and Applications

• Superior Visual Acuity: By preserving edge sharpness while smoothing continuous gradients, Cubic Convolution produces the most visually appealing results for continuous optical imagery, high-resolution aerial photography, orthophotography, and satellite panchromatic bands.
• Computational Burden: Highly compute-intensive. Evaluating 16 cell lookups and solving multiple cubic polynomial equations per output pixel requires roughly 4 to 8 times more processing cycles than Bilinear Interpolation and 16+ times more than Nearest Neighbor.
• The Overshoot and Undershoot (Ringing) Phenomenon: Because the cubic weighting function contains negative coefficients, the interpolated output value can mathematically exceed the maximum or fall below the minimum value present among the 16 input cells.
  • Near sharp contrast boundaries (e.g., a dark water body adjacent to a bright sandy beach, or a bright building roof against dark asphalt), Cubic Convolution can produce localized artificial "halos" or ringing artifacts.
  • In 8-bit imagery (valid range 0 to 255), calculated values can dip below 0 (causing negative values that clip to 0) or exceed 255 (clipping to 255).
  • In continuous elevation models, cubic convolution can introduce artificial depression pits or false ridge crests near cliff faces.

---

5. Majority / Mode Resampling and Aggregation Methods

When converting a categorical raster from a fine resolution to a coarser resolution (downsampling or aggregating—such as converting a 1-meter drone land cover classification into a 30-meter regional model), standard Nearest Neighbor can randomly pick an anomalous or unrepresentative single pixel, while Bilinear Interpolation mathematically corrupts the classes.

To address this specific challenge, GIS software provides Majority (Mode) resampling:

• Mechanism: The algorithm evaluates all input cells within a specified neighborhood kernel (e.g., $3 \times 3$, $4 \times 4$, or the exact bounding footprint of the coarser output cell) and determines the statistical mode—the most frequently occurring class value.
• Tie-Breaking: If two classes tie for the majority count, software typically applies a systematic fallback rule (e.g., selecting the minimum class value, or prioritizing the center pixel).
• Primary Application: Aggregating discrete land cover, wetland classifications, or zoning layers to coarser spatial resolutions. It filters out isolated, single-pixel noise (speckle) while preserving dominant categorical patches.

Summary of Resampling Techniques

---

6. Common GISP Exam Traps & Pitfalls

[!CAUTION] Exam Trap 10.2.1: Resampling Method Selection by Data Scale. This is among the most frequently tested concepts on the GISP exam. If an exam question describes reprojecting a Land Use/Land Cover raster, soil map, or timber classification layer, the ONLY correct answer is Nearest Neighbor (or Majority if downsampling). Any choice recommending Bilinear Interpolation or Cubic Convolution for discrete data is completely false, as averaging class code numbers destroys the categorical integrity of the dataset.

[!CAUTION] Exam Trap 10.2.2: The Overshoot Property of Cubic Convolution. Exam questions often test why an analyst discovered negative surface reflectance values (e.g., $-4$) or elevation pits below sea level in a dataset that originally had a minimum value of zero. The cause is Cubic Convolution overshoot/undershoot (ringing) occurring near high-contrast edges due to the negative sidelobes of the cubic weighting function. Bilinear interpolation can never produce a value outside the local min/max range.

[!CAUTION] Exam Trap 10.2.3: Forward Mapping vs. Backward Mapping. A common conceptual distractor suggests that GIS software projects each input pixel forward into the output coordinate system. If software used forward mapping, the output raster would be filled with unassigned void holes and overlapping multi-pixel collisions. Modern GIS software universally employs backward (inverse) mapping, stepping through each output cell center and interpolating backward from the source raster space.