Automating Shapefile Cleanup with Python

Legacy shapefiles arrive with missing .prj files, self-intersecting rings, and null geometries that crash the first spatial join downstream — this guide builds one reusable function that enforces a coordinate reference system, repairs topology, and prunes bad records before ingestion. It is for anyone scripting a batch pipeline over folders of untrusted vendor data with GeoPandas and Shapely. It sits under Shapefile & GeoJSON Parsing in Geospatial Data Ingestion & Processing Workflows.

Why This Approach / What Goes Wrong

The naive cleanup — open the file in desktop GIS, run "repair geometry", export — is manual, unrepeatable, and invisible to version control. A scripted GeoPandas function makes the whole operation deterministic: the same dirty input always produces the same clean output, and every transformation is auditable in a diff. Three specific failures are what a good cleanup routine has to defend against, and each has a subtle wrong answer.

• A missing CRS is not the same as a wrong CRS. When gdf.crs is None, the shapefile shipped without a .prj sidecar. The trap is reaching for to_crs() to "fix" it — but to_crs() reprojects coordinates and needs a source CRS to do so, and on None it raises. You must first declare the true source CRS with set_crs() (metadata only, no coordinate change), and only then reproject. Confusing the two silently shifts every coordinate.
• make_valid can change a geometry's type. Repairing a self-intersecting Polygon frequently yields a MultiPolygon, or even a GeometryCollection mixing lines and polygons. A pipeline that assumes single-part polygons downstream will break on the repaired output unless you explode or filter by type. The deeper mechanics of this live in fixing self-intersecting polygons programmatically.
• The ESRI Shapefile format silently truncates field names to 10 characters. Write a column called population_density back to .shp and it becomes populatio_ with no error. If two columns collide after truncation, one is lost. This is a format limitation, not a bug you can code around — you rename before writing, or you write to a format without the limit.
• The attribute table is dBASE, and dBASE has almost no type system. A .dbf column is character, numeric, date, or logical — nothing else. There is no timestamp, so a datetime64 column loses its time component on write; there is no 64-bit integer, so an identifier above roughly two billion is either widened to a float and loses its last digits or is rejected; and there is no NULL for numerics, so a missing value is written as zero and becomes indistinguishable from a genuine zero. A cleanup routine that only touches geometry hands the next stage attributes that are quietly wrong in a way no validity check will catch.
• A Shapefile layer holds one geometry type, and repair can violate that. The .shp header declares a single shape type for the whole file. Feed to_file a frame that mixes Polygon and MultiPolygon — exactly what step 2 produces — and the driver promotes or rejects depending on the GDAL build, so the safe move is to decide the output type yourself rather than let the writer decide for you.

Prerequisites

• geopandas>=0.14 — the DataFrame wrapper, file I/O, and vectorized is_valid/is_empty accessors
• shapely>=2.0 — provides make_valid and the GEOS 3.8+ topology engine it delegates to
• pyproj>=3.4 — resolves CRS definitions and datum transforms under to_crs()

conda install -c conda-forge "geopandas=0.14.*" "shapely=2.0.*" "pyproj=3.4.*"

Install from conda-forge, not pip, so the GDAL, GEOS, and PROJ C libraries underneath these bindings stay ABI-compatible — the mismatched-DLL failures that plague pip installs are the same ones covered across Shapefile & GeoJSON Parsing.

Step-by-Step Implementation

1. Load the file and gate on the CRS. Read the shapefile, then branch on whether a CRS is present. Assign a known source CRS with set_crs() when it is missing, and reproject with to_crs() when it is present but wrong. The example targets EPSG:25832 (ETRS89 / UTM 33N), a metric projected grid appropriate for area and distance work — never use Web Mercator (EPSG:3857) for measurement.

import warnings
import geopandas as gpd

def load_and_align_crs(input_path: str, source_crs: str, target_crs: str) -> gpd.GeoDataFrame:
    """Load a shapefile and guarantee it lands in target_crs.

    source_crs is the CRS the data is ACTUALLY in — used only when the
    file has no .prj. It never reprojects on its own.
    """
    parcels = gpd.read_file(input_path)
    target_epsg = int(target_crs.split(":")[1])

    if parcels.crs is None:
        warnings.warn(
            f"No .prj found. Declaring source CRS {source_crs} as metadata. "
            "Verify this matches the true projection before trusting coordinates."
        )
        parcels = parcels.set_crs(source_crs)      # declare — no coordinate change

    if parcels.crs.to_epsg() != target_epsg:
        parcels = parcels.to_crs(target_crs)       # reproject — coordinates change

    return parcels

Keep set_crs and to_crs conceptually separate: the first fixes metadata, the second transforms coordinates. Getting this wrong is the most common CRS error on this site — the full treatment is in Coordinate Reference System Transformations and Coordinate Systems with PyProj.

2. Repair invalid topology in place. Build a boolean mask of invalid rows and apply make_valid only to those — running it over already-valid geometries wastes time and can needlessly reorder vertices. make_valid resolves self-intersections and bad ring orientation without discarding data.

from shapely import make_valid   # shapely 2.x top-level import

def repair_topology(parcels: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    invalid_mask = ~parcels.geometry.is_valid
    n_invalid = int(invalid_mask.sum())
    if n_invalid:
        print(f"Repairing {n_invalid} invalid geometries...")
        parcels.loc[invalid_mask, "geometry"] = parcels.loc[
            invalid_mask, "geometry"
        ].apply(make_valid)
    return parcels

3. Prune null and empty features. A null geometry (no .shp record) and an empty geometry (a valid but zero-area/zero-length shape) are different states, and both break renderers and joins. Drop the nulls with dropna, then filter the empties with the vectorized is_empty predicate.

def prune_features(parcels: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    before = len(parcels)
    parcels = parcels.dropna(subset=["geometry"])
    parcels = parcels[~parcels.geometry.is_empty]
    print(f"Dropped {before - len(parcels)} null/empty features")
    return parcels

4. Normalize field names, then export. Rename any column longer than 10 characters before writing to .shp, otherwise ESRI truncation mangles it silently. If you control the downstream stack, prefer writing to a format without the 10-character ceiling — see GeoParquet vs Shapefile for storage.

def export_clean(parcels: gpd.GeoDataFrame, output_path: str) -> None:
    rename_map = {
        "population_density": "pop_dens",
        "land_use_classification": "landuse",
    }
    parcels = parcels.rename(columns=rename_map)

    long_names = [c for c in parcels.columns if c != "geometry" and len(c) > 10]
    if long_names:
        warnings.warn(f"Fields will be truncated by ESRI Shapefile: {long_names}")

    parcels.to_file(output_path, driver="ESRI Shapefile")
    print(f"Cleanup complete. {len(parcels)} features saved to {output_path}")

5. Settle the geometry type before the writer does. After make_valid a polygon layer typically contains a mixture of Polygon and MultiPolygon, and repairs on badly-noded input can leave a GeometryCollection holding stray lines or points. Decide explicitly: promote everything to multi-part, which keeps every feature and every attribute row intact, or explode to single parts, which changes the row count and duplicates attributes across the pieces. Promotion is the right default for a cleanup pass, because it is the only option that preserves the one-row-per-record relationship the vendor's identifiers assume.

from shapely.geometry import MultiPolygon
from shapely.geometry.base import BaseMultipartGeometry

def homogenise_geometry(parcels: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Force a single, writable geometry type: MultiPolygon for every row."""
    collections = parcels.geom_type == "GeometryCollection"
    if collections.any():
        # Keep only the polygonal members; lines and points are repair debris
        parcels.loc[collections, "geometry"] = parcels.loc[
            collections, "geometry"
        ].apply(lambda g: MultiPolygon([p for p in g.geoms
                                        if p.geom_type == "Polygon"]))

    parcels["geometry"] = parcels.geometry.apply(
        lambda g: g if isinstance(g, BaseMultipartGeometry) else MultiPolygon([g])
    )
    assert parcels.geom_type.eq("MultiPolygon").all()
    return parcels

6. Coerce the attribute table into what dBASE can actually store. This is the step most cleanup scripts skip, and the one that produces the subtlest downstream bugs. Convert timestamps to dates, widen or stringify oversized integers, map booleans to a two-value integer flag, replace vendor sentinel values with genuine nulls, and strip the trailing spaces that fixed-width .dbf character fields carry.

import numpy as np
import pandas as pd

SENTINELS = ["-9999", "-999", "N/A", "NULL", "", "  "]

def coerce_attributes(parcels: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    for col in parcels.columns:
        if col == "geometry":
            continue
        series = parcels[col]

        if pd.api.types.is_datetime64_any_dtype(series):
            parcels[col] = series.dt.date          # dBASE has no time component
        elif pd.api.types.is_bool_dtype(series):
            parcels[col] = series.astype("int32")  # logical -> 0/1 numeric
        elif pd.api.types.is_integer_dtype(series):
            if series.abs().max() > np.iinfo(np.int32).max:
                parcels[col] = series.astype(str)  # keep every digit; widen later
        elif pd.api.types.is_object_dtype(series):
            parcels[col] = (series.astype("string")
                                  .str.strip()
                                  .replace(SENTINELS, pd.NA))
    return parcels

Coercing before the write means the rules are visible in code and in a diff, rather than being applied invisibly by whichever GDAL version the machine happens to have. Where the downstream consumer is under your control, the better answer is not to coerce at all — write GeoPackage or GeoParquet, which carry real types, and reserve the Shapefile export for the one consumer that demands it.

7. Compose the full pipeline and record what changed. Chain the stages into one entry point you can call over a folder of vendor deliveries. Return a per-file record of the repairs, because "the cleanup ran" is not the same claim as "nothing needed cleaning" — a delivery whose invalid-geometry count jumps from 4 to 4,000 is a supplier problem, and you only see it if the numbers are written down.

from pathlib import Path
import pandas as pd

def clean_shapefile(input_path: str, output_path: str,
                    source_crs: str = "EPSG:25832",
                    target_crs: str = "EPSG:25832") -> dict:
    parcels = load_and_align_crs(input_path, source_crs, target_crs)
    n_read = len(parcels)
    n_invalid = int((~parcels.geometry.is_valid).sum())

    parcels = repair_topology(parcels)
    parcels = prune_features(parcels)
    parcels = homogenise_geometry(parcels)
    parcels = coerce_attributes(parcels)
    export_clean(parcels, output_path)

    return {
        "source": Path(input_path).name,
        "features_in": n_read,
        "features_out": len(parcels),
        "repaired": n_invalid,
        "pruned": n_read - len(parcels),
        "target_epsg": target_crs,
    }

audit = [clean_shapefile(str(shp), f"cleaned/{shp.stem}_clean.shp")
         for shp in sorted(Path("raw_deliveries").glob("*.shp"))]
pd.DataFrame(audit).to_csv("cleaned/_audit.csv", index=False)

Verification

Re-open the written output and assert the three invariants the pipeline promises: the CRS is exactly the target, every geometry is valid, and no null or empty records survived.

import geopandas as gpd

cleaned = gpd.read_file("cleaned/berlin_parcels_clean.shp")

assert cleaned.crs.to_epsg() == 25832, f"Wrong CRS: {cleaned.crs.to_epsg()}"
assert cleaned.geometry.is_valid.all(), "Invalid geometries remain"
assert not cleaned.geometry.is_empty.any(), "Empty geometries remain"
assert cleaned.geometry.notna().all(), "Null geometries remain"

print(f"OK: {len(cleaned)} features, all valid, CRS={cleaned.crs.to_epsg()}")
# OK: 4127 features, all valid, CRS=25832

If all four assertions pass, the file is safe to hand to a spatial join, a PostGIS load, or a web-map tiler.

Edge Cases & Debugging

• to_crs() raises "Cannot transform naive geometries". The file had no .prj and you skipped set_crs(). Declare the true source CRS first; there is no way to reproject a geometry whose starting CRS is unknown.
• make_valid split my polygon into a MultiPolygon. Expected behaviour for self-intersecting input. Append .explode(index_parts=False) after repair if downstream code requires single-part features, or filter with parcels.geom_type == "Polygon".
• GEOSException on a single row. One pathological geometry aborts the batch. Wrap make_valid in try/except, log the failing index to a side CSV, and continue — never let one bad record sink the whole delivery.
• Files over ~500 MB exhaust RAM. GeoPandas loads the entire table into memory. Stream feature-by-feature with a fiona iterator, or move the job to dask-geopandas for out-of-core, partitioned execution.
• Reprojection looks shifted by tens of metres. A datum mismatch, not a bug — an ED50 or NAD27 source reprojected as though it were WGS84/ETRS89. Confirm the true source datum before declaring it with set_crs(); the axis-order and datum traps are detailed in Coordinate Reference System Transformations.
• unable to open .shx halts the batch on one delivery. The offset index was lost in a partial copy. Set os.environ["SHAPE_RESTORE_SHX"] = "YES" once at the top of the batch script so GDAL rebuilds it from the .shp for every file, rather than fixing deliveries one at a time.
• The deliveries arrive as ZIP archives. Do not unpack them into a scratch directory you then have to clean up. GDAL reads inside archives directly through its virtual filesystem: gpd.read_file("/vsizip/raw_deliveries/parcels_2026Q1.zip/parcels.shp"), and /vsizip//vsicurl/https://... reads one straight off a web server without landing it on disk at all.
• Identifiers lose their last digits after the round trip. A 64-bit parcel key went through a .dbf numeric field and came back as a float. Cast oversized integers to string before writing, or move the output to GeoPackage.
• Accented place names come back as ü. The .dbf was written in a legacy codepage with no .cpg sidecar. Pass encoding="cp1252" on read for that delivery — the full diagnosis is in Shapefile & GeoJSON Parsing.
• A delivery over ~500 MB exhausts RAM mid-batch. GeoPandas materialises the whole table. Convert those files with ogr2ogr first, or reproject and repair them in streamed batches — the pattern is in Reprojecting Large Datasets Without Memory Errors.

Frequently Asked Questions

Should I clean the Shapefile in place or convert first and clean afterwards? Convert first whenever the downstream consumer will accept it. Every repair written back to .shp has to survive the 10-character field ceiling, the dBASE type system and the single-geometry-type header, so a good half of this guide is working around the output format rather than fixing the data. Reading the dirty Shapefile and writing a clean GeoPackage or GeoParquet deletes those constraints outright — see GeoParquet vs Shapefile for Storage. Keep the Shapefile export only for the specific consumer that insists on one.

How do I know the source CRS when there is no .prj? You do not, and guessing is what causes the tens-of-metres offsets. Ask the supplier, then corroborate: reproject a sample to EPSG:4326 under each candidate CRS and see which one lands the features on the ground they describe. A national grid mistaken for its neighbouring zone puts the data a few hundred kilometres out, which is obvious on a map and invisible in an assertion. Record the answer in a per-supplier configuration file so the next delivery does not restart the investigation.

Is it safe to run make_valid over every feature instead of masking the invalid ones? It is safe but wasteful, and it is not a no-op. make_valid rebuilds the geometry through the noding engine, so already-valid rings can come back with reordered vertices and, on some GEOS versions, a different ring orientation — enough to make a byte-level diff between two runs useless. Masking on ~is_valid keeps untouched features genuinely untouched, which is what makes the audit trail meaningful.

What should happen to features the pipeline cannot repair? Quarantine, never drop. Write the failing rows to a side file with their original index and the exception text, keep the batch running, and count them in the audit record. Silently discarding a handful of features per delivery is how a parcel layer ends up 0.3% short of the register with nobody able to say when it happened.

Does this pipeline scale to a nightly job over hundreds of deliveries? Yes, with two changes. Parallelise across files rather than within them — each delivery is independent, so a ProcessPoolExecutor over the file list is the whole story — and make the run idempotent by writing each output to a temporary name and moving it into place only on success, so an interrupted run leaves no half-written .shp for the next stage to pick up.