Connecting GeoPandas to PostGIS with SQLAlchemy

Moving a layer from a file into PostGIS should be a single call, yet to_postgis fails or silently corrupts geometry more often than it works — usually because the connection object, the geometry type registration, or the CRS is wrong. This guide wires GeoPandas to PostGIS the modern way — a SQLAlchemy 2.0 Engine plus GeoAlchemy2 — and writes, indexes, and reads a GeoDataFrame without losing coordinate reference system metadata. It is for anyone moving from file-based workflows to a spatial database. It sits under PostGIS Integration with Python in Spatial Analysis & Advanced Query Techniques.

GeoDataFrame–PostGIS round trip through SQLAlchemy and GeoAlchemy2 The write path sends a GeoDataFrame through a SQLAlchemy Engine to GeoAlchemy2, which serialises each geometry to Extended Well-Known Binary with the SRID taken from the frame's CRS, producing a typed PostGIS geometry column with a GiST index. The read path reverses it: read_postgis pulls rows back through GeoAlchemy2, which decodes the EWKB into shapely geometries and restores the CRS. GeoDataFrame shapely + .crs EPSG:4326 PostGIS table geometry, SRID 4326 + GiST index Write · to_postgis() SQLAlchemy Engine GeoAlchemy2 encode → EWKB SRID column stores the EPSG code read from .crs read_postgis · decode EWKB shapely geometries, CRS restored Read · read_postgis()
The write path hands GeoPandas a SQLAlchemy Engine; GeoAlchemy2 serialises each geometry to EWKB with the SRID taken from .crs. read_postgis reverses the trip — decoding EWKB back to shapely geometries and restoring the CRS.

Why This Approach / What Goes Wrong

GeoDataFrame.to_postgis and geopandas.read_postgis are thin wrappers over SQLAlchemy and GeoAlchemy2. Three assumptions have to hold, and each breaks in a distinct way when it doesn't.

They need a SQLAlchemy Engine, not a raw DBAPI connection. Passing a bare psycopg connection object raises an AttributeError deep inside the ORM layer, or — worse on some versions — writes the frame but skips geometry handling entirely. The Engine is what lets GeoPandas open its own transactional connection and hand GeoAlchemy2 a place to register types.

GeoAlchemy2 must be importable in the active environment. GeoPandas does not encode geometry itself; it delegates to GeoAlchemy2's Geometry type, which serialises each shapely object to Extended Well-Known Binary (EWKB) with the SRID embedded, and registers a result processor that decodes EWKB back on read. If the package is missing, the geometry column is written as opaque WKT text and comes back as plain Python strings instead of shapely geometries — no error, just a broken frame.

The CRS becomes the SRID. GeoPandas reads GeoDataFrame.crs, resolves it to an EPSG code, and stamps that as the column's SRID. If .crs is None, the table is created with SRID 0, and every later spatial predicate that compares against a real SRID fails with an Operation on mixed SRID geometries error. This is the single most common failure, and it is why a CRS guard belongs in the write path — the same discipline that matters everywhere you touch projections in Python.

There is a fourth assumption that only bites on the second write. The geometry column has a type modifier, and appends must match it. When GeoPandas creates the table it inspects the frame and stamps the column as, say, geometry(Polygon,4326). Load a later batch that contains a single MultiPolygon — a parcel that happens to be split by a road, an island in an administrative unit — and PostgreSQL rejects the whole insert with Geometry type (MultiPolygon) does not match column type (Polygon). Nothing about the first write warned you that the schema had been narrowed by whatever happened to be in the first batch. The durable fix is to normalise geometry type on the way in, so the column is created as the widest type you will ever store.

from shapely.geometry import MultiPolygon

# Promote every single-part polygon so the column is created as MULTIPOLYGON
zoning["geometry"] = [
    geom if geom.geom_type == "MultiPolygon" else MultiPolygon([geom])
    for geom in zoning.geometry
]
print(sorted(zoning.geom_type.unique()))   # ['MultiPolygon']

Mixed points-and-polygons in one frame are a different problem: PostGIS will accept them under a generic geometry column, but no index or predicate then behaves the way a reader expects. Split them into separate tables instead — the reasoning is worked through in Handling Mixed Geometry Types in a GeoDataFrame.

Get the engine, the extension, the CRS, and the geometry type right once, and the round trip is reliable and repeatable.

Three preconditions of a to_postgis write and how each one fails A vertical decision chain beginning at a to_postgis call. The first check asks whether the connection object is a SQLAlchemy Engine; answering no branches right to an AttributeError raised deep in the ORM layer, because a raw driver connection offers nowhere to register types. The second check asks whether GeoAlchemy2 is importable in the active environment; answering no branches to geometry written as opaque text that reads back as Python strings with no error raised. The third check asks whether the GeoDataFrame carries a CRS; answering no branches to a column created with SRID zero, which makes every later predicate fail with a mixed-SRID error. Only when all three answer yes does the write produce a typed geometry column whose SRID is taken from the frame's CRS. Three preconditions, three distinct failures gdf.to_postgis(name, con) con is a SQLAlchemy Engine? geoalchemy2 importable in this interpreter? gdf.crs resolves to an EPSG code? Typed geometry column, SRID from .crs round trip returns shapely geometries with the CRS intact yes yes yes no no no AttributeError inside the ORM layer a raw driver connection has nowhere to register types Geometry stored as opaque text no error raised — it reads back as Python str Column created with SRID 0 every later predicate fails on mixed SRIDs Two of the three failures are silent — only the first one raises at write time.
Each precondition fails in its own way: only the missing Engine raises immediately — a missing GeoAlchemy2 or a missing CRS corrupts the table quietly.

Prerequisites

conda install -c conda-forge "geopandas=0.14.*" "sqlalchemy=2.0.*" "geoalchemy2=0.15.*" "psycopg=3.1.*"

Keep geoalchemy2 in the same environment your script runs in — a common cause of the "geometry comes back as text" symptom is a second interpreter that has GeoPandas but not GeoAlchemy2.

Step-by-Step Implementation

1. Create the engine and enable PostGIS. The engine is created once and reused; it owns a connection pool, so you do not open a new socket per query.

import geopandas as gpd
from sqlalchemy import create_engine, text

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

with engine.begin() as conn:                     # begin() commits on exit
    conn.execute(text("CREATE EXTENSION IF NOT EXISTS postgis"))

Use engine.begin() for anything that writes — it opens a transaction and commits when the block exits cleanly (or rolls back on exception). Use engine.connect() for read-only queries.

2. Load a layer and guarantee it has a CRS. Fail loud rather than write SRID 0.

zoning = gpd.read_file("zoning.gpkg")

if zoning.crs is None:
    raise ValueError("Source has no CRS — assign one before writing to PostGIS")

zoning = zoning.to_crs(epsg=4326)                # store in one known SRID

Storing in EPSG:4326 (lon/lat) is a fine default for a shared table that many clients read. It is not the CRS to compute areas or distances in — for metric work, reproject to an appropriate projected system such as the correct UTM zone for your data before measuring, or do the measurement in SQL with ST_Area(geography(geometry)).

3. Write the frame to a table. GeoPandas creates the geometry column with the matching SRID and a typed geometry definition.

zoning.to_postgis("zoning", engine, if_exists="replace", index=False)

For large layers, stream the insert and skip the row-by-row overhead of building one giant statement:

zoning.to_postgis(
    "zoning", engine,
    if_exists="replace", index=False,
    chunksize=10_000,                            # flush every 10k rows
)

4. Add a spatial index and refresh planner statistics. Without a GiST index PostGIS scans every row; without ANALYZE the planner has no statistics and may ignore the index it has. Both matter — see Spatial Indexing in PostGIS with GiST for how to confirm the index is actually used.

with engine.begin() as conn:
    conn.execute(text(
        "CREATE INDEX IF NOT EXISTS zoning_gix ON zoning USING GIST (geometry)"
    ))
    conn.execute(text("ANALYZE zoning"))

5. Read back with a server-side spatial filter. Push the bounding box into SQL so PostGIS uses the index and returns only the rows you need, instead of pulling the whole table into Python and filtering there.

bbox_sql = """
SELECT zone_code, geometry
FROM zoning
WHERE geometry && ST_MakeEnvelope(%(xmin)s, %(ymin)s, %(xmax)s, %(ymax)s, 4326)
"""
window = gpd.read_postgis(
    bbox_sql, engine, geom_col="geometry",
    params={"xmin": 7.6, "ymin": 45.0, "xmax": 7.8, "ymax": 45.1},
)

The && operator is the bounding-box overlap test that the GiST index accelerates; pair it with ST_Intersects when you need an exact-geometry filter as well. Always pass literals through params rather than f-strings — the same parameterisation that protects against SQL injection also lets the driver adapt Python types correctly.

6. Pin the column definition explicitly when the table is long-lived. Letting to_postgis infer the type is fine for a scratch table and risky for one other services read. Pass a GeoAlchemy2 Geometry through dtype and the column is created exactly as you intend, regardless of what the first batch happened to contain.

from geoalchemy2 import Geometry

zoning.to_postgis(
    "zoning", engine,
    if_exists="replace", index=False,
    dtype={"geometry": Geometry("MULTIPOLYGON", srid=4326)},
)

The same argument accepts Geometry("GEOMETRY", srid=4326) when the column genuinely must hold several types, and Geography("POINT", srid=4326) when you want spheroidal metres from ST_Distance without choosing a projection.

7. Append new batches without dropping the table. if_exists="replace" is destructive in ways that are easy to forget: it drops the table, and with it the GiST index, the ANALYZE statistics, any constraints, and every permission grant you added. For incremental loads, create the table once and append.

from sqlalchemy import text

new_parcels = gpd.read_file("parcels_2026Q2.gpkg").to_crs(epsg=4326)

with engine.begin() as conn:
    conn.execute(text("DELETE FROM zoning WHERE vintage = :v"), {"v": "2026Q2"})

new_parcels.to_postgis("zoning", engine, if_exists="append", index=False, chunksize=10_000)

Deleting the batch's own rows first makes the load idempotent — rerun it after a failure and you get the same table, not duplicates. Wrap both statements in one engine.begin() block if a half-applied reload would be worse than a failed one; the delete and the insert then commit or roll back together.

Two details bite on real deployments. If the table lives outside public, pass the schema explicitly (to_postgis("zoning", engine, schema="planning")) rather than relying on search_path, which differs between your session and the service account's. And keep column names lowercase: PostgreSQL folds unquoted identifiers to lowercase, while GeoPandas quotes them, so a GeoDataFrame column called ZoneCode becomes a case-sensitive "ZoneCode" that every later hand-written query must quote too.

Client-side filtering versus a server-side bounding-box predicate Two lanes share the same source table and the same 318-row answer. In the upper lane no spatial predicate is sent, so the whole 1.2 million row table is serialised across the connection and filtered afterwards in the Python process; the transfer channel is drawn completely full. In the lower lane the double-ampersand envelope test runs inside PostGIS against the GiST index, so only 318 rows are serialised and the same channel is drawn almost empty. The result set is identical in both lanes. Where the bounding box gets applied decides what crosses the connection Filter after the fetch, in pandas zoning table 1.2M rows every geometry serialised · channel full filtered in Python after the transfer 318 rows Filter in SQL with && ST_MakeEnvelope zoning table 1.2M rows 318 rows · only the map window already filtered GiST pruned the rest 318 rows Identical result set — the difference is entirely in how much of the table is serialised and shipped. Push the envelope into the WHERE clause and the index decides what leaves the database.
Both lanes return the same 318 features; only the lane that sends && to the server lets the GiST index prune rows before they are serialised.

Verification

Confirm the SRID survived the round trip and the geometry decoded to real shapely objects rather than strings.

from sqlalchemy import text

with engine.connect() as conn:
    srid = conn.execute(
        text("SELECT Find_SRID('public', 'zoning', 'geometry')")
    ).scalar()
print("Stored SRID:", srid)                       # Stored SRID: 4326
assert srid == 4326

print(type(window.geometry.iloc[0]))              # <class 'shapely.geometry.polygon.Polygon'>
assert window.crs.to_epsg() == 4326
print(f"Read back {len(window)} features")        # Read back 318 features

A Find_SRID result of 0 means the GeoDataFrame had no CRS at write time — go back to step 2. A geometry column whose type() is str means GeoAlchemy2 was not importable in the environment that wrote the table.

Confirm the column type and the row count too, since both can be wrong while the SRID is right. geometry_columns is the catalogue view PostGIS maintains for exactly this purpose, and it reports the type modifier the column was actually created with.

with engine.connect() as conn:
    meta = conn.execute(text("""
        SELECT type, coord_dimension, srid
        FROM geometry_columns
        WHERE f_table_name = 'zoning' AND f_geometry_column = 'geometry'
    """)).one()
    total = conn.execute(text("SELECT count(*) FROM zoning")).scalar()

print(meta, total)          # ('MULTIPOLYGON', 2, 4326) 1204118
assert total == len(zoning), "Row count differs — a chunk failed to commit"

A row count short of the frame's length after a chunked write means a chunk raised and the transaction rolled back partially — check the server log for the geometry-type or SRID mismatch that caused it, rather than re-running the load blind.

Edge Cases & Debugging

Frequently Asked Questions

Should I store the table in EPSG:4326 or in a projected CRS? Store in whichever CRS most clients read, which for a shared table is usually EPSG:4326, and do metric work with an explicit cast or transform at query time — ST_Area(geography(geometry)) for true square metres, or ST_Transform into a local projected SRID for repeated metric analysis. Storing in a projected CRS is the better default only when nearly every query is metric and confined to one region, because then you avoid transforming on every call.

Do I need GeoAlchemy2 if I never import it in my own code? Yes. GeoPandas delegates geometry encoding and decoding to it, so it must be importable in the interpreter that runs the write and the read even though your script never names it. Its absence produces no error — just a column of text where geometries should be.

Is to_postgis fast enough for a multi-million-row load? For a nightly batch, yes, with chunksize set and the spatial index created afterwards. For a bulk initial load of tens of millions of rows, a COPY-based path — ogr2ogr or writing the frame to CSV/WKB and issuing COPY — is substantially faster, because to_postgis goes through parameterised INSERT statements. Use the fast path once, then keep to_postgis for incremental updates.

Why did my query get slower after a big append? The planner is working from stale statistics. ANALYZE after any load that materially changes row counts or the spatial extent, and remember that a replace write drops the index entirely. Spatial Indexing in PostGIS with GiST covers how to confirm the index is being used rather than assuming it.

Can I read a PostGIS table straight into Dask or DuckDB instead of GeoPandas? Yes, and it is the right move once the result set stops fitting in memory. DuckDB Spatial Analytics can attach a PostgreSQL database and query it in place, and Dask can read partitioned ranges in parallel. Both still benefit from pushing the spatial predicate into the PostGIS query rather than pulling rows first.

How do I keep database credentials out of the connection string? Build the URL from environment variables, or point SQLAlchemy at a ~/.pgpass entry by omitting the password from the URL entirely. Never commit a URL containing a password; the engine is created once at startup, so reading two environment variables there costs nothing.

Related