Changing grid resolution and filling missing values
2026-02-27
When the European Space Agency launched Sentinel-2 in 2015, it offered something Landsat users had never had: 10-metre resolution imagery in the visible and near-infrared bands, freely available every five days at any point on Earth. Research groups that had built workflows around Landsat’s 30-metre pixels immediately faced a practical problem: their existing DEMs, land-cover maps, and ancillary datasets were at 30 metres. To combine them with Sentinel-2 data, someone had to resample something — either the Sentinel image to a coarser grid or the existing datasets to a finer one. The choice of resampling method matters more than most users realise: nearest-neighbour resampling is fast and preserves discrete class boundaries, but produces blocky artefacts in continuous data; bilinear interpolation is smooth but blurs edges and can create fractional class values from categorical data.
Resampling — changing the grid spacing of a raster — is one of the most common operations in spatial data processing, and it is unavoidable whenever data from different sources must be combined in a common grid. It is also a more mathematically interesting problem than it appears. Estimating the value at a new grid position from surrounding known values is an interpolation problem; the method chosen controls how smooth or sharp the result is, whether edges are preserved or blurred, and how well the output honours the values at original sample points. Nearest-neighbour, bilinear, and bicubic methods represent a progression of polynomial order with corresponding trade-offs in smoothness, accuracy, and computational cost. This model derives each method from first principles and shows when each is appropriate.
How do you convert a 30-meter resolution satellite image to 10-meter resolution?
Resampling changes the grid spacing of raster data: - Upsampling (zoom in): Fewer pixels → more pixels (decrease cell size) - Downsampling (zoom out): More pixels → fewer pixels (increase cell size) - Reprojection: Change coordinate system (requires resampling)
Interpolation estimates values between measured points: - DEM creation: Points from LiDAR → continuous elevation surface - Climate mapping: Weather stations → temperature grid - Missing data: Cloud gaps in satellite imagery → fill with neighbors
The mathematical question: Given values at known grid positions, how do we estimate values at new positions?
Trade-offs: - Speed vs. quality: Fast methods produce blocky results - Smoothness vs. fidelity: Smooth methods may blur sharp features - Artifacts: Poor methods create visual or analytical problems
Input raster: Grid with cell size \Delta x_1, values z_{i,j}
Output raster: New grid with cell size \Delta x_2
Challenge: Output pixel centers don’t align with input pixels.
Example: Convert 30m DEM to 10m DEM - Input: 100×100 grid (30m cells) - Output: 300×300 grid (10m cells) - Must estimate elevation at 90,000 new positions from 10,000 known values
1. Nearest Neighbor
Use value from closest input pixel.
z_{\text{new}}(x, y) = z_{i,j} \text{ where } (i,j) = \text{nearest to } (x,y)
Pros: Fast, preserves original values (good for categorical data)
Cons: Blocky appearance, discontinuous
2. Bilinear Interpolation
Weighted average of 4 nearest pixels.
z(x, y) = (1-u)(1-v)z_{0,0} + u(1-v)z_{1,0} + (1-u)vz_{0,1} + uvz_{1,1}
Where u, v are fractional positions within cell.
Pros: Smooth, continuous
Cons: Blurs edges, creates new values not in original data
3. Bicubic Interpolation
Weighted average of 16 nearest pixels using cubic polynomials.
Pros: Very smooth, preserves gradients
Cons: Slower, can overshoot (create values outside original range)
Input pixel (i, j) with value z_{i,j}
Pixel center in world coordinates:
x = x_0 + i \cdot \Delta x y = y_0 - j \cdot \Delta y
(Note: y decreases with row index—image convention)
Query point (x_q, y_q): Where do we want to estimate z?
Find containing cell:
i_q = \lfloor (x_q - x_0) / \Delta x \rfloor j_q = \lfloor (y_q - y_0) / \Delta y \rfloor
Fractional position within cell:
u = \frac{x_q - x_0}{\Delta x} - i_q v = \frac{y_q - y_0}{\Delta y} - j_q
Where $0 u, v < 1$.
Algorithm:
i_nearest = round((xq - x0) / dx)
j_nearest = round((yq - y0) / dy)
z_interp = z[i_nearest, j_nearest]Equivalent: Round u, v to 0 or 1.
Decision boundary: Each output pixel assigned to nearest input pixel.
Result: Piecewise constant (stair-step appearance).
Four corners of cell:
Linear interpolation in x:
z_0 = (1-u)z_{0,0} + u z_{1,0} (bottom edge) z_1 = (1-u)z_{0,1} + u z_{1,1} (top edge)
Linear interpolation in y:
z(u, v) = (1-v)z_0 + v z_1
Expanded form:
z = (1-u)(1-v)z_{0,0} + u(1-v)z_{1,0} + (1-u)v z_{0,1} + uv z_{1,1}
Weights sum to 1:
(1-u)(1-v) + u(1-v) + (1-u)v + uv = 1 ✓
Properties: - Continuous (no jumps) - C^0 continuous (values match, derivatives don’t) - Exact at grid points: z(0,0) = z_{0,0}, etc.
Uses 16 pixels (4×4 neighborhood centered on query point).
Cubic convolution kernel:
w(x) = \begin{cases} (a+2)|x|^3 - (a+3)|x|^2 + 1 & \text{if } |x| \leq 1 \\ a|x|^3 - 5a|x|^2 + 8a|x| - 4a & \text{if } 1 < |x| < 2 \\ 0 & \text{if } |x| \geq 2 \end{cases}
Where a = -0.5 (standard choice) or a = -0.75 (sharper).
Interpolation:
z(u, v) = \sum_{m=0}^{3} \sum_{n=0}^{3} z_{m-1, n-1} \cdot w(m - u) \cdot w(n - v)
Properties: - C^1 continuous (smooth derivatives) - Preserves local shape better than bilinear - Can overshoot (produce values < min or > max of neighbors)
Problem: Resample this 2×2 elevation grid to find elevation at (1.3, 0.7) using bilinear interpolation.
Input grid (1m cells):
x=0 x=1 x=2
y=0 100 110
y=1 105 120Grid origin at (0, 0), cell size = 1m.
Query point: (x_q, y_q) = (1.3, 0.7)
Step 1: Find containing cell
i_q = \lfloor 1.3 / 1 \rfloor = 1 j_q = \lfloor 0.7 / 1 \rfloor = 0
Cell corners at: (1,0), (2,0), (1,1), (2,1)
Step 2: Fractional position
u = 1.3 - 1 = 0.3 v = 0.7 - 0 = 0.7
Step 3: Get corner values
Actually, our grid only goes to x=1. Let me recalculate with correct bounds.
Corrected: Query (1.3, 0.7) is outside grid. Let’s use (0.3, 0.7) instead.
Step 1 (corrected):
i_q = \lfloor 0.3 / 1 \rfloor = 0 j_q = \lfloor 0.7 / 1 \rfloor = 0
Step 2:
u = 0.3 - 0 = 0.3 v = 0.7 - 0 = 0.7
Step 3: Corner values
Step 4: Bilinear formula
z = (1-0.3)(1-0.7) \cdot 100 + (0.3)(1-0.7) \cdot 110 \quad + (1-0.3)(0.7) \cdot 105 + (0.3)(0.7) \cdot 120
= 0.7 \times 0.3 \times 100 + 0.3 \times 0.3 \times 110 \quad + 0.7 \times 0.7 \times 105 + 0.3 \times 0.7 \times 120
= 21 + 9.9 + 51.45 + 25.2 = 107.55 \text{ m}
Verification: - Nearest to (0,1) at 105m (weight 0.49) - Also weighted toward 110 and 120 - Result 107.55 is reasonable weighted average
Below is an interactive resampling demonstration.
<label>
Method:
<select id="interp-method">
<option value="nearest">Nearest Neighbor</option>
<option value="bilinear" selected>Bilinear</option>
<option value="bicubic">Bicubic</option>
</select>
</label>
<label>
Output resolution:
<input type="range" id="resolution-slider" min="5" max="30" step="5" value="10">
<span id="resolution-value">10</span>×
</label>
<label>
Show input grid:
<input type="checkbox" id="show-input" checked>
</label>
<div class="resample-info">
<p><strong>Input cells:</strong> 5×5 = 25</p>
<p><strong>Output cells:</strong> <span id="output-cells">--</span></p>
<p><strong>Sampling ratio:</strong> <span id="sample-ratio">--</span>×</p>
</div><canvas id="resample-canvas" width="600" height="300" style="border: 1px solid #ddd;"></canvas>Try this: - Nearest neighbor: Blocky, stair-step appearance (preserves exact values) - Bilinear: Smooth, continuous (but blurs the peak) - Bicubic: Very smooth (may overshoot slightly) - Increase resolution: More output pixels reveal interpolation differences - Color scale: Blue = low elevation, Red = high elevation - Notice: Nearest neighbor preserves sharp gradients, bilinear smooths them
Key insight: Method choice depends on data type—categorical (land cover) needs nearest neighbor, continuous (elevation) benefits from bilinear/bicubic.
Nearest Neighbor: - Categorical data: Land cover, soil types, zoning codes - Preserving exact values: Don’t want interpolated values (e.g., “class 2.3”) - Speed critical: Fastest method - Example: Reproject land cover map—don’t create new land cover classes
Bilinear: - Continuous data: Elevation, temperature, NDVI - General purpose: Good balance of speed and quality - Moderate smoothing acceptable: Slight blur is OK - Example: Resize satellite imagery for display
Bicubic: - High-quality visualization: Print maps, publications - Preserving gradients: Slope/aspect calculation from DEM - Smooth interpolation needed: Climate surfaces, bathymetry - Example: Upsample DEM for detailed terrain visualization
Problem: Downsampling without filtering causes aliasing (high-frequency artifacts).
Example: 100m DEM → 1000m DEM by sampling every 10th pixel - Misses peaks/valleys between sampled points - Creates artificial patterns
Solution: Average all input pixels that fall within output pixel (anti-aliasing filter).
For each output pixel:
Find all input pixels within its extent
z_output = mean(z_input)This is area-weighted averaging, not interpolation.
Changing coordinate systems requires resampling:
1. Create output grid in new projection
2. For each output pixel:
a. Transform pixel center to input projection
b. Interpolate value from input grid
c. Assign to output pixelChallenge: Output pixel may correspond to irregular shape in input space.
Common practice: Use bilinear for most data, nearest neighbor for categorical.
Bicubic can create values outside input range.
Example: - Input range: [0, 255] (8-bit image) - Bicubic output: Values like -5 or 270
Problem: Invalid for physical quantities (negative reflectance, temperatures > known max)
Solution: Clamp output to valid range: value = max(min_valid, min(max_valid, value))
Near grid boundaries: Insufficient neighbors for bilinear/bicubic.
Options: 1. Fallback to nearest neighbor at edges 2. Pad with replicated edge values 3. Assume zeros beyond boundary 4. Reflect boundary (mirror)
Choice depends on data semantics.
Many interpolation methods assume: - Rectangular grid - Uniform spacing - Aligned with axes
Real data may be: - Irregular point clouds (LiDAR) - Rotated grids - Variable resolution
Solution: Use scattered data interpolation (IDW, kriging) or triangulation.
Bilinear/bicubic blur sharp features.
Problem for: - Coastlines (land/water boundary) - Faults (abrupt elevation change) - Roads (linear features)
Solution: - Use nearest neighbor for boundaries - Edge-preserving filters - Or accept smoothing as necessary artifact
For irregular point data:
z(x, y) = \frac{\sum_i w_i z_i}{\sum_i w_i}
Where:
w_i = \frac{1}{d_i^p}
d_i = distance from (x,y) to point i
p = power parameter (typically 2)
Nearby points have more influence (high weight).
Geostatistical interpolation using spatial covariance.
Steps: 1. Fit variogram (model of spatial correlation vs. distance) 2. Solve system of equations for weights 3. Compute weighted sum (optimal in least-squares sense)
Advantage: Provides uncertainty estimates.
Disadvantage: Computationally expensive, requires expertise.
Minimizes bending energy (creates smoothest surface).
Used for: - Climate data interpolation - Surface fitting with constraints - Morphing between shapes
Disadvantage: Global method (slow for large datasets).
Points: (x_0, y_0) and (x_1, y_1)
Interpolate at x:
y = y_0 + \frac{x - x_0}{x_1 - x_0}(y_1 - y_0)
Equivalent form:
y = (1 - t) y_0 + t y_1
Where t = \frac{x - x_0}{x_1 - x_0} is normalized position.
Properties: - Linear in x - Exact at endpoints: y(x_0) = y_0, y(x_1) = y_1 - t = 0 → y = y_0; t = 1 → y = y_1
Apply linear interpolation twice:
Result: Bilinear = product of two linear interpolations.
Not the same as 2D linear (plane through 3 points)—bilinear is a hyperbolic paraboloid (saddle shape).