🧭 Concept Spatial statistics Maps & visualisation Spatial Statistics · step 2 of 13 Updated

Kernel Density Explained: What a Heatmap Actually Shows

Problem statement

Someone hands you a heatmap of incidents and asks where the problem is. The map has three glowing red areas, so you point at them.

Then you make the same map with a different bandwidth and there are twelve. Or one.

 bandwidth   distinct areas above half the peak   share of map
      0.7 m                                  12          0.4%
      1.5 m                                   3          1.1%
      5.0 m                                   3          4.9%
     14.0 m                                   2         23.3%
     25.0 m                                   1         59.2%

Same 600 points. Same code. The number of hotspots on the map is a parameter you chose, not a property of the data — and on most heatmaps nobody records what it was.

Kernel density estimation is the right tool for turning points into a continuous surface. It just has one knob that decides the answer, and understanding that knob is the whole subject.

Quick answer

KDE replaces every point with a small hill and adds the hills up:

from sklearn.neighbors import KernelDensity
import numpy as np

xy = np.column_stack([points.geometry.x, points.geometry.y])
kde = KernelDensity(bandwidth=500, kernel="gaussian").fit(xy)   # bandwidth in map units

gx, gy = np.mgrid[minx:maxx:250j, miny:maxy:250j]
density = np.exp(kde.score_samples(np.column_stack([gx.ravel(), gy.ravel()])))
surface = density.reshape(gx.shape)

Three parameters, in decreasing order of importance:

Parameter What it controls How much it matters
bandwidth how wide each hill is decides the answer
cell size how finely the surface is sampled cosmetic, if smaller than the bandwidth
kernel shape the profile of each hill almost irrelevant

The bandwidth must be in the same units as your coordinates, which means a projected CRS. On EPSG:4326 a bandwidth of 500 means 500 degrees.

Individual points each replaced by a small hill, the hills summed into a continuous density surface.
One hill per point, added together. The width of the hill is the only real choice.

Step-by-step solution

1. Understand what the values mean

A KDE surface is a density: events per unit area, not a count. Integrate it over the whole plane and it comes to 1 (or to n, if you scale by the number of points).

That has two practical consequences people trip over:

  • A cell value is not a count. "0.00088" is not 0.00088 incidents. It is a density, and its units are per square metre.
  • Changing the bandwidth changes every value. The peak density above went from 0.00298 at bandwidth 1.5 to 0.00013 at bandwidth 25 — a factor of 23, on identical data. Absolute KDE values are only comparable between maps made with the same bandwidth.

Which is why almost every heatmap is displayed on a relative scale, and why comparing two heatmaps by colour is meaningless unless both were made the same way.

2. Choose the bandwidth from the process, not from a rule

Automatic rules exist. scipy's default is Scott's rule:

from scipy.stats import gaussian_kde
kde = gaussian_kde(xy.T)
print(f"factor {kde.factor:.4f} -> {kde.factor * xy.std():.2f} map units")
factor 0.3443 -> 7.86 map units

Scott's and Silverman's rules both optimise for a single smooth bump. Real spatial data is clustered and bounded, so they systematically over-smooth it — they are a starting point, not an answer.

The defensible choice comes from the phenomenon:

  • Burglary is a street-level process. A bandwidth of 100–200 m matches how offenders move.
  • Air pollution disperses over kilometres. A 200 m bandwidth would be noise.
  • Retail catchments are walking distance. 400–800 m.

If you cannot name the distance the process operates over, you cannot defend the bandwidth — and then you must show several.

3. Recognise under- and over-smoothing

The measured behaviour of the surface as the bandwidth grows:

 bandwidth   blobs above half the peak   area above half the peak
      0.7                          12                       0.4%
      1.5                           3                       1.1%
      3.0                           3                       2.5%
      5.0                           3                       4.9%
      8.0                           3                       7.5%
     14.0                           2                      23.3%
     25.0                           1                      59.2%

The data really has three clusters. Below 1.5 the surface fragments — it is drawing individual points, not structure. Above 8 it merges — two genuine clusters become one. Between 1.5 and 8 the answer is stable, which is exactly the sign you want.

A bandwidth range over which the answer does not change is the best evidence you can offer. Report it.

4. Keep the cell size well below the bandwidth

The cell size only controls how finely you sample a surface the bandwidth has already smoothed. As a rule, make cells a quarter of the bandwidth or smaller and stop thinking about it.

Cells larger than the bandwidth reintroduce exactly the aggregation artefacts KDE was supposed to avoid — you are back to the modifiable areal unit problem, with extra steps.

5. Ignore the kernel shape

Gaussian, Epanechnikov, quartic, triangular — for a given bandwidth the differences are visually negligible and statistically minor. Gaussian has infinite support (every point influences every cell, which is slower); Epanechnikov and quartic are compact and faster.

Choose on performance. Never present a kernel choice as an analytical decision.

Three bandwidth regimes: fragmented at 0.7, stable at 1.5 to 8 showing three clusters, and merged at 25.
The stable middle band is the finding. Its width is what you should report.

Code examples

Example 1 — a KDE surface with the parameters made explicit 

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.transform import from_origin
from sklearn.neighbors import KernelDensity


def kde_surface(points, bandwidth, *, cell=None, kernel="gaussian", pad=None):
    """Return a density surface plus the affine transform that georeferences it."""
    if points.crs is None or points.crs.is_geographic:
        raise ValueError("KDE needs a projected CRS — the bandwidth is in map units")

    cell = cell or bandwidth / 4          # sample well below the smoothing scale
    pad = pad if pad is not None else bandwidth * 3   # room for the tails

    minx, miny, maxx, maxy = points.total_bounds
    minx, miny, maxx, maxy = minx - pad, miny - pad, maxx + pad, maxy + pad
    ncols = int(np.ceil((maxx - minx) / cell))
    nrows = int(np.ceil((maxy - miny) / cell))

    xs = minx + (np.arange(ncols) + 0.5) * cell
    ys = maxy - (np.arange(nrows) + 0.5) * cell         # north-up, row 0 at the top
    gx, gy = np.meshgrid(xs, ys)

    xy = np.column_stack([points.geometry.x, points.geometry.y])
    model = KernelDensity(bandwidth=bandwidth, kernel=kernel).fit(xy)
    log_density = model.score_samples(np.column_stack([gx.ravel(), gy.ravel()]))

    surface = np.exp(log_density).reshape(nrows, ncols) * len(points)   # -> events / m²
    transform = from_origin(minx, maxy, cell, cell)
    print(f"{nrows}x{ncols} cells of {cell:.0f} m, bandwidth {bandwidth} m, "
          f"peak {surface.max():.3e} events/m²")
    return surface, transform


surface, transform = kde_surface(incidents, bandwidth=300)
440x512 cells of 75 m, bandwidth 300 m, peak 4.812e-06 events/m²

Two details that matter. Multiplying by len(points) converts a probability density into events per square metre, which is a quantity a reader can reason about. And the padding stops the surface being clipped at the edge of the data, where the tails of the kernels still carry weight.

Example 2 — testing whether your conclusion survives the bandwidth 

from scipy.ndimage import label


def bandwidth_sweep(points, bandwidths, threshold=0.5):
    """Count distinct high-density areas at each bandwidth."""
    rows = []
    for bw in bandwidths:
        surface, _ = kde_surface(points, bw, cell=bw / 4)
        hot = surface > threshold * surface.max()
        rows.append({
            "bandwidth": bw,
            "blobs": label(hot)[1],
            "area_share": f"{hot.mean():.1%}",
            "peak": f"{surface.max():.2e}",
        })
    return pd.DataFrame(rows)


print(bandwidth_sweep(incidents, [0.7, 1.5, 3, 5, 8, 14, 25]).to_string(index=False))
 bandwidth  blobs area_share     peak
       0.7     12       0.4% 2.98e-03
       1.5      3       1.1% 1.83e-03
       3.0      3       2.5% 8.80e-04
       5.0      3       4.9% 5.11e-04
       8.0      3       7.5% 3.60e-04
      14.0      2      23.3% 1.44e-04
      25.0      1      59.2% 1.30e-04

"Three hotspots, stable from 1.5 to 8" is a finding. "Three hotspots" is a bandwidth.

Notice the peak column falling by a factor of 23 across the sweep. That is why two heatmaps made with different bandwidths cannot be compared by colour, even when both use the same colour ramp.

Example 3 — density of a rate, not of raw counts

Raw density maps mostly show where people are. If you want the surface to mean something other than population, divide two densities:

def relative_risk(cases, controls, bandwidth, **kwargs):
    """Density of cases divided by density of the population at risk."""
    case_surface, transform = kde_surface(cases, bandwidth, **kwargs)
    control_surface, _ = kde_surface(controls, bandwidth, **kwargs)

    # suppress the ratio where there is essentially no population — it explodes there
    floor = np.percentile(control_surface, 20)
    ratio = np.where(control_surface > floor, case_surface / control_surface, np.nan)

    print(f"relative risk: median {np.nanmedian(ratio):.2f}, "
          f"max {np.nanmax(ratio):.2f}, {np.isnan(ratio).mean():.0%} suppressed")
    return ratio, transform


risk, transform = relative_risk(incidents, households, bandwidth=300)
relative risk: median 1.03, max 3.41, 20% suppressed

Both surfaces must use the same bandwidth and the same grid, or the ratio is comparing differently-smoothed things. The floor is not optional: dividing by a near-zero density produces enormous values in exactly the empty areas where they mean nothing, which is the single commonest way a relative-risk map goes wrong.

Explanation

Why the bandwidth changes the peak value

Each point contributes a fixed total mass of 1, spread over the kernel. Widen the kernel and the same mass covers more area, so the height everywhere falls.

That is why the peak fell from 2.98e-03 to 1.30e-04 across the sweep: the total under the surface is constant, and the bandwidth decides whether it is concentrated or spread. Absolute density values are meaningless without the bandwidth attached.

Why the edge of the study area is biased low

A point on the boundary has half its kernel outside the study area. That mass is not redistributed — it is simply gone from the visible surface. The result is a density that falls off near the edge, whether or not the process does.

This is why coastal cities always appear to have a "quiet" waterfront on an uncorrected heatmap. Options, in increasing order of effort: pad the extent and crop for display, use data from beyond the boundary, or apply an explicit edge correction that reweights boundary kernels.

Why KDE is not a hotspot test

A KDE surface always has a maximum. Run it on 600 uniformly random points and you will get peaks, valleys and something that looks exactly like structure — because random points are not evenly spaced.

KDE describes; it does not test. To claim a concentration is more than chance you need a statistical test — Getis-Ord Gi* on aggregated units, or a point-pattern method like Ripley's K.

Present a heatmap as a description of where events are, never as evidence that the pattern is significant.

A point near the study-area boundary with half its kernel falling outside, producing an artificially low density at the edge.
The mass outside the boundary is lost, not redistributed. Every uncorrected KDE fades at its edges.

Why a projected CRS is non-negotiable

The bandwidth is a distance in coordinate units. In EPSG:4326 those units are degrees, and a degree of longitude is 111 km at the equator and 71 km in Manchester.

So a KDE on unprojected coordinates uses a kernel that is circular in degree-space and increasingly elliptical on the ground — over-smoothing east-west relative to north-south, by a factor that grows with latitude. The map still renders. It is just wrong in a way that varies across itself. See projected vs geographic CRS.

Edge cases or notes

  • KDE on unprojected coordinates is always wrong. Reproject first; the distortion is latitude-dependent and invisible.
  • Duplicate points at the same coordinate (geocoded to a centroid, for instance) create a spike that is an artefact of geocoding, not of the process. See geocoding returns wrong coordinates.
  • Weighted KDE is available via sample_weight in scikit-learn — use it when each record represents several events.
  • Gaussian kernels have infinite support, so every point contributes to every cell. On large datasets a compact kernel is much faster with no visible difference.
  • Do not compare heatmaps from different bandwidths by colour. The scales differ by orders of magnitude.
  • Adaptive bandwidth (wider where points are sparse) is available in some packages and helps with strongly varying density, at the cost of a surface whose smoothing varies across itself.
  • The grid extent should exceed the data extent by two or three bandwidths, or you clip the tails and worsen the edge effect.

FAQ

What does a value on a KDE surface mean?

A density — events per unit area, not a count. Multiply by the number of points to get events per square metre, and always state the bandwidth alongside.

How do I choose the bandwidth?

From the distance the process operates over: street-level crime, a few hundred metres; dispersal, kilometres. Automatic rules over-smooth clustered data. If you cannot justify one value, show several.

Does the kernel shape matter?

Barely. Gaussian, quartic and Epanechnikov produce visually similar surfaces at the same bandwidth. Choose on speed.

Why do my hotspots disappear when I change the bandwidth?

Because the number of hotspots is a function of the bandwidth. Sweep a range and report the band over which the answer is stable.

Why does the density drop at the edge of my study area?

Edge effect — points near the boundary have part of their kernel outside it, and that mass is lost. Pad the extent, include data from beyond the boundary, or apply an edge correction.

Can I compare two heatmaps?

Only if they were made with the same bandwidth, cell size and extent. Peak density varies by orders of magnitude with bandwidth alone.

Is a KDE hotspot statistically significant?

No. KDE always produces peaks, including on uniformly random points. For significance use Getis-Ord Gi* or a point-pattern test.