Loading GeoDataFrames into PostGIS with GeoAlchemy2

One-sentence answer: assert the CRS, promote geometry types to the column’s declared type, stream the rows through COPY into an unlogged staging table as WKB hex, then do one insert-select with ON CONFLICT — and commit the manifest row in the same transaction.

Context

Getting vector data into PostGIS is the least interesting part of a farm platform and one of the easiest to get quietly wrong. The three failures are all silent until much later: a layer loaded with the wrong SRID (every subsequent spatial join returns nothing), a Polygon where MultiPolygon was declared (rows rejected in a batch that was reported as successful), and a load committed without a manifest row (a retry produces a second copy of a boundary, and every area figure doubles).

This guide implements the loading side of PostGIS schema design for farm data.

Row-by-row inserts versus COPY into staging The upper path shows a GeoDataFrame inserted row group by row group through the driver, taking about ninety seconds for one hundred thousand geometries with no upsert available. The lower path shows the same data written as WKB hex through COPY into an unlogged staging table in about four seconds, followed by one insert-select with an on-conflict clause in about three seconds. Default writer — many round trips GeoDataFrame 100,000 geometries batched INSERTs via the driver ≈ 90 s · no upsert clause available COPY into staging, then one statement WKB hex buffer in memory, streamed COPY → unlogged staging ≈ 4 s INSERT … SELECT … ON CONFLICT ≈ 3 s · idempotent ≈ 7 s total

Prerequisites

Beyond the parent topic’s stack: psycopg 3.2.* (for its copy support), SQLAlchemy 2.0.* and GeoAlchemy2 0.15.*. The target table and its constraints must already exist — this is a loading guide, not a migration one.

Step-by-step

1. Validate in Python first. CRS, geometry validity and type, before a byte crosses the wire.

One load, one transaction Four steps: assert the coordinate reference system, geometry validity and type in Python; stream the rows as well-known-binary hex into an unlogged staging table with COPY; build the geometry and apply the conflict rule in a single insert-select; and commit the load together with its manifest row. Validate in Python CRS, validity, type COPY to staging WKB hex, unlogged INSERT … SELECT build geometry, set SRID ON CONFLICT DO NOTHING Commit load and manifest row in one transaction Validating in Python reports every bad row at once; a database constraint reports one and rolls the transaction back.

2. COPY into an unlogged staging table — unlogged because the data is about to be copied into a logged table anyway, and skipping the write-ahead log roughly halves the time.

3. One insert-select that constructs the geometry, sets the SRID and applies the conflict rule.

4. Commit the manifest row in the same transaction.

5. Read a sample back and compare.

PYTHON 
import io

import geopandas as gpd
from shapely.geometry import MultiPolygon, Polygon
from sqlalchemy import create_engine, text


def prepare(gdf: gpd.GeoDataFrame, epsg: int = 4326) -> gpd.GeoDataFrame:
    """Assert and normalise before the database ever sees the data."""
    assert gdf.crs is not None, "GeoDataFrame has no CRS — refusing to guess"
    if gdf.crs.to_epsg() != epsg:
        gdf = gdf.to_crs(epsg=epsg)

    invalid = ~gdf.geometry.is_valid
    assert not invalid.any(), (
        f"{int(invalid.sum())} invalid geometr(ies); repair with make_valid before loading — "
        "see the invalid-geometry guide")
    assert not gdf.geometry.is_empty.any(), "empty geometries in the layer"

    # Promote to the multi form so a MultiPolygon column accepts every row.
    gdf = gdf.copy()
    gdf["geometry"] = gdf.geometry.apply(
        lambda g: MultiPolygon([g]) if isinstance(g, Polygon) else g)
    kinds = set(gdf.geom_type)
    assert kinds == {"MultiPolygon"}, f"mixed geometry types after promotion: {kinds}"
    return gdf


def load_boundaries(gdf: gpd.GeoDataFrame, dsn: str, source: str, checksum: str) -> int:
    """Bulk load field boundaries and record the load in the manifest, atomically."""
    gdf = prepare(gdf)

    buf = io.StringIO()
    frame = gdf.assign(wkb=gdf.geometry.to_wkb(hex=True))[["field_id", "valid_from", "wkb"]]
    frame.to_csv(buf, index=False, header=False)
    buf.seek(0)

    engine = create_engine(dsn, future=True)
    with engine.begin() as conn:                       # one transaction for load + manifest
        raw = conn.connection.driver_connection
        with raw.cursor() as cur:
            cur.execute("""
                CREATE UNLOGGED TABLE IF NOT EXISTS stage_boundary
                    (field_id bigint, valid_from timestamptz, wkb text)
            """)
            cur.execute("TRUNCATE stage_boundary")
            with cur.copy("COPY stage_boundary (field_id, valid_from, wkb) "
                          "FROM STDIN WITH (FORMAT csv)") as cp:
                cp.write(buf.read())

        inserted = conn.execute(text("""
            INSERT INTO field_boundary (field_id, geom, source, valid_from)
            SELECT s.field_id,
                   ST_SetSRID(ST_GeomFromWKB(decode(s.wkb, 'hex')), 4326),
                   :source,
                   s.valid_from
            FROM stage_boundary s
            ON CONFLICT DO NOTHING
        """), {"source": source}).rowcount

        conn.execute(text("""
            INSERT INTO ingest_manifest (source, external_id, checksum)
            VALUES (:source, :ext, :sum)
            ON CONFLICT (source, external_id, checksum) DO NOTHING
        """), {"source": source, "ext": f"boundaries:{source}", "sum": checksum})
    return inserted

Inline verification — read a sample back and compare geometry and SRID, not just row counts:

PYTHON 
engine = create_engine(DSN, future=True)
back = gpd.read_postgis(
    "SELECT field_id, geom FROM field_boundary WHERE source = %(s)s AND valid_to IS NULL",
    engine, geom_col="geom", params={"s": source})

assert len(back) == len(gdf), f"loaded {len(gdf)} rows, read back {len(back)}"
assert back.crs.to_epsg() == 4326, f"round-tripped CRS is {back.crs} — SRID was not set"

a = gdf.set_index("field_id").to_crs(32615).geometry.area
b = back.set_index("field_id").to_crs(32615).geometry.area
rel = ((a - b).abs() / a).max()
assert rel < 1e-9, f"geometry changed on the round trip (max relative area error {rel:.2e})"
print(f"{len(back)} boundaries round-tripped, max relative area error {rel:.2e}")

Settings worth being deliberate about

Setting Default here Why it matters
Canonical SRID 4326 One SRID for every geometry column in the schema. Mixed SRIDs make spatial joins return zero rows with no error, which is the most expensive silent failure in the whole stack
Geometry encoding WKB hex Avoids every CSV quoting problem that WKT introduces with embedded commas and line wraps in long coordinate lists
Staging table UNLOGGED Skips the write-ahead log for data that is about to be copied into a logged table anyway — commonly half the load time
Batch size ~50,000 rows Larger batches inflate the staging table and the transaction’s lock footprint; smaller ones lengthen the load without reducing risk
Geometry type MultiPolygon Promote on the way in. Farm boundaries mix single and multi forms constantly, and relaxing the column type instead removes a useful constraint
Conflict rule DO NOTHING on the natural key What makes a re-run harmless — provided the unique index it needs actually exists

Do the validation in Python, before the wire. A failed constraint reports one offending row at a time and rolls the transaction back; a GeoDataFrame assertion reports how many rows are wrong, which ones, and why, in a single message.

Gotchas and edge cases

  • ST_SetSRID does not reproject — it relabels. Applying it to coordinates that are actually in UTM tells PostGIS they are degrees, and everything afterwards is wrong by hundreds of kilometres with no error anywhere. Reproject in Python with a known CRS, then set the SRID to match.
Four loading mistakes that do not raise an error A table of four bulk-loading mistakes, each of which produces no exception: relabelling a geometry's SRID instead of reprojecting it, using an on-conflict clause with no unique index to conflict on, inserting a plain polygon into a multi-polygon column, and skipping ANALYZE after a load. Mistake Symptom Guard ST_SetSRID on unprojected metres relabels, does not reproject Geometry lands in the wrong hemisphere Reproject in Python first ON CONFLICT with no unique index the clause is a no-op Re-runs quietly duplicate rows Assert the constraint exists Polygon into a MultiPolygon column mixed geometry types Batch rejected, load reported as done Promote to multi on the way in No ANALYZE after the load statistics from an empty table First spatial query scans everything ANALYZE the partition

  • ON CONFLICT DO NOTHING needs a constraint to conflict on. Without a unique index covering the natural key, it is a no-op and a re-run inserts duplicates cheerfully. Confirm the constraint exists before relying on the clause; this is where the boundary exclusion constraint from the parent topic earns its place.

  • Unlogged staging tables do not survive a crash. That is the point — but it also means a transaction that spans a server restart loses the staging data and the load has to start again. Keep each load in one transaction and one process.

  • CSV quoting and WKB hex are a good pair. Hex encoding avoids embedded commas, quotes and newlines entirely, which is what makes the COPY path safe for arbitrary geometries. Passing WKT through CSV works until a coordinate list wraps a line.

  • Very large geometries hit COPY line limits. A single field boundary with a hundred thousand vertices produces an enormous hex string. Simplify machine-generated boundaries before loading — the tolerance should come from the positional error budget in RTK GPS accuracy and positional error budgets, not from taste.

  • ANALYZE after the load. A fresh table with no statistics gets a plan built for an empty relation, which turns the first spatial query into a sequential scan.

This guide is part of PostGIS Schema Design for Farm Data — see there for the tables, constraints and manifest this load writes into.

Parent section PostGIS Schema Design for Farm Data Related Spatial Indexing and Partitioning Yield Points in PostGIS