Calculate Mean Centre and Standard Distance in PyQGIS

Descriptive statistics for numbers — mean, median, standard deviation — have direct spatial equivalents. The mean centre is the average location of a set of points. The standard distance is how far, typically, points lie from that centre. The standard deviational ellipse adds direction: is the spread stretched along a valley, a coastline or a motorway? Together they compress thousands of points into a few numbers and shapes that are easy to compare between groups and between years — how the centre of a city's population has drifted, or whether one species ranges more widely than another.

This recipe belongs to Spatial Statistics & Pattern Analysis. It computes weighted and unweighted mean centres with Processing, the median centre and standard distance with NumPy, builds a standard deviational ellipse, and compares these summaries across groups and time.

Centre, circle and ellipseA scatter of points elongated from lower left to upper right. The mean centre is marked at the average location. A standard distance circle around it has a radius equal to the root mean square distance of points from the centre. A standard deviational ellipse is also drawn, rotated along the direction of greatest spread and narrower across it, capturing the elongation that the circle misses.Where, how spread, and in which directionmean centreaverage locationstandard distancecircle, no directiondeviational ellipsespread + orientation

Prerequisites

  • QGIS 3.34 LTR or newer, or the QGIS 4 series. NumPy ships with QGIS.
  • A point layer in a projected CRS. A mean of latitudes and longitudes is not a meaningful centre over large areas, and distances in degrees are meaningless.
  • Polygon data converted to representative points first — centroids or points on surface — if you want to summarise polygons.

Mean centre with Processing

The mean coordinates algorithm computes the average x and y of a layer, optionally weighted by a field and optionally grouped by a category.

import processing
from qgis.core import QgsProject

households = QgsProject.instance().mapLayersByName("households")[0]

centre = processing.run("native:meancoordinates", {
    "INPUT": households, "WEIGHT": None, "UID": None,
    "OUTPUT": "memory:mean_centre"})["OUTPUT"]
weighted = processing.run("native:meancoordinates", {
    "INPUT": households, "WEIGHT": "persons", "UID": None,
    "OUTPUT": "memory:mean_centre_persons"})["OUTPUT"]

for layer in (centre, weighted):
    f = next(layer.getFeatures())
    print(layer.name(), round(f["MEAN_X"]), round(f["MEAN_Y"]))
QgsProject.instance().addMapLayers([centre, weighted])

Breakdown: The output is a point layer with MEAN_X and MEAN_Y fields as well as the geometry. Weighting by persons per household moves the centre towards large households — the population centre rather than the dwelling centre. With a UID field, the algorithm produces one centre per group, which is the basis for the comparisons later on. The mean centre is sensitive to outliers: a few distant points pull it noticeably, which is one reason to compute the median centre as well.

Median centre

The median centre is the location that minimises the total straight-line distance to all points. It resists outliers in the way the ordinary median does, and has a practical meaning: it is the best single site for a facility serving everyone, if travel is by straight line.

Mean versus median centre with outliersA compact group of points and three distant outliers to the right. The mean centre is pulled noticeably towards the outliers. The median centre, which minimises total distance, stays within the main group. Weiszfeld's algorithm finds the median centre by repeatedly re-weighting points by the inverse of their distance to the current estimate.Outliers move the mean, not the medianmedian centremean centreoutliersWeiszfeld iteration: re-weight by 1 / distance until it settles

import numpy as np

def median_centre(xy, weights=None, tol=0.01, max_iter=200):
    w = np.ones(len(xy)) if weights is None else np.asarray(weights, float)
    est = np.average(xy, axis=0, weights=w)
    for _ in range(max_iter):
        d = np.linalg.norm(xy - est, axis=1)
        d = np.where(d < 1e-9, 1e-9, d)
        k = w / d
        new = (xy * k[:, None]).sum(axis=0) / k.sum()
        if np.linalg.norm(new - est) < tol:
            return new
        est = new
    return est

xy = np.array([[f.geometry().asPoint().x(), f.geometry().asPoint().y()]
               for f in households.getFeatures()])
persons = np.array([f["persons"] or 0 for f in households.getFeatures()], float)
med = median_centre(xy, persons)
mean = np.average(xy, axis=0, weights=persons)
print("median", med.round(), "mean", mean.round(), "gap", round(float(np.linalg.norm(med - mean))), "m")

Breakdown: Weiszfeld's algorithm starts from the mean and repeatedly recomputes a weighted average in which each point's weight is divided by its distance from the current estimate, so nearby points count more each round. It converges quickly for typical data; the tolerance is in metres. The guard against zero distance prevents a division by zero when the estimate lands exactly on a point. A large gap between mean and median centre is itself a finding: it means the distribution is lopsided, with a long tail in the direction of the mean.

Standard distance

The standard distance is the spatial standard deviation: the root of the mean squared distance from the mean centre. Drawn as a circle, it encloses roughly two thirds of the points for a compact, round distribution.

from qgis.core import QgsGeometry, QgsPointXY, QgsVectorLayer, QgsFeature

d2 = ((xy - mean) ** 2).sum(axis=1)
sd = float(np.sqrt(np.average(d2, weights=persons)))
print(f"standard distance: {sd:.0f} m")

circle = QgsVectorLayer(f"Polygon?crs={households.crs().authid()}&field=sd_m:double",
                        "standard distance", "memory")
f = QgsFeature(circle.fields())
f.setAttributes([sd])
f.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(*mean)).buffer(sd, 64))
circle.dataProvider().addFeatures([f])
QgsProject.instance().addMapLayer(circle)

Breakdown: Weighting the squared distances by the same weights as the centre keeps the two statistics consistent. Buffering the centre point by the standard distance with 64 segments produces a smooth circle. The number is directly comparable between groups measured in the same CRS: a standard distance of 2.1 km for one group and 4.8 km for another means the second is more than twice as spread out. It ignores direction, which the ellipse addresses.

Standard deviational ellipse

Spread is often directional — settlement along a river, incidents along a road corridor. The standard deviational ellipse finds the axis of greatest spread and measures the spread along and across it.

Ellipse from the covariance matrixThe ellipse comes from the covariance matrix of the point coordinates around the mean centre. Its eigenvectors give the directions of greatest and least spread, and the square roots of its eigenvalues give the standard deviations along each. The ratio of the axes measures elongation; the angle of the major axis gives orientation.Covariance → axes and anglecovariancevar(x) cov(x,y)cov(x,y) var(y)eigen-decompositionvectors → directionsvalues → spreadsellipsemajor, minor axisangle, elongation√2 × axis lengths is the common convention for display

import math

def deviational_ellipse(xy, w, centre, scale=math.sqrt(2)):
    d = xy - centre
    cov = np.cov(d.T, aweights=w, bias=True)
    vals, vecs = np.linalg.eigh(cov)          # ascending
    minor, major = np.sqrt(vals) * scale
    angle = math.degrees(math.atan2(vecs[1, 1], vecs[0, 1]))
    return major, minor, angle

major, minor, angle = deviational_ellipse(xy, persons, mean)
print(f"major {major:.0f} m, minor {minor:.0f} m, angle {angle:.0f}°, "
      f"elongation {major / minor:.2f}")

circle_pts = [QgsPointXY(math.cos(t) * major, math.sin(t) * minor)
              for t in np.linspace(0, 2 * math.pi, 73)]
ellipse = QgsGeometry.fromPolygonXY([circle_pts])
ellipse.rotate(-angle, QgsPointXY(0, 0))
ellipse.translate(mean[0], mean[1])

Breakdown: The weighted covariance matrix of coordinates around the centre holds everything needed: its eigenvectors are the ellipse's axes and the square roots of its eigenvalues are the standard deviations along them. Multiplying by √2 is a common display convention that makes the ellipse enclose a larger share of points; whatever you choose, use the same factor for every group you compare. The angle is measured anticlockwise from the x-axis; QgsGeometry.rotate rotates clockwise, hence the sign. An elongation near 1 means no preferred direction; values above 2 indicate strong directionality.

Compare groups and years

The summaries become informative when compared: the population centre in 2015 and 2025, the spread of two species, the orientation of incidents by time of day.

from collections import defaultdict

groups = defaultdict(list)
for f in households.getFeatures():
    p = f.geometry().asPoint()
    groups[f["census_year"]].append((p.x(), p.y(), f["persons"] or 0))

summary = {}
for year, rows in sorted(groups.items()):
    a = np.array(rows)
    pts, wts = a[:, :2], a[:, 2]
    c = np.average(pts, axis=0, weights=wts)
    sdist = math.sqrt(np.average(((pts - c) ** 2).sum(axis=1), weights=wts))
    summary[year] = (c, sdist)
    print(year, c.round(), round(sdist))

years = sorted(summary)
for a, b in zip(years, years[1:]):
    shift = np.linalg.norm(summary[b][0] - summary[a][0])
    print(f"{a}→{b}: centre moved {shift:.0f} m, spread changed "
          f"{summary[b][1] - summary[a][1]:+.0f} m")

Breakdown: Computing each group's centre and standard distance with the same weights and CRS makes them directly comparable. The movement of the centre between periods — its distance and, with atan2, its direction — is a compact way to describe drift such as suburban growth. Drawing successive centres as points joined in order with native:pointstopath produces the familiar "centre of population" track. Small groups give noisy centres; report the count with each.

Draw the centre's track over time

The comparison above prints numbers; a map of the centre moving from year to year communicates the same result instantly. Writing each period's centre as a point with its year, then joining the points in order, produces the classic "centre of population" track.

from qgis.core import QgsVectorLayer, QgsFeature, QgsGeometry, QgsPointXY

track = QgsVectorLayer(f"Point?crs={households.crs().authid()}&field=year:integer"
                       "&field=sd_m:double", "population centre by year", "memory")
rows = []
for year in years:
    c, sdist = summary[year]
    f = QgsFeature(track.fields())
    f.setAttributes([int(year), round(sdist, 1)])
    f.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(*c)))
    rows.append(f)
track.dataProvider().addFeatures(rows)

path = processing.run("native:pointstopath", {
    "INPUT": track, "ORDER_EXPRESSION": '"year"', "GROUP_EXPRESSION": "",
    "OUTPUT": "memory:centre_path"})["OUTPUT"]
QgsProject.instance().addMapLayers([track, path])

Breakdown: Storing the standard distance with each centre lets you scale the point symbols by spread, so the map shows both where the centre went and how the distribution widened or narrowed. native:pointstopath orders points by the expression and joins them into one line; with an arrow symbol layer on the line, as in arrow symbol layers, direction is clear at a glance. Label points with the year and keep the basemap quiet so the track stands out.

QGIS version compatibility

native:meancoordinates has the same parameters on QGIS 3.34 LTR, 3.40 LTR and QGIS 4 (earlier 3.x releases used qgis:meancoordinates). The NumPy functions used — average, cov with aweights, linalg.eigh — are available in every NumPy version shipped with current QGIS. QgsGeometry.rotate and translate work in place on all releases.

Troubleshooting

  • The centre is in the sea. The distribution is crescent-shaped or bimodal; the mean centre need not lie inside it. Report the median centre, or analyse the parts separately.
  • Distances are tiny decimals. The layer is geographic; reproject first.
  • The ellipse points the wrong way. The rotation sign is reversed; QgsGeometry.rotate uses clockwise degrees.
  • Weighted results equal unweighted. The weight field is mostly NULL, which becomes 0 or 1 depending on how you read it; check it first.

Conclusion

Compute the mean centre with native:meancoordinates, add the median centre for robustness to outliers, measure spread with the standard distance and its direction with a covariance-based ellipse, keep weights and CRS identical across groups, and compare groups and years by centre movement and change in spread.

Frequently Asked Questions

Can I compute these for polygons? Use centroids or points on surface, weighted by area or population, as inputs.

What share of points does the standard distance circle contain? About 63–68 % for a roughly circular normal distribution; less for skewed or clustered data.

Is the median centre the best facility location? Only for straight-line distance. For road travel, use network analysis such as an origin–destination matrix.

Does QGIS have a deviational ellipse algorithm? Not in core Processing; the NumPy approach above is short and transparent.