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

Problem statement #

The cell has been running for two hours and the progress bar says 14%.

gdf["area_m2"] = gdf.apply(lambda row: row.geometry.area, axis=1)

Four million rows. It will finish in about fourteen hours, and this does the same thing in under a second:

gdf["area_m2"] = gdf.geometry.area

The gap is not 10× or 100×. On a wide frame it is routinely a thousand times, and the reason is not that apply is a slow loop — it is that apply(axis=1) does something expensive on every single row that has nothing to do with your function.

Knowing what that is tells you which apply calls to fix first and what to replace them with.

Quick answer #

# 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.

Step-by-step solution #

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.041 s
narrow: apply(axis=1)                    14.882 s
wide (41 cols): geometry.area             0.041 s
wide (41 cols): apply(axis=1)            61.204 s

The vectorised call is unaffected by column count. apply(axis=1) is four times slower on a frame with four times the columns, 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.

So the cost is not the loop — it is 200,000 × 41 unnecessary value copies plus 200,000 Series constructions.

2. Find the apply calls that matter #

import re
from pathlib import Path

def find_applies(root="src"):
    hits = []
    for path in Path(root).rglob("*.py"):
        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}:{n}", line.strip()[:88]))
    for kind, where, code in sorted(hits, key=lambda h: {"iterrows": 0, "axis=1": 1}.get(h[0], 2)):
        mark = {"iterrows": "✗✗", "axis=1": "✗ ", "apply": "· "}[kind]
        print(f"  {mark} {where:<34} {code}")
    return hits

find_applies()

  ✗✗ src/clean.py:88                  for idx, row in gdf.iterrows():
  ✗  src/pipeline.py:42               gdf["area"] = gdf.apply(lambda r: r.geometry.area, axis=1)
  ✗  src/pipeline.py:61               gdf["zone"] = gdf.apply(find_zone, axis=1)
  ·  src/clean.py:12                  gdf["ok"] = gdf.geometry.apply(is_acceptable)

The two-tier marking matters. iterrows and axis=1 are almost always fixable; a plain .apply on a GeoSeries with a genuinely custom function is only 3× slower than vectorised and 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:   412.9 s
sjoin:                1.9 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):           61.2 s
zip over arrays:          1.8 s
itertuples:               3.4 s

Thirty-four times faster 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.09 s

680× faster than the apply. 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 #

# ❌ 61 s
gdf["result"] = gdf.apply(lambda r: custom(r.geometry), axis=1)

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

Same number of calls to custom, 68× faster, because the frame-level apply was building 200,000 Series that custom never looked at.

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%|██████████| 4012884/4012884 [06:12<00:00, 10784.11it/s]

The rate is the useful number: at 10,784 rows per second you can predict the finish and decide whether to optimise further. See how to add a progress bar and ETA.

Code examples #

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

import time
import numpy as np
import geopandas as gpd

# ── before: three apply(axis=1) calls, ~11 minutes 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, ~1.4 seconds ────────────────────────────────────
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"]:
    same = (a[col].fillna(-1) == b[col].fillna(-1)).all()
    print(f"  {col:<16} identical: {same}")

apply(axis=1) x6         664.20 s
vectorised                 1.42 s
  area_m2          identical: True
  compactness      identical: True
  zone_name        identical: True
  size_class       identical: True

468× faster with identical output. The per-column equality 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 (typically 60x 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:   412.8 s
6 workers:                   78.4 s   (5.3x)

5.3× 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 468× 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.

Explanation #

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 the measurement that surprises people: the cost scales with column count, not just row count. The same operation on a 41-column frame is four times slower than on an 11-column frame, even though the function uses one column in both cases.

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 120× 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 3× 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). Vectorising alone gives perhaps 100×; adding the index gives another 100×; sjoin gives both in one line. 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.

Edge cases or notes #

• apply(axis=1) cost scales with column count. Dropping unused columns helps even before vectorising.
• iterrows has the same problem and is worse. itertuples is roughly 30× faster.
• gdf.geometry.apply(f) is ~60× 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) is only about 3× slower than a vectorised call 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 is typically 30× faster than apply(axis=1) with the function unchanged.

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 468× vectorisation beats a 6× 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.