PostGIS Integration with Python: Server-Side Spatial SQL from GeoPandas

When a spatial dataset outgrows what fits comfortably in memory, the analysis moves into the database. PostGIS — the spatial extension to PostgreSQL — runs predicates, joins, and aggregations server-side over indexed geometry, returning only the rows you need to Python. This stage of Spatial Analysis & Advanced Query Techniques covers the round trip between a GeoPandas GeoDataFrame and PostGIS, and sits beside the two other out-of-core engines in this section: DuckDB Spatial Analytics for embedded, file-based analysis, and Scaling with Dask-GeoPandas for distributed, partitioned processing. Reach for PostGIS when the data must be shared, transactional, and indexed for many concurrent clients rather than analysed once and discarded.

GeoPandas to PostGIS round trip GeoPandas writes a GeoDataFrame to a PostGIS table through GeoAlchemy2; PostGIS applies a GiST index and runs server-side spatial SQL; only the result rows return to GeoPandas. GeoPandas GeoDataFrame in-memory PostGIS GiST index ST_Intersects ST_DWithin server-side SQL Result set filtered rows back to Python to_postgis read_postgis
PostGIS pushes the spatial work to an indexed database; Python sends geometry in and pulls only matching rows back.

Architecture & Data Structures

PostGIS stores geometry in one of two column types, and the choice governs how every predicate behaves. A geometry column holds planar coordinates and is tagged with an SRID — the database's name for an EPSG code — so all measurement happens in that projection's units. A geography column holds longitude/latitude on the WGS84 spheroid, and its functions return true metres and compute great-circle distances directly. Use geometry in a projected SRID when most work is metric and confined to a region; use geography when data spans continents or you want spheroidal distance without picking a projection. Every geometry also carries a type modifier (point, linestring, polygon, and their multi- variants) and an optional Z or M dimension.

The Python bridge is GeoPandas' to_postgis/read_postgis, which lean on SQLAlchemy for the connection and GeoAlchemy2 for geometry type handling. SQLAlchemy owns the driver, the connection pool, and transaction boundaries; GeoAlchemy2 registers the Geometry type so pandas' writer knows to emit WKB and PostGIS knows to decode it. You rarely call GeoAlchemy2 directly, but it must be importable in the environment or to_postgis cannot serialise the geometry column.

Spatial predicates (ST_Intersects, ST_DWithin, ST_Contains) run in SQL and are accelerated by a GiST index over the geometry column. GiST is a generalised search tree; for geometry it stores each row's bounding box and lets the planner discard, in one indexed descent, every row whose box cannot satisfy the predicate. That bounding-box test is the cheap filter stage; PostGIS then runs the exact GEOS predicate only on the surviving candidates — the refine stage. Without the index, the filter stage degrades to a sequential scan over every row.

import geopandas as gpd
from sqlalchemy import create_engine

engine = create_engine("postgresql+psycopg://gis:gis@localhost:5432/gisdb")

# Read with a spatial predicate executed server-side
sql = """
SELECT parcel_id, land_use, geom
FROM parcels
WHERE ST_DWithin(geom::geography, ST_SetSRID(ST_MakePoint(%(lon)s, %(lat)s), 4326)::geography, 500)
"""
nearby = gpd.read_postgis(sql, engine, geom_col="geom", params={"lon": 7.69, "lat": 45.07})

The win is that parcels may hold tens of millions of rows; only the few hundred within 500 m cross the wire. read_postgis reconstructs the CRS from the geometry column's SRID, so nearby returns as a fully-typed GeoDataFrame ready for the operations in Mastering Core Geospatial Python Libraries.

Underneath the type name, a PostGIS geometry is a gserialized blob: a short header carrying the SRID, the type modifier, flags for Z/M dimensions and a cached bounding box, followed by the raw coordinate array. Two facts about that layout drive most performance surprises. First, the cached box is what && and the GiST index read, which is why a bounding-box test costs four float comparisons regardless of how many vertices the geometry has. Second, PostgreSQL TOASTs any row wider than roughly two kilobytes — it compresses the value and moves it to a side table — so a coastline polygon with 80 000 vertices is decompressed and reassembled every time a query reads that column. A query that only needs parcel_id and land_use but writes SELECT * pays that detoasting cost on every row, which is frequently the difference between a two-second and a two-minute response on a table of complex polygons.

The database also keeps a small catalogue you will end up querying during debugging. spatial_ref_sys is the table of known SRIDs with their PROJ and WKT definitions; geometry_columns is a view that reports, per table, the geometry column name, its SRID, its declared type and its dimensionality, derived from the column's type modifier. When to_postgis creates a table it writes a typmod-constrained column (geometry(MultiPolygon, 4326)), which is what makes geometry_columns accurate and what rejects a later insert of a mismatched type or SRID at write time rather than at query time. A column declared as bare geometry accepts anything, reports SRID 0 in the view, and defers every mistake to the first predicate that fails.

-- What the server thinks it is storing, straight from the catalogue
SELECT f_table_name, f_geometry_column, coord_dimension, srid, type
FROM geometry_columns
WHERE f_table_schema = 'public';
-- census_blocks | geometry | 2 | 4326 | MULTIPOLYGON

Environment Configuration & Dependency Resolution

conda install -c conda-forge "geopandas>=1.0" "sqlalchemy>=2.0" "geoalchemy2>=0.15" "psycopg>=3.1"
# PostGIS server (Docker is the simplest dev setup):
#   docker run -e POSTGRES_PASSWORD=gis -p 5432:5432 postgis/postgis:16-3.4

Use psycopg (version 3) with the postgresql+psycopg:// URL prefix; the older psycopg2 still works via postgresql+psycopg2:// but is on a separate maintenance track, and mixing the two URL schemes against one codebase is a common source of confusion. Keep the client's PostGIS-facing libraries aligned with the server: a client built against a much older GEOS can serialise WKB the server reads fine, but predicate results are defined by the server's GEOS, so treat the server version as the source of truth for query semantics.

Enable the extension once per database before any spatial DDL runs. This is a privileged operation, so do it as a superuser or an owner with the right grant:

from sqlalchemy import create_engine, text

engine = create_engine("postgresql+psycopg://gis:gis@localhost:5432/gisdb")
with engine.begin() as conn:
    conn.execute(text("CREATE EXTENSION IF NOT EXISTS postgis"))
    version = conn.execute(text("SELECT postgis_full_version()")).scalar()
    print(version)  # confirms GEOS, PROJ, and GDAL versions the server links against

Version differences that change the answer

Four version boundaries decide whether a query you copied from a colleague works at all:

Keep the client and the server on the same major PostGIS line where you can. A newer client library against an older server is usually fine because the wire format is EWKB and stable, but the semantics of a predicate — how a marginally invalid polygon is treated, whether an overlay throws — come from the server's GEOS. Reproduce production bugs against the production server version, not against a locally installed Docker tag that happens to be two minor versions ahead.

Vectorized Operations & Core Workflow

The standard loop: write a GeoDataFrame to a table, index it, then query with spatial SQL. Reproject to a single, deliberate storage CRS before writing so the SRID stamped on the column is correct from the start. Connection, chunking, and write-mode details are covered in depth in Connecting GeoPandas to PostGIS with SQLAlchemy.

import geopandas as gpd
from sqlalchemy import create_engine, text

engine = create_engine("postgresql+psycopg://gis:gis@localhost:5432/gisdb")

census_blocks = gpd.read_file("census_blocks.gpkg").to_crs(epsg=4326)

# chunksize streams the write so a large frame never buffers as one giant INSERT
census_blocks.to_postgis(
    "census_blocks", engine, if_exists="replace", index=False, chunksize=10_000
)

with engine.begin() as conn:
    conn.execute(text("CREATE INDEX ON census_blocks USING GIST (geometry)"))
    conn.execute(text("ANALYZE census_blocks"))  # refresh planner statistics

Two details make or break this write. if_exists="replace" drops and recreates the table, which also drops any index and constraints — recreate the GiST index after every bulk replace, not once at setup. And ANALYZE is not optional: the planner chooses an index scan over a sequential scan based on row-count and selectivity statistics, and immediately after a load those statistics are stale, so the first query can run the slow plan until statistics catch up.

When to_postgis stops scaling

to_postgis builds parameterised multi-row INSERT statements. That is convenient and safe, but each row carries per-value protocol overhead, and past roughly a million geometries the write becomes the slowest step in the pipeline — commonly ten to twenty times slower than the same data going in through COPY. PostgreSQL's COPY path skips statement parsing and per-row parameter binding entirely, and psycopg 3 exposes it directly. Feed it hex-encoded EWKB and let the server cast on ingest:

import geopandas as gpd
import shapely
from sqlalchemy import create_engine, text

engine = create_engine("postgresql+psycopg://gis:gis@localhost:5432/gisdb")
buildings = gpd.read_parquet("buildings.parquet").to_crs(epsg=4326)

# to_wkb(hex=True, include_srid=True) produces EWKB the server casts for free
wkb_hex = shapely.to_wkb(buildings.geometry.values, hex=True, include_srid=True)

with engine.begin() as conn:
    conn.execute(text("""
        CREATE TABLE IF NOT EXISTS buildings_stage (
            building_id bigint,
            height_m    double precision,
            geometry    geometry(Polygon, 4326)
        )
    """))
    raw = conn.connection.driver_connection          # the underlying psycopg 3 connection
    with raw.cursor().copy(
        "COPY buildings_stage (building_id, height_m, geometry) FROM STDIN"
    ) as copy:
        for bid, height, geom in zip(buildings.building_id, buildings.height_m, wkb_hex):
            copy.write_row((int(bid), float(height), geom))

Note int(bid) and float(height): pandas hands you numpy scalars, and psycopg 3 refuses to adapt them. This is the single most common failure when moving a working psycopg2 loader to psycopg 3.

The staging table in that snippet is doing a second job. Loading into *_stage and then swapping it into place is how you refresh a table that live clients are querying without a window where the data is half-written or the index is missing:

with engine.begin() as conn:
    conn.execute(text("CREATE INDEX ON buildings_stage USING GIST (geometry)"))
    conn.execute(text("ANALYZE buildings_stage"))
    # The swap itself is a metadata-only operation inside one transaction
    conn.execute(text("ALTER TABLE IF EXISTS buildings RENAME TO buildings_old"))
    conn.execute(text("ALTER TABLE buildings_stage RENAME TO buildings"))
    conn.execute(text("DROP TABLE IF EXISTS buildings_old"))

Because the index and ANALYZE happen on the staging table before the rename, readers never see the table without its index, and the first query after the swap already has fresh statistics. Compare that with if_exists="replace", where the table exists but unindexed for the whole duration of the load — on a large layer that is minutes of every query falling back to a sequential scan.

Geometry / Data Processing Details

Server-side joins are where PostGIS earns its place. A spatial join that would materialise a huge intermediate in memory runs as an indexed nested loop in the database: for each row on the left, GiST returns the handful of right-hand candidates whose bounding boxes overlap, and the exact predicate runs only on those.

Indexed spatial join: filter, then refine Each left-table row probes the GiST bounding-box index on the right table, which cheaply returns a small candidate set — the filter stage. The exact GEOS predicate ST_Contains then runs only on those candidates to keep the true matches — the refine stage. Without the index, the planner falls back to a sequential scan that runs the exact predicate on every left-by-right pair. Indexed spatial join: filter, then refine probe filter refine Left table census_blocks each row → GiST index R-tree of bounding boxes Candidates boxes overlap a handful of rows Matches predicate holds true hits only Filter stage cheap bounding-box test, index-only — discards non-overlapping rows Refine stage exact ST_Contains Without a GiST index the planner falls back to a sequential scan — the exact predicate runs on every left × right pair (N × M tests), with no candidate set to shrink the work
The GiST index does the cheap work first: one indexed descent per left row yields a small candidate set, and the exact GEOS predicate runs only on those — instead of every row pair.
import geopandas as gpd
from sqlalchemy import create_engine

engine = create_engine("postgresql+psycopg://gis:gis@localhost:5432/gisdb")

# Count sensors per census block entirely in SQL
join_sql = """
SELECT b.block_id, b.geometry, COUNT(s.sensor_id) AS sensor_count
FROM census_blocks b
LEFT JOIN sensors s ON ST_Contains(b.geometry, s.geom)
GROUP BY b.block_id, b.geometry
"""
result = gpd.read_postgis(join_sql, engine, geom_col="geometry")

The equivalent operation on in-memory data is covered in Spatial Joins & Merging; PostGIS is the answer when "in memory" stops being an option.

Two patterns repay the effort of learning them. For nearest-neighbour queries, the index-backed <-> distance operator combined with ORDER BY ... LIMIT k performs a true KNN search that walks the GiST tree instead of scanning — the server-side analogue of the approaches in Nearest Neighbor & KD-Tree Search:

-- Five nearest hydrants to a given fire, index-assisted via <->
SELECT h.hydrant_id, h.geom <-> f.geom AS dist
FROM hydrants h, fires f
WHERE f.fire_id = 42
ORDER BY h.geom <-> f.geom
LIMIT 5;

And for joins against very large or complex polygons — coastlines, admin boundaries with millions of vertices — pre-splitting the polygons with ST_Subdivide shrinks each row's bounding box so the GiST filter stage rejects far more candidates, often turning a slow join fast without changing the result.

The mechanism is worth understanding because it explains a class of "the index is there but the query is still slow" reports. Consider a single watershed polygon with 120 000 vertices that snakes across a whole province. Its bounding box covers the entire province, so every sensor in the province survives the filter stage and goes to the exact predicate — the index prunes nothing, and GEOS then walks 120 000 vertices per candidate point. Subdividing that one row into a few hundred pieces of at most 512 vertices each gives every piece a small, tight box; now the filter stage discards almost everything, and the exact test runs against a fragment instead of the whole monster:

-- Split any polygon over 512 vertices into tiles; keep the parent id
CREATE TABLE watersheds_split AS
SELECT watershed_id, name, ST_Subdivide(geom, 512) AS geom
FROM watersheds;

CREATE INDEX ON watersheds_split USING GIST (geom);
ANALYZE watersheds_split;

-- The join now touches fragments, so DISTINCT restores one row per parent
SELECT DISTINCT w.watershed_id, s.sensor_id
FROM watersheds_split w
JOIN sensors s ON ST_Intersects(w.geom, s.geom);

Two caveats. ST_Subdivide multiplies row count — a table of 4 000 watersheds can become 300 000 fragments — so aggregate back to the parent id with DISTINCT or a GROUP BY, and keep the split table as a derived artefact rather than the system of record. And the vertex ceiling is a tuning knob, not a constant: values between 100 and 1 000 are the usual range, with smaller values buying tighter boxes at the cost of more rows.

The complementary trick is prepared geometry. When one geometry is tested against many others in the same query, PostGIS builds an internal edge index for the repeated side and caches it for the duration of the statement. That happens automatically inside a join, which is why a single JOIN ... ON ST_Contains(...) beats issuing one query per polygon from a Python loop by an order of magnitude even when both use the same index — the loop throws away the prepared geometry on every round trip, and pays network latency on top.

Finally, some analyses that people export to Python belong in SQL. ST_ClusterDBSCAN is a window function that assigns cluster ids in a single pass, so density clustering over a few million points can run beside the data instead of streaming it out:

-- Cluster crash points: 50 m radius, at least 5 points to seed a cluster
SELECT crash_id,
       ST_ClusterDBSCAN(geom, eps := 50, minpoints := 5) OVER () AS cluster_id
FROM crashes
WHERE crash_date >= DATE '2025-01-01';

The distance is in the column's units, so this only means 50 metres on a projected column. The trade-off against the scikit-learn path — parameter sweeps, HDBSCAN, silhouette scoring — is laid out in Spatial Clustering Algorithms; the rule of thumb is that SQL wins when the answer is a label you want to store back, and Python wins when you are still choosing the parameters.

CRS Alignment & Projection Pipeline

Every PostGIS geometry carries an SRID, and predicates between mismatched SRIDs raise an error rather than silently misalign — a feature, not a bug. Decide the storage CRS deliberately: EPSG:4326 for general storage and geography-based distance, or a projected SRID (the local UTM zone, never Web Mercator EPSG:3857) when most work is metric. Reproject in Python via Coordinate Systems with PyProj before writing, or with ST_Transform in SQL when the source rows already live in the database.

from sqlalchemy import create_engine, text

engine = create_engine("postgresql+psycopg://gis:gis@localhost:5432/gisdb")

# Distance in metres: cast 4326 geometry to geography, which measures on the spheroid
metric_sql = """
SELECT a.hydrant_id
FROM hydrants a, fires f
WHERE ST_DWithin(a.geom::geography, f.geom::geography, 300)
"""
# Avoid ST_Distance / ST_DWithin on raw 4326 geometry — it returns degrees, not metres.

Two axis-order traps recur. PostGIS constructors like ST_MakePoint(x, y) and WKB always take longitude first, so data ingested latitude-first lands transposed with no error raised — normalise axis order at ingest, the same always_xy discipline you apply in PyProj. And ST_Transform re-projects but does not reproject the concept of distance: after transforming to a projected SRID, measurements are back in that projection's units, so choose the target zone from the data's centroid to keep distortion small across the extent.

A subtler mismatch bites teams that reproject in both places. ST_Transform runs against the server's PROJ installation and its spatial_ref_sys table, while to_crs runs against the client's PROJ. If a datum shift needs a transformation grid — NADCON/NTv2 files for North American or Australian datums, for example — and the grid is installed in the Python environment but not in the container running PostgreSQL, the two paths silently disagree by up to a couple of metres. That is invisible in a web map and fatal in a cadastral join. Decide which side owns reprojection, and if the answer is the server, install the grids there; the client-side mechanics are covered in Datum Shifts & Grid Files in PyProj.

You can compare the two directly, which is the fastest way to prove a suspicion:

import geopandas as gpd
from sqlalchemy import create_engine, text

engine = create_engine("postgresql+psycopg://gis:gis@localhost:5432/gisdb")

with engine.connect() as conn:
    server_xy = conn.execute(text(
        "SELECT ST_X(p), ST_Y(p) FROM ("
        "  SELECT ST_Transform(ST_SetSRID(ST_MakePoint(-93.10, 44.95), 4326), 26915) AS p"
        ") t"
    )).one()

client = gpd.GeoSeries.from_xy([-93.10], [44.95], crs="EPSG:4326").to_crs(26915)
print(server_xy)                             # (491234.87, 4977512.44)
print(client.x.iloc[0], client.y.iloc[0])    # 491234.87 4977512.44
# A disagreement beyond a few centimetres means the two PROJ installs
# picked different transformation pipelines — usually a missing grid file.

If your work uses a local grid that is not in the EPSG registry — a municipal system, a legacy mine grid — insert it into spatial_ref_sys with an SRID in the user range (900000 and above by convention) and its full PROJ string, then reference it like any other SRID. Keep that insert in the same migration that creates the tables using it, or a restore into a fresh database will fail on every ST_Transform call.

Latitude-first ingest versus longitude-first ingest Two identical world graticules holding the same hydrant record. On the left the coordinates were passed latitude first, so the point plots near the equator off East Africa instead of at its dashed target on the 45 degree north line, and a dashed connector shows how far it slid. On the right the same numbers passed longitude first land on the target. No error is raised in either case, which is why the transposition has to be caught at ingest. PostGIS constructors are longitude-first the same hydrant at 7.69 E, 45.07 N, ingested two ways — neither raises an error ST_MakePoint(lat, lon) 45° N 0° intended lands here ST_MakePoint(lon, lat) 45° N 0° on target row stored silently ~4000 km off, WKB and all normalise axis order once, at ingest
A transposed pair is not a syntax error — it is a valid geometry in the wrong hemisphere, so the only place to catch it is a bounds assertion before the write.

Production Export & Integration

Two of those bullets deserve code. Streaming a large read is the one people discover too late, usually after a worker is OOM-killed halfway through a nightly export. By default the driver buffers the entire result set in client memory before read_postgis sees a single row; a server-side cursor changes that to batches:

import geopandas as gpd
import pandas as pd
from sqlalchemy import create_engine, text

engine = create_engine("postgresql+psycopg://gis:gis@localhost:5432/gisdb")

sql = text("SELECT parcel_id, land_use, geom FROM parcels WHERE land_use = 'residential'")

frames = []
with engine.connect().execution_options(stream_results=True, max_row_buffer=50_000) as conn:
    for chunk in gpd.read_postgis(sql, conn, geom_col="geom", chunksize=50_000):
        # Reduce each chunk immediately — never accumulate the whole table
        frames.append(chunk.dissolve(by="land_use", aggfunc="sum"))

summary = gpd.GeoDataFrame(pd.concat(frames), crs=frames[0].crs)

And the tile path is worth adopting whole rather than reinventing: ST_AsMVT with ST_AsMVTGeom produces a Mapbox Vector Tile binary directly from a query, so a FastAPI or Flask route can serve tiles from live data without a build step.

-- One tile of parcels at z/x/y, clipped and simplified for the tile grid
WITH bounds AS (SELECT ST_TileEnvelope(:z, :x, :y) AS geom),
     src AS (
       SELECT p.parcel_id, p.land_use,
              ST_AsMVTGeom(ST_Transform(p.geom, 3857), bounds.geom, 4096, 64, true) AS geom
       FROM parcels p, bounds
       WHERE p.geom && ST_Transform(bounds.geom, 4326)
     )
SELECT ST_AsMVT(src, 'parcels', 4096, 'geom') FROM src;

EPSG:3857 appears here for a legitimate reason — it is the tile grid's own coordinate system, and the geometry is being drawn, not measured. That is the one context in which Web Mercator is the right answer; the moment a distance or an area is involved, go back to a projected CRS chosen for the data's extent. When tiles are static rather than live, the pre-baked path through Vector Tile Pipelines with PMTiles costs less to serve.

Where PostGIS stops being the right tool

PostGIS is a transactional row store, and that shapes its ceiling. Point queries, radius searches, and joins that touch a small slice of a big table stay fast essentially forever, because the index bounds the work. Analytical scans that touch every row — "sum this attribute over 200 million rows grouped by region" — are where a columnar engine pulls ahead by an order of magnitude, because it reads only the columns involved and stores them compressed. That is the argument for keeping DuckDB Spatial Analytics in the toolkit alongside PostGIS rather than instead of it, and the detailed comparison lives in DuckDB Spatial vs PostGIS for Analytics. The other ceiling is write throughput on a single primary: heavy concurrent ingest plus heavy analytical reads on one instance is the classic reason a nightly job and a live map start timing each other out. Move the reads to a replica before you start tuning the queries.

Windows / Platform Edge Cases & Debugging

Frequently Asked Questions

Should I store geometry in EPSG:4326 or in a projected SRID? Store in 4326 when the table is a shared source of truth read by many clients — web maps, tile servers and exports all expect it, and ::geography gives you true metres on demand for distance and area. Store in a projected SRID when the table exists to serve one region's metric analysis and those casts would run on every query, because casting to geography per row costs real time at scale. What you should not do is pick Web Mercator for storage on the theory that it is "already metric": its scale factor grows with latitude, so measurements are wrong by roughly 40% at 45° and by a factor of two at 60°.

Is it faster to do the spatial join in PostGIS or in GeoPandas? If the data is already in the database and the result is a small fraction of it, PostGIS wins outright, because the index bounds the work and only matching rows cross the network. If the data is already in memory and the join is the whole job, GeoPandas is competitive and often faster, since it skips serialisation entirely. The expensive answer is the hybrid nobody intended: pulling two full tables into Python to join them, which pays the network cost of a database and the memory cost of an in-memory engine at the same time. The in-memory mechanics are covered in Spatial Joins & Merging.

Why does my query still do a sequential scan when a GiST index exists? Three causes, in descending order of frequency: ANALYZE has not run since the load, so the planner has no selectivity estimate; the predicate wraps the indexed column in a function, which makes the stored envelopes unusable; or the query genuinely matches most of the table, in which case a scan is the correct plan and the planner is right. Diagnose with EXPLAIN (ANALYZE, BUFFERS) and work through Spatial Indexing in PostGIS with GiST.

Can I skip SQLAlchemy and use psycopg directly? For reads, yes — gpd.read_postgis accepts anything with a DBAPI cursor, and a raw psycopg connection works fine. For writes it is not optional: to_postgis needs a SQLAlchemy connectable so GeoAlchemy2 can register the geometry type, and passing a bare connection either raises inside the ORM layer or writes geometry as text. Details are in Connecting GeoPandas to PostGIS with SQLAlchemy.

How big can a PostGIS table get before I need to partition it? Row count matters less than working-set size. A hundred million points with a GiST index that fits in RAM answers window queries in milliseconds; ten million complex polygons whose index does not fit will thrash. Partitioning by region or by time helps when queries always carry a partition key, because the planner can prune whole partitions before touching an index — but it adds real operational complexity, so reach for ST_Subdivide, a covering index, and more RAM first.

Does PostGIS handle raster data well enough to skip Python? It can store and tile rasters through postgis_raster, and ST_Clip plus ST_SummaryStats will produce zonal statistics server-side. In practice the Python path is faster to develop, faster to run on anything cloud-hosted, and easier to parallelise, because a Cloud-Optimized GeoTIFF supports windowed reads without a database at all. Keep vectors in PostGIS, keep rasters in object storage, and join them with Zonal Statistics & Raster Sampling.

Related