Create a Kernel Density Heatmap Raster in PyQGIS

A map of ten thousand points says little beyond "there are many". A density surface turns the same points into a continuous field: high where points crowd together, low where they thin out, smooth in between. Kernel density estimation does this by placing a small hill — the kernel — over every point and adding the hills up. The result is a raster you can classify, contour, sample and compare, which a renderer-only heatmap cannot give you.

This recipe belongs to Spatial Statistics & Pattern Analysis. It runs qgis:heatmapkerneldensityestimation, explains how radius and pixel size shape the result, adds weights, chooses kernels and output units, and styles and normalises the surface for honest comparison. For the on-the-fly display version, see the heatmap renderer.

Kernels add up to a surfaceAlong a line of five points, each point contributes a smooth kernel shaped like a hill with a width equal to the search radius. Where points are close together their kernels overlap and sum to a tall peak; an isolated point produces a single low hill. The density surface is the sum of all kernels.Each point adds a hill; the hills add upclose points: tall peakisolated: low hillblue: individual kernels · amber: their sum, the density

Prerequisites

  • QGIS 3.34 LTR or newer, or the QGIS 4 series.
  • A point layer in a projected CRS with metre units; radius and pixel size are in layer units.
  • An idea of the scale of the pattern you care about — a city block, a neighbourhood, a region — because that decides the radius.

Run the algorithm

The algorithm needs a radius, a pixel size and an output path. Everything else has workable defaults to start from.

import processing
from qgis.core import QgsProject, QgsRasterLayer

collisions = QgsProject.instance().mapLayersByName("road_collisions")[0]
out = processing.run("qgis:heatmapkerneldensityestimation", {
    "INPUT": collisions,
    "RADIUS": 400,             # metres
    "RADIUS_FIELD": None,
    "PIXEL_SIZE": 20,          # metres
    "WEIGHT_FIELD": None,
    "KERNEL": 0,               # 0 = quartic
    "DECAY": 0,
    "OUTPUT_VALUE": 0,         # 0 = raw, 1 = scaled
    "OUTPUT": "/data/results/collision_density.tif",
})["OUTPUT"]

density = QgsRasterLayer(out, "collision density")
stats = density.dataProvider().bandStatistics(1)
print(density.width(), "x", density.height(), "pixels; max", round(stats.maximumValue, 2))
QgsProject.instance().addMapLayer(density)

Breakdown: The output is a single-band GeoTIFF covering the points' extent plus the radius. Each pixel's value is the sum of kernel contributions from all points within the radius of the pixel centre. With raw output, values are in "kernel units" — not directly points per square kilometre — so use them for relative comparison within one surface, or switch to scaled output and normalise as described below. A 20 m pixel over a city produces a modest raster; halving the pixel size quadruples it.

Choose radius and pixel size

The radius is the single most important choice. It sets the scale at which the surface sees structure: small radii show individual hot corners, large radii show broad districts.

Radius sets the scale of what you seeThe same collisions with three radii. 100 metres: many small spots at individual junctions, noisy and fragmented. 400 metres: corridors and clusters along main roads, the scale of the question. 1500 metres: a single broad blob over the city centre that hides the road structure. Pixel size should be several times smaller than the radius.Too small, about right, too large100 m: junction noise400 m: road corridors1500 m: one blobpixel size: about a fifth of the radius or smaller

There is no universally correct radius; there is a radius that matches the question. For junction safety, a radius of 50–100 m keeps junctions separate. For neighbourhood-level patterns, 300–800 m. A practical starting point is a few times the median nearest-neighbour distance — the nearest neighbour analysis gives it — and then compare a smaller and a larger value.

for radius in (150, 400, 1000):
    path = f"/data/results/collision_density_r{radius}.tif"
    processing.run("qgis:heatmapkerneldensityestimation", {
        "INPUT": collisions, "RADIUS": radius, "PIXEL_SIZE": max(radius / 10, 10),
        "KERNEL": 0, "OUTPUT_VALUE": 1, "OUTPUT": path})
    QgsProject.instance().addMapLayer(QgsRasterLayer(path, f"density r={radius}"))

Breakdown: Generating a few radii side by side makes the scale dependence visible and prevents the first surface you produce from being taken as the truth. Tying pixel size to radius keeps each surface smooth without wasting pixels: a pixel much larger than a fifth of the radius makes the hills blocky, while one much smaller adds file size without adding information.

Weight points

Not every point counts the same. A collision with fatalities matters more than a minor one; a shop with 200 staff more than one with two. A weight field multiplies each point's kernel.

weighted = processing.run("qgis:heatmapkerneldensityestimation", {
    "INPUT": collisions, "RADIUS": 400, "PIXEL_SIZE": 20,
    "WEIGHT_FIELD": "severity_weight",     # e.g. 1 minor, 3 serious, 10 fatal
    "KERNEL": 0, "OUTPUT_VALUE": 1,
    "OUTPUT": "/data/results/collision_density_weighted.tif"})["OUTPUT"]

Breakdown: Weights are multiplied into the kernel heights, so a weight of 10 is equivalent to ten points at the same place. NULL weights count as zero, which silently drops points — check the weight field for missing values first. Weighting by a severity score is a judgement call; document the weights next to the map, because a reader cannot recover them from the surface. RADIUS_FIELD similarly lets each point have its own radius, useful when points represent features of different reach.

Sample the density at locations

A density surface becomes an attribute when you read it at points of interest: the collision density at each school entrance, the crime density at each bus stop. Sampling turns the raster into numbers that can be ranked, joined and reported.

schools = QgsProject.instance().mapLayersByName("school_entrances")[0]
sampled = processing.run("native:rastersampling", {
    "INPUT": schools, "RASTERCOPY": weighted,
    "COLUMN_PREFIX": "kde_", "OUTPUT": "memory:schools_kde"})["OUTPUT"]

ranked = sorted(sampled.getFeatures(), key=lambda f: f["kde_1"] or 0, reverse=True)
for f in ranked[:10]:
    print(f"{f['school_name']:<30} {f['kde_1'] * 1e6:8.1f} weighted collisions / km²")
QgsProject.instance().addMapLayer(sampled)

Breakdown: native:rastersampling adds one field per band, named with the prefix and band number, holding the pixel value under each point. Multiplying scaled values by a million converts them to events per square kilometre, which is the figure to show in a ranking table. Points outside the raster's extent get NULL, hence the or 0 in the sort key. Because the surface already integrates everything within the radius, a single sample summarises the surroundings of each location — a quicker and smoother measure than counting points in a buffer, and directly comparable between locations. The general technique is covered in sampling raster values at points.

Choose kernel and output values

The kernel shape matters much less than the radius, but the output scaling matters for interpretation.

Kernel shapes and output scalingKernel shapes include quartic, the default, a smooth dome; triangular, a cone; uniform, a flat disk; triweight, a steeper dome; and Epanechnikov. All fall to zero at the radius. Output values can be raw, the summed kernel heights, or scaled, which normalises each kernel to integrate to one so the surface estimates points per unit area.Shape matters little; scaling matters moreKERNEL0 quartic (smooth, default)1 triangular · 2 uniform3 triweight · 4 Epanechnikovall reach zero at the radiusOUTPUT_VALUE0 raw: summed heights1 scaled: density per areascaled → comparable surfaces

The quartic default produces smooth, natural-looking surfaces and is the right choice almost always. A uniform kernel counts points within the radius and gives a flat-topped, blocky result; it is useful when you want "number of events within 400 m" literally. Scaled output divides each kernel by its integral, so pixel values estimate events per square map unit — per square metre for a metric CRS. Multiply by 1,000,000 for events per square kilometre.

Normalise for comparison

Two density surfaces can only be compared if they are on the same scale. Comparing this year's collisions with last year's, or one city with another, needs scaled output and a common unit.

import numpy as np
from osgeo import gdal

def per_km2(path):
    ds = gdal.Open(path)
    arr = ds.GetRasterBand(1).ReadAsArray().astype("float64") * 1e6
    return arr, ds

this_year, ds = per_km2("/data/results/collisions_2026_scaled.tif")
last_year, _ = per_km2("/data/results/collisions_2025_scaled.tif")
change = this_year - last_year
print(f"peak density {this_year.max():.1f} /km²; largest increase {change.max():.1f} /km²")

Breakdown: Scaled values multiplied by a million are events per square kilometre, a unit people understand. Differencing two surfaces needs identical extents and pixel grids — run both with the same PIXEL_SIZE and clip both to a common extent, or the arrays will not align; resampling and aligning rasters covers the general case. Writing change back to a GeoTIFF with the same geotransform gives a map of where density rose and fell.

Style the surface so it reads correctly

A density raster's look is part of its message. A sequential colour ramp, transparent low values and a stretch that does not let one extreme pixel flatten everything else make the difference between a useful map and a misleading one.

from qgis.core import (QgsSingleBandPseudoColorRenderer, QgsColorRampShader,
                       QgsRasterShader, QgsStyle)

ramp = QgsStyle.defaultStyle().colorRamp("Inferno")
p95 = float(np.percentile(this_year[this_year > 0], 95)) / 1e6
shader_fn = QgsColorRampShader(0, p95, ramp)
shader_fn.classifyColorRamp(8)
items = shader_fn.colorRampItemList()
lowest = items[0].color
lowest.setAlpha(0)                              # lowest class transparent
items[0].color = lowest
shader_fn.setColorRampItemList(items)
shader = QgsRasterShader()
shader.setRasterShaderFunction(shader_fn)
density.setRenderer(QgsSingleBandPseudoColorRenderer(density.dataProvider(), 1, shader))
density.triggerRepaint()

Breakdown: Stretching to the 95th percentile of non-zero values rather than the maximum stops a single extreme pixel from compressing the rest of the surface into the bottom of the ramp. Making the lowest class transparent lets the basemap show through where density is negligible. A perceptually uniform ramp such as Inferno or Viridis keeps equal steps in value looking like equal steps in colour. More on raster styling in applying a colour ramp to a raster.

QGIS version compatibility

qgis:heatmapkerneldensityestimation has the same parameters on QGIS 3.34 LTR, 3.40 LTR and QGIS 4. Kernel and output-value enums are numbered identically. QgsColorRampShader.classifyColorRamp with a class count has been available throughout QGIS 3.

Troubleshooting

  • The raster is enormous or the run takes forever. Pixel size is far too small for the extent; check width × height before running.
  • The surface is one smooth blob. The radius is too large for the pattern.
  • All values are tiny decimals. You chose scaled output in a metric CRS; multiply by 1e6 for per km².
  • Weighted output looks like unweighted. The weight field is text or mostly NULL.

Conclusion

Run kernel density with a radius chosen for the question and a pixel size about a fifth of it, compare a few radii before trusting one, weight points deliberately, use scaled output for anything you compare, and style with a sequential ramp stretched to a high percentile with transparent low values.

Frequently Asked Questions

Is this the same as the heatmap renderer? The renderer computes a similar surface on the fly for display only. The algorithm produces a raster you can analyse, export and compare.

Does density respect roads or barriers? No, kernels are circular in straight-line distance. Network kernel density needs specialised tools.

Can I get a hot spot significance map from this? Not directly; density shows where events concentrate, not whether that is significant. Use Getis-Ord Gi* for that.

How do I contour the surface? Run gdal:contour on the output, as shown in creating contours from a DEM.