How to Reduce GIS File Size in Python Without Wrecking the Data
On this page
Problem statement
The layer works but nobody can use it. A large shapefile will not open in a browser, will not fit in the email, takes minutes to load in QGIS, and costs real money to move between a laptop and a server every day.
The size usually comes from a small number of causes:
- vertex density far beyond the accuracy of the source
- coordinate precision stored to more digits than the survey supports
- 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
Every lever on this page is measured on one real dataset: HM Land Registry's INSPIRE Index Polygons for Manchester — 121,064 freehold parcels that mostly tile the city like a jigsaw, which is exactly the kind of data where careless simplification shows. (Mostly: registered titles overlap, by about 10 km² in total, and that turns out to matter.)
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 — which pays off mainly in text formats
- simplify geometry with a tolerance you can defend, and coverage-aware if the polygons share edges
- 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 = ["INSPIREID", "geometry"]
slim = gdf[keep].copy()
metric = slim if slim.crs.is_projected else slim.to_crs(slim.estimate_utm_crs())
metric["geometry"] = set_precision(metric.geometry.values, grid_size=0.01) # 1 cm
metric.to_parquet("data/out/parcels.parquet", compression="zstd")
metric.to_file("data/out/parcels.gpkg", driver="GPKG")
print(len(metric), "parcels ->", round(Path("data/out/parcels.parquet").stat().st_size / 1e6, 1), "MB of Parquet,",
round(Path("data/out/parcels.gpkg").stat().st_size / 1e6, 1), "MB of GeoPackage")121064 parcels -> 12.9 MB of Parquet, 60.8 MB of GeoPackage
Dropping columns and changing to GeoParquet are lossless for the geometry and did most of the work here; GeoPackage is the better interchange format but was not the smaller file. Simplification is the only step that changes shape, so it comes last and with a measured tolerance.
Where the bytes are
Columns, format, precision, then simplification
Measure where the size actually is
Guessing wastes effort. Split the file into geometry and attributes and look.
import shapely
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")
vertices = shapely.get_num_coordinates(gdf.geometry.values)
print(f"\nfeatures : {len(gdf):,}")
print(f"vertices : {vertices.sum():,} (mean {vertices.mean():.0f}, max {vertices.max():,})")
print(f"columns : {len(gdf.columns)}")
usage = gdf.drop(columns="geometry").memory_usage(deep=True, index=False).sort_values(ascending=False)
print("\nattribute columns in memory (MB):", (usage / 1e6).round(2).to_dict()).cpg 0.0 MB
.dbf 20.1 MB
.prj 0.0 MB
.shp 46.1 MB
.shx 1.0 MB
features : 121,064
vertices : 2,457,901 (mean 20, max 4,745)
columns : 7
attribute columns in memory (MB): {'gml_id': 5.98, 'BEGINLIFES': 4.03, 'VALIDFROM': 3.96, 'INSPIREID': 0.48, 'NATIONALCA': 0.48, 'LABEL': 0.48}
Here the geometry is the bigger part — 46 MB of .shp against 20 MB of .dbf — but the attribute table is still nearly a third of the file, and it is Land Registry's identifiers and dates, most of which a given job never uses. That third costs nothing to drop, and no amount of simplification would touch it.
Drop columns — lossless and free
USED = ["INSPIREID", "geometry"]
print("dropping:", [c for c in gdf.columns if c not in USED])
slim = gdf[USED].copy()
slim.to_file("data/out/parcels_slim.shp")
print(f"{sum(f.stat().st_size for f in Path('data/out').glob('parcels_slim.*')) / 1e6:.1f} MB as a shapefile "
f"with one attribute, down from {sum(f.stat().st_size for f in Path('data/raw').glob('parcels.*')) / 1e6:.1f} MB")dropping: ['gml_id', 'LABEL', 'NATIONALCA', 'VALIDFROM', 'BEGINLIFES']
48.3 MB as a shapefile with one attribute, down from 67.2 MB
Change the format
The container matters more than most people expect.
import shutil
out = Path("data/out")
shutil.rmtree(out / "formats", ignore_errors=True)
sizes = {}
writers = {
"shapefile": lambda p: slim.to_file(p.with_suffix(".shp")),
"geopackage": lambda p: slim.to_file(p.with_suffix(".gpkg"), driver="GPKG"),
"geojson": lambda p: slim.to_file(p.with_suffix(".geojson"), driver="GeoJSON"),
"flatgeobuf": lambda p: slim.to_file(p.with_suffix(".fgb"), driver="FlatGeobuf"),
"geoparquet (snappy)": lambda p: slim.to_parquet(p.with_suffix(".snappy.parquet")),
"geoparquet (zstd)": lambda p: slim.to_parquet(p.with_suffix(".zstd.parquet"), compression="zstd"),
}
for name, write in writers.items():
stem = out / "formats" / "parcels"
stem.parent.mkdir(exist_ok=True)
before = set(stem.parent.iterdir())
write(stem)
sizes[name] = sum(f.stat().st_size for f in set(stem.parent.iterdir()) - before)
base = sizes["shapefile"]
for name, n in sorted(sizes.items(), key=lambda kv: kv[1]):
print(f"{name:20} {n / 1e6:7.1f} MB {n / base:5.0%} of shapefile")geoparquet (zstd) 14.7 MB 31% of shapefile
geoparquet (snappy) 20.4 MB 42% of shapefile
shapefile 48.3 MB 100% of shapefile
flatgeobuf 52.3 MB 108% of shapefile
geopackage 59.7 MB 124% of shapefile
geojson 126.4 MB 262% of shapefile
GeoParquet with zstd was the smallest by a wide margin — under a third of the shapefile — and GeoJSON the largest at over two and a half times. GeoPackage came out larger than the shapefile, as did FlatGeobuf: both build a spatial index into the file, and a GeoPackage stores every geometry as WKB with a header. It remains the best interchange choice, but not a way to save space.
Reduce coordinate precision
These parcels are recorded to the centimetre, and a double stores far more digits than that. Snapping to a coarser grid removes digits — but whether that shrinks the file depends on the format:
from shapely import set_precision
moves = {}
for grid in (0.01, 0.1, 1.0):
snapped = slim.copy()
snapped["geometry"] = set_precision(slim.geometry.values, grid_size=grid)
snapped.to_file(f"data/out/prec_{grid}.gpkg", driver="GPKG")
snapped.to_parquet(f"data/out/prec_{grid}.parquet", compression="zstd")
snapped.to_file(f"data/out/prec_{grid}.geojson", driver="GeoJSON")
moved = moves[grid] = snapped.geometry.hausdorff_distance(slim.geometry)
mb = lambda suffix: Path(f"data/out/prec_{grid}.{suffix}").stat().st_size / 1e6
print(f"grid {grid:>5} m: GeoPackage {mb('gpkg'):5.1f} MB, Parquet {mb('parquet'):5.1f} MB, "
f"GeoJSON {mb('geojson'):6.1f} MB | moved > 2 grid cells: {int((moved > 2 * grid).sum()):>4} | "
f"empty: {int(snapped.is_empty.sum())}")
worst = moves[0.01].idxmax()
at_1cm = set_precision(slim.geometry[worst], grid_size=0.01)
print(f"at 1 cm the largest vertex move was {moves[0.01].max():.0f} m, on a parcel whose area went "
f"{slim.area[worst]:,.1f} -> {at_1cm.area:,.1f} m²")grid 0.01 m: GeoPackage 60.8 MB, Parquet 12.9 MB, GeoJSON 125.7 MB | moved > 2 grid cells: 916 | empty: 0
grid 0.1 m: GeoPackage 60.5 MB, Parquet 11.7 MB, GeoJSON 127.7 MB | moved > 2 grid cells: 992 | empty: 10
grid 1.0 m: GeoPackage 56.2 MB, Parquet 7.3 MB, GeoJSON 67.6 MB | moved > 2 grid cells: 716 | empty: 98
at 1 cm the largest vertex move was 194 m, on a parcel whose area went 1,003,318.5 -> 1,003,318.7 m²
A GeoPackage stores every coordinate as an 8-byte double whatever its value, so snapping barely changed its size until the grid was a metre. The compressed Parquet shrank as the grid coarsened, because repeated low digits compress well, and GeoJSON — text — halved at 1 m. set_precision does more than round, though: it re-nodes each polygon, and at the source's own 1 cm it still changed 916 parcels. The largest "move", nearly 200 m on a 100-hectare parcel whose area did not change, was a zero-width spike in the source outline that snapping removed. At coarser grids small parcels collapsed to empty — 98 of them at 1 m. Pick a grid no coarser than the data's real accuracy, and check what it changed.
For text formats written through GDAL you can also cap the digits directly:
slim.to_crs(4326).to_file("data/out/parcels_4326_p6.geojson", driver="GeoJSON", COORDINATE_PRECISION=6)
slim.to_crs(4326).to_file("data/out/parcels_4326.geojson", driver="GeoJSON")
print({p.name: round(p.stat().st_size / 1e6, 1) for p in sorted(Path("data/out").glob("parcels_4326*.geojson"))}){'parcels_4326.geojson': 122.4, 'parcels_4326_p6.geojson': 78.1}
Six decimal places of latitude is about 0.1 m, still ten times coarser than these parcels' centimetres — so this is a trade, not a free lunch.
Simplify geometry — the only lossy step
before = shapely.get_num_coordinates(slim.geometry.values).sum()
total_area = slim.area.sum()
footprint = shapely.union_all(slim.geometry.values).area
for tol in (0.5, 1, 2, 5):
simplified = slim.geometry.simplify(tol, preserve_topology=True)
v = shapely.get_num_coordinates(simplified.values).sum()
overlap = simplified.area.sum() - shapely.union_all(simplified.values).area
print(f"tol {tol:>3} m -> {v:>9,} vertices ({v / before:5.1%}), area {simplified.area.sum() / total_area - 1:+.3%}, "
f"new overlaps {(overlap - (total_area - footprint)) / 1e4:7.1f} ha")tol 0.5 m -> 885,210 vertices (36.0%), area -0.018%, new overlaps 4.1 ha
tol 1 m -> 813,615 vertices (33.1%), area -0.060%, new overlaps 8.4 ha
tol 2 m -> 752,655 vertices (30.6%), area -0.196%, new overlaps 17.8 ha
tol 5 m -> 662,926 vertices (27.0%), area -1.686%, new overlaps 27.9 ha
Two thirds of the vertices went at half a metre with the total area barely moving, which says how over-digitised the outlines are. But preserve_topology=True stops a polygon collapsing or self-intersecting; it does not keep shared boundaries between neighbours aligned. Parcels simplified independently each get their own version of a shared boundary, and at 2 m that created 17.8 ha of new overlaps. For a coverage, simplify the shared edges once:
tol = 2
independent = slim.geometry.simplify(tol, preserve_topology=True)
coverage = gpd.GeoSeries(shapely.coverage_simplify(slim.geometry.values, tol), crs=slim.crs)
for label, result in (("simplify, per parcel", independent), ("coverage_simplify", coverage)):
overlap = result.area.sum() - shapely.union_all(result.values).area
print(f"{label:22} {shapely.get_num_coordinates(result.values).sum():>9,} vertices, "
f"overlapping area {overlap / 1e4:6.2f} ha (source: {(total_area - footprint) / 1e4:.2f} ha)")simplify, per parcel 752,655 vertices, overlapping area 1018.48 ha (source: 1000.71 ha)
coverage_simplify 1,159,255 vertices, overlapping area 1000.49 ha (source: 1000.71 ha)
shapely.coverage_simplify (Shapely 2.1 with GEOS 3.12) simplifies each shared edge once, so neighbours stayed neighbours: the overlapping area stayed at the source's own 1,000 ha instead of growing by 18 ha. It kept more vertices, because every point where three parcels meet has to stay. The catch is its assumption: it expects a valid coverage, and shapely.coverage_is_valid() said this one is not — the registered titles overlap. On that input the verification below found 140 parcels moved further than the tolerance, the worst by 24 m on a parcel of 7 m². Clean the coverage first, or check what the simplifier did where it overlapped.
Verify before you ship
def size_report(original, reduced, label=""):
print(f"-- {label}")
print(f"features : {len(original):,} -> {len(reduced):,}")
print(f"total area : {original.area.sum() / 1e6:,.3f} -> {reduced.area.sum() / 1e6:,.3f} km²")
print(f"max shift : {reduced.hausdorff_distance(original.geometry).max():.2f} m")
print(f"all valid : {bool(reduced.is_valid.all())} | empty: {int(reduced.is_empty.sum())}")
size_report(slim, coverage, "coverage_simplify @ 2 m")-- coverage_simplify @ 2 m
features : 121,064 -> 121,064
total area : 119.920 -> 119.890 km²
max shift : 24.47 m
all valid : True | empty: 0
Feature count, total area, maximum vertex displacement, validity. If any of those surprises you, the tolerance is too aggressive.
Shrinking Manchester's parcel shapefile
The data used on this page
HM Land Registry's INSPIRE Index Polygons for Manchester City Council (© Crown copyright and database rights; Open Government Licence), converted once from the published GML into a shapefile with all its attributes, as a colleague's download might arrive. Reading the GML with DOWNLOAD_SCHEMA=NO matters: without it, GDAL spent over two minutes trying to fetch an XML schema from a Land Registry server that is not public.
from pathlib import Path
import geopandas as gpd
import requests
for folder in ("data/raw", "data/out"):
Path(folder).mkdir(parents=True, exist_ok=True)
if not Path("data/raw/parcels.shp").exists():
if not Path("data/Manchester.zip").exists():
Path("data/Manchester.zip").write_bytes(requests.Session().get(
"https://use-land-property-data.service.gov.uk/datasets/inspire/download/Manchester_City_Council.zip",
timeout=600).content)
raw = gpd.read_file("/vsizip/data/Manchester.zip/Land_Registry_Cadastral_Parcels.gml", DOWNLOAD_SCHEMA="NO")
raw.to_file("data/raw/parcels.shp")
gdf = gpd.read_file("data/raw/parcels.shp")
print(len(gdf), "parcels,", list(gdf.columns), gdf.crs)121064 parcels, ['gml_id', 'INSPIREID', 'LABEL', 'NATIONALCA', 'VALIDFROM', 'BEGINLIFES', 'geometry'] EPSG:27700
Example 1: A reduction pipeline with a measured report
def reduce_size(src, dest, keep_columns=None, grid_m=0.01, simplify_m=None):
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:
gdf = gdf[[c for c in keep_columns if c in gdf.columns] + [gdf.geometry.name]]
metric = gdf if gdf.crs.is_projected else gdf.to_crs(gdf.estimate_utm_crs())
area_before = metric.area.sum()
geoms = metric.geometry.values
if grid_m:
geoms = shapely.set_precision(geoms, grid_size=grid_m)
if simplify_m:
geoms = shapely.coverage_simplify(geoms, simplify_m)
metric = metric.set_geometry(gpd.GeoSeries(shapely.make_valid(geoms), crs=metric.crs, index=metric.index))
metric = metric[metric.geometry.notna() & ~metric.geometry.is_empty]
out = metric.to_crs(gdf.crs)
if dest.suffix == ".parquet":
out.to_parquet(dest, compression="zstd")
else:
out.to_file(dest, driver="GPKG")
report.update({"out_features": len(out), "out_bytes": dest.stat().st_size,
"ratio": round(dest.stat().st_size / report["in_bytes"], 3),
"area_change_pct": round((metric.area.sum() / area_before - 1) * 100, 4)})
return out, report
_, report = reduce_size("data/raw/parcels.shp", "data/out/parcels_reduced.parquet",
keep_columns=["INSPIREID"], grid_m=0.01, simplify_m=2)
for k, v in report.items():
print(f"{k:16} {v}")in_features 121064
in_columns 7
in_bytes 67199280
out_features 121064
out_bytes 6216584
ratio 0.093
area_change_pct -0.0224
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):
tol = tolerance_for_scale(scale)
v = shapely.get_num_coordinates(shapely.coverage_simplify(slim.geometry.values, tol)).sum()
print(f"1:{scale:>7,} -> tolerance {tol:6.2f} m -> {v / before:6.1%} of the vertices")1: 1,000 -> tolerance 0.25 m -> 74.9% of the vertices
1: 10,000 -> tolerance 2.50 m -> 45.9% of the vertices
1: 50,000 -> tolerance 12.50 m -> 39.5% of the vertices
1:250,000 -> tolerance 62.50 m -> 38.6% of the vertices
Detail finer than a quarter of a millimetre on the final map cannot be seen; removing it costs nothing visible. The returns flatten fast: past 1:10,000 the vertices that remain are the ones where parcels meet, which a coverage simplifier cannot remove — and at 1:250,000 the parcels themselves are smaller than a pixel.
Example 3: Build multiple products from one master
PRODUCTS = {
"analysis": {"tol": None, "fmt": "parquet"},
"desktop": {"tol": 1, "fmt": "gpkg"},
"web": {"tol": 5, "fmt": "fgb"},
}
for name, spec in PRODUCTS.items():
g = slim.copy()
if spec["tol"]:
g["geometry"] = shapely.coverage_simplify(g.geometry.values, spec["tol"])
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:6.1f} MB ({str(spec['tol']) + ' m' if spec['tol'] else 'full detail'})")analysis 14.7 MB (full detail)
desktop 39.0 MB (1 m)
web 29.8 MB (5 m)
Keep the full-fidelity master and derive the small ones. Never simplify in place. Note that the desktop GeoPackage, simplified at 1 m, was still larger than the unsimplified Parquet — the container matters as much as the tolerance.
Example 4: Vector tiles when a single file will not do
tippecanoe -o data/out/parcels.pmtiles -Z 8 -z 14 \
--drop-densest-as-needed --extend-zooms-if-still-dropping \
--force data/out/parcels_4326_p6.geojsonTiling changes the problem rather than the file: the client fetches only the tiles it needs, at a resolution appropriate to the zoom. tippecanoe was not installed where this page ran, so the command above is the documented form rather than a measured one; building vector tiles in Python covers it end to end.
Where the bytes in a vector file go
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.
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 handful of identifiers and dates repeated for every parcel made up nearly a third of the file. GeoParquet is columnar with real compression and was by far the smallest; GeoPackage (SQLite) and FlatGeobuf both carry a spatial index and were larger than the shapefile.
Coordinate precision is nearly lossless and widely misunderstood. A double holds about 15 significant digits whatever you store in it, so in a binary format snapping to a grid changes the values but not the bytes — it pays off in text formats and under compression. 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.
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 it, and overlaps appear between them. For a dataset that forms a clean coverage, coverage_simplify treats shared edges once; for one that only nearly does, like these titles, it needs the coverage repaired first. 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.
Simplification that breaks shared boundaries
- Simplify in a projected CRS: a tolerance in degrees varies with latitude and is meaningless as a distance. Reproject, simplify, reproject back.
preserve_topologyis per feature: it does not preserve shared boundaries between features. Useshapely.coverage_simplify(ortopojson'stoposimplify) for coverages.set_precisioncan produce empty geometries: very thin slivers collapse. Checkis_emptyafter snapping and decide whether that is acceptable.- Shapefile has a 2 GB limit: both the
.shpand.dbfare 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
.dbfdwarfs the.shp, no geometry work will help. Drop columns and use categories.
Internal links
- How to Simplify Geometry in Python with GeoPandas and Shapely
- How to Read and Write GeoPackage Files in Python
- How to Speed Up GeoPandas: Tips for Large Datasets
- Fixing Memory Errors in GeoPandas When Working with Large Files
- How to Batch Export a GeoDataFrame to Multiple Formats in Python
- How to Process a Very Large GeoPackage in Chunks with Python
FAQ
What is the smallest format for vector data?
GeoParquet with zstd compression, by a large margin: 31% of the shapefile for Manchester's parcels, with snappy-compressed Parquet second at 42%. GeoPackage and FlatGeobuf were both larger than the shapefile, because they carry a spatial index; GeoPackage is still 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 — 17.8 ha of new overlaps at 2 m on Manchester's parcels. Use a coverage-aware tool such as shapely.coverage_simplify() or topojson's toposimplify(), and check the input is a valid coverage first.
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.