How to Reduce GIS File Size in Python Without Wrecking the Data

Problem statement

The layer works but nobody can use it. A 1.4 GB shapefile will not open in a browser, will not fit in the email, takes four minutes to load in QGIS, and costs real money to move between a laptop and a server every day.

>>> import geopandas as gpd
>>> gdf = gpd.read_file("data/raw/parcels.shp")
>>> len(gdf), len(gdf.columns)
(842_119, 47)
>>> gdf.geometry.apply(lambda g: len(g.exterior.coords) if g.geom_type == "Polygon" else 0).mean()
612.4

Six hundred vertices per parcel, 47 columns of which analysis uses six, and a format that stores everything as fixed-width text. Every one of those is a lever, and each has a different cost in fidelity.

The size usually comes from a small number of causes:

vertex density far beyond the accuracy of the source
coordinate precision stored to 15 decimal places when the survey is accurate to a centimetre
columns nobody uses, carried from the original download
a format with no compression (shapefile) or a verbose one (GeoJSON)
duplicate or overlapping features, and multipart bundles exploded earlier in the pipeline
attributes stored as wide text where a code or category would do

Quick answer

Attack size in this order — cheapest fidelity cost first:

drop columns and rows you do not need
change format: GeoPackage or GeoParquet instead of shapefile or GeoJSON
reduce coordinate precision to match the data's real accuracy
simplify geometry with a tolerance you can defend
only then consider tiling or splitting the dataset

import geopandas as gpd
from shapely import set_precision

gdf = gpd.read_file("data/raw/parcels.shp")

keep = ["parcel_id", "class", "area_m2", "updated", "geometry"]
slim = gdf[keep].copy()

# 1 cm grid in a metric CRS — beyond survey accuracy anyway
metric = slim.to_crs(slim.estimate_utm_crs())
metric["geometry"] = set_precision(metric.geometry.values, grid_size=0.01)
slim = metric.to_crs(slim.crs)

slim.to_parquet("data/out/parcels.parquet", compression="zstd")
slim.to_file("data/out/parcels.gpkg", driver="GPKG")

Dropping columns and changing format are lossless for the geometry and typically account for most of the reduction. Simplification is the only step that changes shape, so it comes last and with a measured tolerance.

Where the bytes are

Bar chart of file size after each reduction step from raw shapefile to simplified Parquet. — Format and columns do most of the work before a single vertex is touched.

Step-by-step solution

Vertical steps: measure, drop columns, change format, reduce precision, simplify, verify. — Lossless steps first; the lossy one last, and measured.

Measure where the size actually is

Guessing wastes effort. Split the file into geometry and attributes and look.

from pathlib import Path
import geopandas as gpd
import numpy as np

path = Path("data/raw/parcels.shp")
for f in sorted(path.parent.glob(f"{path.stem}.*")):
    print(f"{f.suffix:6} {f.stat().st_size / 1e6:8.1f} MB")

gdf = gpd.read_file(path)
coords = gdf.geometry.apply(lambda g: len(g.exterior.coords) if g.geom_type == "Polygon" else 0)
print(f"\nfeatures     : {len(gdf):,}")
print(f"vertices     : {coords.sum():,} (mean {coords.mean():.0f}, max {coords.max():,})")
print(f"columns      : {len(gdf.columns)}")

usage = gdf.drop(columns="geometry").memory_usage(deep=True).sort_values(ascending=False)
print("\nheaviest attribute columns (MB):")
print((usage / 1e6).head(8).round(2))

A .dbf bigger than the .shp means the attributes are the problem, and no amount of simplification will help.

Drop columns and rows — lossless and free

USED = ["parcel_id", "class", "area_m2", "updated", "geometry"]
print("dropping:", [c for c in gdf.columns if c not in USED])
slim = gdf[USED].copy()

# rows too, if the deliverable has a scope
slim = slim[slim["class"].isin(["residential", "commercial"])]

Then shrink what remains. Categorical dtypes and narrower numerics cost nothing in fidelity:

slim["class"] = slim["class"].astype("category")
slim["area_m2"] = slim["area_m2"].astype("float32")
slim["parcel_id"] = slim["parcel_id"].astype("string")
print(slim.memory_usage(deep=True).sum() / 1e6, "MB in memory")

Change the format

The container matters more than most people expect.

from pathlib import Path
import geopandas as gpd

out = Path("data/out"); out.mkdir(parents=True, exist_ok=True)
sizes = {}

slim.to_file(out / "parcels.shp", driver="ESRI Shapefile")
sizes["shapefile"] = sum(f.stat().st_size for f in out.glob("parcels.*")
                         if f.suffix in {".shp", ".shx", ".dbf", ".prj", ".cpg"})

slim.to_file(out / "parcels.gpkg", driver="GPKG")
sizes["geopackage"] = (out / "parcels.gpkg").stat().st_size

slim.to_file(out / "parcels.geojson", driver="GeoJSON")
sizes["geojson"] = (out / "parcels.geojson").stat().st_size

slim.to_file(out / "parcels.fgb", driver="FlatGeobuf")
sizes["flatgeobuf"] = (out / "parcels.fgb").stat().st_size

slim.to_parquet(out / "parcels.parquet", compression="zstd")
sizes["geoparquet"] = (out / "parcels.parquet").stat().st_size

base = sizes["shapefile"]
for name, n in sorted(sizes.items(), key=lambda kv: kv[1]):
    print(f"{name:12} {n/1e6:8.1f} MB   {n/base:5.0%} of shapefile")

GeoParquet with zstd is usually the smallest by a wide margin, GeoPackage is the best interchange choice, and GeoJSON is the largest — useful to know before emailing one.

Reduce coordinate precision

Fifteen decimal places of longitude is a picometre. Storing it is pure waste, and it also makes geometry operations slower.

from shapely import set_precision
import geopandas as gpd

metric = slim.to_crs(slim.estimate_utm_crs())

for grid in (0.001, 0.01, 0.1, 1.0):
    snapped = metric.copy()
    snapped["geometry"] = set_precision(metric.geometry.values, grid_size=grid)
    snapped.to_file(f"/tmp/prec_{grid}.gpkg", driver="GPKG")
    size = Path(f"/tmp/prec_{grid}.gpkg").stat().st_size / 1e6
    moved = snapped.geometry.hausdorff_distance(metric.geometry).max()
    print(f"grid {grid:>6} m → {size:6.1f} MB, max vertex shift {moved:.3f} m")

set_precision snaps coordinates to a grid and re-nodes the result, so it also removes the micro-self-intersections that break overlays. Pick a grid coarser than your data's real accuracy and no coarser.

For formats written through GDAL you can also cap the digits directly:

slim.to_file("data/out/parcels.geojson", driver="GeoJSON",
             COORDINATE_PRECISION=6)          # ~0.1 m in degrees

Simplify geometry — the only lossy step

import geopandas as gpd

metric = slim.to_crs(slim.estimate_utm_crs())
before_vertices = metric.geometry.apply(lambda g: len(g.exterior.coords)
                                        if g.geom_type == "Polygon" else 0).sum()

for tol in (0.5, 1, 2, 5, 10):
    simplified = metric.geometry.simplify(tol, preserve_topology=True)
    v = simplified.apply(lambda g: len(g.exterior.coords) if g.geom_type == "Polygon" else 0).sum()
    area_change = (simplified.area.sum() - metric.geometry.area.sum()) / metric.geometry.area.sum()
    print(f"tol {tol:>5} m → {v:>10,} vertices ({v/before_vertices:5.1%}), "
          f"area change {area_change:+.3%}")

preserve_topology=True stops a polygon collapsing or self-intersecting, but it does not keep shared boundaries between neighbours aligned — two adjacent parcels simplified independently will develop slivers and gaps between them. For a coverage, use a topology-aware simplifier:

import topojson as tp

topo = tp.Topology(slim.to_crs(slim.estimate_utm_crs()), prequantize=False)
simplified = topo.toposimplify(5).to_gdf()      # shared edges stay shared
simplified.crs = slim.estimate_utm_crs()

Verify before you ship

def size_report(original, reduced, label=""):
    om = original.to_crs(original.estimate_utm_crs())
    rm = reduced.to_crs(original.estimate_utm_crs())
    print(f"── {label}")
    print(f"features   : {len(original):,} → {len(reduced):,}")
    print(f"total area : {om.geometry.area.sum()/1e6:,.2f} → {rm.geometry.area.sum()/1e6:,.2f} km²")
    print(f"area change: {(rm.geometry.area.sum()/om.geometry.area.sum() - 1):+.4%}")
    print(f"max shift  : {rm.geometry.hausdorff_distance(om.geometry).max():.2f} m")
    print(f"all valid  : {rm.geometry.is_valid.all()}")

size_report(gdf, simplified, "after simplify @ 5 m")

Feature count, total area, maximum vertex displacement, validity. If any of those surprises you, the tolerance is too aggressive.

Code examples

Example 1: a reduction pipeline with a measured report

from pathlib import Path
import geopandas as gpd
from shapely import set_precision

def reduce_size(src, dest, keep_columns=None, grid_m=0.01, simplify_m=None,
                driver="GPKG", compression="zstd"):
    src, dest = Path(src), Path(dest)
    gdf = gpd.read_file(src)
    report = {"in_features": len(gdf), "in_columns": len(gdf.columns),
              "in_bytes": sum(f.stat().st_size for f in src.parent.glob(f"{src.stem}.*"))}

    if keep_columns:
        cols = [c for c in keep_columns if c in gdf.columns]
        if gdf.geometry.name not in cols:
            cols.append(gdf.geometry.name)
        gdf = gdf[cols]

    metric_crs = gdf.estimate_utm_crs()
    metric = gdf.to_crs(metric_crs)
    area_before = metric.geometry.area.sum()

    if grid_m:
        metric["geometry"] = set_precision(metric.geometry.values, grid_size=grid_m)
    if simplify_m:
        metric["geometry"] = metric.geometry.simplify(simplify_m, preserve_topology=True)

    metric["geometry"] = metric.geometry.make_valid()
    metric = metric[metric.geometry.notna() & ~metric.geometry.is_empty]

    out = metric.to_crs(gdf.crs)
    dest.parent.mkdir(parents=True, exist_ok=True)
    if dest.suffix == ".parquet":
        out.to_parquet(dest, compression=compression)
    else:
        out.to_file(dest, driver=driver)

    report.update({
        "out_features": len(out),
        "out_columns": len(out.columns),
        "out_bytes": dest.stat().st_size,
        "ratio": round(dest.stat().st_size / report["in_bytes"], 3),
        "area_change_pct": round((metric.geometry.area.sum() / area_before - 1) * 100, 4),
    })
    return out, report

_, report = reduce_size("data/raw/parcels.shp", "data/out/parcels.parquet",
                        keep_columns=["parcel_id", "class", "area_m2"],
                        grid_m=0.01, simplify_m=2)
for k, v in report.items():
    print(f"{k:16} {v}")

Example 2: pick a simplification tolerance from the map scale

Tolerance is not arbitrary — it follows from how the data will be displayed.

def tolerance_for_scale(scale_denominator: int, pixel_mm: float = 0.25) -> float:
    """Ground distance, in metres, of one screen-detectable unit at a given scale."""
    return scale_denominator * (pixel_mm / 1000)

for scale in (1_000, 10_000, 50_000, 250_000, 1_000_000):
    print(f"1:{scale:>9,}  →  simplify tolerance ≈ {tolerance_for_scale(scale):8.1f} m")

Detail finer than a quarter of a millimetre on the final map cannot be seen; removing it costs nothing visible and can remove most of the vertices.

Example 3: build multiple products from one master

from pathlib import Path
import geopandas as gpd

master = gpd.read_file("data/clean/parcels.gpkg")
metric = master.to_crs(master.estimate_utm_crs())

PRODUCTS = {
    "analysis":  {"tol": None, "cols": None,                          "fmt": "parquet"},
    "desktop":   {"tol": 1,    "cols": ["parcel_id", "class"],        "fmt": "gpkg"},
    "web":       {"tol": 10,   "cols": ["parcel_id", "class"],        "fmt": "fgb"},
}

for name, spec in PRODUCTS.items():
    g = metric.copy()
    if spec["tol"]:
        g["geometry"] = g.geometry.simplify(spec["tol"], preserve_topology=True)
    if spec["cols"]:
        g = g[spec["cols"] + ["geometry"]]
    g = g.to_crs(master.crs)

    dest = Path(f"data/out/parcels_{name}.{spec['fmt']}")
    if spec["fmt"] == "parquet":
        g.to_parquet(dest, compression="zstd")
    else:
        g.to_file(dest, driver={"gpkg": "GPKG", "fgb": "FlatGeobuf"}[spec["fmt"]])
    print(f"{name:9} {dest.stat().st_size/1e6:7.1f} MB")

Keep the full-fidelity master and derive the small ones. Never simplify in place.

Example 4: vector tiles when a single file will not do

import subprocess
from pathlib import Path

gdf.to_file("data/out/parcels.geojson", driver="GeoJSON", COORDINATE_PRECISION=6)

subprocess.run([
    "tippecanoe", "-o", "data/out/parcels.pmtiles",
    "-Z", "8", "-z", "14",
    "--drop-densest-as-needed", "--extend-zooms-if-still-dropping",
    "--force", "data/out/parcels.geojson",
], check=True)

print(Path("data/out/parcels.pmtiles").stat().st_size / 1e6, "MB")

Tiling changes the problem rather than the file: the client fetches only the tiles it needs, at a resolution appropriate to the zoom.

Explanation

A vector file's size is roughly the number of vertices times the bytes per coordinate, plus the attribute table, plus whatever overhead the format imposes. Each of those three terms has its own lever, and they differ enormously in what they cost you.

Grid comparing shapefile, GeoPackage, GeoJSON, FlatGeobuf and GeoParquet on size, compression and use. — The same data in five containers — with what each one is actually for.

Dropping columns and switching format are lossless: the geometry is untouched and the remaining attributes are identical. They are also where most of the easy win lives. A shapefile stores attributes in a dBase table with fixed-width fields and no compression, so a 47-column table of mostly-empty text can dwarf the geometry. GeoPackage is SQLite, GeoParquet is columnar with real compression, and FlatGeobuf is a compact binary layout with a spatial index — all of them dramatically smaller for the same content.

Coordinate precision is nearly lossless and widely misunderstood. A double holds about 15 significant digits, so a longitude written in full implies sub-atomic precision. Snapping to a grid coarser than the data's true accuracy discards only noise, and set_precision re-nodes the geometry while it does so, which incidentally clears up the micro-intersections that cause TopologyException in overlays.

Simplification is the only genuinely lossy step, and it is where care is required. Douglas–Peucker with preserve_topology=True keeps each geometry valid, but it processes features independently — so two parcels sharing a boundary get two different simplified versions of that boundary, and slivers appear between them. For any dataset that forms a coverage, a topology-preserving simplifier that treats shared edges once is the right tool. Choosing the tolerance from the intended display scale, rather than by trial and error, turns an aesthetic judgement into a defensible one.

Finally, keep the master. Every reduction is a product built for an audience — analysis, desktop, web — and the moment you simplify in place you have lost the ability to make a different product later.

Edge cases or notes

Simplify in a projected CRS: A tolerance in degrees varies with latitude and is meaningless as a distance. Reproject, simplify, reproject back.
preserve_topology is per feature: It does not preserve shared boundaries between features. Use topojson/toposimplify for coverages.
set_precision can produce empty geometries: Very thin slivers collapse. Check is_empty after snapping and decide whether that is acceptable.
Shapefile has a 2 GB limit: Both the .shp and .dbf are capped. Hitting it is a sign to change format, not to simplify harder.
GeoJSON is text: It compresses well over HTTP with gzip or brotli, so the on-disk size overstates the transfer cost.
Parquet needs a reader: GeoParquet is superb for pipelines and poor for handing to someone with an old desktop GIS. Match the format to the recipient.
Attributes may be the whole problem: If the .dbf dwarfs the .shp, no geometry work will help. Drop columns and use categories.

Internal links

FAQ

What is the smallest format for vector data?

GeoParquet with zstd compression, usually by a large margin, followed by FlatGeobuf. GeoPackage is the best compromise when the recipient needs to open it in a desktop GIS.

How much can I simplify without damaging the data?

Derive the tolerance from the display scale: roughly scale_denominator × 0.00025 metres. Detail finer than that is invisible on the final map, so removing it costs nothing a viewer can see.

Why do gaps appear between polygons after simplifying?

Because simplify() processes each geometry independently, so a shared boundary is simplified twice, differently. Use a topology-aware tool such as topojson.Topology(...).toposimplify() for coverages.

Does reducing coordinate precision lose real information?

Only if you snap coarser than the data's true accuracy. Survey data accurate to a centimetre loses nothing on a 1 cm grid, and you shed the meaningless digits that made every coordinate 15 characters long.

My shapefile's `.dbf` is larger than the `.shp` — what do I do?

The attributes are the problem. Drop unused columns, convert repeated strings to categories, and move to GeoPackage or Parquet, which compress attribute data properly.

Should I simplify before or after reprojecting?

Simplify in a projected CRS with metre units, then reproject the result if you need a different output CRS. A tolerance expressed in degrees is not a distance.

When should I switch to vector tiles?

When the map is interactive and the dataset is too large to send in one piece at any tolerance. Tools like tippecanoe build zoom-dependent tiles so the client only downloads what is in view.

How to Reduce GIS File Size in Python Without Wrecking the Data #

Problem statement #

Quick answer #

Where the bytes are #

Step-by-step solution #

Measure where the size actually is #

Drop columns and rows — lossless and free #

Change the format #

Reduce coordinate precision #

Simplify geometry — the only lossy step #

Verify before you ship #

Code examples #

Example 1: a reduction pipeline with a measured report #

Example 2: pick a simplification tolerance from the map scale #

Example 3: build multiple products from one master #

Example 4: vector tiles when a single file will not do #

Explanation #

Edge cases or notes #

Internal links #

FAQ #

What is the smallest format for vector data? #

How much can I simplify without damaging the data? #

Why do gaps appear between polygons after simplifying? #

Does reducing coordinate precision lose real information? #

My shapefile's .dbf is larger than the .shp — what do I do? #

Should I simplify before or after reprojecting? #

When should I switch to vector tiles? #