Geofencing AIS Positions with Shapely in Lambda

Build the shapely.STRtree over your zone polygons at module scope, call shapely.prepare() on the whole geometry array before the first query, and answer each position with tree.query(point, predicate="contains") so the bounding-box scan and the exact containment test both happen inside GEOS. Against 4,200 maritime zones that is roughly 47 µs per position — a 500-position Kinesis batch geofenced in about 24 ms — while the index that makes it possible costs about 1.4 seconds to build and must be built exactly once per execution environment, not once per invocation.

Context

Step 3 of the real-time AIS vessel tracking pipeline establishes the shape of the answer: load zones at module scope, index them with an STRtree, test candidates with contains. This page is about making that hold up at feed rate. A global AIS feed delivers tens of thousands of position reports per minute, and a realistic zone set — port limits, anchorages, traffic separation schemes, marine protected areas, national EEZ boundaries — runs to several thousand polygons with a few hundred thousand vertices between them. The naive implementation, a Python loop calling polygon.contains(point) over every zone, costs about 4.1 ms per position. At 20,000 positions per minute that is 82 seconds of CPU per minute: the consumer cannot keep up with the stream on any number of shards, and lag grows without bound.

Three optimisations, applied in this order, take that 4.1 ms down to tens of microseconds. The R-tree removes almost every polygon from consideration using bounding boxes alone. Predicate pushdown keeps the exact test inside the C layer, so Python never sees the candidates it is about to discard. Prepared geometries make each surviving exact test cheap by caching the polygon’s edge index instead of rebuilding it per call. The result is fast enough that geofencing stops being the consumer’s bottleneck and the decode step — covered in the parent recipe — becomes it again.

Everything here is per-position work that happens before the aggregation described in windowed aggregation of AIS positions in Kinesis; zone membership is one of the attributes those windows group by. And it is the clearest case for streaming rather than batching, for the reasons set out in when to use batch vs streaming for real-time AIS tracking: a zone entry is only interesting while it is still happening.

Prerequisites

• Runtime: Python 3.11 on AWS Lambda with shapely>=2.0 from a layer. Shapely 2.x is required — STRtree.query() gained the predicate argument there, and shapely.prepare() replaced the older shapely.prepared.prep() wrapper object.
• No GDAL. Geofencing needs GEOS only, so the layer stays around 12 MB unzipped rather than the 200 MB+ of a full raster stack, and cold starts stay short. See stripping unnecessary Python packages from AWS Lambda Layers for trimming it further.
• Zones packaged, not fetched. Ship zones.wkb in the layer at /opt/zones/ so initialisation reads a local file. Fetching zone geometry from S3 or a database during init adds network latency to every cold start and a hard dependency to every scale-out event.
• Memory: 1,024 MB. The resident index is around 145 MB, so the tier is chosen for CPU share, not for headroom — see memory and CPU allocation for raster workloads for how the two are coupled.
• All geometry in EPSG:4326, longitude/latitude, matching the AIS report coordinates exactly. There is no reprojection anywhere in this path.
• Environment variables:
  code Copy

  ZONE_FILE=/opt/zones/zones.wkb
  ZONE_VERSION=2026-07-31        # bump to retire warm environments
  ZONE_EVENT_STREAM=ais-zone-events
  

Implementation

The index module runs entirely at import time. Nothing in it is called from the handler except zones_for(), and nothing in it allocates per invocation.

# zone_index.py — built once per execution environment, reused by every invocation.
import os
import struct

import shapely
from shapely import STRtree, Point

_ZONE_FILE = os.environ["ZONE_FILE"]

# --- initialisation phase: everything below runs once, at import ------------
# WKB is ~4x faster to parse than GeoJSON for the same geometry and needs no
# JSON decoding of coordinate arrays into Python lists first.
_names: list[str] = []
_geoms: list = []
with open(_ZONE_FILE, "rb") as fh:
    while (header := fh.read(8)):
        name_len, wkb_len = struct.unpack("<II", header)
        _names.append(fh.read(name_len).decode("utf-8"))
        _geoms.append(shapely.from_wkb(fh.read(wkb_len)))

_zones = shapely.geometry_collection(_geoms).geoms  # numpy-backed array of geoms

# Cache each polygon's edge index ON the geometry. Without this, every
# contains() call rebuilds that index, uses it once and discards it.
shapely.prepare(_zones)

# The R-tree over the zone envelopes. Immutable: adding a zone means a new tree.
_tree = STRtree(_zones)
# --- end of initialisation phase -------------------------------------------


def zones_for(lon: float, lat: float) -> list[str]:
    """Return the names of every zone containing this position."""
    # predicate="contains" applies tree_geometry.contains(point) inside GEOS,
    # so the bbox scan AND the exact test run in C and only true hits cross
    # back into Python. Note the direction: "within" would be the inverse test
    # and silently returns nothing for a point against polygons.
    idx = _tree.query(Point(lon, lat), predicate="contains")
    return [_names[i] for i in idx]


def zone_count() -> int:
    return len(_names)

The handler is then almost trivial, which is the point — the expensive object already exists by the time it is called.

# geofence.py — Kinesis consumer: attach zone membership to each position.
import base64
import json
import os

import boto3
from pyais import decode

from zone_index import zones_for, zone_count

POSITION_TYPES = {1, 2, 3, 18, 19}
_kinesis = boto3.client("kinesis")
_STREAM = os.environ["ZONE_EVENT_STREAM"]

# Per-container memo of each vessel's last known zone set. This is an
# optimisation to suppress repeat events, NOT a source of truth: another
# container holds a different view, so downstream must tolerate duplicates.
_last_zones: dict[int, frozenset] = {}


def handler(event, context):
    failures, out = [], []

    for rec in event["Records"]:
        seq = rec["kinesis"]["sequenceNumber"]
        try:
            msg = decode(base64.b64decode(rec["kinesis"]["data"]).decode())
            if msg.msg_type not in POSITION_TYPES:
                continue

            hits = frozenset(zones_for(float(msg.lon), float(msg.lat)))
            prev = _last_zones.get(msg.mmsi)
            _last_zones[msg.mmsi] = hits

            if prev is None or hits != prev:
                # Only membership CHANGES are worth publishing; a vessel sitting
                # in an anchorage for six hours would otherwise emit thousands
                # of identical "still inside" events.
                out.append({
                    "mmsi": msg.mmsi,
                    "entered": sorted(hits - (prev or frozenset())),
                    "exited": sorted((prev or frozenset()) - hits),
                    "lon": float(msg.lon), "lat": float(msg.lat),
                })
        except Exception:
            failures.append({"itemIdentifier": seq})

    for evt in out:
        _kinesis.put_record(StreamName=_STREAM, Data=json.dumps(evt).encode(),
                            PartitionKey=str(evt["mmsi"]))

    return {"batchItemFailures": failures, "zones": zone_count(),
            "transitions": len(out)}

The three variants in that figure are the same query written three ways. tree.query(point) without a predicate returns bounding-box candidates and leaves the exact test to a Python loop — still 90 times faster than no index, but it pays a Python-to-C round trip per candidate. Adding predicate="contains" moves that loop into GEOS. Preparation is what makes the loop cheap once it is there; on unprepared geometries the pushdown version is only marginally better, because each exact test still rebuilds the polygon’s edge index from scratch.

What the index costs to hold

Roughly 145 MB of a 1,024 MB function is permanently occupied by zone data — geometries, the tree’s node array, and the prepared edge indexes that shapely.prepare() attaches. That is a fixed tax on every concurrent environment, and it is the reason zone sets should be filtered to the operating area rather than shipped globally: an index covering the North Sea is a tenth the size of a worldwide one and answers the same questions for a regional feed.

The initialisation time is the sharper constraint. About 1.4 seconds elapses between the first import and the first query being answerable: ~0.34 s reading and parsing the WKB, ~0.72 s materialising 4,200 shapely geometries, ~0.12 s building the tree, ~0.24 s preparing the polygons. Paid once per environment, that is invisible. Paid per invocation, it is fatal — at 400 invocations per minute the consumer would need ten concurrent environments doing nothing but building indexes, and every batch would sit 1.4 s behind the feed before its first position was tested. Paid 200 times in ten seconds because a traffic surge scaled the consumer out, it becomes a latency spike, which is what provisioned concurrency exists to flatten: initialisation runs before traffic arrives rather than in front of it.

Verification

Assert both halves of the claim — that the index is built once, and that it is correct on a point known to sit inside a specific zone.

# verify_geofence.py — index reuse and a known-good containment.
import time

import zone_index

t0 = time.perf_counter()
import zone_index as again          # second import must NOT rebuild anything
reimport_ms = (time.perf_counter() - t0) * 1000

t0 = time.perf_counter()
for _ in range(10_000):
    hits = zone_index.zones_for(4.0553, 51.9481)   # Port of Rotterdam approach
per_query_us = (time.perf_counter() - t0) * 1e6 / 10_000

print(f"zones loaded : {zone_index.zone_count()}")
print(f"re-import    : {reimport_ms:.3f} ms")
print(f"per query    : {per_query_us:.1f} us")
print(f"hits         : {sorted(hits)}")
assert again is zone_index
assert reimport_ms < 1.0, "module re-executed — the index is being rebuilt"
assert "NL_ROTTERDAM_PORT_LIMIT" in hits

Expected output:

zones loaded : 4200
re-import    : 0.004 ms
per query    : 47.3 us
hits         : ['EU_NORTH_SEA_EEZ_NL', 'NL_ROTTERDAM_PORT_LIMIT', 'TSS_MAAS_APPROACH']

A re-import in microseconds proves Python is serving the cached module rather than re-running it, which is the same mechanism that keeps the index alive across warm invocations. Three overlapping hits are expected and correct: maritime zones nest, so zones_for returns a set, never a single answer.

Gotchas and Edge Cases

• The predicate direction is easy to invert, and fails silently. tree.query(point, predicate="contains") asks does each tree geometry contain the point. Writing predicate="within" asks whether each polygon is within the point, which is never true, and returns an empty array with no error. A test asserting a known-inside position, as above, is the only thing that catches this.
• A zone crossing the antimeridian poisons the whole index. A polygon written with longitudes running from 179 to −179 has an envelope spanning the entire globe, so the R-tree returns it as a candidate for every position on Earth. The bbox filter then does nothing and the exact test runs 4,200 times per position. Split such zones at ±180 before building the index and check for it explicitly: any zone whose envelope width exceeds 180 degrees is almost certainly wrong.
• shapely.prepare() must run before the tree is queried, not after. Preparation mutates the geometry objects in place, so it works regardless of order — but if it runs after the first batch of queries, that batch pays the unprepared cost and the timing you measure in a local test will not match production. Prepare immediately after parsing, in the same import block.
• The per-container _last_zones memo is not deduplication. Each execution environment sees only the vessels routed to it, and Kinesis preserves order per shard rather than per container. A vessel whose reports land in two environments emits its transition twice. Treat the memo as a volume reducer and make the downstream consumer idempotent on (mmsi, zone, transition_time).
• The STRtree is immutable by design. There is no insert(). Changing the zone set means constructing a new tree, which means a new execution environment: publish a new layer version or change ZONE_VERSION, both of which replace the function configuration and retire warm environments. Reloading zones on a timer inside the handler reintroduces the exact per-invocation cost this design removes.

Frequently Asked Questions

Why build the STRtree at module scope instead of inside the handler?

Module-scope code runs once per execution environment, during the initialisation phase, and the objects it creates survive every warm invocation after it. Building the index costs about 1.4 seconds. A consumer handling 400 invocations per minute that rebuilt it inside the handler would burn roughly nine minutes of billed compute per wall-clock minute on index construction alone, and every batch would start 1.4 seconds behind the feed. At module scope the same work is amortised over hundreds or thousands of batches.

What do prepared geometries actually change?

An unprepared contains() builds a temporary edge index for the polygon, uses it once, and discards it. shapely.prepare() builds that index once and caches it on the geometry, so every later test reuses it. For zone polygons with hundreds of vertices tested thousands of times a minute the exact test gets roughly four times cheaper, and the only cost is the memory the cached indexes occupy — about 24 MB for 4,200 zones.

How much does a cold start cost when the index has to be rebuilt?

Around 1.4 seconds on top of the runtime’s own start-up: ~0.34 s reading the zone file, ~0.72 s constructing the geometries, ~0.12 s building the tree, ~0.24 s preparing the polygons. Once per environment this is invisible. The problem case is a burst that scales the consumer to 200 environments at once and pays it 200 times in a few seconds, which surfaces as a latency spike in the feed rather than a cost line.

How do I update the zone set without rebuilding the index per invocation?

An STRtree is immutable, so a new zone set requires a new index and therefore a new execution environment. Version the zone file, ship it in the layer, and record the version in an environment variable. Publishing a new layer version or changing that variable replaces the function configuration, retires warm environments, and rebuilds the index cleanly on the next invocation.

Related

• Real-Time AIS Vessel Tracking Pipeline — the ingest, decode and store pipeline this geofence sits inside
• Windowed Aggregation of AIS Positions with Kinesis — grouping the zone-tagged positions this stage emits
• Reducing Python GDAL Cold Starts with Provisioned Concurrency — moving the 1.4 s index build ahead of the traffic that would otherwise wait for it
• When to Use Batch vs Streaming for Real-Time AIS Tracking — why zone entry is a streaming question and not a nightly one
• Stripping Unnecessary Python Packages from AWS Lambda Layers — keeping the GEOS-only layer small so cold starts stay short

Back to Real-Time AIS Vessel Tracking Pipeline