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 #

Check the three in this order:

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  ", surface.max())

CRS          EPSG:4326 | geographic: True
extent       [-2. 53. -2. 54.]
bandwidth    300
cell         75
bandwidth/cell = 4.0
peak value   -1.2e+02

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

Step-by-step solution #

1. Compare the bandwidth to the extent #

The single most useful diagnostic is one division:

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 9,000 · bandwidth 300 · ratio 30.0

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", surface.min(), "max", surface.max())

min -302.7 max -121.9

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(shape)

4. Check the cell-to-bandwidth ratio #

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

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.

cases, _ = density_surface(incidents, bandwidth=300)
controls, transform = density_surface(households, bandwidth=300)

floor = np.percentile(controls, 20)
risk = np.where(controls > floor, cases / controls, np.nan)
print(f"raw peak at {np.unravel_index(cases.argmax(), cases.shape)}")
print(f"risk peak at {np.unravel_index(np.nanargmax(risk), risk.shape)}")

raw peak at (188, 241)
risk peak at (96, 402)

Two different places. The raw peak is the city centre, because that is where people are. The risk peak is where events are common relative to the population — which is almost always the question that was actually asked.

Code examples #

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


diagnose_kde(incidents, surface, bandwidth=300, cell=75)

extent 9,000 · bandwidth 300 · cell 75 · ratio 30.0 · peak 7.225e-04
  ✓ parameters are in a sane relationship

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

extent 0 · bandwidth 300 · cell 75 · ratio 0.0 · peak -1.219e+02
  ✗ geographic CRS (WGS 84) — bandwidth is in DEGREES
  ✗ extent/bandwidth = 0.0 — kernel covers the study area
  ✗ negative values (min -302.7) — np.exp() is missing

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

from scipy.ndimage import label
import pandas as pd


def stability_sweep(points, bandwidths, threshold=0.5):
    rows = []
    for bw in bandwidths:
        surface, _ = density_surface(points, bw, cell=bw / 4)
        hot = surface > threshold * surface.max()
        rows.append({"bandwidth": bw, "blobs": label(hot)[1],
                     "area": f"{hot.mean():.1%}"})
    frame = pd.DataFrame(rows)

    counts = frame["blobs"].tolist()
    stable = [(counts[i], frame["bandwidth"][i]) for i in range(len(counts))]
    from itertools import groupby
    runs = [(k, [b for _, b in g]) for k, g in groupby(stable, key=lambda t: t[0])]
    longest = max(runs, key=lambda r: len(r[1]))
    print(frame.to_string(index=False))
    print(f"\nstable at {longest[0]} blobs for bandwidth "
          f"{min(longest[1])}–{max(longest[1])}")
    return frame


stability_sweep(incidents, [0.7, 1.5, 3, 5, 8, 14, 25])

 bandwidth  blobs  area
       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%

stable at 3 blobs for bandwidth 1.5–8

The plateau is the finding. Twelve blobs at 0.7 is the estimator drawing individual points; one blob at 25 is two genuine clusters merged. Report the plateau and the bandwidth range that produces it.

Example 3 — the duplicate-coordinate spike #

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

(384766.2, 398772.1)    147
(385010.9, 398201.4)      6
(384880.6, 398702.4)      3
147 coordinates with more than one point

One coordinate with 147 points is a geocoding artefact — every unmatched address fell back to the same city centroid, as in geocoding returns wrong coordinates. KDE faithfully renders it as an enormous spike that dominates the colour ramp.

Two honest responses:

# a) drop the artefact if you know what it is
clean = points[~points.geometry.apply(lambda g: (round(g.x, 1), round(g.y, 1)))
                .eq((384766.2, 398772.1))]

# b) collapse genuine repeats to weighted points
grouped = points.groupby([points.geometry.x, points.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=points.crs,
)
surface, transform = density_surface(weighted, bandwidth=300, weights=weighted["weight"])

(b) is right when repeats are real — several incidents at one address. (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.

Explanation #

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.6× in Manchester, 2× at 60°N.

The map renders and looks fine. It has simply been smoothed 60% 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.

Edge cases or notes #

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