Your Heatmap Looks Wrong: KDE Bandwidth and Cell Size

Problem statement

The heatmap renders. It just does not look like anything:

• One enormous orange blob covering most of the study area.
• A field of tiny dots, one per point, which is the scatter plot you were trying to replace.
• A uniform wash of a single colour, with no structure at all.
• A hotspot in a park, where nothing happens and nobody lives.
• Values that are negative, which no density can be.

None of these raise an exception. Kernel density estimation always succeeds; it produces a surface for any parameters you hand it, including parameters that make the surface meaningless.

Almost all of it comes down to three numbers being in the wrong relationship to each other: the bandwidth, the cell size, and the units your coordinates are in.

Quick answer

Here is a typical broken attempt on Leeds' 9,229 crimes (see The data used on this page) — the points left in the CRS they arrived in, and the exponential forgotten — followed by the three checks, in order:

points = incidents.to_crs("EPSG:4326")                  # as the data arrived
bandwidth, cell = 300, 75
xy = np.column_stack([points.geometry.x, points.geometry.y])
minx, miny, maxx, maxy = points.total_bounds
gx, gy = np.meshgrid(np.arange(minx, maxx + cell, cell), np.arange(miny, maxy + cell, cell))
grid_points = np.column_stack([gx.ravel(), gy.ravel()])
model = KernelDensity(bandwidth=bandwidth).fit(xy)
surface = model.score_samples(grid_points)              # np.exp() missing

print("CRS         ", points.crs, "| geographic:", points.crs.is_geographic)
print("extent      ", points.total_bounds.round(0))
print("bandwidth   ", bandwidth)
print("cell        ", cell, f"(bandwidth/cell = {bandwidth / cell:.1f}, want >= 4)")
print("peak value  ", round(surface.max(), 6))

CRS          EPSG:4326 | geographic: True
extent       [-2. 54. -1. 54.]
bandwidth    300
cell         75 (bandwidth/cell = 4.0, want >= 4)
peak value   -13.245442

Two faults visible at once: a geographic CRS with a bandwidth of 300 (degrees), and a negative peak (the np.exp was skipped).

Bandwidth against extent, CRS, exponential, cell ratio

1. Compare the bandwidth to the extent

The single most useful diagnostic is one division:

points = incidents                                      # projected: EPSG:27700
minx, miny, maxx, maxy = points.total_bounds
extent = max(maxx - minx, maxy - miny)
print(f"extent {extent:,.0f} · bandwidth {bandwidth} · ratio {extent / bandwidth:.1f}")

extent 30,707 · bandwidth 300 · ratio 102.4

A ratio of roughly 10 to 100 produces a readable surface. Outside that:

• Ratio below 5 — the kernel is a large fraction of the study area. Everything merges into one blob.
• Ratio above 500 — the kernel is a speck. You get one dot per point.

If your extent is in degrees and shows as 0.15, the CRS is the problem, not the bandwidth.

2. Rule out the CRS

if points.crs.is_geographic:
    raise ValueError(
        f"KDE on {points.crs.name}: bandwidth {bandwidth} means {bandwidth} DEGREES"
    )

This is the cause of both the "one giant blob" and the "uniform wash" symptoms, and it is invisible in the output. A bandwidth of 300 degrees smooths the entire planet into a constant; a bandwidth of 0.003 degrees is about 200 m at UK latitudes and roughly 330 m east-west at the equator — so the same number means different distances in different places on the same map.

Reproject first, always. Choose a projected CRS for your area.

3. Check for the missing exponential

print("min", round(surface.min(), 6), "max", round(surface.max(), 6))

min -13.307697 max -13.245442

Densities cannot be negative. KernelDensity.score_samples() returns the logarithm of the density, because real density values underflow float64. Forgetting np.exp gives a map with roughly the right shape and entirely wrong values — the ramp still renders, so nothing looks broken.

surface = np.exp(model.score_samples(grid_points)).reshape(gx.shape)

4. Check the cell-to-bandwidth ratio

print(f"bandwidth/cell = {bandwidth / cell:.1f}")

bandwidth/cell = 4.0

Below about 4, the surface is visibly blocky and you have quietly reintroduced the grid artefacts KDE exists to avoid — see the modifiable areal unit problem. Above about 20, you are paying for resolution the smoothing has already removed.

Cell size is the one parameter that is genuinely cosmetic, provided it stays well under the bandwidth.

5. Ask what the denominator is

A hotspot in an empty park is usually not an error in the code. It is a density of events where you wanted a density of risk.

residents = zones.assign(geometry=zones.representative_point())     # each LSOA's population at its centre
extent = zones.total_bounds                                          # one grid for both surfaces
cases, transform = density_surface(incidents, bandwidth=1000, bounds=extent)
controls, _ = density_surface(residents, bandwidth=1000, weights=residents["population"], bounds=extent)
peaks = {"raw density": cases.argmax()}
for name, floor in (("risk, floor at the 20th percentile", np.percentile(controls, 20)),
                    ("risk, floor 1,000 residents/km²", 1000 / 1e6)):
    risk = np.where(controls > floor, cases / controls, np.nan)      # crimes per resident
    peaks[name] = np.nanargmax(risk)
for name, index in peaks.items():
    row, col = np.unravel_index(index, cases.shape)
    x, y = transform * (col + 0.5, row + 0.5)
    print(f"{name:<35} peak at ({x:,.0f}, {y:,.0f})  {cases[row, col] * 1e6:6.1f} crimes/km²  "
          f"{controls[row, col] * 1e6:7,.0f} residents/km²")

raw density                         peak at (430,345, 433,800)   204.8 crimes/km²    5,526 residents/km²
risk, floor at the 20th percentile  peak at (449,845, 431,800)     0.0 crimes/km²        0 residents/km²
risk, floor 1,000 residents/km²     peak at (430,095, 433,300)   195.6 crimes/km²    4,894 residents/km²

Three results, and the middle one is the trap. With the floor at the 20th percentile of the control surface, the "risk" peak is in empty countryside at the edge of the padded grid: the percentile of a mostly-empty grid is essentially zero, so a sliver of density divided by a sliver of population wins. With a floor in real units — 1,000 residents per km² — the risk peak is back in the city centre, 560 m from the raw peak.

That last result is not a failure of the method; it says residents are the wrong denominator here. A city centre's crime happens to the people who visit it, and few of them live there. The denominator has to be the population at risk — footfall, shop floor space, visitors — and choosing it is the actual question, not a detail of the code.

A diagnostic that names the fault in a heatmap

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 heatmaps are of the crime points, incidents; the population in zones is the denominator for relative risk. density_surface is the function from the heatmap guide.

import io
import time

from rasterio.transform import from_origin
from sklearn.neighbors import KernelDensity

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")


def density_surface(points, bandwidth, *, cell=None, weights=None, pad_factor=3, bounds=None):
    """Gaussian KDE in events per square map unit on a padded grid, with its transform.
    Pass the same bounds= to surfaces you intend to divide or compare, so their grids match."""
    cell = cell or bandwidth / 4
    pad = bandwidth * pad_factor
    base = points.total_bounds if bounds is None else np.asarray(bounds)
    minx, miny, maxx, maxy = base + np.array([-pad, -pad, pad, pad])
    xs = minx + (np.arange(int(np.ceil((maxx - minx) / cell))) + 0.5) * cell
    ys = maxy - (np.arange(int(np.ceil((maxy - miny) / cell))) + 0.5) * cell
    gx, gy = np.meshgrid(xs, ys)
    xy = np.column_stack([points.geometry.x, points.geometry.y])
    model = KernelDensity(bandwidth=bandwidth).fit(xy, sample_weight=weights)
    total = len(points) if weights is None else float(np.sum(weights))
    log_density = model.score_samples(np.column_stack([gx.ravel(), gy.ravel()]))
    return np.exp(log_density).reshape(gx.shape) * total, from_origin(minx, maxy, cell, cell)

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

Example 1 — a diagnostic that names the fault

import numpy as np


def diagnose_kde(points, surface, bandwidth, cell):
    """Print every check that distinguishes a good density surface from a broken one."""
    problems = []

    if points.crs is None:
        problems.append("no CRS set — bandwidth units are unknown")
    elif points.crs.is_geographic:
        problems.append(f"geographic CRS ({points.crs.name}) — bandwidth is in DEGREES")

    minx, miny, maxx, maxy = points.total_bounds
    extent = max(maxx - minx, maxy - miny)
    ratio = extent / bandwidth
    if ratio < 5:
        problems.append(f"extent/bandwidth = {ratio:.1f} — kernel covers the study area")
    elif ratio > 500:
        problems.append(f"extent/bandwidth = {ratio:.1f} — kernel is a speck; one dot per point")

    if bandwidth / cell < 4:
        problems.append(f"bandwidth/cell = {bandwidth / cell:.1f} — surface will look blocky")

    if surface.min() < 0:
        problems.append(f"negative values (min {surface.min():.1f}) — np.exp() is missing")

    finite = surface[np.isfinite(surface)]
    if finite.size and finite.max() > 0:
        spread = finite.max() / max(np.percentile(finite, 50), 1e-12)
        if spread < 2:
            problems.append(f"peak is only {spread:.1f}x the median — surface is nearly flat")

    dupes = points.geometry.duplicated().sum()
    if dupes > len(points) * 0.02:
        problems.append(f"{dupes} duplicate coordinates ({dupes / len(points):.0%}) — spikes likely")

    print(f"extent {extent:,.0f} · bandwidth {bandwidth} · cell {cell} · "
          f"ratio {ratio:.1f} · peak {finite.max():.3e}")
    for problem in problems:
        print(f"  ✗ {problem}")
    if not problems:
        print("  ✓ parameters are in a sane relationship")
    return problems


good, _ = density_surface(incidents, 300, cell=75)
diagnose_kde(incidents, good, bandwidth=300, cell=75)

extent 30,707 · bandwidth 300 · cell 75 · ratio 102.4 · peak 9.075e-04
  ✗ 5226 duplicate coordinates (57%) — spikes likely

Run it on the broken case from the quick answer and it names all of them:

diagnose_kde(incidents.to_crs("EPSG:4326"), model.score_samples(grid_points), bandwidth=300, cell=75)

extent 0 · bandwidth 300 · cell 75 · ratio 0.0 · peak -1.325e+01
  ✗ geographic CRS (WGS 84) — bandwidth is in DEGREES
  ✗ extent/bandwidth = 0.0 — kernel covers the study area
  ✗ negative values (min -13.3) — np.exp() is missing
  ✗ 5226 duplicate coordinates (57%) — spikes likely

Example 2 — finding the bandwidth range where the answer is stable

from itertools import groupby

import pandas as pd
from scipy.ndimage import label


def stability_sweep(points, bandwidths, thresholds=(0.5, 0.1)):
    """Count separate hot areas above a fraction of the peak, at each bandwidth and threshold."""
    rows = []
    for bw in bandwidths:
        surface, _ = density_surface(points, bw, cell=bw / 4)
        row = {"bandwidth": bw}
        for t in thresholds:
            hot = surface > t * surface.max()
            row[f"blobs>{t:g}"] = label(hot)[1]
            row[f"area>{t:g}"] = f"{hot.mean():.1%}"
        rows.append(row)
    frame = pd.DataFrame(rows)
    print(frame.to_string(index=False))
    for t in thresholds:
        counts = frame[f"blobs>{t:g}"].tolist()
        runs = [(k, len(list(g))) for k, g in groupby(counts)]
        k, n = max(runs, key=lambda r: r[1])
        print(f"above {t:g} of the peak: longest run {n} bandwidths at {k} blob(s)")
    return frame


stability_sweep(incidents, [150, 300, 600, 1000, 1500, 2500])

 bandwidth  blobs>0.5 area>0.5  blobs>0.1 area>0.1
       150          1     0.0%         30     0.8%
       300          1     0.1%          6     1.9%
       600          1     0.3%          5     7.2%
      1000          1     1.1%          2    11.7%
      1500          1     2.0%          2    13.2%
      2500          1     3.4%          1    14.7%
above 0.5 of the peak: longest run 6 bandwidths at 1 blob(s)
above 0.1 of the peak: longest run 2 bandwidths at 2 blob(s)

Read the two thresholds separately. Above half the peak there is one blob at every bandwidth — perfectly "stable", and meaningless: on Leeds crime the city centre is so much denser than anywhere else that nothing else ever reaches half its height. Above a tenth of the peak the surface has structure — 30 blobs at 150 m (the estimator drawing individual busy streets), 5 or 6 between 300 and 600 m, 2 from 1,000 m — but no long plateau. A threshold relative to the peak is only meaningful when no single peak dominates; report the bandwidth, the threshold and the counts together, and when they never settle, say that.

Example 3 — the duplicate-coordinate spike

counts = incidents.geometry.apply(lambda g: (round(g.x, 1), round(g.y, 1))).value_counts()
print(counts.head(3))
print(f"{(counts > 1).sum()} coordinates with more than one point")

geometry
(430269.0, 433400.0)    69
(430415.0, 433639.0)    45
(428400.0, 429312.9)    41
Name: count, dtype: int64
1866 coordinates with more than one point

69 crimes at one coordinate, 45 at the next. These are not errors: police.uk moves every crime to an anonymised "on or near" point, so a busy street's offences share one. Geocoding produces the same pattern as an artefact — unmatched addresses falling back to one city centroid, as in geocoding returns wrong coordinates. Either way KDE renders it faithfully as a spike.

Two honest responses:

# a) drop the artefact if you know what it is
top = counts.index[0]
clean = incidents[~incidents.geometry.apply(lambda g: (round(g.x, 1), round(g.y, 1))).isin([top])]

# b) collapse genuine repeats to weighted points
grouped = incidents.groupby([incidents.geometry.x, incidents.geometry.y]).size()
weighted = gpd.GeoDataFrame(
    {"weight": grouped.values},
    geometry=gpd.points_from_xy(grouped.index.get_level_values(0),
                                grouped.index.get_level_values(1)),
    crs=incidents.crs,
)
surface, transform = density_surface(weighted, bandwidth=300, weights=weighted["weight"])
print(f"{len(incidents):,} points collapsed to {len(weighted):,} weighted locations; "
      f"peak {surface.max():.3e} (unweighted duplicates give the same surface: {np.isclose(surface.max(), good.max())})")

9,229 points collapsed to 4,003 weighted locations; peak 9.075e-04 (unweighted duplicates give the same surface: True)

(b) is right when repeats are real — several incidents at one address — and it changes nothing about the surface: less than half as many locations, the same peak. It is a speed-up, not a correction. (a) is right when they are a geocoding failure. Deciding which requires knowing where the data came from, which no parameter can do for you.

Plausible wrong maps and hotspots in parks

Why a geographic CRS produces a plausible wrong map

If the bandwidth is far too large in degrees, you get a flat wash and notice. The dangerous case is a bandwidth that is roughly right — say 0.003 degrees.

A degree of latitude is about 111 km everywhere. A degree of longitude is 111 km at the equator and about 69 km at 51°N. So a circular kernel in degree-space is an ellipse on the ground, stretched east-west by 1/cos(latitude) — 1.7× in Leeds, 2× at 60°N.

The map renders and looks fine. It has simply been smoothed 70% more in one direction than the other, and the distortion changes across the map. Nothing in the output reveals it.

Why negative values do not raise

score_samples returns log density for a numerical reason: a density of 1e-320 underflows to exactly zero in float64, and any ratio or comparison involving it breaks. Its logarithm, −737, is an unremarkable number.

Matplotlib maps any float range onto a colour ramp, so a surface of log densities produces a map with the correct shape — the peaks are still peaks — and completely wrong values. It is only obviously wrong if you look at the numbers, which is why the diagnostic checks surface.min() < 0 first.

Why "the hotspot is in a park" is usually not a bug

A raw KDE of events shows where events are, and events happen where people are. On any city-wide map the peak is the busiest place, which is a fact about footfall, not about risk.

The correction is a relative risk surface: divide the case density by a control density built the same way — population, households, footfall, whatever the appropriate denominator is. Both surfaces must use the same bandwidth and grid, or you are comparing differently smoothed quantities.

And the ratio must be floored where the denominator is near zero, or empty countryside produces enormous risk values from two events and almost no population.

Why the cell size gets blamed for bandwidth problems

Blockiness looks like a resolution problem, so the instinct is to shrink the cells. That fixes blockiness — and blockiness only.

Everything else on the symptom list is the bandwidth. Halving the cell size on a surface that is one giant blob produces a smoother giant blob, four times slower. Check the extent-to-bandwidth ratio before touching the cell size.

Rule-of-thumb bandwidths, all-NaN surfaces and log ramps

• Rule-of-thumb bandwidths over-smooth spatial data. Scott's and Silverman's rules assume a single unimodal bump; clustered, bounded point patterns are neither.
• gaussian_filter takes sigma in cells. If you switched to the binned approximation, sigma = bandwidth / cell — passing the bandwidth directly smooths by a factor of cell too much.
• A surface that is entirely NaN usually means the grid and the points are in different CRSs, so no point is near any cell.
• Memory scales with the grid, not the points. nrows * ncols * 8 bytes before the estimator allocates anything.
• Clipping before estimating creates a false edge. Clip the finished surface for display instead.
• Log-scale colour ramps hide bandwidth problems. They make any surface look structured. Use a linear ramp while diagnosing.
• Compare like with like. Two heatmaps are only comparable if bandwidth, cell size and extent all match.

Internal links

• Kernel density explained — what the bandwidth is actually doing
• How to make a point density heatmap in Python — the implementation with the guardrails in place
• Choose a projected CRS for your area — fixing the commonest cause
• Geocoding returns wrong or missing coordinates — the source of duplicate-coordinate spikes
• The modifiable areal unit problem explained — what an oversized cell reintroduces
• How to find hotspots with Getis-Ord Gi* in Python — when you need a test rather than a picture
• Projected vs geographic CRS explained — why degrees break distance parameters
• Raster and vector do not line up in Python — when the surface is right but misplaced

FAQ

Why is my heatmap one giant blob?

The bandwidth is too large relative to the extent — usually because the CRS is geographic and the bandwidth is being read as degrees. Check extent / bandwidth; you want roughly 10 to 100.

Why does my heatmap look like the scatter plot I was replacing?

The bandwidth is far too small. Each point is being drawn as its own peak. Raise it to the distance the process actually operates over.

Why are my density values negative?

score_samples returns log density. Wrap it in np.exp(). The map still renders without it, which is why this is easy to miss.

Why is the hotspot somewhere nothing happens?

You are mapping event density, not risk. Divide by a control density — population or households — built with the same bandwidth and grid, and floor the ratio where the denominator is near zero.

Why is my surface blocky?

The cells are too large relative to the bandwidth. Use $\text{cell} \le \text{bandwidth}/4$. This is the only symptom that cell size causes.

Why does the density fade at the edges of my study area?

Edge effect — kernels of boundary points extend outside the grid and that mass is lost. Pad the grid by about three bandwidths and crop for display.

How do I know I picked the right bandwidth?

Sweep a range and look for a plateau where the number of distinct hotspots does not change. Report the plateau, not a single value.