Tracing PostGIS Query Spans in Spatial ETL

Wrap every PostGIS query span — ST_Intersects joins, ST_MakeValid repairs, and spatial index scans — in a child OpenTelemetry span that records db.statement, the geometry count, and the returned row count, so a slow spatial join stops being a mystery in the load stage and becomes a labelled bar in the trace waterfall. Once each query carries its timing and cardinality, correlating a 4-second ST_Intersects span with a sequential scan over a table missing its GiST index takes seconds, because the span tells you both the SQL text and the geometry volume it processed.

When to Use This Pattern#

Add explicit query spans when the database is a first-class actor in the pipeline’s latency and correctness, which for spatial ETL is nearly always. Typical triggers:

• A load or enrichment task runs spatial predicates (ST_Intersects, ST_DWithin, ST_Contains) and you cannot tell which query dominates task latency.
• You suspect a spatial join is doing a sequential scan because a GiST index is missing or not being used.
• You want database timing nested under the flow trace produced by OpenTelemetry tracing for spatial tasks rather than living in disconnected database logs.
• Geometry repair with ST_MakeValid runs inline and you need to attribute its cost separately from the join itself.

Why Spatial Queries Deserve Their Own Spans#

Generic database auto-instrumentation will happily wrap cur.execute() and record a span with a duration, but it treats a spatial join the same as a primary-key lookup. That is inadequate for spatial ETL, because the cost of a PostGIS predicate is dominated by geometry complexity and index behaviour, neither of which a plain duration reveals. A one-second ST_Intersects might be entirely reasonable over a million-vertex coastline, or a scandal over ten simple points against an unindexed table — the number that distinguishes them is the geometry count set against the row count, not the wall-clock time alone. By hand-authoring the span you get to attach geo.input_geom_count, geo.output_row_count, and the operand SRID, turning a bare timing into a diagnosable record. When you later sort the trace backend’s spans by duration, the spatial spans carry enough context to triage without opening psql, and the ones that need a plan inspection are already labelled with the SQL that produced them.

Complete Working Example#

The helper below opens a child span per query, records a truncated db.statement, the input geometry count, and the output row count, and marks the span as errored on failure. It wraps both a raw psycopg cursor and a SQLAlchemy execution so either access layer produces the same span shape.

from contextlib import contextmanager
from typing import Iterator, Any

import psycopg
from sqlalchemy import create_engine, text
from sqlalchemy.engine import Engine

from opentelemetry import trace
from opentelemetry.trace import Status, StatusCode

tracer = trace.get_tracer("geo.postgis")


@contextmanager
def query_span(name: str, statement: str, geom_count: int | None = None) -> Iterator[Any]:
    """Open a child span describing one PostGIS query."""
    with tracer.start_as_current_span(name) as span:
        span.set_attribute("db.system", "postgresql")
        span.set_attribute("db.statement", statement.strip()[:1024])  # bound cardinality
        if geom_count is not None:
            span.set_attribute("geo.input_geom_count", geom_count)
        try:
            yield span
        except Exception as exc:  # noqa: BLE001 - re-raised after recording
            span.set_status(Status(StatusCode.ERROR, str(exc)[:200]))
            span.record_exception(exc)
            raise


def spatial_join_psycopg(dsn: str, input_geoms: list[str]) -> int:
    """Find admin regions each input point intersects, traced as one span."""
    sql = """
        SELECT p.geom_wkt, a.region_id
        FROM   unnest(%(wkts)s::text[]) AS p(geom_wkt)
        JOIN   admin_boundaries a
          ON   ST_Intersects(a.geom, ST_GeomFromText(p.geom_wkt, 4326))
    """
    with psycopg.connect(dsn) as conn, conn.cursor() as cur:
        with query_span("postgis.st_intersects_join", sql, geom_count=len(input_geoms)) as span:
            cur.execute(sql, {"wkts": input_geoms})
            rows = cur.fetchall()
            span.set_attribute("geo.output_row_count", len(rows))
            span.set_attribute("geo.epsg", 4326)
            return len(rows)


def make_valid_sqlalchemy(engine: Engine, table: str) -> int:
    """Repair invalid geometries in place, traced with before/after counts."""
    count_sql = f"SELECT count(*) FROM {table} WHERE NOT ST_IsValid(geom)"
    repair_sql = f"UPDATE {table} SET geom = ST_MakeValid(geom) WHERE NOT ST_IsValid(geom)"
    with engine.begin() as conn:
        invalid = conn.execute(text(count_sql)).scalar_one()
        with query_span("postgis.st_makevalid", repair_sql, geom_count=invalid) as span:
            result = conn.execute(text(repair_sql))
            span.set_attribute("geo.repaired_count", result.rowcount)
            return result.rowcount

Each query becomes a self-describing span: postgis.st_intersects_join carries the input geometry count and output row count, while postgis.st_makevalid records how many geometries it repaired. When the join span is slow, the ratio of geo.output_row_count to geo.input_geom_count next to a long duration is the first signal of a scan problem.

The query_span helper is deliberately access-layer agnostic. Whether the query runs through a raw psycopg cursor for a bulk COPY-adjacent workload or through a SQLAlchemy Connection in a service that already owns an engine, the span shape is identical: same db.system, same truncated db.statement, same geometry-count attributes. That uniformity matters when you query the trace backend, because a single filter such as db.statement =~ "ST_Intersects" returns every spatial join across the codebase regardless of how it was executed. The ST_MakeValid example also shows a useful pattern for repair operations: capture the invalid count before the update in the same transaction, so the span records how much work the repair actually did rather than how many rows the table holds. Wrapping the count and the update in one engine.begin() block keeps that measurement consistent even under concurrent writers.

Parameter & Option Reference#

Verification & Testing#

First confirm the span attributes are attached, then verify the underlying query actually uses the GiST index rather than scanning. An in-memory exporter checks the span; EXPLAIN ANALYZE checks the plan.

def test_join_span_records_counts(exporter, dsn: str) -> None:
    spatial_join_psycopg(dsn, ["POINT(12.5 41.9)", "POINT(2.3 48.8)"])
    span = next(s for s in exporter.get_finished_spans()
                if s.name == "postgis.st_intersects_join")
    assert span.attributes["geo.input_geom_count"] == 2
    assert "ST_Intersects" in span.attributes["db.statement"]
    assert span.attributes["geo.epsg"] == 4326

To prove a slow join is an indexing problem, run the query plan in psql for the exact statement the span recorded — because db.statement holds the parameterised SQL, you can copy it straight from the trace into an EXPLAIN. A Seq Scan on admin_boundaries under the intersect is the smoking gun; an Index Scan using ... gist is what you want.

EXPLAIN (ANALYZE, BUFFERS)
SELECT a.region_id
FROM   admin_boundaries a
JOIN   input_points p
  ON   ST_Intersects(a.geom, p.geom);

-- Fix: create the GiST index the planner needs, then re-EXPLAIN.
CREATE INDEX IF NOT EXISTS admin_boundaries_geom_gist
  ON admin_boundaries USING GIST (geom);
ANALYZE admin_boundaries;

The BUFFERS option is worth including because it exposes how many pages the query touched. A join that reads far more buffers than the result size would suggest is scanning data it should have skipped — a second, independent signal of a missing or unused index that complements the raw timing on the span. Confirm from the shell that the index now exists and is the right type:

psql "$DSN" -c "\di+ admin_boundaries_geom_gist"

After creating the index, the postgis.st_intersects_join span duration should drop sharply on the next run, and the trace waterfall in your backend will show the load stage shrinking accordingly. Capture a before-and-after pair of traces for the same input so the improvement is concrete: the span attributes hold constant (same geo.input_geom_count, same geo.output_row_count) while only the duration changes, which is exactly the evidence you want when justifying the index in a review. Keeping both traces also protects against a regression where a later schema migration drops or invalidates the index and the sequential scan quietly returns.

It is worth watching the plan node type as well as the raw timing. A Bitmap Heap Scan fed by a Bitmap Index Scan on the GiST index is the healthy shape for a selective spatial join returning many rows, whereas a Nested Loop with an inner Index Scan suits a small, highly selective probe. If EXPLAIN ANALYZE shows the planner estimating one row but the join actually returns hundreds of thousands, stale statistics are misleading it — another reason to run ANALYZE after any bulk load. Recording geo.output_row_count on the span gives you the real cardinality to compare against the planner’s estimate without re-running the query.

Common Pitfalls#

• SRID mismatch disables the index. If admin_boundaries.geom is SRID 3857 and the input is 4326, ST_Intersects either errors or forces an implicit transform that bypasses the GiST index. Record geo.epsg on the span so a mismatch is visible; align SRIDs before the join, echoing the CRS discipline in state management in geospatial flows.
• Putting full result geometries in db.statement. Interpolating returned WKT into the attribute explodes storage and can leak large payloads. Keep db.statement to the parameterised SQL text and store counts, not geometries.
• Tracing the connect, not the query. Wrapping psycopg.connect() in the span instead of cur.execute() measures pool checkout, not spatial work. Open the span tightly around the execute-and-fetch so the duration reflects the predicate.
• Missing ANALYZE after index creation. A fresh GiST index the planner has no statistics for may still be ignored; the span stays slow and the pitfall looks like the index “not working.” Run ANALYZE, then re-check the plan. This connects the load stage back to the flow trace assembled via propagating trace context across Prefect tasks.
• Wrapping a lazy result set. With psycopg server-side cursors or a lazily evaluated SQLAlchemy result, the query may not actually execute until rows are consumed outside the span’s with block, so the span records near-zero duration while the real work happens later. Force materialisation — fetchall(), or iterate the result — inside the span so the timing reflects execution, not just statement preparation.
• Untraced implicit transforms. A predicate like ST_Intersects(a.geom, ST_Transform(p.geom, 4326)) may reproject millions of geometries on the fly, and the cost hides inside a single join span. Record the operand SRIDs as attributes and prefer aligning SRIDs at load time so the span reflects a clean index-backed join rather than a hidden per-row transform.
• One giant span for a multi-statement transaction. Bundling a whole load transaction — insert, index refresh, and analyze — into a single span hides which statement is slow. Open a child span per meaningful statement so the waterfall attributes cost precisely; the small extra span count is well worth the diagnostic clarity when a nightly load regresses.

Related: This page complements OpenTelemetry tracing for spatial tasks, shares its flow trace with propagating trace context across Prefect tasks, feeds slow-query signals to Grafana dashboards for GIS workflows, and pairs with structured logging for geospatial flows. The PostGIS spatial indexing documentation explains GiST index behaviour in depth.

← Back to OpenTelemetry Tracing for Spatial Tasks