Interpolation, density and hexagon grids

Sixty noise meters give sixty readings, and a map needs a value at every one of the town’s 4,800 hectare cells. Estimating values between readings is interpolation, and each estimate depends on how much the nearest reading counts: at one place, the same four readings give 60.90 dB or 67.74 dB depending on a single setting.

This lesson builds such estimates with inverse distance weighting, tests them against readings they did not use, and sets out the methods beside it: krigingkriging: Interpolation that takes its weights from a model of the variogram, so that the readings' own spatial structure decides how much each one counts, and that gives a standard error for every estimate. for values with an error, kernel densitykernel density estimation: Turning points that are events into a smooth surface of events per unit area, by spreading each point over a kernel of a chosen radius, the bandwidth, and adding the kernels. for events, and grids of squares or hexagons for counts.

Inverse distance weighting trusts the nearest readings most

Inverse distance weightinginverse distance weighting (IDW): An interpolation method that estimates the value at a place as a weighted mean of known samples, each weighted by one over its distance raised to a power, so that nearer samples count for more. More in “What GIS is and the questions it answers” estimates the value at a place as a weighted mean of the readings, each weighted by one over its distance to a power p:

In words: the nearer a reading, the more it counts, and the larger p, the faster that influence fades. At a reading’s own place the estimate is the reading itself. Step through one estimate at T, 64 m from a 70 dB reading and farther from three quieter ones:

The program repeats it for four powers. Press Run, then move T or change the readings:

"""Inverse distance weighting at one place, T, from four noise readings: the distances, the weights for powers 0 to 3,
and the estimate each power gives. At a reading's own place IDW returns that reading."""
import math

READINGS = {"A": (0, 0, 70), "B": (200, 0, 52), "C": (0, 200, 55), "D": (200, 200, 50)}  # x m, y m, dB
T = (50, 40)


def idw(x, y, power):
    weights, total = {}, 0.0
    for name, (sx, sy, value) in READINGS.items():
        d = math.dist((x, y), (sx, sy))
        if d == 0:
            return value, {name: 1.0}  # a reading's own place: exact
        weights[name] = 1 / d**power
        total += weights[name]
    shares = {k: w / total for k, w in weights.items()}
    return sum(shares[k] * READINGS[k][2] for k in shares), shares


for name, (sx, sy, value) in READINGS.items():
    print(f"{name}: {value} dB, {math.dist(T, (sx, sy)):.1f} m from T")
for power in (0, 1, 2, 3):
    estimate, shares = idw(*T, power)
    print(f"power {power}: shares " + ", ".join(f"{k} {s:.3f}" for k, s in shares.items()) + f"; estimate {estimate:.2f} dB")
print(f"at A itself: {idw(0, 0, 2)[0]} dB")

A: 70 dB, 64.0 m from T
B: 52 dB, 155.2 m from T
C: 55 dB, 167.6 m from T
D: 50 dB, 219.3 m from T
power 0: shares A 0.250, B 0.250, C 0.250, D 0.250; estimate 56.75 dB
power 1: shares A 0.479, B 0.198, C 0.183, D 0.140; estimate 60.90 dB
power 2: shares A 0.714, B 0.121, C 0.104, D 0.061; estimate 65.04 dB
power 3: shares A 0.869, B 0.061, C 0.048, D 0.022; estimate 67.74 dB
at A itself: 70 dB

With every weight equal, power 0, the estimate is the plain mean, 56.75 dB; power 1 gives 60.90 dB, power 2 65.04 dB and power 3 67.74 dB, closer and closer to A. QGIS’s IDW interpolation takes the power as its distance coefficient P, 2 by default (QGIS), and its introduction to GIS warns that the surface’s maximum and minimum “can only occur at sample data points” (QGIS).

Hidden readings test a surface, and other methods answer other questions

An interpolation cannot be judged by how smooth it looks. Hide some readings, predict them from the others, and measure the errors: that is cross-validationcross-validation: Testing a model or an interpolation on data it was not built from: some samples are held out, predicted from the rest, and the errors at the held-out samples are measured., and leave-one-out does it for each reading in turn. Three other methods answer different questions about points:

• Kriging fits a model, such as an exponential or a Gaussian curve, to the variogramvariogram: A description of how different two measurements are, on average, as a function of the distance between them: half the mean squared difference of the pairs at each distance., half the mean squared difference between readings at each distance apart, and takes its weights from it; it also gives a standard error for every estimate (Matheron, Economic Geology).
• Kernel density spreads each event, such as a crash, over a disc of a chosen radius and adds the discs, giving events per km². The radius acts as a smoothing parameter: a large one gives a smooth, high-bias surface (scikit-learn). The kernel sets how the spread falls off across the disc: QGIS’s heatmap offers the quartic kernel, its default, and triangular, uniform, triweight and Epanechnikov ones (QGIS).
• Hexagon grids count or average points in equal hexagons. Each hexagon has six neighbours at one distance, while squares have neighbours at two distances (H3). H3H3: A global grid that divides the Earth into hexagonal cells, with 12 pentagons, at 16 resolutions, each finer one with about seven times as many cells, and gives every cell an id. fixes such cells for the whole planet, each with an id, at 16 resolutions, each with about seven times as many cells as the next coarser one and exactly 12 pentagons among its hexagons (H3).

Which fits depends on what the points carry and what you need. Answer three questions about your own points:

In short: interpolate measured values, with kriging when every estimate needs an error; turn events into a density; and count or average on a grid when the numbers must be compared or joined.

Choosing IDW's power and radius by cross-validation, measured Pro: opens with a Pro, Premium or Ultimate pass

A variogram of the readings, and what kriging adds Pro: opens with a Pro, Premium or Ultimate pass

Density and grids: radii, kernel shapes, squares, hexagons and H3 Pro: opens with a Pro, Premium or Ultimate pass

Practise: IDW with a search radius

Write IDW with a search radius, then its leave-one-out check. The sample tests use the lesson’s four readings, a reading’s own place, a radius that leaves a place out of reach, and a line of readings scored at two powers:

import math


def idw(samples, x, y, power=2, radius=None):
    """Inverse distance weighting at (x, y) (see the prompt)."""
    return None


def loo_rmse(samples, power=2, radius=None):
    """The leave-one-out RMSE of IDW on the samples (see the prompt)."""
    return None

from interpolate import idw, loo_rmse

FOUR = [(0, 0, 70), (200, 0, 52), (0, 200, 55), (200, 200, 50)]  # the lesson's four readings, dB
LINE = [(0, 0, 10), (100, 0, 20), (200, 0, 30), (300, 0, 40)]


def test_the_lessons_estimate_at_t():
    """gives 65.04 dB at T with power 2 and 60.90 dB with power 1, as the lesson's program printed"""
    assert round(idw(FOUR, 50, 40), 2) == 65.04
    assert round(idw(FOUR, 50, 40, power=1), 2) == 60.9


def test_exact_at_a_sample():
    """returns a reading's own value at its place, without dividing by zero"""
    assert idw(FOUR, 0, 0) == 70


def test_a_search_radius():
    """uses only the samples within the radius, and gives None when none is in reach"""
    assert round(idw(FOUR, 50, 40, radius=160), 2) == 67.38
    assert idw(FOUR, 1_000, 1_000, radius=100) is None


def test_leave_one_out_on_a_line():
    """measures the error of each sample predicted from the others, for two powers"""
    assert loo_rmse(LINE) == 9.653
    assert loo_rmse(LINE, power=1) == 11.912


def test_leave_one_out_with_a_radius():
    """skips the samples that get no estimate, and gives None when none gets one"""
    assert loo_rmse(LINE, radius=150) == 7.071
    assert loo_rmse([(0, 0, 10), (100, 0, 20), (5_000, 0, 30)], radius=500) == 10.0
    assert loo_rmse(FOUR, radius=150) is None

Then test yourself on IDW and on checking a surface, one question at a time, the idea you know least first:

Interview questions

Warm-up (fresher): What does the power in inverse distance weighting control?

It controls how fast a reading’s influence falls with distance: each weight is one over the distance to the power p. With p = 0 every reading counts the same and the estimate is their plain mean; as p grows, the nearest reading dominates. At the lesson’s place T, four readings gave 56.75 dB at power 0, 60.90 dB at power 1, 65.04 dB at power 2 and 67.74 dB at power 3, closer and closer to the 70 dB reading 64 m away. No power is right in general: choose it by hiding readings and measuring the errors, and report it.

Three longer questions take these methods to real decisions about sensors, crashes and hexagons; their first steps and follow-ups are set out in the module’s bank:

Key takeaways

• Inverse distance weighting estimates a place’s value as a weighted mean of readings, weighted by one over distance to a power; at T the power moved the estimate from 56.75 to 67.74 dB.
• An interpolation is judged by the errors at readings it did not use, never by how smooth it looks.
• Kriging adds a standard error from a variogram, kernel density turns events into events per km² at a chosen radius, and hexagon grids give every cell six neighbours at one distance.