Connecting GeoPandas to PostGIS with SQLAlchemy

This guide wires GeoPandas to a PostGIS database the modern way — a SQLAlchemy 2.0 engine plus GeoAlchemy2 — and writes, indexes, and reads a GeoDataFrame without losing CRS 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.

Why This Approach / What Goes Wrong

GeoPandas' to_postgis and read_postgis need a SQLAlchemy Engine, not a raw DBAPI connection — passing a bare psycopg connection raises a type error or silently drops geometry handling. They also rely on GeoAlchemy2 being importable to register the geometry type; without it, the geometry column is written as opaque text and read back as strings. The third trap is the CRS: if the GeoDataFrame has no .crs, the table gets SRID 0 and every later spatial predicate against a real SRID fails on mismatch.

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

Prerequisites

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

Step-by-Step Implementation

1. Create the engine and enable PostGIS.

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:
    conn.execute(text("CREATE EXTENSION IF NOT EXISTS postgis"))

2. Load a layer and guarantee it has a CRS.

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 a known SRID

3. Write it to a table. GeoPandas creates the geometry column with the right SRID.

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

4. Add a spatial index and refresh planner statistics.

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.

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},
)

Verification

Confirm the SRID survived the round trip and the geometry decoded correctly.

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.

Edge Cases & Debugging

Related