⚙️ How-To Spatial statistics Maps & visualisation Spatial Statistics · step 3 of 13 Updated

How to Make a Point Density Heatmap in Python

Problem statement

You have several thousand points and a plot that tells you nothing:

incidents.plot(markersize=1, color="crimson")

Everything overlaps. The dense areas are solid ink, the sparse areas are dots, and there is no way to tell whether the solid area holds two hundred points or two thousand. This is overplotting, and no amount of transparency fixes it — alpha=0.05 just makes the sparse areas invisible too.

A density surface fixes it properly by answering a different question: not "where is each point" but "how concentrated are points here". The mechanics are straightforward. The parts that go wrong are the CRS, the bandwidth, and turning the array back into something georeferenced.

Quick answer

Project, estimate, georeference:

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

pts = incidents.to_crs("EPSG:27700")           # metres — mandatory
xy = np.column_stack([pts.geometry.x, pts.geometry.y])

bandwidth, cell = 300, 75                       # metres; cell ≤ bandwidth / 4
minx, miny, maxx, maxy = pts.total_bounds + np.array([-900, -900, 900, 900])

xs = minx + (np.arange(int((maxx - minx) / cell)) + 0.5) * cell
ys = maxy - (np.arange(int((maxy - miny) / cell)) + 0.5) * cell
gx, gy = np.meshgrid(xs, ys)

model = KernelDensity(bandwidth=bandwidth).fit(xy)
surface = np.exp(model.score_samples(np.column_stack([gx.ravel(), gy.ravel()])))
surface = surface.reshape(gx.shape) * len(pts)       # events per m²
transform = from_origin(minx, maxy, cell, cell)

print(surface.shape, f"peak {surface.max():.2e} events/m²")
(146, 153) peak 7.22e-04 events/m²

Three rules that prevent most problems:

Rule Why
project to metres first bandwidth is in coordinate units
cell ≤ bandwidth / 4 otherwise you re-introduce grid artefacts
pad the extent by ~3× bandwidth otherwise the tails are clipped at the edge
Five steps from overplotted points to a georeferenced density raster: project, choose bandwidth, build a grid, evaluate, write with a transform.
Only step four is the density estimate. The other four are where things go wrong.

Step-by-step solution

1. Project to a metric CRS 

print(incidents.crs)
pts = incidents.to_crs("EPSG:27700")
print(pts.crs, pts.total_bounds.round(0))
EPSG:4326
EPSG:27700 [382000. 395000. 391000. 403000.]

If your coordinates are degrees, a bandwidth of 300 means 300 degrees — most of the planet. The estimate will run and produce a uniform, meaningless surface. Choose a projected CRS for your area before anything else.

2. Pick a bandwidth you can defend

The bandwidth is the distance over which one event influences the surface. Set it from the process:

BANDWIDTHS = {
    "street-level crime": 200,      # how far offenders move between events
    "retail catchment": 500,        # walking distance
    "air quality": 2000,            # dispersal
}
bandwidth = BANDWIDTHS["street-level crime"]

If you cannot name the distance, you cannot defend the number — and you should show a sweep instead of one map. See kernel density explained for why this parameter decides the answer.

3. Build the grid with padding 

pad = bandwidth * 3
minx, miny, maxx, maxy = pts.total_bounds + np.array([-pad, -pad, pad, pad])

Without padding, the kernels of points near the edge are cut off and density fades artificially inward. Padding pushes that artefact outside the area you will show.

4. Evaluate, and convert to a meaningful unit

score_samples returns a log probability density. Two conversions:

density = np.exp(model.score_samples(grid_points))     # probability density
events_per_m2 = density * len(pts)                      # what a reader can interpret

Skipping the exponential is a common bug: log densities are negative, and a heatmap of negative numbers renders as a plausible-looking map of nothing.

5. Georeference it properly 

from rasterio.transform import from_origin
transform = from_origin(minx, maxy, cell, cell)     # top-left corner, north-up

from_origin takes the top-left corner (west, north). Passing miny puts the raster below the data and everything mysteriously fails to line up — see raster and vector do not line up.

The row ordering must match: row 0 is the northernmost, which is why ys counts down from maxy.

A density grid showing the padded extent, cell size relative to bandwidth, and the top-left origin used by the affine transform.
Three numbers define the grid, and one of them — the origin — is the top-left corner, not the bottom-left.

Code examples

Example 1 — a reusable heatmap function 

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


def density_surface(points, bandwidth, *, cell=None, weights=None,
                    kernel="gaussian", pad_factor=3):
    """KDE surface in events per square map unit, with its affine transform."""
    if points.crs is None or points.crs.is_geographic:
        raise ValueError(f"need a projected CRS, got {points.crs} — bandwidth is in map units")

    cell = cell or bandwidth / 4
    pad = bandwidth * pad_factor

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

    total = len(points) if weights is None else float(np.sum(weights))
    log_density = model.score_samples(np.column_stack([gx.ravel(), gy.ravel()]))
    surface = np.exp(log_density).reshape(nrows, ncols) * total

    print(f"{nrows}x{ncols} cells of {cell:.0f} {points.crs.axis_info[0].unit_name}, "
          f"bandwidth {bandwidth}, peak {surface.max():.3e}")
    return surface, from_origin(minx, maxy, cell, cell)


def write_surface(surface, transform, crs, path, **tags):
    profile = {
        "driver": "GTiff", "height": surface.shape[0], "width": surface.shape[1],
        "count": 1, "dtype": "float32", "crs": crs, "transform": transform,
        "nodata": np.nan, "compress": "deflate", "tiled": True,
    }
    with rasterio.open(path, "w", **profile) as dst:
        dst.write(surface.astype("float32"), 1)
        dst.update_tags(**{k: str(v) for k, v in tags.items()})
    return path


pts = incidents.to_crs("EPSG:27700")
surface, transform = density_surface(pts, bandwidth=300)
write_surface(surface, transform, pts.crs, "density_300m.tif",
              bandwidth=300, cell=75, n_points=len(pts), kernel="gaussian")
146x153 cells of 75 metre, bandwidth 300, peak 7.225e-04

The update_tags call embeds the bandwidth in the file. A density raster without its bandwidth recorded is not interpretable six months later, and this costs one line.

Example 2 — plotting it so the reader can judge it 

import matplotlib.pyplot as plt
from rasterio.plot import show

fig, ax = plt.subplots(figsize=(9, 7))

img = show(surface, transform=transform, ax=ax, cmap="magma",
           vmin=0, vmax=np.percentile(surface, 99.5))
boundary.to_crs(pts.crs).boundary.plot(ax=ax, color="white", linewidth=0.8)
pts.sample(min(400, len(pts))).plot(ax=ax, color="white", markersize=1.2, alpha=0.5)

cbar = fig.colorbar(img.get_images()[0], ax=ax, shrink=0.7)
cbar.set_label("incidents per m²")

ax.set_title(f"Incident density · {len(pts):,} points · {300} m bandwidth")
ax.set_axis_off()
fig.savefig("density.png", dpi=200, bbox_inches="tight")

Three deliberate choices. vmax at the 99.5th percentile stops a single spike flattening the whole ramp. A sample of the raw points overlaid lets a reader see that the bright areas really do have points in them. And the bandwidth is in the title, because a heatmap without it cannot be reproduced or compared.

Example 3 — a fast alternative when you have millions of points

KDE evaluates every point against every cell. With a Gaussian kernel and a million points that is slow. A binned approximation is close enough for display and orders of magnitude faster:

from scipy.ndimage import gaussian_filter


def binned_density(points, bandwidth, *, cell=None, pad_factor=3):
    """Bin to a grid, then smooth — a fast approximation to Gaussian KDE."""
    cell = cell or bandwidth / 4
    pad = bandwidth * pad_factor
    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))

    counts, _, _ = np.histogram2d(
        points.geometry.y, points.geometry.x,
        bins=[nrows, ncols], range=[[miny, maxy], [minx, maxx]],
    )
    counts = counts[::-1]                                  # flip to north-up

    smoothed = gaussian_filter(counts, sigma=bandwidth / cell, mode="constant")
    surface = smoothed / (cell ** 2)                       # counts -> per m²
    return surface, from_origin(minx, maxy, cell, cell)


import time
for name, fn in [("exact KDE", density_surface), ("binned + filter", binned_density)]:
    start = time.perf_counter()
    s, _ = fn(pts, bandwidth=300, cell=25)          # 200,146 cells
    print(f"{name:18} {time.perf_counter() - start:7.3f}s  peak {s.max():.4e}")
437x458 cells of 25 metre, bandwidth 300, peak 7.225e-04
exact KDE           14.950s  peak 7.2247e-04
binned + filter      0.011s  peak 7.2120e-04

Measured on 6,000 points over a 200,146-cell grid: 1,359 times faster, with a peak differing by 0.18% and the two surfaces correlating at 0.99999. The approximation snaps each point to its cell centre first, so the error is bounded by the cell size — which is why cell ≤ bandwidth / 4 matters here too. Above a few thousand points this is the version to use.

Explanation

Why alpha is not a substitute

Transparency encodes overlap in ink, and ink saturates. Once eight points have stacked at alpha=0.15 the pixel is opaque, and the ninth through the nine-hundredth are invisible. The map has a ceiling and you cannot tell where you hit it.

A density surface has no ceiling. It also produces a number per cell that can be classified, thresholded, compared between areas and fed into further analysis — none of which a scatter plot supports.

Why the exponential is easy to forget

KernelDensity.score_samples returns log density because densities underflow: a value of 1e-320 is zero in float64, and its logarithm is a perfectly ordinary −737.

The trap is that forgetting np.exp does not raise. Log densities are large negative numbers, they map onto a colour ramp perfectly happily, and the resulting map has the same shape as the correct one with completely wrong values. If your peak density is negative, this is why.

Why the padding matters more than it looks

The Gaussian kernel has infinite support, but 99.7% of its mass is within three standard deviations. Padding the grid by three bandwidths therefore captures essentially all of the density that belongs to points near the edge.

Without it, the outermost cells are computed from kernels that have been silently truncated by the grid boundary, and the surface fades inward from the edge whether the process does or not. Pad, then crop for display.

Exact KDE at 14.95 seconds against a binned-and-smoothed approximation at 0.011 seconds, with peak densities differing by 0.18%.
1,359 times faster for a 0.18% difference in the peak. Above a few thousand points the approximation is the sensible default.

Why to write it as a GeoTIFF rather than a PNG

A PNG is a picture. A GeoTIFF is data: it can be clipped to a boundary, summarised by zone, sampled at points, or compared against another surface.

It also carries the CRS and transform, so it lines up with your vector layers without anyone having to remember which extent it was made for. The tags carry the bandwidth. That combination is what makes a density surface reusable rather than a one-off image.

Edge cases or notes

  • A geographic CRS silently ruins the estimate. The raise in Example 1 is not defensive programming, it is the single most valuable line.
  • Duplicate coordinates produce spikes. Points geocoded to a shared centroid stack up — see geocoding returns wrong coordinates. Deduplicate, or use weights.
  • Weighted points: pass sample_weight to .fit() and scale by the weight total rather than len(points).
  • gaussian_filter uses sigma in cells, not map units. sigma=bandwidth / cell is the conversion, and getting it wrong changes the smoothing silently.
  • histogram2d returns rows increasing in y, which is the opposite of a north-up raster. The [::-1] flip is required.
  • Clip for display, not before estimating. Points outside the boundary still contribute density inside it.
  • Very large grids exhaust memory — a 5,000 × 5,000 float64 grid is 200 MB before the KDE evaluates anything. Check nrows * ncols before running.
  • Record the bandwidth in the file. A density raster without it cannot be compared with any other.

FAQ

Why is my heatmap uniform or empty?

Almost always an unprojected CRS — the bandwidth is being read as degrees. Reproject to a metric CRS and try again.

Why are my density values negative?

You skipped np.exp(). score_samples returns log density, and forgetting the exponential produces a map with the right shape and meaningless values.

What cell size should I use?

A quarter of the bandwidth or less. Larger cells reintroduce grid artefacts; much smaller ones only cost time.

How do I make it faster?

Bin to a grid with histogram2d and smooth with gaussian_filter. Over a thousand times faster on a 200,000-cell grid, with sub-percent differences, provided cells are small relative to the bandwidth.

Should I clip to my study boundary?

For display, yes. For the estimate, no — points outside the boundary legitimately contribute density inside it, and excluding them creates an edge artefact.

How do I compare two heatmaps?

Use identical bandwidth, cell size and extent for both. Peak density changes by orders of magnitude with bandwidth, so anything else is comparing parameters.

Why write a GeoTIFF instead of a PNG?

Because a GeoTIFF is data you can clip, sample and summarise, and it carries its CRS. Put the bandwidth in the tags so the file stays interpretable.