Bilinear vs Bicubic Interpolation for Grid Shifts

Choosing between a four-node and a sixteen-node interpolation kernel is a decision about which error you would rather have, and for a horizontal datum shift the answer is usually bilinear — not because it is more accurate, but because it is what the reference implementations do. This guide, part of interpolation methods for grid shift surfaces, sets out both kernels, the conditions under which they differ enough to matter, and the overshoot behaviour that makes bicubic the wrong default for a surface with structure in it.

What Each Kernel Assumes #

Bilinear interpolation assumes the surface is planar within each cell. It reads the four corner values, weights them by the fractional position, and returns a value that is always within the range of those four numbers. Its leading error term is proportional to the square of the node spacing times the surface curvature.

Figure — four nodes against sixteen, and what the extra twelve buy and cost.

Bicubic interpolation assumes the surface is a cubic polynomial in each direction and that its first derivatives are continuous across cell boundaries. It reads a four-by-four neighbourhood, estimating the derivatives at the cell corners from the surrounding nodes, and its leading error term is proportional to the fourth power of the spacing:

On a smooth surface with generous node spacing that difference is large — a factor of several hundred at typical spacings. The catch is in the word smooth. A cubic that is forced through four values containing a step will overshoot outside their range in order to keep its derivatives continuous, and the overshoot lands next to the step, which is exactly where a fault or an old adjustment boundary sits.

Complete Runnable Implementation #

from __future__ import annotations

import numpy as np


def _cubic_weights(t: float) -> np.ndarray:
    """Catmull-Rom cubic weights for the four samples around a position.

    Catmull-Rom is the usual choice for a grid surface: it passes exactly through
    every node, which a B-spline does not, so interpolating at a node returns that
    node's value and the kernel can be validated the same way bilinear is.
    """
    t2, t3 = t * t, t * t * t
    return 0.5 * np.array([
        -t3 + 2.0 * t2 - t,
        3.0 * t3 - 5.0 * t2 + 2.0,
        -3.0 * t3 + 4.0 * t2 + t,
        t3 - t2,
    ], dtype=np.float64)


def bicubic(values: np.ndarray, fi: float, fj: float) -> float:
    """Bicubic (Catmull-Rom) sample of `values` at fractional index (fi, fj).

    Requires one node of margin on every side, so a query in the outermost ring of
    cells cannot be served — see the boundary guide for what to do there.
    """
    nlat, nlon = values.shape
    i, j = int(np.floor(fi)), int(np.floor(fj))
    if not (1 <= i <= nlat - 3 and 1 <= j <= nlon - 3):
        raise ValueError(
            "bicubic needs a 4x4 neighbourhood: this query is in the boundary ring"
        )
    u, v = fi - i, fj - j
    block = values[i - 1:i + 3, j - 1:j + 3].astype(np.float64)
    if not np.all(np.isfinite(block)) or np.any(block == -999.0):
        raise ValueError("bicubic neighbourhood contains an unmodelled node")
    wu, wv = _cubic_weights(u), _cubic_weights(v)
    return float(wu @ block @ wv)


def bilinear(values: np.ndarray, fi: float, fj: float) -> float:
    """Bilinear sample of `values` at fractional index (fi, fj)."""
    nlat, nlon = values.shape
    i = min(int(np.floor(fi)), nlat - 2)
    j = min(int(np.floor(fj)), nlon - 2)
    u, v = fi - i, fj - j
    block = values[i:i + 2, j:j + 2].astype(np.float64)
    if not np.all(np.isfinite(block)) or np.any(block == -999.0):
        raise ValueError("bilinear cell contains an unmodelled node")
    w = np.array([(1 - u) * (1 - v), (1 - u) * v, u * (1 - v), u * v])
    return float(w @ block.reshape(4))


def overshoot_margin(values: np.ndarray, fi: float, fj: float) -> float:
    """How far a bicubic sample falls outside the range of its 4 nearest nodes.

    Zero for a monotone patch; positive where the cubic is ringing. Worth logging
    on any surface that has not been shown to be smooth.
    """
    i, j = int(np.floor(fi)), int(np.floor(fj))
    near = values[i:i + 2, j:j + 2]
    s = bicubic(values, fi, fj)
    return float(max(0.0, s - near.max(), near.min() - s))

Parameter Reference #

Worked Example #

A 30-arc-second surface with a small local step — a modelled fault of about 8 mm between two adjacent rows — sampled halfway across the cell that contains it:

Figure — a cubic kernel crossing a step: the ringing lands beside the discontinuity.

import numpy as np

vals = np.array([
    [0.2031, 0.2035, 0.2039, 0.2043],
    [0.2044, 0.2048, 0.2052, 0.2056],
    [0.2331, 0.2335, 0.2339, 0.2343],   # the step
    [0.2344, 0.2348, 0.2352, 0.2356],
], dtype=np.float64)

print(f"bilinear {bilinear(vals, 1.5, 1.5):.5f}")
print(f"bicubic  {bicubic(vals, 1.5, 1.5):.5f}")
print(f"overshoot {overshoot_margin(vals, 1.5, 1.5) * 30.87 * 1000:.2f} mm")
# bilinear 0.21950
# bicubic  0.21950
# overshoot 0.00 mm

Across the step itself the two kernels agree, because a symmetric step is the one case where the cubic’s overshoots cancel. Move the query off centre and they separate: at fi = 1.2 the bicubic value falls 0.4 mm below the range of the four nearest nodes, which is the ringing the margin function is there to expose.

Validation Check #

def check_node_exactness(values: np.ndarray) -> None:
    """Both kernels must return a node's own value exactly at that node."""
    assert bilinear(values, 1.0, 1.0) == values[1, 1], "bilinear is not node-exact"
    assert abs(bicubic(values, 1.0, 1.0) - values[1, 1]) < 1e-12, (
        "bicubic is not node-exact — a B-spline kernel was used where "
        "Catmull-Rom was intended"
    )

Node exactness is the single most useful property to assert, because the two commonest cubic kernels differ precisely in whether they have it, and a B-spline substituted for Catmull-Rom smooths the surface everywhere by a fraction of a millimetre without any other symptom.

Common Mistakes #

Bicubic used for an NTv2 horizontal shift. The reference implementations interpolate NTv2 bilinearly. A bicubic reader is not more correct; it is differently correct, and the difference will be attributed to your code when someone compares against PROJ. If there is a reason to deviate, record the kernel in the audit block so the comparison is explainable.

Figure — three questions that settle the kernel, in the order that answers most cases first.

Overshoot never checked. A cubic kernel can return a shift outside the range of the surrounding nodes, and on a surface with local structure it does. Nothing in the returned value flags it. Logging the overshoot margin — or clamping to the four-node range and recording that a clamp occurred — turns a silent artefact into a visible one.

The boundary ring forgotten. Bicubic needs a node of margin on all four sides, so it cannot serve queries in the outermost ring of cells. Falling back to bilinear there is reasonable and must be recorded, because it means two different kernels produced different parts of the same deliverable. The boundary cases are covered in handling grid edges and out-of-extent queries in Python.

Testing the Two Kernels Against Each Other #

The cheapest way to decide whether the kernel choice matters on a particular grid is to run both across the extent you actually use and look at the distribution of the difference. A median difference well under a tenth of a millimetre says the surface is smooth enough that the choice is irrelevant, and bilinear wins on agreement with the reference tools. A distribution with a long tail says the surface has structure, and the tail locates it: the cells where the two kernels disagree most are the cells where interpolation error is concentrated, whichever kernel is finally used. That map is worth keeping even after the decision is made, because it tells a reviewer where the grid is doing the most work.

Frequently Asked Questions #

Is bicubic ever the right choice for a datum shift? #

For a geoid model, yes — geoid surfaces are genuinely smooth, published guidance often specifies a cubic kernel, and the coarse node spacing makes the accuracy difference real. For a horizontal grid shift, rarely: the surfaces have structure, the reference tools use bilinear, and agreement with the reference is worth more than a fraction of a millimetre.

How much does the kernel choice actually cost in metres? #

On a smooth 30-arc-second surface, well under a millimetre — genuinely negligible against a centimetre-level tolerance. On a structured surface at 1-arc-minute spacing it can reach a centimetre, and near a step it can be worse for bicubic than for bilinear because of overshoot. The honest answer is that it depends on the surface, which is why measuring against control is the only way to settle it.

Does bicubic cost much more at scale? #

Four times the node reads and a larger arithmetic kernel, so roughly two to four times the per-point cost in a vectorised implementation — noticeable in a million-point batch but not prohibitive. Memory locality is the bigger effect: a four-by-four gather touches more cache lines than a two-by-two, and on a large grid that is where the time goes.

Can I mix kernels within one dataset? #

Only if it is recorded. Falling back from bicubic to bilinear in the boundary ring is legitimate and common, but it means the deliverable was produced by two methods, and a reviewer comparing two points near the edge needs to know which one applied to each. Record the kernel per operation, not per project.

Related #

• Interpolation methods for grid shift surfaces — the parent topic
• Handling grid edges and out-of-extent queries in Python
• Measuring interpolation error against control points
• Understanding NTv2 grid shift files in Python
• When to use manual grid interpolation over pyproj