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 a handful of glowing areas, so you point at them.

Then you make the same map with a different bandwidth and there are thirty. Or one. On Leeds' 9,229 police-recorded crimes for June 2026 (Example 2):

 bandwidth   distinct areas above a tenth of the peak   share of map
     150 m                                        30           0.8%
     300 m                                         6           1.9%
     600 m                                         5           7.2%
   1,000 m                                         2          11.7%
   2,500 m                                         1          14.7%

Same 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

points = incidents                                  # Leeds crimes, EPSG:27700 metres
minx, miny, maxx, maxy = points.total_bounds
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:

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.

What the values mean, bandwidth from the process, cell size

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. On the Leeds crimes the peak density fell from 1.58 × 10⁻³ per m² at a 150 m bandwidth to 8.23 × 10⁻⁵ at 2,500 m — a factor of 19, on identical data (Example 2). 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.2183 -> 1076.99 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 and street crime are street-level processes: a street or two, a few hundred metres at most.
• Air pollution disperses over kilometres, so a 200 m bandwidth would be mapping noise.
• Retail catchments are walking distance: several hundred metres.

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

Example 2 sweeps the bandwidth on the Leeds crimes and counts the separate areas above a tenth of the peak. Small bandwidths fragment the surface — it draws individual busy streets, not structure. Large ones merge everything into one blob around the centre. In between, look for a range where the count holds.

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.

Leeds crime surfaces and a bandwidth sensitivity test

The data used on this page

Every example runs on real data: the 488 Lower Layer Super Output Areas (LSOAs) of Leeds from the ONS, their Census 2021 population from Nomis, and the street-level crimes West Yorkshire Police published for June 2026 on data.police.uk. The surfaces are of the crime points, incidents; zones supplies the residents for the relative-risk example.

import io
import time

import geopandas as gpd
import numpy as np
import pandas as pd
import requests

ONS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
       "Lower_layer_Super_Output_Areas_December_2021_Boundaries_EW_BGC_V5/FeatureServer/0/query")
features, offset = [], 0
while True:                                   # the service pages at 1,000 features
    page = requests.get(ONS, params={"where": "LSOA21NM LIKE 'Leeds %'", "outFields": "LSOA21CD,LSOA21NM",
                                     "outSR": 27700, "f": "geojson", "resultOffset": offset}, timeout=120).json()
    features += page["features"]
    if not page.get("properties", {}).get("exceededTransferLimit"):
        break
    offset += len(page["features"])
zones = gpd.GeoDataFrame.from_features(features, crs=27700).set_index("LSOA21CD").sort_index()

population = pd.read_csv("https://www.nomisweb.co.uk/api/v01/dataset/NM_2021_1.data.csv?date=latest"
                         f"&geography={','.join(zones.index)}&c2021_restype_3=0&measures=20100"
                         "&select=geography_code,obs_value").set_index("GEOGRAPHY_CODE")["OBS_VALUE"]
zones["population"] = population

west, south, east, north = zones.to_crs(4326).total_bounds
crimes = []
for x0, x1 in zip(np.linspace(west, east, 5)[:-1], np.linspace(west, east, 5)[1:]):
    for y0, y1 in zip(np.linspace(south, north, 5)[:-1], np.linspace(south, north, 5)[1:]):   # ≤10,000 a request
        for attempt in range(5):                        # the API answers bursts with 429
            reply = requests.post("https://data.police.uk/api/crimes-street/all-crime", timeout=120,
                                  data={"poly": f"{y0},{x0}:{y0},{x1}:{y1},{x1}:{y1},{x0}", "date": "2026-06"})
            if reply.status_code != 429:
                break
            time.sleep(2 ** attempt)
        reply.raise_for_status()
        crimes += [{"id": c["id"], "category": c["category"], "lon": float(c["location"]["longitude"]),
                    "lat": float(c["location"]["latitude"])} for c in reply.json()]
crimes = pd.DataFrame(crimes).drop_duplicates("id")
incidents = gpd.GeoDataFrame(crimes, geometry=gpd.points_from_xy(crimes.lon, crimes.lat), crs=4326).to_crs(27700)
incidents = gpd.sjoin(incidents, zones[["geometry"]], predicate="within").rename(columns={"index_right": "lsoa", "LSOA21CD": "lsoa"})
zones["crimes"] = incidents["lsoa"].value_counts().reindex(zones.index, fill_value=0)
zones["rate"] = 1000 * zones["crimes"] / zones["population"]
print(f"{len(zones)} LSOAs, {zones['population'].sum():,} residents, {len(incidents):,} crimes; "
      f"rate median {zones['rate'].median():.1f}, max {zones['rate'].max():.1f} per 1,000")

488 LSOAs, 811,957 residents, 9,229 crimes; rate median 8.4, max 174.4 per 1,000

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, weights=None, bounds=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 if bounds is None else bounds   # share bounds= to align two surfaces
    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, sample_weight=weights)
    log_density = model.score_samples(np.column_stack([gx.ravel(), gy.ravel()]))

    total = len(points) if weights is None else float(np.sum(weights))
    surface = np.exp(log_density).reshape(nrows, ncols) * total        # -> 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)

358x434 cells of 75 m, bandwidth 300 m, peak 9.075e-04 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

import pandas as pd
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, [150, 300, 600, 1000, 1500, 2500], threshold=0.1).to_string(index=False))

691x843 cells of 38 m, bandwidth 150 m, peak 1.580e-03 events/m²
358x434 cells of 75 m, bandwidth 300 m, peak 9.075e-04 events/m²
191x229 cells of 150 m, bandwidth 600 m, peak 3.869e-04 events/m²
124x147 cells of 250 m, bandwidth 1000 m, peak 2.048e-04 events/m²
91x106 cells of 375 m, bandwidth 1500 m, peak 1.360e-04 events/m²
64x74 cells of 625 m, bandwidth 2500 m, peak 8.233e-05 events/m²
 bandwidth  blobs area_share     peak
       150     30       0.8% 1.58e-03
       300      6       1.9% 9.08e-04
       600      5       7.2% 3.87e-04
      1000      2      11.7% 2.05e-04
      1500      2      13.2% 1.36e-04
      2500      1      14.7% 8.23e-05

Thirty areas at 150 m is the estimator drawing individual busy streets; one at 2,500 m is everything merged into the centre. The count holds only briefly — 5 or 6 between 300 and 600 m, 2 from 1,000 to 1,500 m — and that is itself the finding: Leeds crime has one dominant centre and a handful of secondary concentrations, and "how many hotspots" has no single answer. "Six hotspots" is a bandwidth.

Notice the peak column falling by a factor of 19 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, *, weights=None, floor_per_km2=1000, **kwargs):
    """Density of cases divided by density of the population at risk, on one shared grid."""
    bounds = np.vstack([cases.total_bounds, controls.total_bounds])
    bounds = (*bounds[:, :2].min(axis=0), *bounds[:, 2:].max(axis=0))
    case_surface, transform = kde_surface(cases, bandwidth, bounds=bounds, **kwargs)
    control_surface, _ = kde_surface(controls, bandwidth, bounds=bounds, weights=weights, **kwargs)

    # suppress the ratio where there is essentially no population — it explodes there
    floor = floor_per_km2 / 1e6                         # people per m²
    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


residents = zones.assign(geometry=zones.representative_point())   # each LSOA's population at its centre
risk, transform = relative_risk(incidents, residents, bandwidth=1000, weights=residents["population"])

127x150 cells of 250 m, bandwidth 1000 m, peak 2.055e-04 events/m²
127x150 cells of 250 m, bandwidth 1000 m, peak 7.369e-03 events/m²
relative risk: median 0.01, max 0.04, 79% suppressed

The ratio is crimes per resident in the month: a median of 0.01 and a maximum of 0.04, or 40 per 1,000, with 79% of the padded grid suppressed because fewer than 1,000 people per km² live there. 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.

Peak heights, edge bias and why KDE is not a test

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 1.58e-03 to 8.23e-05 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.

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 66 km in Leeds.

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.

Duplicate points, weights and adaptive bandwidths

• 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.

Internal links

• How to make a point density heatmap in Python — the practical implementation
• Your heatmap looks wrong: KDE bandwidth and cell size — diagnosing a surface that misleads
• The modifiable areal unit problem explained — the aggregation problem KDE avoids
• How to find hotspots with Getis-Ord Gi* in Python — a test, where KDE is only a description
• How to bin points into hexagons in Python — the discrete alternative
• Projected vs geographic CRS explained — why bandwidth needs metres
• The raster data model explained — what the output surface is
• Choropleth classification explained — the same "breaks decide the story" problem in colour
• Re-identification risk in spatial data explained — why a fine-bandwidth surface still contains the points

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.