How to Speed Up GeoPandas: Tips for Large Datasets

On this page

Problem statement

Leicestershire's 114,187 OpenStreetMap buildings show where GeoPandas time goes. Their areas, computed row by row with iterrows(), took 1.51 s; the vectorised .area took 0.005 s, 298 times faster, with the same answer. Written as GeoJSON they took 1.30 s and 60.4 MB, as Parquet 0.14 s and 11.2 MB, and read back 0.94 s against 0.07 s.

Every tip below is timed on that data — the buildings, the county's land use and its wards. The timings are from one machine; the ratios are what carry over.

Quick answer

The fastest ways to speed up GeoPandas on large datasets are:

read only the columns you need
filter rows as early as possible
use faster formats like Parquet or GeoPackage for intermediate outputs
avoid Python loops and apply() when vectorized operations are available
let the spatial index do the work in joins and filters
reproject only when a task actually requires it
break large workflows into smaller saved steps

Reduce data at each step

Shrinking the data at each step of a GeoPandas workflow.

Eight ways to a faster GeoPandas workflow

Every item removes work rather than doing the same work faster.

1. Identify where GeoPandas is slow

Before changing code, measure which step is slow: file reading, filtering, geometry operations, joins, or export:

timings = {}
with timer(timings, "read shapefile from zip"):
    buildings = gpd.read_file(f"{OSM}gis_osm_buildings_a_free_1.shp")
with timer(timings, "filter to houses"):
    houses = buildings[buildings["type"] == "house"]
with timer(timings, "reproject to EPSG:27700"):
    houses = houses.to_crs(27700)
for step, seconds in timings.items():
    print(f"{step:26} {seconds:6.2f} s")
print(f"{len(buildings):,} buildings, {len(houses):,} houses | "
      f"{buildings.memory_usage(deep=True).sum() / 1e6:.0f} MB in memory")

read shapefile from zip      0.31 s
filter to houses             0.00 s
reproject to EPSG:27700      0.14 s
114,187 buildings, 39,780 houses | 8 MB in memory

This helps you avoid optimizing the wrong part of the workflow.

2. Reduce the amount of data loaded into memory

Read only the columns you need

for label, kwargs in (("all columns", {}), ("columns=['osm_id', 'type']", {"columns": ["osm_id", "type"]}),
                      ("read_geometry=False", {"columns": ["type"], "read_geometry": False})):
    t0 = time.perf_counter()
    frame = gpd.read_file("data/buildings.gpkg", **kwargs)
    print(f"{label:28} {time.perf_counter() - t0:5.2f} s, {frame.memory_usage(deep=True).sum() / 1e6:5.1f} MB")

all columns                   0.25 s,   7.6 MB
columns=['osm_id', 'type']    0.22 s,   4.3 MB
read_geometry=False           0.07 s,   1.4 MB

The OSM building layer has only a handful of columns, so dropping them saves little; skipping the geometry entirely, when you only need attributes, saves much more. Wide tables with many text fields gain the most from columns=.

Filter features early — in the read itself

t0 = time.perf_counter()
houses_only = gpd.read_file("data/buildings.gpkg", where="type = 'house'")
where_s = time.perf_counter() - t0
t0 = time.perf_counter()
area = wards.loc[wards["WD24NM"].str.startswith("Loughborough")].to_crs(4326).total_bounds
local = gpd.read_file("data/buildings.gpkg", bbox=tuple(area))
bbox_s = time.perf_counter() - t0
print(f"where=: {len(houses_only):,} houses in {where_s:.2f} s | bbox=: {len(local):,} buildings around Loughborough in {bbox_s:.2f} s")

where=: 39,780 houses in 0.10 s | bbox=: 15,289 buildings around Loughborough in 0.06 s

where= and bbox= are passed to GDAL, so the rows you do not want are never turned into Python objects. A where on a column must include that column in columns= if you also pass one, or GDAL returns no rows.

3. Use faster file formats for intermediate results

File format has a large effect on runtime:

sample = buildings.to_crs(27700)
for suffix, write, read in (("shp", lambda p: sample.to_file(p), gpd.read_file),
                            ("geojson", lambda p: sample.to_file(p), gpd.read_file),
                            ("gpkg", lambda p: sample.to_file(p), gpd.read_file),
                            ("parquet", lambda p: sample.to_parquet(p), gpd.read_parquet)):
    path = Path(f"output/buildings.{suffix}")
    t0 = time.perf_counter(); write(path); write_s = time.perf_counter() - t0
    t0 = time.perf_counter(); read(path); read_s = time.perf_counter() - t0
    files = [path] if suffix != "shp" else [f for f in path.parent.glob("buildings.*") if f.suffix in (".shp", ".shx", ".dbf", ".prj", ".cpg")]
    size = sum(f.stat().st_size for f in files)
    print(f"{suffix:8} write {write_s:5.2f} s  read {read_s:5.2f} s  {size / 1e6:6.1f} MB")

shp      write  0.58 s  read  0.27 s    57.0 MB
geojson  write  1.30 s  read  0.94 s    60.4 MB
gpkg     write  0.38 s  read  0.27 s    31.2 MB
parquet  write  0.14 s  read  0.07 s    11.2 MB

Parquet was the fastest to write and read and among the smallest, and GeoJSON the slowest. Save cleaned intermediate results in Parquet so you do not repeat expensive steps; keep GeoPackage or GeoJSON for delivery to other tools.

4. Improve geometry operation performance

Simplify geometries when acceptable

Complex geometries make overlays, joins, clipping and reprojection slower. The wards are the complex layer here:

wards_bng = wards.to_crs(27700)
simple = wards_bng.assign(geometry=wards_bng.simplify(10, preserve_topology=True))
for label, layer in (("full-resolution wards", wards_bng), ("simplified 10 m", simple)):
    t0 = time.perf_counter()
    rows = len(gpd.sjoin(sample, layer, predicate="within"))
    print(f"{label:22} {shapely.get_num_coordinates(layer.geometry.values).sum():>8,} vertices  "
          f"join {time.perf_counter() - t0:.2f} s, {rows:,} rows")

full-resolution wards   214,358 vertices  join 0.15 s, 113,710 rows
simplified 10 m          23,109 vertices  join 0.09 s, 113,681 rows

Use this carefully. Simplified boundaries can change spatial results near edges — the row counts above differ for exactly that reason.

Check invalid geometries before overlay or joins

invalid = ~buildings.geometry.is_valid
print(int(invalid.sum()), "invalid buildings")
buildings.loc[invalid, "geometry"] = buildings.loc[invalid, "geometry"].make_valid()

0 invalid buildings

Repair with make_valid() rather than buffer(0), which can silently drop part of a self-intersecting polygon.

Reproject only when needed

t0 = time.perf_counter()
_ = buildings.to_crs(27700)
print(f"one reprojection of {len(buildings):,} buildings: {time.perf_counter() - t0:.2f} s")
if sample.crs.to_epsg() != 27700:
    sample = sample.to_crs(27700)          # the check costs nothing; the call does

one reprojection of 114,187 buildings: 0.25 s

5. Speed up spatial joins and spatial filters

Let the spatial index do it

fresh = sample.copy()
print("index before first use:", fresh.has_sindex)
t0 = time.perf_counter(); _ = fresh.sindex; build_s = time.perf_counter() - t0
print(f"built in {build_s:.3f} s, now: {fresh.has_sindex}")

index before first use: False
built in 0.038 s, now: True

GeoPandas builds the index the first time a join or query needs it, and caches it on the frame. The expensive mistake is not a missing index but a loop that never uses one — step 6.

Limit candidate features before joining

minx, miny, maxx, maxy = wards_bng[wards_bng["WD24NM"].str.startswith("Loughborough")].total_bounds
candidates = sample.cx[minx:maxx, miny:maxy]
t0 = time.perf_counter()
joined = gpd.sjoin(candidates, wards_bng, predicate="within")
print(f".cx kept {len(candidates):,} of {len(sample):,}; join {time.perf_counter() - t0:.2f} s")

.cx kept 15,319 of 114,187; join 0.02 s

.cx is a bounding-box prefilter, not an exact spatial match, so the real join still runs after it — on a fraction of the rows.

6. Replace slow row-by-row patterns

t0 = time.perf_counter()
areas = [row.geometry.area for _, row in sample.iterrows()]
loop_s = time.perf_counter() - t0
t0 = time.perf_counter()
vector = sample.geometry.area
vector_s = time.perf_counter() - t0
print(f"iterrows: {loop_s:.2f} s | .area: {vector_s:.4f} s | {loop_s / vector_s:,.0f}× faster, same answer: "
      f"{bool((abs(vector.values - areas) < 1e-6).all())}")

iterrows: 1.51 s | .area: 0.0051 s | 298× faster, same answer: True

For attribute logic, use pandas vectorized operations instead of apply(): sample["is_large"] = sample.area > 1000.

7. Manage memory in long workflows

Drop unused columns and temporary outputs as soon as possible, and if the whole workflow does not fit in memory, process pieces and save each:

for number, start in enumerate(range(0, len(sample), 50_000), 1):
    sample.iloc[start:start + 50_000].to_parquet(f"output/buildings_part{number}.parquet")
parts = sorted(Path("output").glob("buildings_part*.parquet"))
print(len(parts), "parts,", sum(len(gpd.read_parquet(p)) for p in parts), "rows in total")

3 parts, 114187 rows in total

GeoPandas is an in-memory tool. For very large pipelines, writing intermediate files is often better than keeping everything in one session.

8. Know when GeoPandas is not the right tool for the full job

GeoPandas is strong for analysis, QA, preprocessing and export. For millions of features or repeated heavy joins, move the largest filtering and join steps into PostGIS or DuckDB with its spatial extension — see running spatial SQL from Python — and use GeoPandas for the smaller final steps.

Leicestershire buildings and land use, processed fast

The data used on this page

Leicestershire's buildings and land use from Geofabrik's OpenStreetMap extract of 1 January 2026 (© OpenStreetMap contributors, ODbL), and the wards of Leicester and the county's districts from the ONS December 2024 boundaries (OGL). timer() records how long a block takes.

import time
from contextlib import contextmanager
from pathlib import Path

import geopandas as gpd
import requests
import shapely

Path("data").mkdir(exist_ok=True)
Path("output").mkdir(exist_ok=True)
if not Path("data/leicestershire.shp.zip").exists():
    Path("data/leicestershire.shp.zip").write_bytes(requests.get(
        "https://download.geofabrik.de/europe/united-kingdom/england/leicestershire-260101-free.shp.zip", timeout=600).content)
LADS = ["E06000016", "E07000129", "E07000130", "E07000131", "E07000132", "E07000133", "E07000134", "E07000135"]
if not Path("data/wards.geojson").exists():
    Path("data/wards.geojson").write_bytes(requests.get(
        "https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
        "Wards_December_2024_Boundaries_UK_BFC/FeatureServer/0/query",
        params={"where": "LAD24CD IN ('" + "','".join(LADS) + "')", "outFields": "WD24NM,LAD24NM",
                "outSR": 27700, "f": "geojson"}, timeout=180).content)
OSM = "/vsizip/data/leicestershire.shp.zip/"
wards = gpd.read_file("data/wards.geojson")
if not Path("data/buildings.gpkg").exists():
    gpd.read_file(OSM + "gis_osm_buildings_a_free_1.shp").to_file("data/buildings.gpkg")


@contextmanager
def timer(store, label):
    start = time.perf_counter()
    yield
    store[label] = time.perf_counter() - start

print(len(wards), "wards;", "buildings.gpkg", round(Path("data/buildings.gpkg").stat().st_size / 1e6), "MB")

161 wards; buildings.gpkg 31 MB

Example 1: the whole workflow, fast

t0 = time.perf_counter()
houses = gpd.read_file("data/buildings.gpkg", columns=["osm_id", "type"], where="type = 'house'")   # read less
houses = houses.to_crs(27700)                                                                       # reproject once
houses = houses[~houses.is_empty]
result = gpd.sjoin(houses, wards[["WD24NM", "geometry"]], predicate="within", how="inner")           # indexed join
per_ward = result.groupby("WD24NM").size().sort_values(ascending=False)                             # vectorised
result[["osm_id", "WD24NM", "geometry"]].to_parquet("output/houses_by_ward.parquet")                 # fast format
print(f"{len(houses):,} houses -> {len(per_ward)} wards in {time.perf_counter() - t0:.2f} s; top: {per_ward.head(2).to_dict()}")

39,780 houses -> 147 wards in 0.26 s; top: {'Hinckley De Montfort': 2649, 'Burbage Sketchley and Stretton': 2476}

Example 2: the relative cost of common operations

landuse = gpd.read_file(OSM + "gis_osm_landuse_a_free_1.shp", where="fclass = 'residential'").to_crs(27700)
costs = {}
with timer(costs, "attribute filter"):
    _ = sample[sample["type"] == "house"]
with timer(costs, "buffer 10 m"):
    _ = sample.buffer(10)
with timer(costs, "reproject"):
    _ = sample.to_crs(4326)
with timer(costs, "spatial join"):
    _ = gpd.sjoin(sample, wards_bng, predicate="within")
with timer(costs, "overlay with land use"):
    _ = gpd.overlay(sample.head(20_000), landuse, how="intersection")
for operation, seconds in sorted(costs.items(), key=lambda item: item[1]):
    print(f"{operation:22} {seconds:7.3f} s")

attribute filter         0.003 s
spatial join             0.134 s
overlay with land use    0.249 s
reproject                0.315 s
buffer 10 m              1.134 s

Buffering every building was the slowest step and reprojecting them the second; the overlay, run on a sixth of the buildings, came third, and the indexed join and the attribute filter barely registered. Order a pipeline cheapest-first — filter, then join, then the geometry work — and the expensive steps get smaller inputs.

Why filtering first speeds up every later step

Order your pipeline cheapest-first and the expensive steps get smaller inputs.

Early filtering has a large effect because every later step runs on fewer rows: reducing a layer to the houses, or to one town's bounding box, before a join also reduces index building, geometry comparisons and output size.

Disk speed, geometry complexity and Parquet's limits

Disk speed matters: SSD vs network storage can change file read and write performance significantly.
Geometry complexity matters: 10,000 very complex polygons can be slower than 100,000 simple points.
Parquet is not a universal interchange format: it is most useful when your workflow stays mainly in Python or other modern data tools.

Internal links

FAQ

Why is GeoPandas so slow with large datasets?

Usually because too much data is loaded at once, geometries are complex, file formats are inefficient, or the workflow uses loops instead of vectorized operations.

Does GeoPandas use a spatial index automatically?

GeoPandas uses spatial indexing in many operations such as spatial joins, but behavior depends on your environment and installed backends. It is still useful to check has_sindex and reduce candidate features before the operation.

Which file format is faster than Shapefile for GeoPandas workflows?

GeoPackage is usually better than Shapefile, and Parquet is often faster for repeated Python-based analysis workflows.

Should I use `apply()` for GeoPandas performance-sensitive tasks?

Usually no. Prefer vectorized GeoPandas geometry methods and pandas column operations when possible.

When should I switch from GeoPandas to PostGIS or another database?

If your workflow involves millions of features, repeated heavy joins, or memory failures, a spatial database is often a better place for the largest filtering and join steps.

Can I read only a spatial subset of a file instead of loading it all?

Yes. gpd.read_file() accepts a bbox= tuple or a mask= geometry to read only features in that extent, and reading from GeoParquet supports similar pushdown filtering. This avoids loading the whole layer into memory before filtering.

Does GeoPandas use multiple CPU cores for spatial operations?

Core GeoPandas geometry operations are mostly single-threaded. For parallelism across cores or larger-than-memory data, look at dask-geopandas, which partitions a GeoDataFrame and runs operations like spatial joins in parallel.

Will dropping the geometry column speed up pure attribute work?

Yes. If a step only touches attributes, convert to a plain DataFrame with pd.DataFrame(gdf.drop(columns="geometry")) so pandas does not carry the geometry array along. Reattach geometry afterward only when you need spatial operations again.

How can I make a heavy spatial join faster than `sjoin` alone?

Filter both layers to a shared bounding box with .cx[] first, ensure the right-hand layer has a built spatial index via .sindex, and choose the most specific predicate= (such as within instead of intersects). For truly large joins, push the join into DuckDB or PostGIS and use GeoPandas for the final steps.