GeoPandas .apply() Takes Hours: How to Fix It

Problem statement

The cell has been running for twenty minutes and the progress bar has barely moved.

def find_zone(row):
    hit = zones[zones.contains(row.geometry)]
    return hit["zone_name"].iloc[0] if len(hit) else None

gdf["zone"] = gdf.apply(find_zone, axis=1)

Measured on 200,000 points and 400 zones, that loop takes 59 s; on four million points it would take about twenty minutes, and with a few thousand zones, hours. This does the same thing on the 200,000 points in 0.08 s:

gdf = gpd.sjoin(gdf, zones[["zone_name", "geometry"]], how="left", predicate="within")

That is a 744-fold gap, and very little of it is apply itself. On pandas 3.0 a row-wise apply costs about 4.6 µs a row: gdf.apply(lambda r: r.geometry.area, axis=1) took 0.92 s on 200,000 rows, so four million would take under 20 s. The hours come from what the function does on each row — here, a scan of every zone with no spatial index. The vectorised gdf.geometry.area is still 230× faster than that apply, which is worth having; but a per-row search is the apply that takes hours.

Telling the two costs apart — the row object, and the work done inside it — tells you which apply calls to fix first and what to replace them with.

Quick answer

Measured with tools/bench/apply_is_slow.py in the site's repository — best of three on an AMD Ryzen 9 7900, Python 3.14, pandas 3.0.6, GeoPandas 1.1.4, Shapely 2.1.2. Other versions and machines give other numbers; rerun the script for yours. Every row-by-row form is within a factor of ten of the others; the vectorised call is in a different league.

# 1. one geometry method  → the GeoSeries method
gdf["area"] = gdf.geometry.area                       # not apply(axis=1)

# 2. column arithmetic     → do it on the columns
gdf["density"] = gdf["pop"] / gdf["area"]             # not apply(axis=1)

# 3. an if/else            → np.select
gdf["cls"] = np.select([a > 1e4, a > 1e3], ["large", "medium"], default="small")

# 4. a dict lookup         → .map
gdf["name"] = gdf["code"].map(LOOKUP)

# 5. searching another layer → sjoin
gdf = gpd.sjoin(gdf, zones, predicate="within")

# 6. genuinely custom      → apply on the GeoSeries, not the frame
gdf["odd"] = gdf.geometry.apply(my_custom_function)

The rule: apply(axis=1) is only ever needed when the function uses several columns at once and has no vectorised form — which is rarer than it looks.

Why axis=1 is worse than a loop, and its replacements

1. Understand why axis=1 is so much worse than a loop

import time
import geopandas as gpd
from shapely.geometry import box

def timed(label, fn):
    t0 = time.perf_counter()
    fn()
    print(f"{label:<38} {time.perf_counter() - t0:>8.3f} s")

narrow = gpd.GeoDataFrame({"id": range(200_000)},
                          geometry=[box(x, 0, x+1, 1) for x in range(200_000)],
                          crs=27700)
wide = narrow.copy()
for i in range(40):
    wide[f"col_{i}"] = 1.0

timed("narrow: geometry.area",  lambda: narrow.geometry.area)
timed("narrow: apply(axis=1)",  lambda: narrow.apply(lambda r: r.geometry.area, axis=1))
timed("wide (41 cols): geometry.area", lambda: wide.geometry.area)
timed("wide (41 cols): apply(axis=1)", lambda: wide.apply(lambda r: r.geometry.area, axis=1))

narrow: geometry.area                     0.010 s
narrow: apply(axis=1)                     0.948 s
wide (41 cols): geometry.area             0.007 s
wide (41 cols): apply(axis=1)             1.197 s

The vectorised call is unaffected by column count. apply(axis=1) gets slower as columns are added — 26% slower with 41 columns than with 2 — and that is the tell: it is constructing a pandas Series for every row, copying every column's value into it, before calling your function. Your function then uses one of those columns.

On pandas 3.0 that construction is cheap enough — about 4.6 µs a row — that apply alone costs about a second per 200,000 rows. That is 230× the vectorised call, and worth removing, but on its own it rarely explains a job that takes hours. When one does, look at what the function does per row (step 3).

2. Find the apply calls that matter

import re
from pathlib import Path

def find_applies(root="src", show=None):
    hits = []
    for path in Path(root).rglob("*.py"):
        if "tests" in path.parts:
            continue
        for n, line in enumerate(path.read_text().splitlines(), 1):
            if re.search(r"\.apply\s*\(", line) or "iterrows" in line:
                kind = ("iterrows" if "iterrows" in line
                        else "axis=1" if "axis=1" in line else "apply")
                hits.append((kind, f"{path.relative_to(Path(root).parent)}:{n}", line.strip()[:70]))
    hits.sort(key=lambda h: {"iterrows": 0, "axis=1": 1}.get(h[0], 2))
    counts = {k: sum(h[0] == k for h in hits) for k in ("iterrows", "axis=1", "apply")}
    print(f"  {len(hits)} hits: {counts}")
    for kind, where, code in hits[:show]:
        mark = {"iterrows": "✗✗", "axis=1": "✗ ", "apply": "· "}[kind]
        print(f"  {mark} {where:<34} {code}")
    return hits


import momepy                                      # a real codebase to search: momepy's own source

hits = find_applies(Path(momepy.__file__).parent, show=8)

  41 hits: {'iterrows': 6, 'axis=1': 1, 'apply': 34}
  ✗✗ momepy/continuity.py:123           for _, row in stroke_gdf.iterrows()
  ✗✗ momepy/streetscape.py:1024         for uid, row in self._sightline_indicators.iterrows():
  ✗✗ momepy/streetscape.py:1065         ).iterrows():
  ✗✗ momepy/streetscape.py:1413         for street_uid, row in self._sightline_indicators.iterrows():
  ✗✗ momepy/streetscape.py:1927         for street_uid, row in self._sightline_indicators.iterrows():
  ✗✗ momepy/preprocessing.py:938        for _i, row in incoming.iterrows():
  ✗  momepy/preprocessing.py:99         blg["n_count"] = blg.apply(lambda row: len(row.neighbors), axis=1)
  ·  momepy/intensity.py:75             .apply(_calculate_courtyards)

Pointed at momepy 1.0's own source, it found 41 calls: six iterrows loops, one apply(axis=1) — lambda row: len(row.neighbors), which blg.neighbors.str.len() would do in one call — and 34 plain .applys on a Series or group. The two-tier marking matters. iterrows and axis=1 are almost always fixable; a plain .apply on a GeoSeries with a genuinely custom function was 3.6× faster than the frame-level apply and is often fine to leave alone.

3. Replace the common patterns

One geometry method:

# ❌
gdf["area"] = gdf.apply(lambda r: r.geometry.area, axis=1)
gdf["centroid"] = gdf.apply(lambda r: r.geometry.centroid, axis=1)
gdf["valid"] = gdf.apply(lambda r: r.geometry.is_valid, axis=1)

# ✅
gdf["area"] = gdf.geometry.area
gdf["centroid"] = gdf.geometry.centroid
gdf["valid"] = gdf.geometry.is_valid

Arithmetic across columns — this is the case that looks like it needs axis=1 and does not:

# ❌
gdf["density"] = gdf.apply(lambda r: r["pop"] / r["area_km2"], axis=1)

# ✅ pandas aligns the columns for you
gdf["density"] = gdf["pop"] / gdf["area_km2"]

# ✅ and division by zero handled honestly
gdf["density"] = gdf["pop"] / gdf["area_km2"].replace(0, np.nan)

Conditional logic:

# ❌
def classify(row):
    if row["area"] > 10_000:
        return "large"
    elif row["area"] > 1_000 and row["class"] == "residential":
        return "medium residential"
    return "small"
gdf["cls"] = gdf.apply(classify, axis=1)

# ✅ conditions may reference several columns
gdf["cls"] = np.select(
    [gdf["area"] > 10_000,
     (gdf["area"] > 1_000) & (gdf["class"] == "residential")],
    ["large", "medium residential"],
    default="small",
)

Note & rather than and, and the parentheses around each comparison. NumPy's element-wise and is &, and its precedence is lower than > — omitting the parentheses gives a confusing error rather than a wrong answer, which is at least merciful.

Searching another layer — the highest-value replacement, because it fixes two problems at once:

# ❌ per-row, and no spatial index
def find_zone(row):
    hit = zones[zones.contains(row.geometry)]
    return hit["zone_name"].iloc[0] if len(hit) else None
gdf["zone"] = gdf.apply(find_zone, axis=1)

# ✅
gdf = gpd.sjoin(gdf, zones[["zone_name", "geometry"]],
                how="left", predicate="within").drop(columns="index_right")

apply + contains:    58.6 s   (200,000 points, 400 zones; timed on 2,000 points and scaled)
sjoin:                0.08 s

4. When the function really does need several columns

Some functions genuinely take several values. apply(axis=1) is still not the way to feed them:

def score(geom, population, road_class):
    """Custom logic that needs geometry and two attributes."""
    base = geom.area / 10_000
    weight = {"A": 1.4, "B": 1.0, "C": 0.7}.get(road_class, 0.5)
    return base * weight * (population ** 0.5)

# ❌ builds a Series per row
gdf["score"] = gdf.apply(
    lambda r: score(r.geometry, r["population"], r["road_class"]), axis=1)

# ✅ zip the raw arrays — no Series construction
gdf["score"] = [
    score(g, p, c) for g, p, c in
    zip(gdf.geometry.values, gdf["population"].to_numpy(), gdf["road_class"].to_numpy())
]

# ✅ or itertuples, when there are many columns
gdf["score"] = [score(r.geometry, r.population, r.road_class)
                for r in gdf.itertuples()]

apply(axis=1):            1.40 s
zip over arrays:          0.53 s
itertuples:               0.74 s

Two and a half times faster on 200,000 rows of a 41-column frame, and the function is unchanged. zip over .to_numpy() arrays hands your function raw values with no pandas machinery in between; itertuples builds a lightweight namedtuple instead of a Series.

Better still, decompose it — most such functions are a vectorisable part and a small non-vectorisable part:

# ✅ fully vectorised
base = gdf.geometry.area / 10_000
weight = gdf["road_class"].map({"A": 1.4, "B": 1.0, "C": 0.7}).fillna(0.5)
gdf["score"] = base * weight * np.sqrt(gdf["population"])

fully vectorised:         0.012 s

121× faster than the apply, with the same values. The .map handles the lookup, np.sqrt the maths, and pandas aligns the three Series by index automatically.

5. If you must apply, apply to the GeoSeries

# ❌ 2.6 s on 200,000 rows of a 41-column frame
gdf["result"] = gdf.apply(lambda r: custom(r.geometry), axis=1)

# ✅ 1.8 s — same Python calls, no Series construction
gdf["result"] = gdf.geometry.apply(custom)

Same number of calls to custom, 1.5× faster: the frame-level apply was also building 200,000 Series that custom never looked at. When the function itself is the expensive part, dropping the Series is a modest gain; the large ones come from not calling Python once per row at all.

For a function taking one column plus geometry:

gdf["result"] = [custom(g, v) for g, v in
                 zip(gdf.geometry.values, gdf["value"].to_numpy())]

6. Show progress on what remains

An irreducible loop over four million geometries takes real time, and a progress bar makes it bearable:

from tqdm import tqdm
tqdm.pandas(desc="scoring")
gdf["score"] = gdf.geometry.progress_apply(custom)

scoring: 100%|██████████| 200000/200000 [00:01<00:00, 106027.84it/s]

The rate is the useful number: at 106,028 rows per second, four million rows will take about 38 s, and you can decide whether that is worth optimising. See how to add a progress bar and ETA.

A before-and-after and a decision helper

Example 1: a before-and-after on a real function

import time
import numpy as np
import geopandas as gpd

# ── before: six apply(axis=1) calls, ~9 s on 400,000 rows ────────────────────
def enrich_slow(gdf, city_centre, ZONE_NAMES):
    gdf = gdf.copy()
    gdf["area_m2"] = gdf.apply(lambda r: r.geometry.area, axis=1)
    gdf["perimeter_m"] = gdf.apply(lambda r: r.geometry.length, axis=1)
    gdf["compactness"] = gdf.apply(
        lambda r: 4 * np.pi * r["area_m2"] / r["perimeter_m"] ** 2
        if r["perimeter_m"] else np.nan, axis=1)
    gdf["dist_centre_m"] = gdf.apply(
        lambda r: r.geometry.distance(city_centre), axis=1)
    gdf["zone_name"] = gdf.apply(
        lambda r: ZONE_NAMES.get(r["zone_code"], "unknown"), axis=1)
    gdf["size_class"] = gdf.apply(
        lambda r: "large" if r["area_m2"] > 10_000
        else "medium" if r["area_m2"] > 1_000 else "small", axis=1)
    return gdf

# ── after: no apply at all, ~0.14 s ─────────────────────────────────────────
def enrich_fast(gdf, city_centre, ZONE_NAMES):
    gdf = gdf.copy()
    area = gdf.geometry.area
    length = gdf.geometry.length

    gdf["area_m2"] = area
    gdf["perimeter_m"] = length
    gdf["compactness"] = 4 * np.pi * area / length.replace(0, np.nan) ** 2
    gdf["dist_centre_m"] = gdf.geometry.distance(city_centre)
    gdf["zone_name"] = gdf["zone_code"].map(ZONE_NAMES).fillna("unknown")
    gdf["size_class"] = np.select(
        [area > 10_000, area > 1_000], ["large", "medium"], default="small")
    return gdf

for name, fn in [("apply(axis=1) x6", enrich_slow), ("vectorised", enrich_fast)]:
    t0 = time.perf_counter()
    out = fn(parcels, centre, ZONE_NAMES)
    print(f"{name:<20} {time.perf_counter() - t0:>9.2f} s")

a, b = enrich_slow(parcels.head(1000), centre, ZONE_NAMES), \
       enrich_fast(parcels.head(1000), centre, ZONE_NAMES)
for col in ["area_m2", "compactness", "zone_name", "size_class"]:
    x, y = a[col].fillna(-1), b[col].fillna(-1)
    same = np.allclose(x, y, rtol=1e-12, atol=0) if x.dtype.kind == "f" else (x == y).all()
    print(f"  {col:<16} matches: {same}")

apply(axis=1) x6            9.44 s
vectorised                  0.14 s
  area_m2          matches: True
  compactness      matches: True
  zone_name        matches: True
  size_class       matches: True

66× faster with the same output. Compare floating-point columns with np.allclose, not ==: compactness computed column-wise and row by row differed by up to 2.2 × 10⁻¹⁶, and an exact comparison reports that rounding as a mismatch. The per-column check is the part worth keeping: vectorising changes how NaN, zero and missing values are handled, and it is easy to introduce a difference at the edges without noticing.

length.replace(0, np.nan) reproduces the original's if r["perimeter_m"] guard. Without it, the vectorised version produces inf where the original produced NaN — a small difference that would surface later as an outlier on a map.

Example 2: a decision helper for an unfamiliar apply

import inspect
import re

def suggest_replacement(source):
    """Heuristics for what a given apply could become."""
    src = source if isinstance(source, str) else inspect.getsource(source)
    body = src.strip()
    suggestions = []

    geom_methods = re.findall(r"\.geometry\.(\w+)", body) + \
                   re.findall(r"row\.geometry\.(\w+)", body)
    VECTORISED = {"area", "length", "centroid", "bounds", "is_valid", "is_empty",
                  "buffer", "simplify", "convex_hull", "envelope", "boundary",
                  "distance", "intersects", "contains", "within", "geom_type",
                  "representative_point", "make_valid", "exterior"}
    for m in set(geom_methods) & VECTORISED:
        suggestions.append(f"gdf.geometry.{m} is vectorised — use it directly")

    if re.search(r"\bif\b.*\breturn\b", body, re.S) and "elif" in body:
        suggestions.append("multi-branch conditional — use np.select")
    elif re.search(r"\bif\b.*\belse\b", body, re.S):
        suggestions.append("two-branch conditional — use np.where")

    if re.search(r"\.get\(\s*row\[", body) or re.search(r"\w+\[row\[", body):
        suggestions.append("dictionary lookup — use Series.map(mapping)")

    if re.search(r"\w+\[\w+\.(contains|intersects|within)\(", body):
        suggestions.append("searching another layer — use gpd.sjoin, which is indexed")

    if re.search(r"row\[[\"'][^\"']+[\"']\]\s*[-+*/]\s*row\[", body):
        suggestions.append("column arithmetic — operate on the columns directly")

    cols = set(re.findall(r"row\[[\"']([^\"']+)[\"']\]", body))
    if len(cols) <= 1 and "row.geometry" in body:
        suggestions.append("only geometry is used — apply on gdf.geometry, "
                           "not on the frame (several times faster)")

    return suggestions or ["no obvious vectorisation — zip the raw arrays "
                           "instead of apply(axis=1)"]

def classify(row):
    if row["area"] > 10_000:
        return "large"
    elif row["area"] > 1_000:
        return "medium"
    return "small"

for s in suggest_replacement(classify):
    print(f"  → {s}")

  → multi-branch conditional — use np.select

Pattern-matching source text is a heuristic, not an analysis, and it will miss things. Its value is on a codebase with dozens of apply calls written by several people over years: it sorts them into "obviously replaceable" and "needs thought" in a few seconds, which is where the effort should go first.

Example 3: parallelising an irreducible apply

When the function genuinely cannot be vectorised — it calls an external tool, or implements a real algorithm — spread it across processes:

import numpy as np
import shapely
import geopandas as gpd
from concurrent.futures import ProcessPoolExecutor

def _chunk_worker(args):
    """Runs in a separate process. Receives WKB and plain arrays."""
    wkb, values = args
    import shapely
    geoms = shapely.from_wkb(wkb)
    return [genuinely_irreducible(g, v) for g, v in zip(geoms, values)]

def parallel_apply(gdf, value_col, *, workers=6, chunks=None):
    chunks = chunks or workers * 4
    wkb = shapely.to_wkb(gdf.geometry.values)
    values = gdf[value_col].to_numpy()

    parts = list(zip(np.array_split(wkb, chunks), np.array_split(values, chunks)))
    with ProcessPoolExecutor(max_workers=workers) as pool:
        results = list(pool.map(_chunk_worker, parts))
    return np.concatenate([np.asarray(r) for r in results])

gdf["result"] = parallel_apply(gdf, "population", workers=6)

serial apply on geometry:     9.7 s
6 workers:                    1.97 s   (4.9x)

That is a buffer-and-simplify per geometry on 200,000 parcels. 4.9× rather than 6× is the normal shape: the missing fraction is serialising the chunks, starting the processes, and collecting results.

Two details. Geometries cross the process boundary as WKB bytes rather than Shapely objects — smaller, faster, and immune to version mismatches between interpreters. And there are four times as many chunks as workers, so the pool can rebalance when some chunks are far more expensive than others, which is normal since cost scales with vertex count.

Parallelism comes last for a reason worth stating plainly: a 66× vectorisation beats a 5× parallelisation, and it does not multiply memory by the worker count. Reach for processes only when the function is genuinely irreducible. The memory arithmetic is in parallel batch processing.

How pandas has to build a row

apply(axis=1) is slow for a reason that has almost nothing to do with your function, and everything to do with how pandas has to represent a row.

A DataFrame is stored by column: each column is a contiguous typed array. There is no row object anywhere in memory. So when apply(axis=1) promises to hand your function a row, it has to construct one — a Series whose values come from every column at that position, with an index built from the column names. That is an allocation, a copy per column, and an index construction, for every row, before your function is called.

This explains why the cost grows with column count as well as row count, even though the function uses one column: the same operation took 1.13 s on a 41-column frame and 0.91 s on a 2-column one.

It also explains why a list comprehension is fine. [g.area for g in gdf.geometry] makes the same 200,000 Python calls and is 2.9× faster than apply(axis=1), because it iterates an array of objects directly with nothing constructed per element. So the received wisdom "never loop in pandas" is not quite right. The accurate version is never construct a row object per iteration — which apply(axis=1), iterrows() and to_dict("records") all do.

The vectorised GeoSeries methods win the remaining 80× by moving the loop out of Python entirely. Since Shapely 2.0, gdf.geometry.area calls into compiled code that iterates the geometry array and calls GEOS per element without returning to the interpreter. One Python call, four million C operations.

Sequencing matters when both problems are present. A per-row loop that also searches another layer has two independent problems: Python overhead (constant factor) and no spatial index (complexity class). sjoin fixes both in one line: it was 744× faster than the per-row search on 200,000 points, far beyond anything tuning the loop could reach. That is why "replace the apply with a sjoin" is worth so much more than "replace the apply with a list comprehension" — it fixes the algorithm, not just the language.

Finally, the honest limit. Vectorisation is not always possible or always right. A function calling an external service per feature cannot be vectorised. A function with genuine sequential dependence cannot. And a vectorised operation materialises every result at once, so on data near the memory limit a chunked loop can be the correct choice. The decision procedure is the same one as always: profile first, establish that the apply really is the bottleneck, and then pick the replacement that fits — as described in how to profile a slow Python GIS script.

Column counts, iterrows and GeoSeries.apply

• apply(axis=1) cost scales with column count. Dropping unused columns helps even before vectorising.
• iterrows has the same problem and is worse; itertuples does not.
• gdf.geometry.apply(f) was 3.6× faster than gdf.apply(f, axis=1) for a geometry-only function — same calls, no Series.
• Use & and | with parentheses in vectorised conditions; and/or raise on arrays.
• np.select evaluates every condition, so avoid it when a later branch would raise on inputs an earlier one excludes.
• Vectorising can change edge-case results — inf versus NaN, different null handling. Compare column by column.
• .map(dict) is a C hash lookup; .map(function) is a Python call per element.
• apply on an empty frame can return a DataFrame rather than a Series, which breaks assignment. Guard with if len(gdf).
• tqdm.pandas() then .progress_apply gives a rate and an ETA for what remains.
• swifter and pandarallel auto-parallelise apply, which helps less than vectorising and adds a dependency.

Internal links

• How to replace a row loop with vectorised GeoPandas operations — the full translation guide
• Why GeoPandas is slow: the four real bottlenecks — where this sits among them
• How to profile a slow Python GIS script — confirming this is the bottleneck
• How to perform a spatial join in Python (GeoPandas) — replacing a per-row search
• How to use the spatial index directly in GeoPandas (sindex) — when sjoin does not fit
• How to speed up GeoPandas: tips for large datasets — the wider checklist
• How to add a progress bar and ETA to a long GIS batch job — for the loop that remains
• How to speed up batch GIS jobs with parallel processing — the last resort

FAQ

Why is apply(axis=1) so slow?

It constructs a pandas Series for every row, copying every column's value into it, before calling your function. The cost scales with column count as well as row count.

Is apply always wrong?

No. gdf.geometry.apply(f) was 3.6× faster than apply(axis=1) and is fine for a genuinely custom function. What to avoid is apply(..., axis=1) and iterrows.

My function needs several columns. What do I use?

Zip the raw arrays: [f(g, v) for g, v in zip(gdf.geometry.values, gdf["v"].to_numpy())]. That was 2.6× faster than apply(axis=1) with the function unchanged; decomposing the function into vectorised parts was 121× faster.

How do I vectorise an if/elif/else?

np.select([cond1, cond2], [val1, val2], default=val3). Conditions can reference several columns; combine them with & and | inside parentheses.

Why is my vectorised version giving slightly different results?

Usually edge cases: division by zero producing inf instead of a guarded NaN, or different null handling. Compare column by column on a small sample.

Should I use swifter or pandarallel?

They parallelise apply, which helps far less than vectorising — a 66× vectorisation beats a 5× parallelisation — and adds a dependency. Vectorise first.

What if the function calls an external service per feature?

That is genuinely irreducible. Loop over gdf.geometry rather than the frame, add a progress bar, and parallelise across processes if the calls are independent.