⚙️ How-To PostGIS & databases Batch processing Spatial SQL & PostGIS · step 6 of 15 Updated

How to Load a Folder of Shapefiles into PostGIS with Python

Problem statement

You have 340 shapefiles from a data supplier and a PostGIS database. The loop looks trivial:

for path in Path("data").glob("*.shp"):
    gpd.read_file(path).to_postgis(path.stem, engine, if_exists="replace")

It runs for two hours and then you discover what it actually produced:

  • Some tables are in EPSG:4326 and some in EPSG:27700, because the source files were mixed.
  • Half the tables have no spatial index, so every subsequent query scans.
  • One file was Cheshire East.shp, giving a table called Cheshire East that needs quoting forever after.
  • The MultiPolygon/Polygon split means one table rejects half its own data on the next append.
  • The job died at file 212 and rerunning it starts from zero.
  • Column names were truncated to ten characters by the shapefile format long before Python saw them.

None of these are hard problems individually. The point is that a folder load has to solve all of them at once, and the one-line version solves none.

Quick answer 

from pathlib import Path
import geopandas as gpd
from sqlalchemy import create_engine, text
import re

engine = create_engine("postgresql+psycopg://user@localhost/gis")
SRID = 27700

def safe_name(path):
    name = re.sub(r"[^a-z0-9_]+", "_", path.stem.lower()).strip("_")
    return f"t_{name}" if not name[0].isalpha() else name

for path in sorted(Path("data").glob("*.shp")):
    table = safe_name(path)
    gdf = gpd.read_file(path)
    if gdf.crs is None:
        print(f"  ✗ {path.name}: no CRS, skipped"); continue
    gdf = gdf.to_crs(SRID)
    gdf.columns = [c.lower() for c in gdf.columns]
    gdf.to_postgis(table, engine, if_exists="replace", index=False, chunksize=10_000)
    with engine.begin() as con:
        con.execute(text(f'CREATE INDEX ON "{table}" USING GIST (geometry)'))
        con.execute(text(f'ANALYZE "{table}"'))
    print(f"  ✓ {table:<28} {len(gdf):>8,} rows")
Vertical steps from discovery through CRS check, name sanitising, load, index and verify.
Six steps. The one-line version does step four only.
Step Why it cannot be skipped
audit before loading mixed CRS and geometry types are the norm, not the exception
sanitise table names spaces and capitals mean quoting every query forever
unify the CRS otherwise no two tables can be joined
load in chunks a single 4-million-row insert is one long transaction
index and ANALYZE an unindexed table makes every later query a scan
verify counts a partial load looks exactly like a complete one

Step-by-step solution

1. Audit the folder before writing anything 

from pathlib import Path
import geopandas as gpd, pandas as pd, pyogrio

def audit(folder, pattern="*.shp"):
    rows = []
    for path in sorted(Path(folder).glob(pattern)):
        try:
            info = pyogrio.read_info(path)
            rows.append({
                "file": path.name,
                "rows": info["features"],
                "crs": info["crs"],
                "geom_type": info["geometry_type"],
                "fields": len(info["fields"]),
                "mb": round(sum(p.stat().st_size for p in path.parent.glob(path.stem + ".*")) / 1e6, 1),
            })
        except Exception as exc:
            rows.append({"file": path.name, "error": str(exc)[:60]})
    df = pd.DataFrame(rows)
    print(df.groupby("crs", dropna=False)["rows"].agg(["count", "sum"]))
    print(df["geom_type"].value_counts().to_dict())
    return df

df = audit("data/")
                count      sum
crs
EPSG:27700        318  8412993
EPSG:4326          21   104882
None                1        0
{'MultiPolygon': 294, 'Polygon': 41, 'Point': 5}

pyogrio.read_info reads the header only, so auditing 340 shapefiles takes seconds rather than loading gigabytes. Three findings are already visible: twenty-one files need reprojecting, one has no CRS at all, and the geometry type is split between Polygon and MultiPolygon.

2. Sanitise table names

PostgreSQL folds unquoted identifiers to lower case. A table created as "Cheshire East" can only ever be referenced as "Cheshire East", quotes included, by everyone, forever.

import re

def safe_name(stem, prefix="", max_len=63):
    name = re.sub(r"[^a-z0-9_]+", "_", stem.lower()).strip("_")
    name = re.sub(r"_+", "_", name)
    if not name or not name[0].isalpha():
        name = "t_" + name
    return (prefix + name)[:max_len]        # 63 bytes is Postgres's limit

for stem in ["Cheshire East", "2026-boundaries", "Ward Areas (final)"]:
    print(f"{stem!r:<24} → {safe_name(stem)}")
'Cheshire East'          → cheshire_east
'2026-boundaries'        → t_2026_boundaries
'Ward Areas (final)'     → ward_areas_final

Watch for collisions: Ward-Areas and Ward Areas both sanitise to ward_areas. Check before loading rather than discovering it when the second file silently replaces the first.

from collections import Counter
names = [safe_name(p.stem) for p in Path("data").glob("*.shp")]
clashes = {n: c for n, c in Counter(names).items() if c > 1}
if clashes:
    raise SystemExit(f"table name collisions: {clashes}")

3. Unify the CRS and the geometry type 

TARGET_SRID = 27700

def prepare(gdf, srid=TARGET_SRID, promote_multi=True):
    if gdf.crs is None:
        raise ValueError("no CRS — cannot load safely")
    if gdf.crs.to_epsg() != srid:
        gdf = gdf.to_crs(srid)

    gdf = gdf[~gdf.geometry.isna() & ~gdf.geometry.is_empty].copy()

    if promote_multi:
        from shapely.geometry import MultiPolygon, MultiLineString
        types = set(gdf.geom_type)
        if types <= {"Polygon", "MultiPolygon"} and len(types) > 1:
            gdf["geometry"] = [
                g if g.geom_type == "MultiPolygon" else MultiPolygon([g])
                for g in gdf.geometry
            ]
        elif types <= {"LineString", "MultiLineString"} and len(types) > 1:
            gdf["geometry"] = [
                g if g.geom_type == "MultiLineString" else MultiLineString([g])
                for g in gdf.geometry
            ]
    gdf.columns = [c.lower() for c in gdf.columns]
    return gdf

The mixed-type promotion matters because to_postgis creates a typed geometry column from whatever it finds first. A file whose first row is a Polygon gets a geometry(Polygon, 27700) column, and the first MultiPolygon row then fails:

Geometry type (MultiPolygon) does not match column type (Polygon)

Promoting the other way round — a single Polygon into a geometry(MultiPolygon) column — is accepted by PostGIS, which promotes it silently. Only MultiPolygon into geometry(Polygon) is rejected. Promoting everything to multi is therefore the safe direction. See PostGIS write fails on SRID or geometry type.

4. Load, then index, then analyse — in that order 

from sqlalchemy import text

def load_one(gdf, table, engine, srid=TARGET_SRID, chunksize=10_000):
    gdf.to_postgis(table, engine, if_exists="replace",
                   index=False, chunksize=chunksize)
    with engine.begin() as con:
        con.execute(text(f'CREATE INDEX "{table}_geom_idx" ON "{table}" USING GIST (geometry)'))
        con.execute(text(f'ANALYZE "{table}"'))
    return table

The order is not stylistic. if_exists="replace" drops the table and every index on it, so an index created beforehand is discarded — and even if it survived, maintaining it through the insert would cost more than building it once at the end. ANALYZE gives the planner the statistics it needs to actually choose the new index; without it, a fresh index is often ignored. See PostGIS spatial indexes explained.

chunksize bounds how much is in flight at once. Without it, a 4-million-row frame becomes one enormous statement and one very long transaction.

5. Make the job resumable

A folder load is a batch job, and it should behave like one:

def already_loaded(engine, table):
    with engine.begin() as con:
        return con.execute(text(
            "SELECT to_regclass(:t) IS NOT NULL"), {"t": f"public.{table}"}).scalar()

Skipping existing tables turns a two-hour job that died at file 212 into a twenty-minute resume. The general pattern is in how to build a resumable batch GIS job.

6. Verify 

def verify(engine, table, expected_rows, srid=TARGET_SRID):
    with engine.begin() as con:
        got = con.execute(text(f'SELECT COUNT(*) FROM "{table}"')).scalar()
        srids = con.execute(text(
            f'SELECT DISTINCT ST_SRID(geometry) FROM "{table}"')).scalars().all()
        invalid = con.execute(text(
            f'SELECT COUNT(*) FROM "{table}" WHERE NOT ST_IsValid(geometry)')).scalar()
        has_idx = con.execute(text(
            "SELECT COUNT(*) FROM pg_indexes WHERE tablename = :t "
            "AND indexdef LIKE '%gist%'"), {"t": table}).scalar()
    problems = []
    if got != expected_rows: problems.append(f"{expected_rows:,} in → {got:,} out")
    if srids != [srid]:      problems.append(f"SRIDs {srids}")
    if not has_idx:          problems.append("no GiST index")
    if invalid:              problems.append(f"{invalid:,} invalid geometries")
    return problems

Row count in versus out is the check that catches a partial load, which otherwise looks exactly like a complete one.

Panels contrasting what a shapefile declares with what a PostGIS table requires.
Every file has an opinion. The table needs one answer.

Code examples

Example 1: the complete loader 

from pathlib import Path
from collections import Counter
import re, time
import geopandas as gpd, pandas as pd, pyogrio
from sqlalchemy import create_engine, text

def load_folder(folder, engine, *, srid=27700, pattern="*.shp",
                prefix="", resume=True, chunksize=10_000):
    paths = sorted(Path(folder).glob(pattern))
    names = [safe_name(p.stem, prefix) for p in paths]
    clashes = {n: c for n, c in Counter(names).items() if c > 1}
    if clashes:
        raise ValueError(f"table name collisions: {clashes}")

    results = []
    for path, table in zip(paths, names):
        t0 = time.perf_counter()
        if resume and already_loaded(engine, table):
            results.append({"file": path.name, "table": table, "status": "skipped"})
            continue
        try:
            expected = pyogrio.read_info(path)["features"]
            gdf = prepare(gpd.read_file(path), srid=srid)
            load_one(gdf, table, engine, srid=srid, chunksize=chunksize)
            problems = verify(engine, table, len(gdf), srid=srid)
            results.append({
                "file": path.name, "table": table, "rows": len(gdf),
                "dropped": expected - len(gdf),
                "status": "ok" if not problems else "warn",
                "problems": "; ".join(problems),
                "seconds": round(time.perf_counter() - t0, 1),
            })
        except Exception as exc:
            results.append({"file": path.name, "table": table,
                            "status": "failed", "problems": str(exc)[:120]})

    df = pd.DataFrame(results)
    marks = {"ok": "✓", "warn": "!", "skipped": "·", "failed": "✗"}
    for r in results:
        print(f"  {marks[r['status']]} {r['table']:<30} "
              f"{r.get('rows', 0):>9,}  {r.get('problems', '')}")
    print(df["status"].value_counts().to_dict())
    return df

engine = create_engine("postgresql+psycopg://user@localhost/gis")
report = load_folder("data/", engine, srid=27700)
report.to_csv("load_report.csv", index=False)
  ✓ cheshire_east                    182,447
  ✓ cheshire_west                    204,118
  ! stockport                        118,204  412 invalid geometries
  ✗ ward_areas_final                       0  no CRS — cannot load safely
  · trafford                               0
{'ok': 316, 'warn': 12, 'skipped': 8, 'failed': 4}

Four statuses, not two. "Warn" means the table loaded but something needs attention; collapsing it into "ok" hides 412 invalid geometries that will break the next spatial join, and collapsing it into "failed" would send someone re-running a load that worked. dropped records the difference between the header count and the loaded count — null and empty geometries removed by prepare, which you want to know about rather than discover later.

Example 2: many files into one table

Often the 340 files are tiles of one dataset, and the useful output is one table:

def append_folder(folder, table, engine, *, srid=27700,
                  source_col="source_file", pattern="*.shp"):
    """Load every file into a single table, tagged with its source."""
    paths = sorted(Path(folder).glob(pattern))
    total = 0
    for i, path in enumerate(paths):
        gdf = prepare(gpd.read_file(path), srid=srid)
        gdf[source_col] = path.stem
        gdf.to_postgis(table, engine,
                       if_exists="replace" if i == 0 else "append",
                       index=False, chunksize=10_000)
        total += len(gdf)
        print(f"  {path.name:<28} +{len(gdf):>8,}  total {total:>10,}")

    with engine.begin() as con:
        con.execute(text(f'CREATE INDEX ON "{table}" USING GIST (geometry)'))
        con.execute(text(f'CREATE INDEX ON "{table}" ({source_col})'))
        con.execute(text(f'ANALYZE "{table}"'))
    return total

append_folder("tiles/", "parcels", engine)

Three things make this work where the naive version fails. replace on the first file and append afterwards creates the schema once — using replace throughout leaves you with only the last file. The source_file column preserves provenance that the merge would otherwise destroy, and makes a single tile reloadable with a DELETE WHERE source_file = …. And the indexes are built once at the end rather than maintained across 340 inserts.

Column mismatches are the remaining risk: appending a frame with different columns raises. Align them first if the supplier's schema drifts between tiles.

Example 3: ogr2ogr when Python is not the right tool

For a straight bulk load with no per-file logic, GDAL's own loader is faster because it streams and never builds a DataFrame:

import subprocess
from pathlib import Path

def ogr_load(path, table, dsn, srid=27700):
    cmd = [
        "ogr2ogr", "-f", "PostgreSQL", f"PG:{dsn}", str(path),
        "-nln", table,
        "-nlt", "PROMOTE_TO_MULTI",       # solves the Polygon/MultiPolygon split
        "-t_srs", f"EPSG:{srid}",
        "-lco", "GEOMETRY_NAME=geom",
        "-lco", "SPATIAL_INDEX=GIST",
        "-lco", "PRECISION=NO",
        "-overwrite",
        "--config", "PG_USE_COPY", "YES",  # COPY instead of INSERT — much faster
    ]
    proc = subprocess.run(cmd, capture_output=True, text=True)
    if proc.returncode != 0:
        raise RuntimeError(proc.stderr.strip()[:300])
    return table

DSN = "host=localhost dbname=gis user=gis"
for path in sorted(Path("data").glob("*.shp")):
    try:
        ogr_load(path, safe_name(path.stem), DSN)
        print(f"  ✓ {path.name}")
    except RuntimeError as exc:
        print(f"  ✗ {path.name}: {exc}")

PG_USE_COPY=YES switches from per-row INSERT to COPY, which is typically three to five times faster on a large load. PROMOTE_TO_MULTI handles the geometry-type split in the loader rather than in your code, and SPATIAL_INDEX=GIST builds the index as part of table creation.

The trade-off is real: ogr2ogr is faster and knows more about formats, but you cannot inspect or clean the data on the way through. Use it when the files are trustworthy, and the Python route when they need cleaning — which, with supplier data, is most of the time.

Explanation

Panels comparing indexing before a bulk load against indexing after it.
An index maintained through four million inserts is rebalanced four million times.

A folder load is a schema negotiation, and that is what makes it harder than it looks. Every shapefile carries its own implicit schema — field names, types, CRS, geometry type — and the database needs one explicit schema per table. Loading one file lets the file decide. Loading 340 means reconciling 340 opinions, and the reconciliation has no default.

The shapefile format makes this worse in specific, well-known ways. Field names are limited to ten characters, so population_density arrives as populatio — truncated by the format long before Python sees it, and unrecoverable without the supplier's documentation. The CRS lives in a separate .prj file that is routinely lost in transit, giving crs = None. There is no distinction between Polygon and MultiPolygon at the file level, so a single file can hold both. And there is no NULL: empty strings and zeros stand in for missing values. Every one of these becomes a database problem the moment you load without checking. GIS vector file formats compared covers why GeoPackage avoids most of them.

Typed geometry columns are where the format's laxity meets the database's strictness. PostGIS's geometry(MultiPolygon, 27700) is a constraint, and constraints reject data. Confirmed against PostGIS 3.4: a single Polygon inserted into a geometry(MultiPolygon, …) column is accepted and promoted, while a MultiPolygon into a geometry(Polygon, …) column is rejected with Geometry type (MultiPolygon) does not match column type (Polygon). The asymmetry is why promoting everything to multi is the safe direction, and why -nlt PROMOTE_TO_MULTI exists.

Indexing after loading rather than during is not a micro-optimisation. A GiST index maintained across four million inserts is rebalanced continuously; built once at the end it is a single bulk operation on data already in place. The difference is routinely an order of magnitude. The same reasoning applies to constraints and foreign keys, and it is why ogr2ogr's SPATIAL_INDEX=GIST creates the index as part of table creation rather than before the copy.

Finally, treat this as a batch job rather than a script. It has the same three requirements as any batch GIS job: per-file error isolation so one bad file does not sink the run, resumability so a failure at file 212 does not discard 211 successes, and a report that distinguishes "loaded", "loaded with problems" and "not loaded". The database part is the easy half.

Edge cases or notes

  • Shapefile field names are truncated to 10 characters by the format. Nothing downstream can recover them; check the supplier's documentation.
  • A missing .prj means crs = None. Do not guess — assigning a CRS relabels coordinates without moving them.
  • Encoding matters. Shapefiles have no reliable encoding declaration; pass encoding="latin-1" or check the .cpg file. See garbled attribute text.
  • if_exists="replace" drops indexes and constraints along with the table.
  • to_postgis names the geometry column geometry, while ogr2ogr defaults to wkb_geometry unless you pass GEOMETRY_NAME=geom. Pick one convention.
  • Reserved words — a table or column called order, user or default needs quoting forever. Sanitise those too.
  • 63 bytes is the identifier limit. Longer names are silently truncated, which can create collisions.
  • Loading into a schema needs schema="raw" on to_postgis and matching qualification in your index DDL.
  • PG_USE_COPY=YES makes ogr2ogr three to five times faster on large loads.
  • Wrap per-file work in its own transaction so a failure rolls back one table, not the whole run.

FAQ

Why does my table name need quotes in every query?

It was created with capitals or spaces. PostgreSQL folds unquoted identifiers to lower case, so "Cheshire East" can only be referenced with quotes. Sanitise names before loading.

Why does the load fail with "Geometry type does not match column type"?

The table has a typed column such as geometry(Polygon, 27700) and the file contains MultiPolygon rows. Promote everything to multi before loading, or use -nlt PROMOTE_TO_MULTI.

Should I create the spatial index before or after loading?

After. if_exists="replace" drops it anyway, and maintaining an index through a bulk insert costs far more than building it once. Follow with ANALYZE.

How do I make the load resumable?

Check whether the target table already exists and skip it. SELECT to_regclass('public.tablename') IS NOT NULL is a cheap test that turns a restart into a resume.

Is ogr2ogr faster than GeoPandas?

Yes, substantially — it streams and can use COPY. Use it when the files need no inspection, and Python when they need cleaning or per-file logic.

What if the files have different CRS?

Reproject each to a single target SRID as it loads. Tables in different SRIDs cannot be joined without a per-query transform, and PostGIS raises rather than returning nothing.

Why are my column names truncated to ten characters?

The shapefile format truncates them on write. The full names never reached Python, and no code can recover them — check the supplier's documentation.