GIS & Geospatial Data (QGIS, GeoPandas, PostGIS) Module 3 – Vector geometry and file formats
Geometry types, WKT, WKB and validity
Write geometries as WKT, WKB and EWKT, find why a polygon is invalid and repair it without losing area, orient its rings, and read a DE-9IM matrix.
What you will learn
- Write the simple-features geometry types as WKT, WKB and PostGIS EWKT, and read the bytes of a WKB point
- Find why a polygon is invalid and repair it while checking that its area survives
- Set the direction of polygon rings for RFC 7946 or a shapefile, and explain why the name right-hand rule is ambiguous
- Describe how two geometries relate with a DE-9IM matrix and the named predicate it matches
Before you start
What you’ll learn · Before you start
What you will learn
- Write the simple-features geometry types as WKT, WKB and PostGIS EWKT, and read the bytes of a WKB point
- Find why a polygon is invalid and repair it while checking that its area survives
- Set the direction of polygon rings for RFC 7946 or a shapefile, and explain why the name right-hand rule is ambiguous
- Describe how two geometries relate with a DE-9IM matrix and the named predicate it matches
Before you start
On this page
A parcel whose ring crosses itself reports an area of 0 m² although it covers 2 m², and any total, overlay or join built on it is wrong, or stops with an error, long after the parcel was drawn. Most such shapes look fine on a map; only the rules of the geometry model tell you they are broken.
This lesson goes one level below the map. You write the seven geometry types as text and as bytes, read a point’s 21 bytes of WKB by hand, check polygons against the five validity rules and repair them without losing area, make rings run the way GeoJSON or a shapefile wants, and describe exactly how two shapes meet with a nine-character matrix. Each claim is shown twice: by code written from the rules, which you can run in your browser, and by GEOS, the engine inside Shapely, PostGIS and QGIS.
Seven types describe almost every feature
The geometry model of almost every GIS is the simple features: The OGC and ISO model of vector geometry: points, line strings and polygons with straight edges between positions, their multi-part forms and the geometry collection, each optionally with Z and M values. model of the OGC’s Simple Feature Access standard: points, lines and polygons with straight edges between their positions. The track’s lesson on vector and raster models used three types and their multi-part forms; the full list has seven:
Every simple-features geometry is one of seven types
Text description of the diagram
Seven small panels in two columns, each with a drawing and the type's name.
Point: a single dot. LineString: a line of three straight segments joining four positions. Polygon: a square area with a smaller square hole inside it, drawn as two rings. MultiPoint: four separate dots. MultiLineString: two separate lines. MultiPolygon: two separate squares. GeometryCollection: a dot, a line and a small square together, parts of different types in one value.
The eighth place holds a note: any of the types may add Z, a height, or M, a measure, as in POINT Z (504000 2005000 560).
Look at the polygon in the figure: one exterior ring and one interior ring, a hole. Each ring is a linear ring: A closed line whose first and last positions are the same and which does not cross itself. A polygon is bounded by one exterior ring and any number of interior rings, its holes., a line that ends where it starts and never crosses itself, which the PostGIS manual defines as “both closed and simple” (PostGIS manual). A multi type holds several parts of one type, such as a ward cut in two by a river, and a GeometryCollection holds parts of any types.
In well-known text (WKT): A plain-text form of a geometry from the OGC Simple Features standard, such as POINT (4000 5000). The same name is also used for a text format of coordinate reference systems. More in “Vector and raster data models” each type’s name is its keyword. Clinic C2 of the track’s town, now placed in
UTM (Universal Transverse Mercator): A system of 60 zones, each 6° of longitude wide and each a transverse Mercator projection with a scale factor of 0.9996 on its central meridian. Mainland India spans zones 42N to 47N. More in “Map projections and their distortions” so that its position is in metres, is POINT (504000 2005000). Ward W02 with a pond is
POLYGON ((…shell…), (…hole…)), the shell first. A third number is a height, Z, and a fourth a measure, M, such as the
distance along a road: POINT Z (504000 2005000 560). A geometry can also be empty, POLYGON EMPTY, which is what a
clip leaves when nothing overlaps.
None of these texts says which coordinate reference system (CRS): The definition that gives coordinates a meaning on the Earth: what the numbers measure, in which units and from which origin, including the datum and, for flat maps, the map projection. More in “What GIS is and the questions it answers” its numbers belong to. PostGIS adds
one in its extended form, EWKT and EWKB: PostGIS's extended forms of WKT and WKB, which carry the SRID inside the value: SRID=32643;POINT(504000 2005000) as text, or a flag in the type code followed by four bytes of SRID in binary.: SRID=32643;POINT(504000 2005000), where the
SRID (spatial reference identifier): A whole number that names a coordinate reference system in a database or file, such as 4326 for WGS 84 latitude and longitude or 32643 for WGS 84 / UTM zone 43N. 32643 names WGS 84 / UTM zone 43N.
WKT is written for people, WKB for programs
Text is easy to read and slow to parse, and decimal digits can round a coordinate. Programs pass geometry as well-known binary (WKB): The binary form of a geometry from the OGC Simple Features standard: one byte for the byte order, four for the type code, then each coordinate as an 8-byte floating-point number. A 2D point takes 21 bytes. instead, which the PostGIS manual calls “a portable, full-precision representation”. Step through the 21 bytes of clinic C2 as a reader takes them apart:
The drawing shows the whole scene; its 8 steps are listed below.
Every step as text
- Twenty-one bytes. Clinic C2, at 504,000 m east and 2,005,000 m north in UTM zone 43N, is 21 bytes of WKB. In hex they read 01 01000000 0000000000C31E41 0000000008983E41.
- The byte order. The first byte says how to read every number after it. 01 means little-endian, least significant byte first, which is what almost every program writes; 00 would mean big-endian.
- The type code. The next four bytes, 01 00 00 00 read little-endian, are the number 1: a Point. 2 is a LineString, 3 a Polygon, 4 to 6 their Multi forms and 7 a GeometryCollection.
- No SRID in plain WKB. Plain WKB carries no coordinate reference system, so the reader goes straight to the coordinates. Which CRS the numbers belong to must travel beside the bytes, in a column or a file header.
- Eight bytes of x. A Point has two coordinates, each an 8-byte floating-point number. 00 00 00 00 00 C3 1E 41 is 504000.0, the easting in metres, exact to the last bit.
- Eight bytes of y. The last eight bytes, 00 00 00 00 08 98 3E 41, are 2005000.0, the northing. 1 + 4 + 8 + 8 makes the 21 bytes, and the value is read.
- PostGIS's EWKB. PostGIS sets a flag in the type, 0x20000000, so its bytes read 01000020, and four more bytes, 837F0000, give SRID 32643 before the coordinates: 25 bytes in all.
- A third number. A point with a height has a third double. ISO WKB writes the type as 1001, E9030000, while PostGIS's EWKB keeps 1 and sets another flag, 0x80000000, so its bytes read 01000080.
The program below writes the same point byte by byte with Python’s struct module, in both byte orders and as
PostGIS’s extended form, and then reads it back. Change x and y and run it again to see which bytes move:
"""Clinic C2 of the town, written as well-known binary (WKB) byte by byte, and read back.
The town lies in UTM zone 43N, so the clinic's position is in metres: 504,000 m east, 2,005,000 m north.
"""
import struct
POINT, LINESTRING, POLYGON = 1, 2, 3 # the type codes of the simple-features standard
EWKB_SRID_FLAG = 0x20000000 # PostGIS's extension: a SRID follows the type
x, y = 504_000.0, 2_005_000.0
def show(label, parts):
"""Prints each field of a WKB value in hex, separated by bars."""
print(f"{label} ({sum(len(p) for p in parts)} bytes)")
print(" " + " | ".join(p.hex().upper() for p in parts))
# Little-endian (byte order 1), as almost every program writes it.
le = [struct.pack("<B", 1), struct.pack("<I", POINT), struct.pack("<d", x), struct.pack("<d", y)]
show("WKB, little-endian", le)
# The same point, big-endian (byte order 0): the same numbers, bytes reversed.
be = [struct.pack(">B", 0), struct.pack(">I", POINT), struct.pack(">d", x), struct.pack(">d", y)]
show("WKB, big-endian", be)
# PostGIS's extended WKB: the type gets a flag, and the SRID 32643 follows it.
ewkb = [struct.pack("<B", 1), struct.pack("<I", POINT | EWKB_SRID_FLAG), struct.pack("<I", 32643), struct.pack("<d", x), struct.pack("<d", y)]
show("EWKB with SRID 32643", ewkb)
# Reading it back: the first byte says how to read the rest.
data = b"".join(le)
order = "<" if data[0] == 1 else ">"
kind, rx, ry = struct.unpack(order + "Idd", data[1:])
print(f"read back: type {kind}, x = {rx:.0f}, y = {ry:.0f}")
# Ward W02, a 2 km square, as WKT and as WKB: one ring of five positions.
ring = [(502_000, 2_004_000), (504_000, 2_004_000), (504_000, 2_006_000), (502_000, 2_006_000), (502_000, 2_004_000)]
wkt = "POLYGON ((" + ", ".join(f"{a} {b}" for a, b in ring) + "))"
wkb = struct.pack("<BII", 1, POLYGON, 1) + struct.pack("<I", len(ring)) + b"".join(struct.pack("<dd", a, b) for a, b in ring)
print(wkt)
print(f"W02 as WKT: {len(wkt)} characters; as WKB: {len(wkb)} bytes")Output
WKB, little-endian (21 bytes) 01 | 01000000 | 0000000000C31E41 | 0000000008983E41 WKB, big-endian (21 bytes) 00 | 00000001 | 411EC30000000000 | 413E980800000000 EWKB with SRID 32643 (25 bytes) 01 | 01000020 | 837F0000 | 0000000000C31E41 | 0000000008983E41 read back: type 1, x = 504000, y = 2005000 POLYGON ((502000 2004000, 504000 2004000, 504000 2006000, 502000 2006000, 502000 2004000)) W02 as WKT: 90 characters; as WKB: 93 bytes
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 wkb_by_hand.py
Runs on this device, in your browser. The first run downloads Python (about 13.5 MB), which is kept for the next runs.
Your run, in this browser
The big-endian version holds the same numbers with each field’s bytes reversed, and the extended form is 25 bytes because the SRID takes four. The last line holds a surprise: ward W02 is 93 bytes of WKB and only 90 characters of WKT. WKB spends 16 bytes on every position whatever its digits, so it is not always smaller; what it buys is exactness, every bit of each double, and no text to parse. GEOS, through Shapely, writes the very same bytes:
"""The same clinic written by Shapely (GEOS): plain WKB, PostGIS's EWKB with a SRID, and a point with a height."""
import shapely
c2 = shapely.from_wkt("POINT (504000 2005000)")
print("WKB: ", shapely.to_wkb(c2, hex=True))
print("EWKB: ", shapely.to_wkb(shapely.set_srid(c2, 32643), hex=True, include_srid=True))
# A third number is a height (Z). ISO WKB adds 1000 to the type code; PostGIS's EWKB sets a flag bit instead.
c2z = shapely.from_wkt("POINT Z (504000 2005000 560)")
print("ISO Z: ", shapely.to_wkb(c2z, hex=True, flavor="iso"))
print("EWKB Z: ", shapely.to_wkb(c2z, hex=True, flavor="extended"))
# Reading WKB back gives the same geometry, and WKT is the readable form of it.
back = shapely.from_wkb(shapely.to_wkb(c2))
print("read back:", back.wkt, "| equal:", back.equals(c2))
print("empty: ", shapely.from_wkt("POLYGON EMPTY").wkt, "| is_empty:", shapely.from_wkt("POLYGON EMPTY").is_empty)
print(f"Shapely {shapely.__version__}, GEOS {shapely.geos_version_string}")Output
WKB: 01010000000000000000C31E410000000008983E41 EWKB: 0101000020837F00000000000000C31E410000000008983E41 ISO Z: 01E90300000000000000C31E410000000008983E410000000000808140 EWKB Z: 01010000800000000000C31E410000000008983E410000000000808140 read back: POINT (504000 2005000) | equal: True empty: POLYGON EMPTY | is_empty: True Shapely 2.1.2, GEOS 3.13.1
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 shapely_wkb.py
For a point with a height, ISO WKB adds 1,000 to the type code, so 1001 is a Point Z, while PostGIS sets a flag bit instead, which is why the two hex strings differ in their second field.
Hex that starts 0101000020 is not plain WKB
Hex copied out of PostGIS often starts 0101000020: the type is 0x20000001, a Point with PostGIS’s SRID flag, followed
by four bytes of SRID. A reader that knows only ISO WKB sees type code 536,870,913 and refuses it. Ask PostGIS for
ST_AsBinary when another program will read the bytes, or a reader that accepts EWKB, such as Shapely’s from_wkb.
Checkpoint
Check: WKT, WKB, EWKT and EWKB · The seven simple-features types
A valid polygon keeps five rules
The validity rules come from the simple-features standard, and the PostGIS manual lists them for polygons in plain words. A polygon is valid 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. when:
- its rings are simple: they do not cross or touch themselves, so there are no cut lines, spikes or loops;
- its rings do not cross each other;
- its rings touch each other, if at all, only at single points;
- its holes lie inside its exterior ring;
- its interior is connected, so no touching splits it into separate parts.
Each ring must also be closed and hold at least four positions. Lines need two distinct points, and points have no rules at all. Here are six polygons that break the rules, each in a different way:
Six polygons, each breaking one validity rule
Text description of the diagram
Six panels in two columns, each naming a broken polygon and how it breaks the rules.
Bow-tie: two triangles meeting at a point, drawn as one ring that crosses itself; a circle marks the crossing at (1, 1). Spike: a square with a dashed line sticking out of its right side to (15, 5), where a circle marks the tip: the ring goes out and comes straight back. Hole outside: a square, and beside it, outside, a small dashed square meant as its hole. Hole crossing: a square with a small dashed rectangle that sits across its right edge, partly inside and partly outside. Too few points: a single line with a dot at each end, a ring of three positions that encloses nothing. Not closed: a square drawn with three solid sides and the last side dashed, because the ring ends at (0, 10) instead of returning to its start.
Finding the bow-tie’s fault is a matter of testing every pair of edges that are not neighbours. Step through a short program that does exactly that on the bow-tie’s ring:
Program · trace/find_crossing.py
# Does the bow-tie cross itself? Test each
# pair of edges that are not neighbours.
ring = [(0, 0), (2, 2), (2, 0), (0, 2), (0, 0)]
edges = list(zip(ring, ring[1:]))
def side(o, a, b):
"""Positive when b is left of o -> a."""
ax, ay = a[0] - o[0], a[1] - o[1]
bx, by = b[0] - o[0], b[1] - o[1]
return ax * by - ay * bx
found = []
last = len(edges) - 1
for i in range(len(edges)):
for j in range(i + 2, len(edges)):
if (i, j) == (0, last):
continue # they meet at (0, 0)
(p, q), (r, s) = edges[i], edges[j]
d1, d2 = side(r, s, p), side(r, s, q)
if d1 * d2 < 0 and side(p, q, r) * side(p, q, s) < 0:
t = d1 / (d1 - d2)
x = p[0] + t * (q[0] - p[0])
y = p[1] + t * (q[1] - p[1])
found.append((x, y))
print(f"edges {i} and {j} cross at {(x, y)}")
print("invalid" if found else "valid")Variables
- just taken by the loop
- changed in this step
- passed by the loop
- read next
- where the loop stopped
- the stop value, never given
Printed so far
Nothing printed yetedges 0 and 2 cross at (1.0, 1.0) invalid
Every step was recorded from a run with Python 3.14.8 on macOS 26 arm64.
The drawing shows the variables when the program ends, and the program's 34 steps are listed below.
Change the program and run it here. Runs on this device, in your browser. The first run downloads Python (about 13.5 MB), which is kept for the next runs.
Your run, in this browser
Every step as text
- Start. No line has run yet, so there are no variables. Python starts with line 3.
- Line 3.
ring = [(0, 0), (2, 2), (2, 0), (0, 2), (0, 0)]makesringa list of 5 items. - Line 4.
edges = list(zip(ring, ring[1:]))setsedgesto[((0, 0), (2, 2)), ((2, 2), (2, 0)), ((2, 0), (0, 2)), ((0, 2), (0, 0))]. - Line 7.
def side(o, a, b):makes the functionside. Its body runs only when it is called. - Line 14.
found = []makesfounda list of 0 items. - Line 15.
last = len(edges) - 1setslastto3. - Line 16.
for i in range(len(edges)):gives its first number:iis0. - Line 17.
for j in range(i + 2, len(edges)):gives its first number:jis2. - Line 18.
(i, j) == (0, last)is false (iis0,jis2andlastis3), so Python skips the block underif (i, j) == (0, last):. - Line 20.
(p, q), (r, s) = edges[i], edges[j]setspto(0, 0), setsqto(2, 2), setsrto(2, 0)and setssto(0, 2). - Line 21.
d1, d2 = side(r, s, p), side(r, s, q)setsd1to4and setsd2to-4. - Line 22.
d1 * d2 < 0 and side(p, q, r) * side(p, q, s) < 0is true (d1is4,d2is-4,pis(0, 0),qis(2, 2),ris(2, 0)andsis(0, 2)), so the block underif d1 * d2 < 0 and side(p, q, r) * side(p, q, s) < 0:runs. - Line 23.
t = d1 / (d1 - d2)setstto0.5. - Line 24.
x = p[0] + t * (q[0] - p[0])setsxto1.0. - Line 25.
y = p[1] + t * (q[1] - p[1])setsyto1.0. - Line 26.
found.append((x, y))adds(1.0, 1.0)to the end offound. Key point Edge 0 runs from (0, 0) to (2, 2) and edge 2 from (2, 0) to (0, 2). Each has the other's ends on opposite sides, so they cross, at (1, 1), and the ring is not simple. - Line 27.
print(f"edges {i} and {j} cross at {(x, y)}")printsedges 0 and 2 cross at (1.0, 1.0). - Line 17.
for j in range(i + 2, len(edges)):gives the next number:jis3. - Line 18.
(i, j) == (0, last)is true (iis0,jis3andlastis3), so the block underif (i, j) == (0, last):runs. - Line 19.
continueskips the rest of this pass and goes back to the first line of the loop. Edges 0 and 3 meet at the ring's start, (0, 0), as the last edge of any ring meets the first. A shared vertex is not a crossing, so this pair is skipped. - Line 17.
for j in range(i + 2, len(edges)):has no number left: a range stops before its stop value,4, so the loop ends. - Line 16.
for i in range(len(edges)):gives the next number:iis1. - Line 17.
for j in range(i + 2, len(edges)):gives its first number:jis3. - Line 18.
(i, j) == (0, last)is false (iis1,jis3andlastis3), so Python skips the block underif (i, j) == (0, last):. - Line 20.
(p, q), (r, s) = edges[i], edges[j]setspto(2, 2)(it was(0, 0)), setsqto(2, 0)(it was(2, 2)), setsrto(0, 2)(it was(2, 0)) and setssto(0, 0)(it was(0, 2)). - Line 21.
d1, d2 = side(r, s, p), side(r, s, q)d1keeps its value,4and setsd2to4(it was-4). - Line 22.
d1 * d2 < 0 and side(p, q, r) * side(p, q, s) < 0is false (d1is4,d2is4,pis(2, 2),qis(2, 0),ris(0, 2)andsis(0, 0)), so Python skips the block underif d1 * d2 < 0 and side(p, q, r) * side(p, q, s) < 0:. - Line 17.
for j in range(i + 2, len(edges)):has no number left: a range stops before its stop value,4, so the loop ends. - Line 16.
for i in range(len(edges)):gives the next number:iis2. - Line 17.
for j in range(i + 2, len(edges)):has nothing to go through, so its body never runs and the loop is over at once. - Line 16.
for i in range(len(edges)):gives the next number:iis3. - Line 17.
for j in range(i + 2, len(edges)):has nothing to go through, so its body never runs and the loop is over at once. - Line 16.
for i in range(len(edges)):has no number left: a range stops before its stop value,4, so the loop ends. - Line 28.
print("invalid" if found else "valid")printsinvalid. That was the last line, so the program ends. One crossing is enough to make the polygon invalid. GEOS finds the same point and reports Self-intersection[1 1].
The same tests, with checks for closure, length, spikes and holes, make a small validity checker. It runs on the six broken polygons and on two that look odd but are valid: one repeats a position, and one has a hole that touches the shell at a single point, which rule 3 allows. The shapes are in a file of their own, which the next programs share:
"""Eight polygons for the validity examples: six broken in a different way each, and two that look odd but are valid.
Each is (name, exterior ring, holes), with positions as (x, y) in metres.
"""
SQUARE = [(0, 0), (10, 0), (10, 10), (0, 10), (0, 0)]
SHAPES = [
("bow-tie", [(0, 0), (2, 2), (2, 0), (0, 2), (0, 0)], []),
("spike", [(0, 0), (10, 0), (10, 5), (15, 5), (10, 5), (10, 10), (0, 10), (0, 0)], []),
("hole outside", SQUARE, [[(12, 2), (14, 2), (14, 4), (12, 4), (12, 2)]]),
("hole crossing", SQUARE, [[(8, 4), (12, 4), (12, 6), (8, 6), (8, 4)]]),
("too few points", [(0, 0), (10, 0), (0, 0)], []),
("not closed", [(0, 0), (10, 0), (10, 10), (0, 10)], []),
("repeated point", [(0, 0), (10, 0), (10, 0), (10, 10), (0, 10), (0, 0)], []),
("hole touches once", SQUARE, [[(0, 5), (5, 3), (5, 7), (0, 5)]]),
]
def wkt(shell, holes):
"""The polygon as well-known text."""
rings = [shell, *holes]
return "POLYGON (" + ", ".join("(" + ", ".join(f"{x} {y}" for x, y in r) + ")" for r in rings) + ")"The checker reads them, says what is wrong with each, and prints the area the usual formula gives:
"""A small validity checker for polygons, written from the rules, run on eight shapes."""
from shapes import SHAPES
def signed_area(ring):
"""Shoelace formula: positive when the ring runs anticlockwise (x east, y north)."""
return sum(x1 * y2 - x2 * y1 for (x1, y1), (x2, y2) in zip(ring, ring[1:])) / 2
def cross(o, a, b):
return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])
def crossing(p, q, r, s):
"""The point where segments pq and rs cross in their interiors, or None."""
d1, d2, d3, d4 = cross(r, s, p), cross(r, s, q), cross(p, q, r), cross(p, q, s)
if d1 * d2 < 0 and d3 * d4 < 0:
t = d1 / (d1 - d2)
return (p[0] + t * (q[0] - p[0]), p[1] + t * (q[1] - p[1]))
return None
def inside(pt, ring):
"""Ray casting: is pt inside the ring (not on it)?"""
x, y = pt
hits = 0
for (x1, y1), (x2, y2) in zip(ring, ring[1:]):
if (y1 > y) != (y2 > y) and x < x1 + (y - y1) * (x2 - x1) / (y2 - y1):
hits += 1
return hits % 2 == 1
def problems(shell, holes):
"""Why the polygon is invalid, or [] when it follows the rules."""
found = []
for ring in [shell, *holes]:
if ring[0] != ring[-1]:
return [f"ring not closed: it starts at {ring[0]} and ends at {ring[-1]}"]
pts = [p for i, p in enumerate(ring) if i == 0 or p != ring[i - 1]] # repeated points are allowed
if len(pts) < 4:
return [f"too few points: a ring needs 4 positions, this one has {len(pts)}"]
segs = list(zip(pts, pts[1:]))
for i, (p, q) in enumerate(segs):
nxt = segs[(i + 1) % len(segs)]
if cross(p, q, nxt[1]) == 0 and (q[0] - p[0]) * (nxt[1][0] - q[0]) + (q[1] - p[1]) * (nxt[1][1] - q[1]) < 0:
found.append(f"spike: the ring goes out to {q} and straight back")
for r, s in segs[i + 2 :]:
x = crossing(p, q, r, s)
if x:
found.append(f"self-intersection: the ring crosses itself at ({x[0]:g}, {x[1]:g})")
for hole in holes:
hits = [crossing(p, q, r, s) for p, q in zip(shell, shell[1:]) for r, s in zip(hole, hole[1:])]
hits = [h for h in hits if h]
if hits:
found.append(f"hole crosses the shell at ({hits[0][0]:g}, {hits[0][1]:g})")
elif not inside(hole[0], shell):
found.append(f"hole lies outside the shell, starting at {hole[0]}")
return found
for name, shell, holes in SHAPES:
issues = problems(shell, holes)
area = abs(signed_area(shell)) - sum(abs(signed_area(h)) for h in holes) if shell[0] == shell[-1] else None
verdict = "valid" if not issues else issues[0]
print(f"{name:18} {verdict}")
if area is not None:
print(f"{'':18} area by the formula: {area:g} m²")
# Repairing the bow-tie by hand: cut it where it crosses itself, into two triangles.
a, b = [(0, 0), (1, 1), (0, 2), (0, 0)], [(1, 1), (2, 2), (2, 0), (1, 1)]
print(f"bow-tie signed area {signed_area(SHAPES[0][1]):g}; cut at (1, 1) into triangles of {signed_area(a):g} and {signed_area(b):g}")
print(f"the two parts cover {abs(signed_area(a)) + abs(signed_area(b)):g} m²")Output
bow-tie self-intersection: the ring crosses itself at (1, 1)
area by the formula: 0 m²
spike spike: the ring goes out to (15, 5) and straight back
area by the formula: 100 m²
hole outside hole lies outside the shell, starting at (12, 2)
area by the formula: 96 m²
hole crossing hole crosses the shell at (10, 4)
area by the formula: 92 m²
too few points too few points: a ring needs 4 positions, this one has 3
area by the formula: 0 m²
not closed ring not closed: it starts at (0, 0) and ends at (0, 10)
repeated point valid
area by the formula: 100 m²
hole touches once valid
area by the formula: 90 m²
bow-tie signed area 0; cut at (1, 1) into triangles of 1 and -1
the two parts cover 2 m²
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 check.py
Runs on this device, in your browser. The first run downloads Python (about 13.5 MB), which is kept for the next runs.
Your run, in this browser
The areas show why validity matters. The bow-tie claims 0 m², because the formula adds one half and subtracts the other; the polygon whose hole lies outside claims 96 m², because the formula subtracts a hole that removes nothing. GEOS, the engine behind Shapely, PostGIS and QGIS, agrees with the checker on every shape and refuses the unclosed ring outright:
"""The same eight polygons checked and repaired by GEOS, through Shapely: its reason, the area it computes, what three
repairs return, and what an overlay does with the invalid bow-tie."""
import shapely
from shapes import SHAPES, wkt
print(f"Shapely {shapely.__version__}, GEOS {shapely.geos_version_string}")
for name, shell, holes in SHAPES:
try:
g = shapely.from_wkt(wkt(shell, holes))
except shapely.errors.GEOSException as e:
print(f"{name}: refused when read: {e}")
continue
print(f"{name}: {shapely.is_valid_reason(g)}, area {g.area:g}")
if not g.is_valid:
line = shapely.make_valid(g)
structure = shapely.make_valid(g, method="structure")
zero = g.buffer(0)
print(f" make_valid: {line.geom_type} {line.area:g} | structure: {structure.geom_type} {structure.area:g} | buffer(0): {zero.geom_type} {zero.area:g}")
bowtie = shapely.from_wkt(wkt(*SHAPES[0][1:]))
print("bow-tie, make_valid:", shapely.make_valid(bowtie).wkt)
print("bow-tie, buffer(0): ", bowtie.buffer(0).wkt)
try:
shapely.intersection(bowtie, shapely.from_wkt(wkt([(0, 0), (2, 0), (2, 2), (0, 2), (0, 0)], [])))
except shapely.errors.GEOSException as e:
print("bow-tie intersected with a square:", e)Output
Shapely 2.1.2, GEOS 3.13.1 bow-tie: Self-intersection[1 1], area 0 make_valid: MultiPolygon 2 | structure: MultiPolygon 2 | buffer(0): Polygon 1 spike: Self-intersection[15 5], area 100 make_valid: GeometryCollection 100 | structure: Polygon 100 | buffer(0): Polygon 100 hole outside: Hole lies outside shell[12 2], area 96 make_valid: MultiPolygon 104 | structure: MultiPolygon 104 | buffer(0): Polygon 100 hole crossing: Self-intersection[10 4], area 92 make_valid: MultiPolygon 100 | structure: Polygon 96 | buffer(0): Polygon 96 too few points: Too few points in geometry component[0 0], area 0 make_valid: LineString 0 | structure: LineString 0 | buffer(0): Polygon 0 not closed: refused when read: IllegalArgumentException: Points of LinearRing do not form a closed linestring repeated point: Valid Geometry, area 100 hole touches once: Valid Geometry, area 90 bow-tie, make_valid: MULTIPOLYGON (((0 2, 1 1, 0 0, 0 2)), ((2 0, 1 1, 2 2, 2 0))) bow-tie, buffer(0): POLYGON ((1 1, 2 2, 2 0, 1 1)) bow-tie intersected with a square: TopologyException: side location conflict at 1 1. This can occur if the input geometry is invalid.
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 geos_check.py
Its last line is the cost of leaving the bow-tie as it is: an intersection with a square stops with a
TopologyException at (1, 1). PostGIS does not check validity when it loads data, as its manual says, so an invalid
shape can sit in a table for years until the first overlay touches it.
Checkpoint
Check: The validity rules of a polygon
One repair keeps the area, the other loses half
GEOS repairs an invalid geometry with make_valid, which Shapely and PostGIS (ST_MakeValid) both call. It has two
methods. The default, linework, keeps every edge and vertex and assigns areas inside and outside by alternation, which
the Shapely documentation warns can give “rather complex results”
(Shapely, make_valid). The structure
method treats exterior rings as area and holes as removed area. An older trick, a buffer of distance 0, also returns a
valid polygon, and the Shapely manual says plainly that with it “your mileage may vary”
(Shapely manual). Here is what each did to the bow-tie:
make_valid keeps both triangles of the bow-tie, buffer(0) drops one
Text description of the diagram
Three panels side by side.
As stored: the bow-tie's single ring crossing itself, with a circle at the crossing. The note says its area is 0 by the formula, although it really covers 2 square metres.
make_valid: the same outline filled as two triangles that meet at the crossing point. The note says 2 triangles, area 2.
buffer(0): the right-hand triangle filled and the left-hand triangle drawn as a dashed grey outline only, because the zero buffer dropped it. The note says 1 triangle, area 1.
| Shape | make_valid | structure | buffer(0) |
|---|---|---|---|
| bow-tie | 2 parts, 2 m² | 2 parts, 2 m² | 1 part, 1 m² |
| spike | polygon and a line, 100 m² | polygon, 100 m² | polygon, 100 m² |
| hole outside | 2 parts, 104 m² | 2 parts, 104 m² | polygon, 100 m² |
| hole crossing | 2 parts, 100 m² | polygon, 96 m² | polygon, 96 m² |
These are the numbers GEOS printed above. The buffer kept one triangle of the bow-tie and silently dropped the other. For the hole that crossed its shell, the methods disagree on what was meant: linework keeps the part of the hole outside the shell as a polygon of its own, 100 m² in all, while structure subtracts only the part inside, 96 m². No method knows the surveyor’s intention, so the repair is a decision you check, not a step you trust.
Key idea
Repair with make_valid, never with buffer(0), and compare the area and the geometry type before and after: a repair that changes either is a decision about what the shape was meant to be.
Checkpoint
Check: Repairing invalid polygons
GeoJSON runs shells anticlockwise, shapefiles clockwise
A ring’s direction is the sign of its area. The shoelace formula, half the sum of x₁y₂ − x₂y₁ over the ring’s edges, is positive when the ring runs anticlockwise with x east and y north, and negative when it runs clockwise. Ward W02 digitised clockwise gives −4,000,000 m²: a 2 km square, run the other way.
Formats disagree about the direction. GeoJSON’s standard says rings “MUST follow the right-hand rule”, which it spells out as “exterior rings are counterclockwise, and holes are clockwise”, and in the same breath that parsers “SHOULD NOT reject” polygons wound otherwise (RFC 7946, section 3.1.6). A shapefile wants the opposite: GDAL’s driver documentation says outer rings “should be oriented clockwise” and inner rings counterclockwise (GDAL, Shapefile driver). QGIS’s “Force right-hand-rule” algorithm also makes exterior rings clockwise, the opposite of RFC 7946’s rule of the same name, and the QGIS 4.2 manual recommends its explicit “Force polygons clockwise” algorithm, added in 4.0, because of this inconsistency (QGIS, vector geometry algorithms). One rule pictures the area on your left as you walk the ring, the other on your right:
The two conventions run the same rings in opposite directions
Text description of the diagram
Two panels, each showing the same square ward with a small square pond in its middle, and arrows on the rings.
RFC 7946: an arrow on the right side of the shell points up and an arrow on the top points left, so the shell runs anticlockwise; the arrow on the left side of the pond points up, so the hole runs clockwise.
Shapefile: an arrow on the top of the shell points right and an arrow on the right side points down, so the shell runs clockwise; the arrow on the right side of the pond points up, so the hole runs anticlockwise.
Turn a ring the other way and its area keeps its size and flips its sign. The program below reads the direction of W02’s rings from the sign, orients them both ways, and checks that the ward’s area does not depend on it:
"""Which way a ring runs, from the sign of its area, and the two conventions that disagree about it."""
def signed_area(ring):
"""Shoelace formula: positive for an anticlockwise ring when x runs east and y north."""
return sum(x1 * y2 - x2 * y1 for (x1, y1), (x2, y2) in zip(ring, ring[1:])) / 2
def direction(ring):
return "anticlockwise" if signed_area(ring) > 0 else "clockwise"
def orient(ring, anticlockwise):
"""The ring itself, or reversed, so that it runs the way asked for."""
return ring if (signed_area(ring) > 0) == anticlockwise else ring[::-1]
# Ward W02 (2 km square) with a pond (a 200 m square hole), in metres of UTM zone 43N.
shell = [(502_000, 2_004_000), (502_000, 2_006_000), (504_000, 2_006_000), (504_000, 2_004_000), (502_000, 2_004_000)]
pond = [(502_900, 2_004_900), (503_100, 2_004_900), (503_100, 2_005_100), (502_900, 2_005_100), (502_900, 2_004_900)]
print(f"as digitised: shell {direction(shell)}, area {signed_area(shell):,.0f}")
print(f" pond {direction(pond)}, area {signed_area(pond):,.0f}")
for name, outer_ac in [("RFC 7946", True), ("Shapefile", False)]:
s, h = orient(shell, outer_ac), orient(pond, not outer_ac)
print(f"{name}: shell {direction(s)}, pond {direction(h)}")
print(f" the shell runs {s[0]} -> {s[1]}")
area = abs(signed_area(shell)) - abs(signed_area(pond))
print(f"ward area {area / 10_000:,.0f} ha, whichever way the rings run")Output
as digitised: shell clockwise, area -4,000,000
pond anticlockwise, area 40,000
RFC 7946: shell anticlockwise, pond clockwise
the shell runs (502000, 2004000) -> (504000, 2004000)
Shapefile: shell clockwise, pond anticlockwise
the shell runs (502000, 2004000) -> (502000, 2006000)
ward area 396 ha, whichever way the rings run
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 winding.py
Runs on this device, in your browser. The first run downloads Python (about 13.5 MB), which is kept for the next runs.
Your run, in this browser
GDAL applies the conventions for you, but only where a format demands one. Written through pyogrio, the same ward digitised both ways came out like this:
"""Ward W02 with its pond, digitised clockwise and then anticlockwise, written by GDAL (through pyogrio) as a
shapefile and as GeoJSON with and without RFC 7946, and read back from the files to see which way each shell runs."""
import json
import struct
import geopandas as gpd
import pyogrio
from shapely.geometry import Polygon
CW_SHELL = [(502_000, 2_004_000), (502_000, 2_006_000), (504_000, 2_006_000), (504_000, 2_004_000), (502_000, 2_004_000)]
CCW_POND = [(502_900, 2_004_900), (503_100, 2_004_900), (503_100, 2_005_100), (502_900, 2_005_100), (502_900, 2_004_900)]
def direction(ring):
area = sum(x1 * y2 - x2 * y1 for (x1, y1), (x2, y2) in zip(ring, ring[1:])) / 2
return "anticlockwise" if area > 0 else "clockwise"
def shp_shell(path):
"""The first ring of the first polygon of a .shp file, read from its bytes (the record starts at byte 100)."""
data = open(path, "rb").read()
nparts, npoints = struct.unpack("<ii", data[144:152])
first, second = struct.unpack("<ii", data[152:160])
start = 152 + 4 * nparts
return [struct.unpack("<dd", data[start + 16 * i : start + 16 * i + 16]) for i in range(second)]
def geojson_shell(path):
return json.load(open(path))["features"][0]["geometry"]["coordinates"][0]
print(f"{'digitised':14}{'shapefile':14}{'GeoJSON':14}RFC 7946")
for shell, pond in [(CW_SHELL, CCW_POND), (CW_SHELL[::-1], CCW_POND[::-1])]:
ward = gpd.GeoDataFrame({"ward": ["W02"]}, geometry=[Polygon(shell, [pond])], crs="EPSG:32643")
ward.to_file("ward.shp")
ward.to_file("plain.geojson", driver="GeoJSON")
ward.to_file("rfc.geojson", driver="GeoJSON", RFC7946="YES")
found = [direction(shp_shell("ward.shp")), direction(geojson_shell("plain.geojson")), direction(geojson_shell("rfc.geojson"))]
print(f"{direction(shell):14}" + "".join(f"{d:14}" for d in found[:2]) + found[2])
print(f"GDAL {pyogrio.__gdal_version_string__} (pyogrio {pyogrio.__version__})")Output
digitised shapefile GeoJSON RFC 7946 clockwise clockwise clockwise anticlockwise anticlockwise clockwise anticlockwise anticlockwise GDAL 3.12.4 (pyogrio 0.13.0)
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 written_by_gdal.py
The shapefile writer makes every shell clockwise. The plain GeoJSON writer keeps whatever direction it was given, so a
clockwise ward produces GeoJSON that breaks RFC 7946, and only the RFC7946=YES option makes it anticlockwise. Shapely
2.1 has orient_polygons for the same job, with anticlockwise shells by default and exterior_cw=True for the
shapefile’s way (Shapely, orient_polygons).
Key idea
Name the direction, never the rule: write “shells anticlockwise, holes clockwise” in a specification, because “right-hand rule” means the opposite in GeoJSON and in QGIS.
Checkpoint
Check: Ring direction in GeoJSON and shapefiles
DE-9IM says exactly how two shapes meet
“Touches”, “within” and “overlaps” sound clear until a clinic stands on a ward boundary. The precise answer comes from DE-9IM: The dimensionally extended nine-intersection model: a 3 by 3 matrix of the dimensions where the interiors, boundaries and exteriors of two geometries meet, written as nine characters such as FF2F11212., the dimensionally extended nine-intersection model. Every geometry has an interior, a boundary and an exterior: a point is its own interior and has no boundary, a line’s boundary is its two ends, and a polygon’s boundary is its rings (PostGIS manual, spatial queries). The model records where the three parts of A meet the three parts of B, as the dimension of each meeting: 2 for an area, 1 for a line, 0 for a point and F for nothing at all. Here are the shapes of the examples:
The shapes whose DE-9IM matrices the programs compute
Text description of the diagram
A map at 1 pixel to 20 metres. Ward W02 is a 2 km square, from 502,000 to 504,000 m east and 2,004,000 to 2,006,000 m north. Ward W03 is the same size directly to its east, so the two share their 2 km edge at 504,000 m east. Ward W07 lies directly below W03, so it meets W02 only at one corner, at 504,000 m east and 2,004,000 m north.
A grey square, the park, sits inside W02 without touching its edges, from 502,500 to 503,500 m east and 2,004,500 to 2,005,500 m north. A dashed square, the flood zone, from 503,000 to 505,000 m east and 2,005,000 to 2,007,000 m north, covers parts of both W02 and W03 and reaches beyond them to the north. A hollow ring marks clinic C2 at 504,000 m east and 2,005,000 m north, exactly on the edge between W02 and W03.
Wards W02 and W03 share their 2 km edge and nothing else. Their matrix, row by row:
| W02 against W03 | Interior | Boundary | Exterior |
|---|---|---|---|
| Interior | F | F | 2 |
| Boundary | F | 1 | 1 |
| Exterior | 2 | 1 | 2 |
Read left to right and top to bottom, that is FF2F11212: the interiors never meet, and the boundaries meet in a line,
the shared edge. The named predicates are patterns over these nine characters, in which T means any of 0, 1 and 2 and
* means anything. PostGIS gives T*F**F*** for within: the interiors meet, and nothing of A is in B’s exterior
(PostGIS, ST_Within). Touches is any of FT*******, F**T***** and
F***T****: the interiors stay apart while a boundary meets something
(PostGIS, ST_Touches). Covered by asks only that nothing of A lies in B’s
exterior (PostGIS, ST_Covers).
This program works out the exact matrix of axis-aligned boxes and points, from the definitions alone, and matches it against those patterns:
"""The DE-9IM matrix of axis-aligned boxes and points, worked out exactly, and the named relationships it gives.
A box is (xmin, ymin, xmax, ymax); a point is (x, y). Rows are A's interior, boundary and exterior; columns are B's.
Each cell is the dimension of where they meet: 2 an area, 1 a line, 0 a point, F nothing.
"""
def edges(b):
x0, y0, x1, y1 = b
return [((x0, y0), (x1, y0)), ((x1, y0), (x1, y1)), ((x0, y1), (x1, y1)), ((x0, y0), (x0, y1))]
def seg_meet(s, t):
"""Dimension of the meeting of two axis-parallel segments: 1, 0 or None."""
(a, b), (c, d) = s, t
lox, hix = max(min(a[0], b[0]), min(c[0], d[0])), min(max(a[0], b[0]), max(c[0], d[0]))
loy, hiy = max(min(a[1], b[1]), min(c[1], d[1])), min(max(a[1], b[1]), max(c[1], d[1]))
if lox > hix or loy > hiy:
return None
return 1 if (hix - lox) + (hiy - loy) > 0 else 0
def seg_in_open(s, b):
"""Does a segment run through the open interior of box b?"""
(a, c), (x0, y0, x1, y1) = s, b
lox, hix = max(min(a[0], c[0]), x0), min(max(a[0], c[0]), x1)
loy, hiy = max(min(a[1], c[1]), y0), min(max(a[1], c[1]), y1)
if a[0] == c[0]:
return x0 < a[0] < x1 and loy < hiy
return y0 < a[1] < y1 and lox < hix
def seg_outside(s, b):
"""Does part of a segment lie outside the closed box b?"""
(a, c), (x0, y0, x1, y1) = s, b
return not all(x0 <= p[0] <= x1 and y0 <= p[1] <= y1 for p in (a, c))
def dim(*ds):
ds = [d for d in ds if d is not None]
return str(max(ds)) if ds else "F"
def relate(a, b):
if len(a) == 2: # a point against a box
x, y = a
x0, y0, x1, y1 = b
inside = x0 < x < x1 and y0 < y < y1
on_edge = x0 <= x <= x1 and y0 <= y <= y1 and not inside
cells = ["0" if inside else "F", "0" if on_edge else "F", "F" if inside or on_edge else "0"]
return "".join(cells) + "FFF" + "212"
if len(b) == 2:
m = relate(b, a)
return "".join(m[i] for i in (0, 3, 6, 1, 4, 7, 2, 5, 8)) # transpose
ii = 2 if max(a[0], b[0]) < min(a[2], b[2]) and max(a[1], b[1]) < min(a[3], b[3]) else None
ib = dim(*(1 if seg_in_open(e, a) else None for e in edges(b)))
ie = "F" if b[0] <= a[0] and b[1] <= a[1] and a[2] <= b[2] and a[3] <= b[3] else "2"
bi = dim(*(1 if seg_in_open(e, b) else None for e in edges(a)))
bb = dim(*(seg_meet(e, f) for e in edges(a) for f in edges(b)))
be = dim(*(1 if seg_outside(e, b) else None for e in edges(a)))
ei = "F" if a[0] <= b[0] and a[1] <= b[1] and b[2] <= a[2] and b[3] <= a[3] else "2"
eb = dim(*(1 if seg_outside(e, a) else None for e in edges(b)))
return dim(ii) + ib + ie + bi + bb + be + ei + eb + "2"
def matches(m, pattern):
"""T: any of 0, 1, 2; F: nothing; *: anything; a digit: that dimension."""
return all(p == "*" or (p == "T" and c != "F") or p == c for c, p in zip(m, pattern))
def named(m):
found = []
if matches(m, "T*F**F***"):
found.append("within")
if m[2] == "F" and m[5] == "F":
found.append("covered by")
if any(matches(m, p) for p in ("FT*******", "F**T*****", "F***T****")):
found.append("touches")
if not matches(m, "FF*FF****"):
found.append("intersects")
return ", ".join(found) or "disjoint"
W02, W03, W07 = (502_000, 2_004_000, 504_000, 2_006_000), (504_000, 2_004_000, 506_000, 2_006_000), (504_000, 2_002_000, 506_000, 2_004_000)
PARK, FLOOD, C2 = (502_500, 2_004_500, 503_500, 2_005_500), (503_000, 2_005_000, 505_000, 2_007_000), (504_000, 2_005_000)
for label, a, b in [("W02, W03", W02, W03), ("W02, W07", W02, W07), ("park, W02", PARK, W02), ("flood zone, W02", FLOOD, W02), ("C2, W02", C2, W02), ("W02, C2", W02, C2)]:
m = relate(a, b)
print(f"{label:16} {m} {named(m)}")Output
W02, W03 FF2F11212 touches, intersects W02, W07 FF2F01212 touches, intersects park, W02 2FF1FF212 within, covered by, intersects flood zone, W02 212101212 intersects C2, W02 F0FFFF212 covered by, touches, intersects W02, C2 FF20F1FF2 touches, intersects
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 boxes.py
Runs on this device, in your browser. The first run downloads Python (about 13.5 MB), which is kept for the next runs.
Your run, in this browser
GEOS computes the same six matrices and the same answers:
"""The same six pairs related by GEOS, through Shapely: the matrix and the named predicates it answers."""
import shapely
from shapely.geometry import Point, box
W02, W03, W07 = box(502_000, 2_004_000, 504_000, 2_006_000), box(504_000, 2_004_000, 506_000, 2_006_000), box(504_000, 2_002_000, 506_000, 2_004_000)
PARK, FLOOD, C2 = box(502_500, 2_004_500, 503_500, 2_005_500), box(503_000, 2_005_000, 505_000, 2_007_000), Point(504_000, 2_005_000)
for label, a, b in [("W02, W03", W02, W03), ("W02, W07", W02, W07), ("park, W02", PARK, W02), ("flood zone, W02", FLOOD, W02), ("C2, W02", C2, W02), ("W02, C2", W02, C2)]:
tests = [("within", a.within(b)), ("covered by", a.covered_by(b)), ("touches", a.touches(b)), ("intersects", a.intersects(b))]
print(f"{label:16} {shapely.relate(a, b)} {', '.join(n for n, t in tests if t) or 'disjoint'}")
print(f"C2 in W02: within {C2.within(W02)}, covered by {C2.covered_by(W02)}; W02 contains C2: {W02.contains(C2)}, covers: {W02.covers(C2)}")
print(f"GEOS {shapely.geos_version_string}")Output
W02, W03 FF2F11212 touches, intersects W02, W07 FF2F01212 touches, intersects park, W02 2FF1FF212 within, covered by, intersects flood zone, W02 212101212 intersects C2, W02 F0FFFF212 covered by, touches, intersects W02, C2 FF20F1FF2 touches, intersects C2 in W02: within False, covered by True; W02 contains C2: False, covers: True GEOS 3.13.1
Recorded with Python 3.14.8 on macOS 26 arm64. To run it yourself: mise exec python@3.14.8 -- python3 geos_relate.py
Clinic C2 is the case to remember. Its matrix against W02 is F0FFFF212: its interior, the point itself, meets W02 only
on the boundary. Within needs the interiors to meet, so C2 is not within W02, and by the same rule not within W03
either. Covered by needs only that nothing of C2 is outside W02, so C2 is covered by both wards.
A point on a boundary falls out of a within count
Counting clinics per ward with within leaves C2 out of both W02 and W03, so the ward totals add up to one fewer than
the clinics. Counting with intersects or covers puts it in both, one too many. Decide the rule for boundary
points, such as the ward with the lower number, and test that the counts add up to the total.
Checkpoint
Check: Reading a DE-9IM matrix
Practise: check, repair and relate geometries
Build three small tools of a geometry checker, with sample tests you can run here: one that orients rings either way, one that names what is wrong with a ring, and one that matches DE-9IM matrices against patterns:
Code lab · 3 steps
Build a geometry checker
Orient rings for GeoJSON or a shapefile, find rings that are open, too short or self-crossing, and test DE-9IM matrices against the named patterns.
- Orient rings for GeoJSON or a shapefile To do
- Name what is wrong with a ring To do
- Match DE-9IM matrices against patterns To do
Exercise · Easy · Python
Orient rings for GeoJSON or a shapefile
Write two functions in rings.py. signed_area(ring) returns the shoelace area of a closed ring, a list of (x, y) positions whose last equals its first: half the sum of x1 * y2 - x2 * y1 over its edges, positive when the ring runs anticlockwise with x east and y north. orient(ring, anticlockwise=True) returns the ring itself when it already runs the way asked for, and the ring reversed otherwise: anticlockwise for a GeoJSON shell, clockwise for a shapefile's.
The sample tests use ward W02 digitised clockwise, its pond digitised anticlockwise, and a triangle.
Starter code · rings.py
def signed_area(ring):
"""The shoelace area of a closed ring: positive when it runs anticlockwise (x east, y north)."""
# Replace this line: add up x1 * y2 - x2 * y1 over the ring's edges, then halve the sum.
return 0.0
def orient(ring, anticlockwise=True):
"""The ring, reversed if needed, so that it runs anticlockwise (or clockwise when anticlockwise is False)."""
# Replace this line: compare the sign of signed_area(ring) with the direction asked for.
return ringThe sample tests · test_rings.py
from rings import orient, signed_area
W02_CW = [(502000, 2004000), (502000, 2006000), (504000, 2006000), (504000, 2004000), (502000, 2004000)]
POND_CCW = [(502900, 2004900), (503100, 2004900), (503100, 2005100), (502900, 2005100), (502900, 2004900)]
TRIANGLE = [(0, 0), (4, 0), (0, 3), (0, 0)]
def test_signed_area():
"""gives the clockwise ward -4,000,000 and the anticlockwise pond 40,000"""
assert signed_area(W02_CW) == -4_000_000
assert signed_area(POND_CCW) == 40_000
assert signed_area(TRIANGLE) == 6
def test_orient_for_geojson():
"""reverses the clockwise ward and keeps the anticlockwise pond for a GeoJSON shell"""
assert orient(W02_CW) == W02_CW[::-1]
assert orient(POND_CCW, True) == POND_CCW
def test_orient_for_a_shapefile():
"""keeps the clockwise ward and reverses the pond for a shapefile's shell"""
assert orient(W02_CW, False) == W02_CW
assert orient(POND_CCW, False) == POND_CCW[::-1]
def test_area_keeps_its_size():
"""keeps the area's size when a ring is turned round"""
assert signed_area(orient(TRIANGLE, False)) == -6
assert abs(signed_area(orient(W02_CW))) == 4_000_000A hint
Pair each position with the next one, zip(ring, ring[1:]), add up x1 * y2 - x2 * y1 and halve it. A ring runs the way asked for when (signed_area(ring) > 0) == anticlockwise; otherwise return ring[::-1].
Results of the sample tests
| Test | Result | Details |
|---|
What your code printed
The sample tests run on this device, in your browser (Pyodide): nothing is sent to mysmartcopilot.com. The first run downloads Python (about 13.5 MB), which is kept for the next runs. A check in your browser is feedback for you, not proof that the code is right for every input.
Exercise · Medium · Python
Name what is wrong with a ring
Write ring_problems(ring) in problems.py. ring is a list of (x, y) positions. Return a list of the problems it has, in this order, or an empty list when it has none:
"not closed"when its last position differs from its first;"too few positions"when it has fewer than 4 positions once a position repeated right after itself is counted once;"crosses itself"when two of its edges that are not neighbours cross each other in their middles, as the bow-tie's edges do at (1, 1). An open ring's edges are its consecutive pairs; a closed ring's first and last edges are neighbours, since they meet at the start.
The sample tests use the lesson's bow-tie and squares, ward W02 digitised with two corners in the wrong order, an open ring, a ring of three positions and a ring with a repeated position, which is valid.
Starter code · problems.py
def ring_problems(ring):
"""The ring's problems, in this order: "not closed", "too few positions", "crosses itself" ([] when none)."""
problems = []
# Replace this line: test that the ring is closed, has 4 or more positions, and has no two edges that are not
# neighbours crossing each other.
return problemsThe sample tests · test_problems.py
from problems import ring_problems
SQUARE = [(0, 0), (10, 0), (10, 10), (0, 10), (0, 0)]
BOWTIE = [(0, 0), (2, 2), (2, 0), (0, 2), (0, 0)]
def test_a_good_ring():
"""finds nothing wrong with a square or a triangle"""
assert ring_problems(SQUARE) == []
assert ring_problems([(0, 0), (4, 0), (0, 3), (0, 0)]) == []
def test_the_bowtie():
"""finds that the bow-tie crosses itself"""
assert ring_problems(BOWTIE) == ["crosses itself"]
def test_a_ward_with_two_corners_swapped():
"""finds the crossing in ward W02 digitised with two corners in the wrong order, and none in the right order"""
swapped = [(502000, 2004000), (504000, 2004000), (502000, 2006000), (504000, 2006000), (502000, 2004000)]
assert ring_problems(swapped) == ["crosses itself"]
assert ring_problems([swapped[0], swapped[1], swapped[3], swapped[2], swapped[0]]) == []
def test_open_rings():
"""reports an open ring, and an open ring that also crosses itself"""
assert ring_problems(SQUARE[:-1]) == ["not closed"]
assert ring_problems(BOWTIE[:-1]) == ["not closed", "crosses itself"]
def test_too_few_positions():
"""reports a ring of three positions, also when one of four repeats"""
assert ring_problems([(0, 0), (10, 0), (0, 0)]) == ["too few positions"]
assert ring_problems([(0, 0), (0, 0), (10, 0), (0, 0)]) == ["too few positions"]
def test_a_repeated_position_is_allowed():
"""accepts a square with one position repeated"""
assert ring_problems([(0, 0), (10, 0), (10, 0), (10, 10), (0, 10), (0, 0)]) == []A hint
Drop each position that equals the one before it first. For the crossings, reuse the lesson's side(o, a, b) test: two edges pq and rs cross in their middles when side(r, s, p) * side(r, s, q) < 0 and side(p, q, r) * side(p, q, s) < 0. Test edge i against the edges from i + 2 on, and skip the pair of the first and last edge of a closed ring.
Results of the sample tests
| Test | Result | Details |
|---|
What your code printed
The sample tests run on this device, in your browser (Pyodide): nothing is sent to mysmartcopilot.com. The first run downloads Python (about 13.5 MB), which is kept for the next runs. A check in your browser is feedback for you, not proof that the code is right for every input.
Exercise · Easy · Python
Match DE-9IM matrices against patterns
Write two functions in de9im.py. matches(matrix, pattern) returns True when a nine-character DE-9IM matrix such as "FF2F11212" fits a pattern such as "T*F**F***": in a pattern, T stands for any of 0, 1 and 2, F for F, * for anything, and a digit for that digit. covered_by(matrix) returns True when the matrix of A against B fits one of "T*F**F***", "*TF**F***", "**FT*F***" and "**F*TF***": nothing of A lies in B's exterior, and A meets B somewhere.
The sample tests use the lesson's matrices: the park inside W02, clinic C2 on W02's edge, W02 against W03 and W07, and the flood zone.
Starter code · de9im.py
def matches(matrix, pattern):
"""True when a DE-9IM matrix such as "FF2F11212" fits a pattern such as "T*F**F***"."""
# Replace this line: compare the two strings character by character (T: 0, 1 or 2; *: anything).
return False
def covered_by(matrix):
"""True when the matrix of A against B fits one of the four covered-by patterns."""
# Replace this line: use matches() with "T*F**F***", "*TF**F***", "**FT*F***" and "**F*TF***".
return FalseThe sample tests · test_de9im.py
from de9im import covered_by, matches
PARK_W02, C2_W02 = "2FF1FF212", "F0FFFF212"
W02_W03, W02_W07, FLOOD_W02 = "FF2F11212", "FF2F01212", "212101212"
def test_within():
"""finds the park within W02 and clinic C2 not within it"""
assert matches(PARK_W02, "T*F**F***") is True
assert matches(C2_W02, "T*F**F***") is False
def test_touches():
"""finds that W02 touches W03 along an edge and W07 at a corner"""
assert matches(W02_W03, "F***T****") is True
assert matches(W02_W07, "F***T****") is True
assert matches(FLOOD_W02, "F***T****") is False
def test_digits():
"""tells a corner from an edge with a digit in the pattern"""
assert matches(W02_W07, "****0****") is True
assert matches(W02_W03, "****0****") is False
def test_covered_by():
"""finds the park and clinic C2 covered by W02, and W03 and the flood zone not"""
assert covered_by(PARK_W02) is True
assert covered_by(C2_W02) is True
assert covered_by(W02_W03) is False
assert covered_by(FLOOD_W02) is FalseA hint
Walk the two strings together with zip(matrix, pattern). A pair fits when the pattern character is *, when it is T and the matrix character is not F, or when the two characters are equal; the whole matrix fits when every pair does.
Results of the sample tests
| Test | Result | Details |
|---|
What your code printed
The sample tests run on this device, in your browser (Pyodide): nothing is sent to mysmartcopilot.com. The first run downloads Python (about 13.5 MB), which is kept for the next runs. A check in your browser is feedback for you, not proof that the code is right for every input.
Then test yourself across the six ideas, from types to matrices; the session starts with the idea your answers so far show to be the weakest:
Practice that adapts to you
Practise this lesson
18 questions about 6 ideas of this lesson. Each answer is explained at once, and what you get right is kept in this browser.
Interview questions
Warm-up (fresher): What is the difference between WKT and WKB, and when would you use each?
Both encode the same simple-features geometry. WKT is text, such as POINT (504000 2005000): easy to read, write by
hand and paste into a query, but slower to parse, and a coordinate written with too few digits is rounded. WKB is
binary: one byte for the byte order, four for the type code and eight for each coordinate, 21 bytes for a 2D point.
It keeps every bit of each double and needs no parsing, so databases and file formats such as GeoPackage and GeoParquet
store it. Neither carries a coordinate reference system; PostGIS’s extended forms, EWKT and EWKB, add the SRID. Use WKT
for people, tests and logs, and WKB between programs.
Two harder questions build on this lesson: a parcel layer whose areas fall short of the district, and a clinic on a ward boundary that a count leaves out. Each opens in the module’s bank with what it tests and a first step:
Model answers open with any pass: ₹79 for 1 week ($2.99), one payment that never renews.
-
Fresher to mid
A district's parcel layer of 40,000 polygons adds up to 3 % less area than the district itself, and some parcels report an area of 0 m². What do you suspect, how do you find the faulty parcels, and how do you repair them without changing the good ones?
Flashcards · 8 cards
Recall it, then check
Say the answer to yourself before you look: recalling it is what makes it stick.
-
Which seven geometry types make up the simple-features model?
Point, LineString and Polygon; their multi forms MultiPoint, MultiLineString and MultiPolygon; and GeometryCollection, which holds parts of any type. Any of them may carry a Z height and an M measure, and any may be empty.
-
How are the 21 bytes of a 2D point laid out in WKB?
One byte for the byte order, four for the type code (1 for a point), then eight bytes, a double, for each coordinate. PostGIS's EWKB adds a flag to the type code and four bytes of SRID, so the same point takes 25.
-
Hex copied out of PostGIS starts 0101000020. Why might a plain WKB reader refuse it?
The type field is 0x20000001: a point with PostGIS's SRID flag, followed by four bytes of SRID. A reader that knows only ISO WKB sees type code 536,870,913. Ask PostGIS for ST_AsBinary, or use a reader that accepts EWKB, such as Shapely's from_wkb.
-
Which five rules does a valid polygon keep?
Its rings do not cross or touch themselves; they do not cross each other; they meet, if at all, only at single points; its holes lie inside the shell; and its interior is in one piece. Every ring is also closed and holds at least four positions.
-
Why repair an invalid polygon with make_valid rather than buffer(0)?
A self-crossing bow-tie of 2 m² reports 0 m² until it is repaired. make_valid kept every edge and returned two parts, 2 m²; buffer(0) silently kept one triangle, 1 m². Check that the area and the kind of geometry survive each repair.
-
Which way should polygon rings run in GeoJSON and in a shapefile?
GeoJSON: shells anticlockwise and holes clockwise, though readers must accept other windings. Shapefile: shells clockwise, holes anticlockwise. Write the direction itself in a specification, since "right-hand rule" means opposite things in RFC 7946 and in QGIS.
-
What does the DE-9IM matrix FF2F11212 say about two shapes?
Each of its nine places gives the dimension where an interior, boundary or exterior of one shape meets one of the other: 2 area, 1 line, 0 point, F nothing. Here the interiors never meet and the boundaries share a line, as for two wards with a common edge.
-
A clinic stands exactly on a ward boundary. Is it within the ward, and is it covered by it?
It is not within: within needs the interiors to meet, and the point meets only the boundary. It is covered by the ward, because no part of it lies outside. Pick one rule for boundary points and check that the ward counts add up to the total.
Key takeaways
- Seven simple-features types describe almost every feature: points, line strings and polygons, their multi forms and the geometry collection, each optionally with Z and M and each possibly empty.
- WKT is text for people and WKB is exact binary for programs, 21 bytes for a 2D point; PostGIS’s EWKT and EWKB add the SRID, and neither plain form carries a coordinate reference system.
- A valid polygon has simple rings that do not cross, touch only at points, keep their holes inside and leave the interior connected; an invalid one gives wrong areas, such as 0 m² for a 2 m² bow-tie, and breaks overlays.
make_validrepairs without dropping area, whilebuffer(0)dropped half of the bow-tie; compare the area and the type before and after every repair.- GeoJSON wants shells anticlockwise and shapefiles want them clockwise, and “right-hand rule” names both, so write the direction itself.
- A DE-9IM matrix records the dimension where each part of two shapes meets, and the named predicates are patterns over it: a point on a boundary is covered by a polygon but not within it.
Terms in this lesson
- DE-9IM
- The dimensionally extended nine-intersection model: a 3 by 3 matrix of the dimensions where the interiors, boundaries and exteriors of two geometries meet, written as nine characters such as FF2F11212.
- EWKT and EWKB
- PostGIS's extended forms of WKT and WKB, which carry the SRID inside the value: SRID=32643;POINT(504000 2005000) as text, or a flag in the type code followed by four bytes of SRID in binary.
- linear ring
- A closed line whose first and last positions are the same and which does not cross itself. A polygon is bounded by one exterior ring and any number of interior rings, its holes.
- simple features
- The OGC and ISO model of vector geometry: points, line strings and polygons with straight edges between positions, their multi-part forms and the geometry collection, each optionally with Z and M values.
- SRID (spatial reference identifier)
- A whole number that names a coordinate reference system in a database or file, such as 4326 for WGS 84 latitude and longitude or 32643 for WGS 84 / UTM zone 43N.
- valid 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.
- well-known binary (WKB)
- The binary form of a geometry from the OGC Simple Features standard: one byte for the byte order, four for the type code, then each coordinate as an 8-byte floating-point number. A 2D point takes 21 bytes.
References
- Simple Feature Access, Part 1: Common Architecture (Open Geospatial Consortium (OGC))
- PostGIS manual: Data management (geometry validity, WKT, WKB, EWKB and EWKT) (PostGIS Project Steering Committee)
- PostGIS manual: Spatial queries (the DE-9IM) (PostGIS Project Steering Committee)
- PostGIS: ST_MakeValid (PostGIS Project Steering Committee)
- PostGIS: ST_Within, ST_Touches and ST_Covers (PostGIS Project Steering Committee)
- Shapely: make_valid (Shapely contributors)
- Shapely: orient_polygons (Shapely contributors)
- The Shapely user manual (buffer and the zero-distance buffer) (Shapely contributors)
- RFC 7946: The GeoJSON Format (Internet Engineering Task Force (IETF))
- GDAL: ESRI Shapefile / DBF driver (Open Source Geospatial Foundation (GDAL))
- QGIS 4.2 user manual: Vector geometry algorithms (QGIS Project)
Related tools
Report a problem with this lesson
Kept only in this browser. Your Learn progress