Digitising, topology and georeferencing

Twelve wards drawn by hand, with four slips: a ring whose corners were clicked out of order, a 3 m gap, a 5 m overlap and a corner clicked twice. Their areas added up to 44,004,000 m² for a town of 48,000,000 m², and QGIS’s gap and overlap checks refused to run at all until the first slip was fixed.

You will digitise a layer so that slips are rare, find the ones that happen, and georeference a scanned map, judging the fit by its residuals. Every count and residual below is what QGIS 4.2.3 computed for the town and a made-up scanned sheet.

Digitising: draw against the neighbours, not beside them

QGIS digitising turns a picture or a survey into features you draw. Start with Layer ▸ Create Layer ▸ New GeoPackage Layer, which asks for the file, a table name, the geometry type, the CRS and the fields (creating layers). Then three settings stop most slips before they happen:

Snappingsnapping: A digitising setting that moves a click onto a nearby vertex or edge of existing features, within a tolerance in screen pixels or in map units, so that shared boundaries match exactly. lands each click exactly on a nearby vertex or edge, tracingtracing: A digitising tool that follows the boundary of existing features between two clicks, adding every vertex along the way, so a shared edge is copied instead of redrawn. follows a neighbour’s boundary between two clicks, and Avoid Overlap on Active Layer cuts any part of a new polygon that would cover an existing one (editing). Topological editing moves shared vertices together, so two wards keep one boundary when you edit it. Digitising with these settings on still needs a check afterwards, as the next section shows.

The attributes need rules too, above all when several people digitise. A form’s widgets, such as a value map, offer only the allowed values; constraints refuse or flag a bad value; and a GeoPackage can store field domainsfield domain: A rule for a field's allowed values stored in the data file itself, such as a list of coded values or a range. GeoPackage supports domains, and QGIS builds form widgets from them., lists or ranges of allowed values, from which QGIS builds the form’s widgets.

Topology errors hide in totals

Topologytopology (in GIS): How features relate to their neighbours, such as polygons that should share boundaries without gaps or overlaps, and lines that should meet at their ends. is how features relate to their neighbours: wards that should tile a town without gaps or overlaps. The edited wards carry four slips:

QGIS finds them with two kinds of check. Check validity tests each feature on its own; the geometry checks of Processing compare neighbours, for gaps, overlaps, repeated vertices and more, and the Topology Checker plugin offers rules such as “Must not have gaps” and “Must not overlap” (Topology Checker). This program sums the areas and runs both kinds:

"""Checking the edited wards with QGIS: their total area, Check validity, and the gap and overlap checks."""
import qgis_start
from edited_wards import edited_rings, layer_of
from qgis.core import QgsProcessingException

processing = qgis_start.processing()
wards = layer_of(edited_rings())
total = sum(f.geometry().area() for f in wards.getFeatures())
print(f"Sum of the 12 ward areas: {total:,.0f} m² (the town is 8 km by 6 km: 48,000,000 m²)")
result = processing.run("native:checkvalidity", {"INPUT_LAYER": wards, "METHOD": 2, "VALID_OUTPUT": "memory:",
                                                 "INVALID_OUTPUT": "memory:", "ERROR_OUTPUT": "memory:"})
print(f"Check validity (GEOS): {result['VALID_COUNT']} valid, {result['INVALID_COUNT']} invalid")
for error in result["ERROR_OUTPUT"].getFeatures():
    print(f"  {error['message']} at {error.geometry().asPoint().x():,.0f} E {error.geometry().asPoint().y():,.0f} N")
for algorithm in ["native:checkgeometrygap", "native:checkgeometryoverlap"]:
    try:
        processing.run(algorithm, {"INPUT": wards, "UNIQUE_ID": "ward", "ERRORS": "memory:", "OUTPUT": "memory:",
                                   "NEIGHBORS": "memory:"} if algorithm.endswith("gap") else
                       {"INPUT": wards, "UNIQUE_ID": "ward", "ERRORS": "memory:", "OUTPUT": "memory:"})
        print(f"{algorithm}: ran")
    except QgsProcessingException as e:
        print(f"{algorithm}: refused: {str(e).split('. ')[0]}.")

Sum of the 12 ward areas: 44,004,000 m² (the town is 8 km by 6 km: 48,000,000 m²)
Check validity (GEOS): 11 valid, 1 invalid
  Self-intersection at 377,000 E 2,049,000 N
native:checkgeometrygap: refused: Feature (3) from “wards_edited” has invalid geometry.
native:checkgeometryoverlap: refused: Feature (3) from “wards_edited” has invalid geometry.

The total was 3,996,000 m² short, almost all of it W03: a bow-tie’s two halves cancel, so a 4 km² ward counted as zero, the invalid geometryvalid geometry: A polygon or multipolygon that follows the simple-features rules: simple rings that do not cross, touch only at points and keep the holes inside the shell, and an interior in one piece. Areas and overlays are reliable only on valid geometry. More in “Geometry types, WKT, WKB and validity” Check validity reported at its crossing. The gap and the overlap could not even be looked for, because Processing refuses invalid input by default. Fix validity first, then topology.

Georeferencing: control points and their residuals

Georeferencing gives a scan coordinates. You pick ground control pointsground control point (GCP): A point identified both on an image or scan and on the ground, with known map coordinates. Georeferencing fits a transformation from image positions to map positions through such points., spots whose map coordinates you know, and QGIS’s georeferencer fits a transformation from the scan’s pixels to the map (georeferencer):

Georeferencing is judged by its residualsresidual: For one control point, the distance between where a fitted transformation puts it and where it really is. The RMSE of the residuals sums up how well the fit matches its points.: how far the fit puts each control point from where it belongs, summed up by the RMSERMSE (root mean square error): The square root of the mean of the squared errors. For positions, the typical size of the distance between where things are recorded and where they really are. More in “Vector and raster data models”. This program fits the lesson’s made-up sheet, 1:25,000 at 300 dots per inch, with Polynomial 1 and with Helmert:

"""The residual of each control point after a Polynomial 1 fit of the scanned sheet, and the RMSE of two fits."""
import math

import qgis_start  # noqa: F401  (starts QGIS without a window)
from qgis.analysis import QgsGcpTransformerInterface as Transformer
from qgis.core import QgsPointXY
from scan import CONTROL, PIXEL, picked

scan_points = picked(CONTROL, seed=7)
for method in [Transformer.TransformMethod.PolynomialOrder1, Transformer.TransformMethod.Helmert]:
    transformer = Transformer.create(method)
    transformer.updateParametersFromGcps([QgsPointXY(*p) for p in scan_points], [QgsPointXY(*p) for p in CONTROL], False)
    residuals = []
    for (x, y), (e, n) in zip(scan_points, CONTROL):
        ok, fitted_e, fitted_n = transformer.transform(x, y)
        residuals.append((fitted_e - e, fitted_n - n))
    rmse = math.sqrt(sum(dx * dx + dy * dy for dx, dy in residuals) / len(residuals))
    print(f"{Transformer.methodToString(method)}: RMSE {rmse:.2f} m = {rmse / PIXEL:.2f} pixels")
    if method == Transformer.TransformMethod.PolynomialOrder1:
        for i, ((e, n), (dx, dy)) in enumerate(zip(CONTROL, residuals), 1):
            print(f"  point {i:>2} at {e:,} E {n:,} N: dE {dx:+.2f} m, dN {dy:+.2f} m, residual {math.hypot(dx, dy):.2f} m")

Polynomial 1: RMSE 1.80 m = 0.85 pixels
  point  1 at 372,000 E 2,050,000 N: dE +1.23 m, dN +0.57 m, residual 1.36 m
  point  2 at 380,000 E 2,050,000 N: dE -1.17 m, dN -0.80 m, residual 1.42 m
  point  3 at 372,000 E 2,044,000 N: dE +0.08 m, dN -1.35 m, residual 1.35 m
  point  4 at 380,000 E 2,044,000 N: dE +1.13 m, dN -0.24 m, residual 1.16 m
  point  5 at 376,000 E 2,050,000 N: dE +2.21 m, dN +0.56 m, residual 2.28 m
  point  6 at 376,000 E 2,044,000 N: dE +1.12 m, dN +2.51 m, residual 2.75 m
  point  7 at 372,000 E 2,047,000 N: dE -1.15 m, dN +0.84 m, residual 1.42 m
  point  8 at 380,000 E 2,047,000 N: dE +0.10 m, dN +0.23 m, residual 0.25 m
  point  9 at 374,000 E 2,046,000 N: dE -1.86 m, dN -2.03 m, residual 2.76 m
  point 10 at 378,000 E 2,048,000 N: dE -1.69 m, dN -0.29 m, residual 1.71 m
Helmert: RMSE 4.88 m = 2.31 pixels

Polynomial 1 fitted the ten points with an RMSE of 1.80 m, 0.85 of a pixel; Helmert, which cannot follow paper that shrank unevenly, 4.88 m. Check points, measured like control points but left out of the fit, judge a fit more fairly. Practise on maps you may use: Survey of India lists its Open Series Maps as PDF files at 1:50,000, at no charge, but the listing states no terms for passing them on (Survey of India), so this lesson uses a sheet of its own.

Repairing the slips, with QGIS's counts before and after each fix Pro: opens with a Pro, Premium or Ultimate pass

Judging a georeference by check points, with five transformations measured Pro: opens with a Pro, Premium or Ultimate pass

Forms that hold: field domains and constraints tested on six records Pro: opens with a Pro, Premium or Ultimate pass

Practise: residuals by hand

Report the residuals of a fit: each point’s distance from its true position, the RMSE in metres and in pixels, and the worst point. The sample tests use the lesson’s own Polynomial 1 fit:

import math


def residual_report(fitted, known, pixel):
    """Each control point's residual in metres, the RMSE in metres and in pixels, and the worst point (from 1)."""
    # Replace these lines: compute every residual first, then the RMSE and the worst point.
    return {"residuals": [], "rmse": 0.0, "rmse_pixels": 0.0, "worst": 1}

from residuals import residual_report

PIXEL = 25_000 * 0.0254 / 300  # 2.1167 m on the ground at 1:25,000 and 300 dots per inch
KNOWN = [(372_000, 2_050_000), (380_000, 2_050_000), (372_000, 2_044_000), (380_000, 2_044_000),
         (376_000, 2_050_000), (376_000, 2_044_000), (372_000, 2_047_000), (380_000, 2_047_000),
         (374_000, 2_046_000), (378_000, 2_048_000)]
SHIFTS = [(1.23, 0.57), (-1.17, -0.80), (0.08, -1.35), (1.13, -0.24), (2.21, 0.56), (1.12, 2.51), (-1.15, 0.84),
          (0.10, 0.23), (-1.86, -2.03), (-1.69, -0.29)]  # dE, dN of the lesson's Polynomial 1 fit, in metres
FITTED = [(e + de, n + dn) for (e, n), (de, dn) in zip(KNOWN, SHIFTS)]


def test_the_lessons_fit():
    """reports the lesson's Polynomial 1 fit: RMSE 1.8 m, 0.85 pixels, worst at point 9"""
    report = residual_report(FITTED, KNOWN, PIXEL)
    assert report["residuals"] == [1.36, 1.42, 1.35, 1.16, 2.28, 2.75, 1.42, 0.25, 2.75, 1.71]
    assert report["rmse"] == 1.8
    assert report["rmse_pixels"] == 0.85
    assert report["worst"] == 9


def test_two_points():
    """squares before averaging: residuals of 3 m and 4 m give an RMSE of 3.54 m, not 3.5"""
    report = residual_report([(3, 0), (0, 4)], [(0, 0), (0, 0)], 2.0)
    assert report == {"residuals": [3.0, 4.0], "rmse": 3.54, "rmse_pixels": 1.77, "worst": 2}


def test_an_exact_fit():
    """gives zero for a fit that passes through its point"""
    assert residual_report([(5.0, 7.0)], [(5.0, 7.0)], 1.0) == {"residuals": [0.0], "rmse": 0.0, "rmse_pixels": 0.0, "worst": 1}

Then test yourself on digitising, topology and residuals:

Interview questions

Warm-up (fresher): What is the difference between an invalid geometry and a topology error in QGIS?

An invalid geometry breaks the rules for one feature on its own, such as a ring that crosses itself; Check validity finds it. A topology error is about neighbours, such as a gap or an overlap between two valid wards, so only checks that compare neighbours find it. Fix validity first, because Processing’s topology checks refuse invalid input by default.

Three questions of this module’s bank start from this lesson:

Key takeaways

• QGIS digitising with snapping, tracing and Avoid Overlap prevents most slips; a check afterwards finds the rest.
• Topology errors hide in totals: one bow-tie made 4 km² vanish, and Processing refused the gap and overlap checks until validity was fixed.
• Georeferencing fits a transformation to control points, and its residuals and RMSE say how well it fits them: 1.80 m for Polynomial 1 on the lesson’s sheet.