Spatial Statistics & Pattern Analysis in PyQGIS

Maps of points invite conclusions. Crimes look clustered around the station, cafés look evenly spaced along the high street, disease cases look concentrated in one district. Some of those impressions are right and some are artefacts of how human eyes read scattered dots, of where people live, or of the scale the map happens to be drawn at. Spatial statistics replaces impressions with measurements: how clustered, where exactly, at what scale, and how surprising compared with chance.

This guide belongs to Spatial Data Processing & Automation. It is written for analysts who work with point events or area counts — public health, policing, retail, ecology, transport — and want to describe and test patterns with the tools QGIS already has, plus a little NumPy where QGIS stops. It explains the main families of methods, the question each one answers, and the choices that decide whether its result means anything.

Questions and the methods that answer themFive questions about a point pattern and the method for each. Is it clustered at all: nearest neighbour index or a Monte Carlo test against random points. Where are the groups: DBSCAN or k-means clustering. How dense is it everywhere: kernel density estimation. Where is it significantly high or low: Getis-Ord Gi star hot spot analysis. Where is it centred and how spread out: mean centre, standard distance and deviational ellipse. Voronoi and Delaunay structures and random sampling support several of these.Start from the questionclustered at all?where are the groups?how dense, everywhere?significantly high/low?nearest neighbour indexDBSCAN · k-meanskernel densityGetis-Ord Gi*supportingVoronoiDelaunayrandom samplescentre and spread

What this guide covers

Before any statistic: three checks

Most misleading spatial statistics fail before the first calculation, at one of three points.

The CRS. Distances, areas and densities must be measured in a projected coordinate system with metre units. Averages of longitudes, nearest distances in degrees and kernel radii of "0.01" are meaningless, and they produce numbers that look plausible enough to go unnoticed. Choosing a projected CRS for analysis covers how to pick one.

Duplicates. Two records of the same event at the same location make any clustering method see a cluster that is not there, and nearest neighbour distances of zero drag the index down. Finding duplicates first is cheap insurance.

The population at risk. Burglaries cluster where houses are; disease cases cluster where people live. A pattern of events that simply mirrors the underlying population is not interesting, and many "hot spots" are population maps in disguise. Wherever possible, analyse rates — events per household, per resident, per kilometre of road — rather than raw counts, or compare the event pattern against the population pattern explicitly.

from qgis.core import QgsProject, QgsWkbTypes

layer = QgsProject.instance().mapLayersByName("incidents")[0]
crs = layer.crs()
print("CRS:", crs.authid(), "| geographic:", crs.isGeographic(),
      "| units:", crs.mapUnits())
print("geometry:", QgsWkbTypes.displayString(layer.wkbType()),
      "| features:", layer.featureCount())

Breakdown: isGeographic() is the one-line guard every statistics script should start with; raising an error when it is true prevents an entire class of silent mistakes. Checking the geometry type catches multipoint layers, which some algorithms treat as single features with many locations.

Global tests: is there a pattern at all?

The first question about a point pattern is whether it differs from randomness. Global tests answer it with one number for the whole study area.

The average nearest neighbour index compares the observed mean distance from each point to its nearest neighbour with the distance expected under complete spatial randomness. Values below 1 indicate clustering, above 1 dispersion. It is fast and built into Processing, but it depends heavily on the study area used to compute the expected distance and it looks at a single scale — the very nearest neighbour.

import processing

r = processing.run("qgis:nearestneighbouranalysis", {
    "INPUT": layer, "OUTPUT_HTML_FILE": "TEMPORARY_OUTPUT"})
print(f"index {r['NN_INDEX']:.2f}, z {r['Z_SCORE']:.1f}, n {r['POINT_COUNT']}")

A Monte Carlo test is the more flexible alternative: place the same number of points at random inside the real study area many times, compute the same statistic each time, and see where the observed value falls among the simulations. It accounts automatically for irregular boundaries and edge effects, and it works for any statistic you can compute. The random points recipe builds one with seeded native:randompointsinpolygons.

Global tests say whether, not where. A significant result is a reason to look further with local methods, not an answer in itself.

Local methods: where is the pattern?

Local methods produce a result for every point or every area, which is what turns a statistic into a map.

Three local views of the same eventsThe same events seen three ways. Clustering assigns each point a group label or marks it as noise. Kernel density produces a continuous raster surface of intensity. Hot spot analysis assigns each grid cell a z-score and a significance class. Clustering answers which points belong together, density answers how intense, and hot spots answer where intensity is unusually high or low.Labels, surface, significanceclusteringlabel per pointnoise allowed (DBSCAN)which belong together?kernel densityvalue per pixelsmooth surfacehow intense, where?hot spotsz-score per cellhot · cold · noneunusually high or low?

Clustering assigns points to groups. DBSCAN grows groups from dense neighbourhoods defined by a distance and a minimum count, and leaves sparse points as noise — ideal for discovering where events concentrate, in shapes of any kind. K-means divides all points into a chosen number of compact groups — ideal for allocation problems such as dividing sites among teams. Parameters matter enormously: the clustering recipe chooses DBSCAN's distance from the k-distance curve and k-means's group count from the spread curve.

Kernel density spreads each event over a smooth hill and sums the hills into a raster. It is the most readable picture of intensity and the easiest to combine with other rasters, but it has no notion of significance — every surface has peaks. Its one crucial parameter is the radius, which sets the scale of the picture.

Hot spot analysis with Getis-Ord Gi* aggregates events to a grid and asks, for each cell, whether the sum over its neighbourhood is higher or lower than chance would produce. It is the method to use when a decision depends on a location being unusually high — targeting patrols, prioritising inspections — because it comes with a z-score and, after correction for multiple tests, a defensible significance level.

These three are complementary rather than competing. A common workflow runs density to see the shape of the pattern, hot spots to decide which concentrations are real, and clustering to extract the events inside them for further work.

Scale is part of the answer

Every local method has a scale parameter: DBSCAN's distance, the kernel radius, the hot spot neighbourhood band, the grid cell size. Change it and the result changes — sometimes a little, sometimes completely. This is not a flaw to be engineered away; patterns genuinely exist at some scales and not others. Shops cluster at the scale of a street and are evenly spread at the scale of a city.

Same data, three scalesResults for the same events at small, medium and large scale parameters. At a small scale, many small tight groups appear, including noise. At a medium scale, a few meaningful concentrations emerge. At a large scale, everything merges into one broad area. Results that persist across a range of scales are the robust findings.Report what survives a change of scalesmall scalemany tiny groupsnoise looks realmedium scalea few concentrationsoften the useful onelarge scaleone broad areastructure lostrerun at 2–3 scales; trust what persists

The practical rule is to run every local method at two or three scales around a sensible starting value, and to treat as findings only the patterns that persist. Starting values come from the data: the median nearest neighbour distance, a few times it for kernel radii, two or three grid cells for hot spot bands. Report the scale with every map, because a reader cannot interpret a hot spot without knowing how big a neighbourhood it summarises.

Structures that support the analysis

Two geometric structures and one technique turn up throughout spatial statistics.

Voronoi polygons give each point the area closer to it than to any other. They allocate locations to their nearest site in a single join, they provide the weights for area-weighted averages of point measurements (the Thiessen method in hydrology), and they show at a glance how evenly a network of sites covers a region.

Delaunay triangulations connect points whose Voronoi cells touch. Their edges are a natural definition of "neighbour" for spatial weights, their triangles are the basis of TIN interpolation, and unusually long edges reveal empty areas between groups of points. The Voronoi and Delaunay recipe builds both and turns edges into a neighbour list.

Random generation supplies the baselines tests compare against and the samples fieldwork depends on. Seeds make both reproducible, which matters as much in analysis as in any other code: a hot spot map that changes every time someone reruns the script cannot support a decision.

Summarising location and spread

Sometimes the question is not about clusters at all but about the distribution as a whole: where is it centred, how spread out is it, is it stretched in a direction, and how have those things changed? Centrographic statistics answer this with a handful of numbers — mean centre, median centre, standard distance, the axes and angle of a deviational ellipse — that are easy to compare between groups and over time.

import numpy as np

pts = np.array([[f.geometry().asPoint().x(), f.geometry().asPoint().y()]
                for f in layer.getFeatures()])
centre = pts.mean(axis=0)
sd = np.sqrt(((pts - centre) ** 2).sum(axis=1).mean())
print("mean centre", centre.round(), "standard distance", round(float(sd)), "m")

Breakdown: Two lines of NumPy give the centre and the spread; the centrography recipe adds weights, the outlier-resistant median centre and the deviational ellipse, and compares them across groups and years. These summaries are especially useful in reports, where a single arrow showing how a population centre moved over a decade communicates more than any table.

Where QGIS stops and Python libraries start

QGIS's Processing toolbox covers nearest neighbour tests, DBSCAN, k-means, kernel density, Voronoi, Delaunay, random points and mean coordinates. It does not include local spatial autocorrelation statistics such as Getis-Ord Gi* or Local Moran's I, Ripley's K function, or HDBSCAN. The recipes here fill those gaps with short NumPy implementations, which have the advantage of being transparent and dependency-free.

For heavier work, the scientific Python ecosystem is the natural extension: scikit-learn for clustering (including HDBSCAN), esda and libpysal from the PySAL project for spatial autocorrelation and weights, pointpats for point pattern statistics. They install into the Python environment QGIS uses, as described in installing Python packages into QGIS, and exchange data with QGIS through NumPy arrays or GeoDataFrames — the bridge covered in PyQGIS and the Python data stack.

The division of labour that works well: QGIS for data preparation, aggregation, grids and maps; NumPy or a library for the statistic; QGIS again for writing results back to layers and styling them honestly.

Mapping results honestly

A statistical result is only as good as the map that communicates it. A few conventions prevent the most common distortions.

  • Show significance as categories — hot, cold, not significant — rather than a continuous ramp of z-scores that invites readers to see hot spots everywhere.
  • Use diverging colours (red and blue) only for results that genuinely have a meaningful midpoint, such as z-scores; use sequential ramps for densities and counts.
  • Stretch density surfaces to a high percentile, not the maximum, so a single extreme pixel does not flatten the rest.
  • State the scale parameter, the study area and the date range on or beside the map.
  • Where results depend on the population at risk, say so — or map rates.

The graduated and categorized renderer guide covers the styling side in detail.

A worked sequence: from events to a decision

The methods combine naturally into a sequence that fits most event analyses. Suppose a city's road safety team wants to choose ten locations for traffic calming from five years of collision records.

  1. Prepare. Reproject collisions to the national grid, remove duplicate records, and filter to the injury collisions that matter for the decision.
  2. Test globally. A nearest neighbour index of 0.55 with a large negative z-score confirms strong clustering; it is worth looking for where.
  3. See the shape. Kernel density at 150 m and 400 m radii shows concentrations at junctions and along two arterial roads.
  4. Test locally. Gi* on a 100 m hexagon grid with a 300 m band, corrected for multiple tests, marks 41 cells as significant hot spots; 33 of them persist at a 450 m band.
  5. Extract and rank. The collisions inside persistent hot spot cells are grouped with DBSCAN at the junction scale, and each group's convex hull and severity-weighted count go into a ranked table.
  6. Report. The map shows hot, cold and not significant cells, the ten chosen sites, and the parameters used.
steps = {
    "density_150": ("qgis:heatmapkerneldensityestimation", {"RADIUS": 150, "PIXEL_SIZE": 15}),
    "density_400": ("qgis:heatmapkerneldensityestimation", {"RADIUS": 400, "PIXEL_SIZE": 40}),
}
for name, (alg, extra) in steps.items():
    processing.run(alg, {"INPUT": layer, "OUTPUT_VALUE": 1, "KERNEL": 0,
                         "OUTPUT": f"/data/results/{name}.tif", **extra})

Breakdown: Writing the parameter sets as data rather than repeating calls makes the multi-scale runs explicit and easy to extend. Each step in the sequence uses a recipe from this guide unchanged; what makes the result defensible is the order — test before mapping, map before deciding — and the scale check in step four. The same sequence works for crime, disease, wildlife sightings or service requests.

Common pitfalls

  • Analysing counts where rates are needed, so hot spots simply mark where people live.
  • Using the bounding box as the study area, which inflates clustering for points that can only occur in part of it.
  • Reading a significant z-score as a strong effect when the index is close to 1 and the sample is huge.
  • Mapping uncorrected significance, so dozens of cells are "hot" by chance.
  • Choosing parameters after looking at results until the map shows what was expected. Fix parameters from the data and the question first.

Key takeaways

  • Work in a projected CRS, remove duplicates, and think about the population at risk before computing anything.
  • Global tests say whether a pattern exists; local methods say where. Use both.
  • Clustering labels points, kernel density draws intensity, and Gi* tests significance — they answer different questions.
  • Every local method has a scale parameter; run at several scales and report what persists.
  • Voronoi cells, Delaunay neighbours and seeded random points are the supporting tools for allocation, weights and baselines.
  • Centre, spread and orientation summarise whole distributions and compare well across groups and time.
  • Map significance as categories and state the parameters, so readers can judge the result.

Frequently Asked Questions

Do I need specialist software for spatial statistics? Not for most work. QGIS plus NumPy covers the common methods; PySAL and scikit-learn extend them inside the same Python environment.

Should I analyse points or aggregate them to areas first? Point methods — nearest neighbour, DBSCAN, kernel density — keep full detail. Area methods — hot spots, rates — need aggregation, which is also how you bring in population denominators. Many analyses use both.

How many points do these methods need? Dozens at minimum for global tests, and enough that most grid cells or neighbourhoods contain several events for local ones. With very few events, describe them rather than test them.

Are results comparable between cities or years? Only if the method, parameters, CRS and study area definitions are identical. Record all of them with every run.

Can I automate these analyses on a schedule? Yes. Every recipe here runs headless; combine them with batch processing to produce monthly hot spot maps or density surfaces automatically.