Spatial Aggregations With Hybrid Expressions

Building spatial aggregations on top of Hybrid Properties for Geometry is what keeps a vessel-tracking API from dragging thousands of ais_positions rows into Python just to compute one convex track or bounding envelope. This page shows how to pair @hybrid_property and hybrid_method with ST_Collect, ST_Union, ST_Extent, and ST_Centroid so every aggregation is compiled into the SQL SELECT list and executed by PostGIS. The scenario throughout is a fleet service that must return, per vessel, a merged track, a centroid, and a bounding box — all in SRID 4326.

Why the naive approach fails

The obvious implementation fetches every position for a vessel and folds the geometries in Python. It is correct but ruinous: it ships the entire point history over the wire, deserialises each WKB blob into a Shapely object, and rebuilds the aggregate that PostGIS could have produced in a single pass.

python
# SLOW: pulls every ais_positions row into Python and aggregates client-side
from geoalchemy2.shape import to_shape
from shapely.ops import unary_union

def vessel_track_bbox(session, mmsi: int):
    rows = session.scalars(
        select(AisPosition).where(AisPosition.mmsi == mmsi)
    ).all()                                   # could be 100k+ WKBElements
    points = [to_shape(r.geom) for r in rows] # 100k Shapely objects in memory
    merged = unary_union(points)              # CPU-bound union in Python
    return merged.bounds                      # (minx, miny, maxx, maxy)

For a vessel with a month of one-second AIS fixes this materialises millions of objects and can exhaust memory before it ever returns the four numbers the caller wanted. The database can compute the same envelope with ST_Extent while touching only the index and the heap once.

Production-ready implementation

The fix is to give each derived geometry a hybrid definition with an @expression companion, so the same attribute name works as a Python fallback on a loaded instance and as pushed-down SQL inside a query. Here is the complete, copy-paste model plus the queries that aggregate in the database.

python
from datetime import datetime
from sqlalchemy import Integer, BigInteger, String, DateTime, ForeignKey, Index, func, select
from sqlalchemy.orm import DeclarativeBase, Mapped, mapped_column, relationship
from sqlalchemy.ext.hybrid import hybrid_property, hybrid_method
from geoalchemy2 import Geometry
from geoalchemy2.shape import to_shape


class Base(DeclarativeBase):
    pass


class Vessel(Base):
    __tablename__ = "vessels"
    id:   Mapped[int] = mapped_column(Integer, primary_key=True)
    mmsi: Mapped[int] = mapped_column(BigInteger, unique=True, index=True)
    name: Mapped[str] = mapped_column(String(120))
    positions: Mapped[list["AisPosition"]] = relationship(back_populates="vessel")


class AisPosition(Base):
    __tablename__ = "ais_positions"
    id:   Mapped[int] = mapped_column(Integer, primary_key=True)
    mmsi: Mapped[int] = mapped_column(BigInteger, ForeignKey("vessels.mmsi"), index=True)
    recorded_at: Mapped[datetime] = mapped_column(DateTime(timezone=True), index=True)
    geom: Mapped[object] = mapped_column(
        Geometry(geometry_type="POINT", srid=4326, spatial_index=False)
    )
    vessel: Mapped["Vessel"] = relationship(back_populates="positions")

    __table_args__ = (
        Index("ais_positions_geom_gist", "geom", postgresql_using="gist"),
    )

    # ── centroid: no arguments, so hybrid_property fits ──────────────────
    @hybrid_property
    def centroid(self):
        # Instance fallback: only valid on a single loaded point
        return to_shape(self.geom).centroid

    @centroid.expression
    def centroid(cls):
        # SQL side: aggregate centroid of all matched points, SRID preserved
        return func.ST_Centroid(func.ST_Collect(cls.geom)).label("centroid")

    # ── bounding envelope of a group of points ───────────────────────────
    @hybrid_property
    def track_extent(self):
        return to_shape(self.geom).bounds

    @track_extent.expression
    def track_extent(cls):
        # ST_Extent returns a box2d; wrap in ST_SetSRID(ST_Envelope(...)) for a
        # real geometry the ORM can decode as a polygon in SRID 4326.
        return func.ST_SetSRID(
            func.ST_Envelope(func.ST_Collect(cls.geom)), 4326
        ).label("track_extent")

    # ── parameterised merge: needs a target SRID, so hybrid_method ───────
    @hybrid_method
    def merged_track(cls, srid: int = 4326):
        # ST_Collect bundles points into a MultiPoint without dissolving;
        # optionally reproject the collection into a projected CRS for length.
        collected = func.ST_Collect(cls.geom)
        if srid != 4326:
            collected = func.ST_Transform(collected, srid)
        return collected.label("merged_track")

With the model in place, the aggregations run entirely in PostGIS. Each query returns one aggregate geometry per vessel via GROUP BY:

python
def fleet_centroids(session):
    """One centroid geometry per vessel, computed in the database."""
    stmt = (
        select(
            AisPosition.mmsi,
            AisPosition.centroid,          # -> ST_Centroid(ST_Collect(geom))
        )
        .group_by(AisPosition.mmsi)
    )
    for mmsi, centroid_wkb in session.execute(stmt):
        pt = to_shape(centroid_wkb)        # WKBElement -> shapely Point (4326)
        yield mmsi, (round(pt.x, 5), round(pt.y, 5))


def fleet_bounding_boxes(session):
    """One envelope polygon per vessel, computed in the database."""
    stmt = (
        select(AisPosition.mmsi, AisPosition.track_extent)
        .group_by(AisPosition.mmsi)
    )
    return {
        mmsi: to_shape(env).bounds
        for mmsi, env in session.execute(stmt)
    }


def merged_track_utm(session, mmsi: int, utm_srid: int = 32631):
    """Merge one vessel's points into a MultiPoint reprojected to UTM 31N."""
    stmt = (
        select(AisPosition.merged_track(utm_srid))
        .where(AisPosition.mmsi == mmsi)
    )
    return to_shape(session.execute(stmt).scalar_one())

ST_Collect is the workhorse for gathering AIS points because it only bundles them into a MultiPoint — it never dissolves geometry, so it stays cheap even for enormous point sets. Reserve ST_Union for the case where dissolved boundaries genuinely matter, such as merging overlapping port_zones polygons into a single coverage footprint:

python
def dissolved_port_coverage(session, srid: int = 4326):
    """True geometric union of overlapping port zones — boundaries dissolved."""
    stmt = select(func.ST_Union(PortZone.geom).label("coverage"))
    return to_shape(session.execute(stmt).scalar_one())

Because UTM zone 31N (EPSG:32631) is a projected, metre-based CRS, the reprojected merged_track can feed a downstream ST_Length for track distance without the degree-based distortion you would get by measuring in raw SRID 4326.

Configuration and tuning knobs

  • work_mem — ST_Collect and ST_Union accumulate geometry in the aggregate state. Set SET LOCAL work_mem = '128MB' for fleet-wide merges so the aggregate does not spill to a temp file.
  • spatial_index and the GiST index — the explicit ais_positions_geom_gist index keeps any WHERE predicate (for example an ST_Intersects against a port zone) index-aware even while the SELECT list aggregates. Aggregation itself is a heap scan, but pre-filtering the group is not.
  • ST_Extent vs ST_Envelope — ST_Extent returns a box2d with no SRID, which the ORM cannot decode as a geometry. Wrap it (or use ST_Envelope(ST_Collect(...)) as above) and re-attach SRID 4326 with ST_SetSRID so to_shape succeeds.
  • GROUP BY cardinality — grouping by mmsi produces one row per vessel; make sure mmsi is indexed so the grouping does not force a full sort of the position table.
  • Statement timeout — a fleet-wide ST_Union over millions of points is analytical, not transactional. Route it to a replica or raise statement_timeout for that session so it does not trip the OLTP budget.

Verification steps

Confirm the aggregation is compiled into SQL rather than executed in Python by printing the compiled statement before it runs:

python
stmt = (
    select(AisPosition.mmsi, AisPosition.centroid)
    .group_by(AisPosition.mmsi)
)
print(stmt.compile(compile_kwargs={"literal_binds": True}))
# Expect ST_Centroid(ST_Collect(ais_positions.geom)) in the SELECT list,
# NOT a bare ais_positions.geom that would force client-side aggregation.

Then confirm the plan and the numeric result at the database:

sql
-- The aggregate should scan the index/heap once, not per group-member
EXPLAIN (ANALYZE, BUFFERS)
SELECT mmsi,
       ST_AsText(ST_Centroid(ST_Collect(geom))) AS centroid
FROM   ais_positions
WHERE  mmsi = 244660000
GROUP  BY mmsi;

-- Sanity-check the envelope SRID survived the round-trip
SELECT ST_SRID(ST_SetSRID(ST_Envelope(ST_Collect(geom)), 4326)) AS srid
FROM   ais_positions
WHERE  mmsi = 244660000;   -- expect 4326

A healthy run shows a single aggregate node in the plan and an SRID of 4326 on the returned envelope. If the compiled SQL still contains a bare geom column with no ST_ wrapper, the @expression companion is missing and the aggregation silently reverted to the Python getter.

Gotchas checklist

  • A hybrid_property without @expression runs in Python. SQLAlchemy quietly uses the instance getter, so the aggregation happens after the rows are fetched — exactly what you were trying to avoid. Always pair the getter with an @expression.
  • ST_Extent has no SRID. It returns a box2d; decoding it as a geometry raises. Re-wrap with ST_Envelope/ST_SetSRID and re-attach SRID 4326.
  • Confusing ST_Collect with ST_Union. ST_Collect bundles (cheap, keeps overlaps); ST_Union dissolves (expensive, merges boundaries). Using ST_Union to gather AIS points wastes CPU for no benefit.
  • Parameterised aggregations need hybrid_method, not hybrid_property. A property cannot take a target SRID or tolerance argument; forcing one leads to closures over module-level state that break under concurrency.
  • Forgetting GROUP BY mixes aggregates and columns. Selecting mmsi alongside an aggregate without grouping raises a grouping error in PostgreSQL — every non-aggregated column in the SELECT must appear in GROUP BY.