Shapely 2: geometry operations in Python

Which of 5,000 households lie within 1 km of one of 8 clinics? Asked one household and one clinic at a time, that is 40,000 distance questions passed one by one from Python to the geometry engine; Shapely 2 takes whole NumPy arrays of geometries and loops over them in compiled code, which is also how GeoPandas works underneath.

Here you build the track’s town as arrays and count the households of 12 wards with one call, predict the sides and the units of a buffer, combine buffers into unions and differences whose areas mean something, repair the input that breaks those operations, and keep results comparable when the engine underneath changes version.

One call tests 5,000 points against 12 wards

Shapely is a Python package “for manipulation and analysis of planar geometric objects”, and its functions are “NumPy ufuncs for operations using arrays of geometries” whose loops “are implemented in C” (Shapely documentation). The geometry work itself is done by GEOSGEOS (Geometry Engine, Open Source): The C and C++ library that computes predicates, buffers, unions and other overlays for Shapely, PostGIS, QGIS and GDAL. Each installation bundles or links its own version, which can change the exact vertices of a result., the engine that PostGIS and QGIS also use. GeoPandas adds tables on top: its documentation says plainly that “geometric operations are performed by shapely” (GeoPandas).

That changes how you write a test. shapely.points(xy) turns a 5,000 × 2 array of coordinates into an array of 5,000 points in one call, and shapely.contains(wards[:, np.newaxis], points[np.newaxis, :]) compares 12 wards with every point at once, by NumPy’s broadcasting, into a table of 12 rows and 5,000 columns. This vectorisedvectorised operation: A function that takes whole arrays, such as a NumPy array of geometries, and loops over them in compiled code, so that one call from Python answers thousands of questions. style is the main change from Shapely 1. The program below makes the town’s households and wards, counts the households of each ward from that table, and checks the counts against plain arithmetic on the grid:

"""The track's synthetic town in metres of UTM zone 43N (EPSG:32643): 12 square wards of 2 km, 5,000 households
from a seeded generator and 8 clinics. Every lesson of the track asks its questions about this town."""
import random

E0, N0 = 500_000, 2_000_000  # the town's south-west corner: easting and northing in metres
WARD = 2_000  # each ward is a 2 km square
COLS, ROWS = 4, 3  # 4 wards from west to east, 3 from south to north

# Metres east and north of the south-west corner. C2 stands exactly on the edge of W02 and W03, on purpose.
CLINICS = {
    "C1": (1_500, 5_200), "C2": (4_000, 5_000), "C3": (6_500, 4_200), "C4": (2_700, 3_100),
    "C5": (5_200, 2_800), "C6": (900, 1_100), "C7": (5_600, 1_200), "C8": (7_400, 900),
}


def ward_bounds():
    """W01 to W12, north-west first, each as (min x, min y, max x, max y)."""
    out = {}
    for row in range(ROWS):
        for col in range(COLS):
            x0, y0 = E0 + col * WARD, N0 + (ROWS - 1 - row) * WARD
            out[f"W{row * COLS + col + 1:02d}"] = (x0, y0, x0 + WARD, y0 + WARD)
    return out


def households(seed=42, n=5_000):
    """(id, x, y) of n households, most of them around three neighbourhood centres."""
    rng = random.Random(seed)
    centres = [(2_600, 4_300, 700), (5_600, 1_500, 900), (6_900, 4_700, 500)]
    homes = []
    while len(homes) < n:
        if rng.random() < 0.7:
            cx, cy, spread = rng.choice(centres)
            # The sum of three random numbers clusters around its middle, like a bell curve.
            x = cx + spread * 2 * (rng.random() + rng.random() + rng.random() - 1.5)
            y = cy + spread * 2 * (rng.random() + rng.random() + rng.random() - 1.5)
        else:
            x, y = rng.uniform(0, COLS * WARD), rng.uniform(0, ROWS * WARD)
        if 0 < x < COLS * WARD and 0 < y < ROWS * WARD:
            homes.append((f"H{len(homes) + 1:04d}", round(E0 + x, 1), round(N0 + y, 1)))
    return homes


def clinics():
    """(id, x, y) of the 8 clinics."""
    return [(k, E0 + x, N0 + y) for k, (x, y) in CLINICS.items()]

"""5,000 households tested against 12 wards in one vectorised call, checked against plain arithmetic, and clinic C2,
which stands on the edge of W02 and W03, under four predicates."""
import numpy as np
import shapely
from town import COLS, E0, N0, ROWS, WARD, clinics, households, ward_bounds

homes = households()
xy = np.array([(x, y) for _, x, y in homes])  # a 5,000 x 2 array of coordinates
points = shapely.points(xy)  # 5,000 Point objects, made in one call
wards = np.array([shapely.box(*b) for b in ward_bounds().values()])
names = list(ward_bounds())
print(f"{points.size} points in a NumPy array of dtype {points.dtype}; the first is {points[0]}")

# One call compares every ward with every household: a 12 x 5,000 table of True and False.
inside = shapely.contains(wards[:, np.newaxis], points[np.newaxis, :])
counts = inside.sum(axis=1)
print(f"contains() answered {inside.size:,} questions; table shape {inside.shape}")
print("households per ward:", ", ".join(f"{w} {n}" for w, n in zip(names, counts)))

# The same counts by arithmetic: a household's column and row of the grid follow from its x and y.
col = ((xy[:, 0] - E0) // WARD).astype(int)
row = ROWS - 1 - ((xy[:, 1] - N0) // WARD).astype(int)
by_arithmetic = np.bincount(row * COLS + col, minlength=COLS * ROWS)
print("same as arithmetic:", counts.tolist() == by_arithmetic.tolist(), "| total", int(counts.sum()))

c2 = shapely.Point(next((x, y) for k, x, y in clinics() if k == "C2"))
print(f"\nC2 at {c2}: intersects within touches covered_by")
for name in ["W02", "W03"]:
    ward = wards[names.index(name)]
    answers = [shapely.intersects(c2, ward), shapely.within(c2, ward), shapely.touches(c2, ward), shapely.covered_by(c2, ward)]
    print(f"  against {name}: " + "  ".join(f"{str(a):>5}" for a in answers))
print("wards that contain C2:", int(shapely.contains(wards, c2).sum()), "| wards C2 intersects:", int(shapely.intersects(wards, c2).sum()))

5000 points in a NumPy array of dtype object; the first is POINT (502076.4 2004143.9)
contains() answered 60,000 questions; table shape (12, 5000)
households per ward: W01 295, W02 696, W03 183, W04 1169, W05 204, W06 416, W07 379, W08 340, W09 126, W10 151, W11 648, W12 393
same as arithmetic: True | total 5000

C2 at POINT (504000 2005000): intersects within touches covered_by
  against W02:  True  False   True   True
  against W03:  True  False   True   True
wards that contain C2: 0 | wards C2 intersects: 2

One call answered 60,000 questions, and the counts are the ones arithmetic gives: a household’s column is its distance east of the town’s corner divided by 2,000 m, rounded down, and likewise its row. Two methods that agree are the best evidence that a spatial count is right. The map shows the result:

The last lines test clinic C2, which stands on the edge of W02 and W03, with four spatial predicatesspatial predicate: A test that answers true or false about how two geometries relate, such as intersects, within, contains, touches or covered by. Joins and selections by location are built on them.. It intersects, touches and is covered by both wards, and is within neither, so a count by contains finds it in no ward at all. The lesson on geometry types and validity explains why with the DE-9IMDE-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. More in “Geometry types, WKT, WKB and validity”; the next lessons turn it into a rule for joins.

A buffer is a polygon in the units of its coordinates

A bufferbuffer: The zone within a given distance of a feature, such as the area within 1 km of a clinic. Counting what lies inside a buffer answers many proximity questions. More in “What GIS is and the questions it answers” of a point is not a circle: it is a polygon whose corners lie on the circle, and quad_segs sets “the number of linear segments in a quarter circle” (Shapely, buffer). The two spellings have different defaults. The function shapely.buffer uses 8 segments a quarter, while the method geometry.buffer uses 16 (Shapely manual). The figure draws clinic C2’s 1 km buffer at four settings:

Each doubling of the sides closes most of the gap that is left. The program measures both defaults, then the share of the circle at six settings:

"""A buffer is a polygon that approximates a circle. The method and the function use different numbers of segments
by default, quad_segs sets how many, and the distance is in the units of the coordinates."""
import math
import warnings

import geopandas as gpd
import numpy as np
import shapely
from pyproj import Transformer
from town import clinics

c2 = shapely.Point(next((x, y) for k, x, y in clinics() if k == "C2"))
circle = math.pi * 1_000**2  # the true area of a 1 km circle, in square metres

method = c2.buffer(1_000)  # the geometry's method: quad_segs=16 by default
function = shapely.buffer(c2, 1_000)  # the vectorised function: quad_segs=8 by default
for label, ring in [("c2.buffer(1000)", method), ("shapely.buffer(c2, 1000)", function)]:
    sides = shapely.get_num_coordinates(ring) - 1  # the ring repeats its first position at the end
    print(f"{label:<25} {sides:>3} sides, {ring.area:,.0f} m², {ring.area / circle:.3%} of the circle")

print("\nquad_segs  sides  share of the circle")
for q in [1, 2, 4, 8, 16, 32]:
    ring = shapely.buffer(c2, 1_000, quad_segs=q)
    print(f"{q:>9}  {4 * q:>5}  {ring.area / circle:.5f}")

# The same 'radius' in degrees of longitude and latitude is not a circle on the ground.
to_lonlat = Transformer.from_crs("EPSG:32643", "EPSG:4326", always_xy=True)
to_utm = Transformer.from_crs("EPSG:4326", "EPSG:32643", always_xy=True)
lon, lat = to_lonlat.transform(c2.x, c2.y)
in_degrees = shapely.buffer(shapely.Point(lon, lat), 0.009, quad_segs=16)
on_ground = shapely.transform(in_degrees, lambda c: np.column_stack(to_utm.transform(c[:, 0], c[:, 1])))
x0, y0, x1, y1 = on_ground.bounds
print(f"\nC2 is at {lon:.4f}° E, {lat:.4f}° N. A buffer of 0.009 degrees, measured in metres:")
print(f"  {x1 - x0:,.0f} m west to east, {y1 - y0:,.0f} m south to north, {on_ground.area / 1e6:.3f} km²")

# GeoPandas knows the CRS of its data, so it can warn about the same mistake.
with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    gpd.GeoSeries([shapely.Point(lon, lat)], crs="EPSG:4326").buffer(0.009)
print("GeoPandas warns:", str(caught[0].message).strip())

c2.buffer(1000)            64 sides, 3,136,548 m², 99.839% of the circle
shapely.buffer(c2, 1000)   32 sides, 3,121,445 m², 99.359% of the circle

quad_segs  sides  share of the circle
        1      4  0.63662
        2      8  0.90032
        4     16  0.97450
        8     32  0.99359
       16     64  0.99839
       32    128  0.99960

C2 is at 75.0378° E, 18.1339° N. A buffer of 0.009 degrees, measured in metres:
  1,904 m west to east, 1,991 m south to north, 2.973 km²
GeoPandas warns: Geometry is in a geographic CRS. Results from 'buffer' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

The method’s buffer has 64 sides and 3,136,548 m², the function’s 32 sides and 3,121,445 m², so two programs that “buffer by 1 km” disagree by 15,103 m² for one clinic. Set quad_segs explicitly in any code whose areas are reported. GeoPandas 1.2 deprecated its own name for the setting, resolution, in favour of quad_segs (GeoPandas changelog).

The distance has no unit of its own. Shapely “does not support coordinate system transformations” and assumes every feature is in “the same Cartesian plane” (the manual), so 1,000 means 1,000 metres only because the town is in UTM zone 43NUTM (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”. On longitude and latitude, a buffer of 0.009 is 0.009 degrees: at C2, 18.1° north, that is 1,904 m from west to east and 1,991 m from south to north, an ellipse about 4 % narrower than it is tall. Shapely cannot see the problem; GeoPandas can, because a GeoSeries carries its CRScoordinate 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”, and the last line of the output is its warning.

Unions count overlaps once, and differences measure gaps

The eight clinics’ buffers add up to 25.092 km², but three pairs overlap, so the area they cover is smaller. shapely.union_all merges an array of geometries in one call, overlaps counted once; unary_union is another name for it (Shapely, union_all). From the union, difference gives what one shape leaves of another, and the symmetric differenceThe parts of two shapes that lie in one of them but not in both: their union minus their intersection. gives what lies in exactly one of the two. Look at the union over the town:

The figure shows five separate patches; the largest is a chain of four overlapping circles in the east, C3, C5, C7 and C8. The program measures them:

"""Eight 1 km clinic buffers combined: the union counts each overlap once, one test against the union replaces
eight, and difference and symmetric difference measure what the buffers leave out of the town."""
import itertools

import numpy as np
import shapely
from town import COLS, E0, N0, ROWS, WARD, clinics, households

names = [k for k, _, _ in clinics()]
sites = shapely.points([(x, y) for _, x, y in clinics()])
rings = shapely.buffer(sites, 1_000, quad_segs=16)  # 8 polygons from one call
cover = shapely.union_all(rings)
print(f"8 buffers: {shapely.area(rings).sum() / 1e6:.3f} km² added up, {cover.area / 1e6:.3f} km² as a union")
print(f"the union is a {cover.geom_type} of {shapely.get_num_geometries(cover)} parts")
for (a, ra), (b, rb) in itertools.combinations(zip(names, rings), 2):
    shared = shapely.intersection(ra, rb).area
    if shared > 0:
        print(f"  {a} and {b} overlap by {shared:,.0f} m²")

homes = shapely.points([(x, y) for _, x, y in households()])
covered = shapely.intersects(cover, homes)  # one polygon against 5,000 points
near = shapely.dwithin(homes[:, np.newaxis], sites[np.newaxis, :], 1_000).any(axis=1)  # true distances
print(f"\nhouseholds inside the union: {int(covered.sum()):,}; within 1,000 m of a clinic: {int(near.sum()):,}")

town = shapely.box(E0, N0, E0 + COLS * WARD, N0 + ROWS * WARD)
uncovered = shapely.difference(town, cover)
outside = shapely.difference(cover, town)
either = shapely.symmetric_difference(town, cover)
print(f"town {town.area / 1e6:.3f} km²; not covered {uncovered.area / 1e6:.3f} km²; buffers beyond the town {outside.area / 1e6:.3f} km²")
print(f"symmetric difference {either.area / 1e6:.3f} km² = the two together")

8 buffers: 25.092 km² added up, 24.695 km² as a union
the union is a MultiPolygon of 5 parts
  C3 and C5 overlap by 34,498 m²
  C5 and C7 overlap by 267,631 m²
  C7 and C8 overlap by 95,095 m²

households inside the union: 2,971; within 1,000 m of a clinic: 2,973
town 48.000 km²; not covered 24.029 km²; buffers beyond the town 0.724 km²
symmetric difference 24.753 km² = the two together

The union is 24.695 km², a MultiPolygon of 5 parts: the overlaps of C3 and C5, C5 and C7, and C7 and C8 add up to the 0.397 km² that the plain sum counted twice. Uncovered town and buffers beyond the town add up to the symmetric difference, 24.029 + 0.724 = 24.753 km².

The household counts hold a lesson of their own. 2,971 households lie inside the union, but 2,973 lie within 1,000 m of a clinic, measured as distances with shapely.dwithin, which the documentation calls “more efficient than computing the distance and comparing the result” (Shapely, dwithin). The two missing households stand in the slivers between a 64-sided polygon and its circle.

Repair the input first, then snap it to a grid

Set operations need valid input, and they are exact to the last bit of every coordinate. Both facts bite in real data. A ring that crosses itself is invalidvalid 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”, and make_valid repairs it (Shapely, make_valid). Two parcels that share an edge, digitised by two people, can miss each other by less than a millimetre, and the union then keeps them apart. A precision gridprecision grid: A grid of fixed spacing, such as 1 cm, onto which coordinates are snapped before and after an operation, so that positions closer than the spacing become one position and nearly meeting edges meet. fixes that: with grid_size, input coordinates are snapped to a grid of that size and so are the results (Shapely, intersection). Follow the lesson’s three shapes through the checks:

The program runs the same steps and adds the opposite case, a parcel that reaches 0.3 mm into its neighbour:

"""Two repairs a pipeline needs: make_valid for a ring that crosses itself, and a precision grid for an edge that was
digitised twice and misses itself by less than a millimetre."""
import shapely
from shapely import Polygon


def describe(g):
    parts = shapely.get_num_geometries(g)
    return f"a {g.geom_type} of {parts} parts" if parts > 1 else f"one {g.geom_type}"


bowtie = Polygon([(502_000, 2_004_000), (502_002, 2_004_002), (502_002, 2_004_000), (502_000, 2_004_002)])
fixed = shapely.make_valid(bowtie)
print(f"bow-tie: {shapely.is_valid_reason(bowtie)}, area {bowtie.area:g} m²")
print(f"make_valid: {describe(fixed)}, area {fixed.area:g} m²")

# Two parcels of 500 m by 500 m that share an edge at x = 502,500. Parcel A was digitised 0.4 mm short of it.
a = Polygon([(502_000, 2_004_000), (502_499.9996, 2_004_000), (502_499.9996, 2_004_500), (502_000, 2_004_500)])
b = Polygon([(502_500, 2_004_000), (503_000, 2_004_000), (503_000, 2_004_500), (502_500, 2_004_500)])
print()
for grid in [None, 0.01]:
    u = shapely.union(a, b, grid_size=grid)
    print(f"union of A and B, grid_size={grid}: {describe(u)}, {u.area:,.1f} m²")

snapped = shapely.set_precision(a, 0.01)  # every position moves to the nearest centimetre, and stays on that grid
print(f"set_precision(A, 0.01): east edge at x = {snapped.bounds[2]:,.2f}, precision {shapely.get_precision(snapped)}")
print(f"then a plain union with B: {describe(shapely.union(snapped, b))}")

# The other way round: parcel C reaches 0.3 mm into B, and their overlap is a sliver.
c = Polygon([(502_000, 2_004_000), (502_500.0003, 2_004_000), (502_500.0003, 2_004_500), (502_000, 2_004_500)])
sliver = shapely.intersection(c, b)
print(f"\nintersection of C and B: {describe(sliver)} of {sliver.area:.2f} m²")
print(f"the same on a 1 cm grid: {shapely.intersection(c, b, grid_size=0.01)}")

bow-tie: Self-intersection[502001 2004001], area 0 m²
make_valid: a MultiPolygon of 2 parts, area 2 m²

union of A and B, grid_size=None: a MultiPolygon of 2 parts, 499,999.8 m²
union of A and B, grid_size=0.01: one Polygon, 500,000.0 m²
set_precision(A, 0.01): east edge at x = 502,500.00, precision 0.01
then a plain union with B: one Polygon

intersection of C and B: one Polygon of 0.15 m²
the same on a 1 cm grid: LINESTRING (502500 2004000, 502500 2004500)

Without a grid the union of A and B is a MultiPolygon of 2 parts, 499,999.8 m², split by a gap no map shows; on a 1 cm grid it is one Polygon of 500,000 m². set_precision snaps a geometry once and keeps the grid, so later operations run “in the precision of the geometry with higher precision” (Shapely, set_precision). The 0.15 m² overlap of C and B shrinks to a line on the grid, because a 0.3 mm strip has no width at 1 cm.

GEOS does the work, so its version is part of the result

GEOS computes every predicate, buffer and union above, and its version travels with each installation of Shapely. The wheels of Shapely 2.1.x bundle GEOS 3.13.1 and those of Shapely 2.2 bundle GEOS 3.14.1 (Shapely 2.x release notes), the browser build of Shapely 2.1.2 in Pyodide 314.0.7 links GEOS 3.12.1 (Pyodide package lock), and GEOS itself has reached 3.15.0 (GEOS releases). Walk through the four:

A new version may return the same shape with its vertices in another order, keep an extra vertex on a straight edge, or round a computed crossing in its last digit. A test that compares WKT text then fails although nothing is wrong. This program shows both kinds of difference and the comparisons that see through them:

"""Compare geometries by shape, not by text. The same ward written from another corner has other WKT, the union of 12
wards keeps vertices that a box of the town does not have, and the GEOS version belongs beside every result."""
import shapely
from shapely import Polygon
from town import COLS, E0, N0, ROWS, WARD, ward_bounds

print(f"Shapely {shapely.__version__} with GEOS {shapely.geos_version_string}")

w02 = shapely.box(*ward_bounds()["W02"])
again = Polygon([(504_000, 2_006_000), (504_000, 2_004_000), (502_000, 2_004_000), (502_000, 2_006_000)])
print(f"\nW02 as box():     {w02.wkt}")
print(f"W02 typed again:  {again.wkt}")
print(f"same WKT: {w02.wkt == again.wkt} | equals: {shapely.equals(w02, again)} | "
      f"same after normalize: {shapely.normalize(w02).wkt == shapely.normalize(again).wkt}")

town = shapely.union_all([shapely.box(*b) for b in ward_bounds().values()])
outline = shapely.box(E0, N0, E0 + COLS * WARD, N0 + ROWS * WARD)
print(f"\nunion of the 12 wards: {shapely.get_num_coordinates(town)} positions; box of the town: {shapely.get_num_coordinates(outline)}")
print(f"equals: {shapely.equals(town, outline)} | equals_exact: {shapely.equals_exact(town, outline, tolerance=0)} | "
      f"area difference: {abs(town.area - outline.area)} m²")
tidy = shapely.simplify(town, 0)  # drops the positions that lie on a straight edge
print(f"after simplify(0): {shapely.get_num_coordinates(tidy)} positions, same normalized WKT: "
      f"{shapely.normalize(tidy).wkt == shapely.normalize(outline).wkt}")

Shapely 2.1.2 with GEOS 3.13.1

W02 as box():     POLYGON ((504000 2004000, 504000 2006000, 502000 2006000, 502000 2004000, 504000 2004000))
W02 typed again:  POLYGON ((504000 2006000, 504000 2004000, 502000 2004000, 502000 2006000, 504000 2006000))
same WKT: False | equals: True | same after normalize: True

union of the 12 wards: 15 positions; box of the town: 5
equals: True | equals_exact: False | area difference: 0.0 m²
after simplify(0): 5 positions, same normalized WKT: True

W02 typed from another corner has other WKT, but equals says the shapes are the same and normalize gives both the same text; the documentation calls normalize “typically useful for testing purposes” (Shapely, normalize). The union of the 12 wards keeps 15 positions where the town’s box has 5, so equals_exact, which requires the coordinates “in the same order” (Shapely, equals_exact), says no; after simplify(0) the normalized texts match.

Practise: three helpers of a geometry toolkit

Build three small helpers in the browser, each with sample tests: one that predicts how much of the circle a buffer covers, one that snaps a ring to a precision grid, and one that compares two rings by shape:

import math


def buffer_share(quad_segs):
    """The area of a buffer polygon with quad_segs sides per quarter circle, as a share of the circle's area."""
    # Replace this line: n = 4 * quad_segs sides, an area of n * r**2 * sin(2 * pi / n) / 2 against pi * r**2.
    return 0.0


def quad_segs_for(share):
    """The smallest quad_segs whose buffer polygon covers at least `share` of the circle."""
    # Replace this line: try 1, 2, 3 … until buffer_share reaches the share.
    return 0

import math

from sides import buffer_share, quad_segs_for


def test_one_segment_is_a_diamond():
    """gives 2 / pi for quad_segs=1, a four-sided diamond"""
    assert math.isclose(buffer_share(1), 2 / math.pi, rel_tol=1e-12)


def test_the_measured_shares():
    """matches the shares Shapely's buffers of C2 measured, to five decimals"""
    measured = {2: 0.90032, 4: 0.9745, 8: 0.99359, 16: 0.99839, 32: 0.9996}
    for q, share in measured.items():
        assert round(buffer_share(q), 5) == share


def test_more_sides_cover_more():
    """grows with every extra segment and stays below the whole circle"""
    shares = [buffer_share(q) for q in range(1, 41)]
    assert all(a < b for a, b in zip(shares, shares[1:]))
    assert shares[-1] < 1


def test_segments_for_a_share():
    """finds the smallest quad_segs for a share: 1 for half, 7 for 99 %, 21 for 99.9 %"""
    assert quad_segs_for(0.5) == 1
    assert quad_segs_for(0.99) == 7
    assert quad_segs_for(0.999) == 21
    assert buffer_share(quad_segs_for(0.995)) >= 0.995 > buffer_share(quad_segs_for(0.995) - 1)

def snap_ring(ring, grid):
    """The ring with every position moved to the nearest multiple of grid, and repeated positions dropped."""
    # Replace this line: snap x and y, round to 9 decimals, skip a position equal to the previous one.
    return list(ring)

from snap import snap_ring

PARCEL_A = [(502000, 2004000), (502499.9996, 2004000), (502499.9996, 2004500), (502000, 2004500), (502000, 2004000)]


def test_the_short_edge_reaches_the_neighbour():
    """moves parcel A's east edge from 502,499.9996 to 502,500 on a 1 cm grid"""
    snapped = snap_ring(PARCEL_A, 0.01)
    assert snapped == [(502000.0, 2004000.0), (502500.0, 2004000.0), (502500.0, 2004500.0), (502000.0, 2004500.0), (502000.0, 2004000.0)]


def test_positions_that_meet_become_one():
    """drops a position 3 mm from the previous one, because both land on the same centimetre"""
    ring = [(0, 0), (10, 0), (10.003, 0.002), (10, 10), (0, 10), (0, 0)]
    assert snap_ring(ring, 0.01) == [(0.0, 0.0), (10.0, 0.0), (10.0, 10.0), (0.0, 10.0), (0.0, 0.0)]


def test_halves_and_coarse_grids():
    """rounds to the nearest multiple of a 1 m grid and keeps the ring closed"""
    ring = [(0.4, 0.6), (5.2, 0.4), (5.6, 4.7), (0.4, 0.6)]
    assert snap_ring(ring, 1.0) == [(0.0, 1.0), (5.0, 0.0), (6.0, 5.0), (0.0, 1.0)]


def test_a_fine_grid_changes_nothing():
    """leaves positions alone when they already sit on the grid"""
    ring = [(1.25, 2.5), (3.75, 2.5), (3.75, 5.0), (1.25, 2.5)]
    assert snap_ring(ring, 0.25) == ring

def same_ring(a, b):
    """True when two closed rings visit the same positions in the same cyclic order, in either direction."""
    # Replace this line: compare a without its last position with each rotation of b, and of b reversed.
    return a == b

from rings import same_ring

W02_BOX = [(504000, 2004000), (504000, 2006000), (502000, 2006000), (502000, 2004000), (504000, 2004000)]
W02_TYPED = [(504000, 2006000), (504000, 2004000), (502000, 2004000), (502000, 2006000), (504000, 2006000)]


def test_another_corner_and_direction():
    """accepts W02 typed from another corner, running the other way"""
    assert same_ring(W02_BOX, W02_TYPED)
    assert same_ring(W02_TYPED, W02_BOX)


def test_every_start_and_direction():
    """accepts every rotation of a ring, forwards and backwards"""
    core = W02_BOX[:-1]
    for i in range(4):
        turned = core[i:] + core[:i]
        assert same_ring(W02_BOX, turned + turned[:1])
        back = turned[::-1]
        assert same_ring(W02_BOX, back + back[:1])


def test_other_shapes():
    """refuses W03, and W02's corners visited in a crossing order"""
    w03 = [(506000, 2004000), (506000, 2006000), (504000, 2006000), (504000, 2004000), (506000, 2004000)]
    crossing = [(504000, 2004000), (502000, 2006000), (504000, 2006000), (502000, 2004000), (504000, 2004000)]
    assert not same_ring(W02_BOX, w03)
    assert not same_ring(W02_BOX, crossing)


def test_an_extra_position_needs_simplify():
    """refuses a ring with an extra position on a straight edge, as the union of two wards has"""
    extra = [(504000, 2004000), (504000, 2005000), (504000, 2006000), (502000, 2006000), (502000, 2004000), (504000, 2004000)]
    assert not same_ring(W02_BOX, extra)

The self-check below covers all five ideas of the lesson and asks first about the one you have answered least well:

Interview questions

Warm-up (fresher): Shapely has no coordinate reference system. What does that mean for point.buffer(1000)?

Shapely works in a flat plane and treats coordinates as plain numbers, so the 1,000 is in whatever units the coordinates use. In metres of a projected CRS such as UTM zone 43N it is a 1 km buffer. On longitude and latitude it is 1,000 degrees, nonsense, and even a small value such as 0.009 makes an ellipse: 1,904 m wide and 1,991 m tall at 18° north. The fix is to transform the data to a projected CRS in metres, buffer there, and transform back only for display. GeoPandas, which does know the CRS, warns when you buffer in a geographic one. Mention too that the buffer is a polygon: 64 sides by default for the method, 32 for shapely.buffer.

Two harder questions continue from here: a coverage count that disagrees with a distance count, and a geometry test that fails after an upgrade. The module’s question bank says, for each, what it tests and how to begin:

Key takeaways

• Shapely 2’s functions take whole NumPy arrays of geometries and loop in C, so one call such as shapely.contains on a 12 × 5,000 table replaces a Python loop of 60,000 steps.
• A buffer is a polygon with its corners on the circle, 64 sides from the method and 32 from the function by default, in the units of the coordinates; buffer only in a projected CRS.
• union_all counts overlaps once, difference and symmetric_difference measure what is left, and a distance test is more exact than a test against a buffer polygon.
• Repair invalid input with make_valid and join edges that nearly meet with grid_size or set_precision, choosing a grid finer than the survey’s accuracy.
• GEOS does the geometry, and its version differs between wheels, browsers and newer releases, so record it and compare results by shape, not by text.