How to Replace a Row Loop with Vectorised GeoPandas Operations #

Problem statement #

This runs for forty minutes:

areas = []
for idx, row in gdf.iterrows():
    areas.append(row.geometry.area)
gdf["area_m2"] = areas

This does the same thing in under a second:

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

Everyone knows the second is better. What is less obvious is which loops have a vectorised equivalent, what to do when the operation seems to need row context, and why some vectorised-looking code is still slow:

gdf["zone"] = gdf.geometry.apply(lambda g: find_zone(g, zones))    # still slow
gdf["dist"] = gdf.apply(lambda r: r.geometry.distance(target), axis=1)  # very slow

Both use apply, which looks vectorised and is not. Vectorising is a translation exercise: for most row loops there is a direct equivalent, and for the rest there is a standard replacement pattern.

Quick answer #

# before
for idx, row in gdf.iterrows():
    gdf.loc[idx, "buffered_area"] = row.geometry.buffer(50).area

# after
gdf["buffered_area"] = gdf.geometry.buffer(50).area

The two things never to do: iterrows() and apply(..., axis=1). Both build a pandas Series per row, copying every column.

Step-by-step solution #

1. Recognise which loops are trivially vectorisable #

Any loop whose body calls one method on row.geometry has a direct equivalent:

# ── measurements ──────────────────────────────────────────────────────────
gdf.geometry.area            gdf.geometry.length         gdf.geometry.bounds
gdf.geometry.is_valid        gdf.geometry.is_empty       gdf.geometry.geom_type
gdf.geometry.has_z           gdf.geometry.is_simple      gdf.geometry.count_coordinates

# ── transformations ───────────────────────────────────────────────────────
gdf.geometry.buffer(50)      gdf.geometry.centroid       gdf.geometry.envelope
gdf.geometry.simplify(1.0)   gdf.geometry.convex_hull    gdf.geometry.boundary
gdf.geometry.make_valid()    gdf.geometry.representative_point()

# ── pairwise, against one geometry or an aligned series ───────────────────
gdf.geometry.distance(target)          gdf.geometry.intersects(target)
gdf.geometry.contains(target)          gdf.geometry.intersection(other_series)

Each returns a Series or GeoSeries aligned to the index, so it assigns straight into a column:

gdf["area_m2"] = gdf.geometry.area
gdf["is_valid"] = gdf.geometry.is_valid
gdf["dist_to_centre"] = gdf.geometry.distance(city_centre)

Since Shapely 2.0 these dispatch to a C loop over the whole array, so the Python overhead is paid once rather than per row.

2. Vectorise conditions with np.where and np.select #

A loop with an if looks unvectorisable and is not:

# ❌
categories = []
for idx, row in gdf.iterrows():
    if row.geometry.area > 10_000:
        categories.append("large")
    elif row.geometry.area > 1_000:
        categories.append("medium")
    else:
        categories.append("small")
gdf["size_class"] = categories

# ✅ two branches
area = gdf.geometry.area
gdf["is_large"] = np.where(area > 10_000, "large", "not large")

# ✅ many branches — conditions are evaluated in order, first match wins
gdf["size_class"] = np.select(
    [area > 10_000, area > 1_000],
    ["large", "medium"],
    default="small",
)

# ✅ or binning, when the branches are ranges of one value
gdf["size_class"] = pd.cut(area, bins=[0, 1_000, 10_000, np.inf],
                           labels=["small", "medium", "large"])

np.select evaluates every condition over the whole array, so it does more total work than a short-circuiting loop — and is still orders of magnitude faster, because the work happens in C.

3. Replace dictionary lookups with map #

# ❌
for idx, row in gdf.iterrows():
    gdf.loc[idx, "zone_name"] = ZONE_NAMES.get(row["zone_code"], "unknown")

# ✅
gdf["zone_name"] = gdf["zone_code"].map(ZONE_NAMES).fillna("unknown")

.map accepts a dict, a Series or a function. With a dict it is a hash lookup per element in C; with a function it is a Python call per element, so prefer the dict.

For a lookup keyed on several columns, merge instead:

gdf = gdf.merge(lookup_df, on=["authority", "zone_code"], how="left")

4. Replace "find the matching feature" loops with a spatial join #

This is the most valuable translation in practice, because it fixes two bottlenecks at once:

# ❌ n × m comparisons, in Python
for idx, row in points.iterrows():
    for _, zone in zones.iterrows():
        if zone.geometry.contains(row.geometry):
            points.loc[idx, "zone"] = zone["name"]
            break

# ✅ indexed and vectorised
points = gpd.sjoin(points, zones[["name", "geometry"]],
                   how="left", predicate="within").drop(columns="index_right")

loop:   2,841 s
sjoin:      1.9 s

Nearest-feature loops have their own vectorised form:

# ✅ nearest, with the distance returned
joined = gpd.sjoin_nearest(points, roads[["road_name", "geometry"]],
                           how="left", max_distance=2_000,
                           distance_col="dist_to_road")

sjoin can return more rows than it received when a feature matches several — a point on a shared boundary matches both zones. That is correct behaviour, and it is the subject of spatial join cardinality explained.

5. Replace accumulation loops with groupby or dissolve #

# ❌
totals = {}
for idx, row in gdf.iterrows():
    totals.setdefault(row["ward"], 0)
    totals[row["ward"]] += row.geometry.area

# ✅
totals = gdf.assign(area=gdf.geometry.area).groupby("ward")["area"].sum()

# ✅ and when the geometry should be merged too
merged = gdf.dissolve(by="ward", aggfunc={"population": "sum", "value": "mean"})

dissolve unions the geometry per group and aggregates the attributes in one pass, which is a loop calling unary_union per group written properly.

6. When a loop is genuinely necessary, loop well #

Some operations have no vectorised form — a per-feature call to an external service, an algorithm with real sequential dependence, a library that takes one geometry at a time.

# ❌ 23 seconds for 200,000 rows
[r.geometry.area for _, r in gdf.iterrows()]

# ⚠ 3.4 seconds — acceptable when there is no alternative
[row.geometry.area for row in gdf.itertuples()]

# ✅ 0.12 seconds — iterate the GeoSeries, not the frame
[g.area for g in gdf.geometry]

# ✅ 0.04 seconds — the vectorised call
gdf.geometry.area

The ranking is stable and worth memorising. itertuples returns a lightweight namedtuple rather than a Series, so it avoids the per-row copying that makes iterrows and apply(axis=1) so expensive.

When you must loop, build a list and assign once rather than writing into the frame per iteration:

# ❌ .loc assignment per row reallocates repeatedly
for row in gdf.itertuples():
    gdf.loc[row.Index, "result"] = expensive(row.geometry)

# ✅
results = [expensive(row.geometry) for row in gdf.itertuples()]
gdf["result"] = results

Code examples #

Example 1: a real pipeline, before and after #

import time
import numpy as np
import geopandas as gpd

# ── before: 4 loops, ~52 minutes on 400,000 rows ────────────────────────────
def process_slow(parcels, zones, roads, city_centre):
    parcels = parcels.copy()
    for idx, row in parcels.iterrows():
        parcels.loc[idx, "area_m2"] = row.geometry.area
        parcels.loc[idx, "perimeter_m"] = row.geometry.length
        parcels.loc[idx, "compactness"] = (
            4 * np.pi * row.geometry.area / row.geometry.length ** 2)
        parcels.loc[idx, "dist_centre"] = row.geometry.distance(city_centre)

    for idx, row in parcels.iterrows():
        for _, z in zones.iterrows():
            if z.geometry.contains(row.geometry.centroid):
                parcels.loc[idx, "zone"] = z["zone_name"]
                break

    for idx, row in parcels.iterrows():
        nearest, best = None, float("inf")
        for _, r in roads.iterrows():
            d = row.geometry.distance(r.geometry)
            if d < best:
                nearest, best = r["road_name"], d
        parcels.loc[idx, "road"] = nearest
        parcels.loc[idx, "dist_road"] = best

    classes = []
    for idx, row in parcels.iterrows():
        a = row["area_m2"]
        classes.append("large" if a > 10_000 else "medium" if a > 1_000 else "small")
    parcels["size_class"] = classes
    return parcels

# ── after: no loops, ~7 seconds ─────────────────────────────────────────────
def process_fast(parcels, zones, roads, city_centre):
    parcels = parcels.copy()

    area = parcels.geometry.area
    length = parcels.geometry.length
    parcels["area_m2"] = area
    parcels["perimeter_m"] = length
    parcels["compactness"] = 4 * np.pi * area / length.replace(0, np.nan) ** 2
    parcels["dist_centre"] = parcels.geometry.distance(city_centre)

    centroids = parcels.set_geometry(parcels.geometry.representative_point())
    zoned = gpd.sjoin(centroids[["geometry"]], zones[["zone_name", "geometry"]],
                      how="left", predicate="within")
    parcels["zone"] = zoned["zone_name"].reindex(parcels.index)

    near = gpd.sjoin_nearest(parcels[["geometry"]], roads[["road_name", "geometry"]],
                             how="left", distance_col="dist_road")
    near = near[~near.index.duplicated(keep="first")]      # ties → keep one
    parcels["road"] = near["road_name"]
    parcels["dist_road"] = near["dist_road"]

    parcels["size_class"] = np.select(
        [area > 10_000, area > 1_000], ["large", "medium"], default="small")
    return parcels

for name, fn in [("loops", process_slow), ("vectorised", process_fast)]:
    t0 = time.perf_counter()
    out = fn(parcels, zones, roads, centre)
    print(f"{name:<12} {time.perf_counter() - t0:>9.1f} s")

loops           3128.4 s
vectorised         6.9 s

Four details in the fast version are worth calling out.

length.replace(0, np.nan) prevents a division-by-zero warning and produces NaN for degenerate geometry rather than inf, which would plot as an extreme value later.

representative_point() rather than centroid, because a centroid can fall outside a concave parcel and land in the wrong zone.

sjoin_nearest can return more rows than it received when distances tie exactly, so deduplicating on the index restores one row per parcel. Dropping that line produces a frame that is silently longer than the input.

And reindex(parcels.index) after the join guards against the join having reordered or dropped rows — assigning a misaligned Series is a subtle way to attach the wrong value to every row.

Example 2: vectorising with Shapely 2.0 functions directly #

For operations without a GeoSeries method, Shapely's array API works on the underlying GeometryArray:

import numpy as np
import shapely
import geopandas as gpd

geoms = gdf.geometry.values          # a GeometryArray, accepted by shapely functions

# vertex counts — no GeoSeries accessor for this
gdf["n_vertices"] = shapely.get_num_coordinates(geoms)

# part counts for multipart geometries
gdf["n_parts"] = shapely.get_num_geometries(geoms)

# interior ring (hole) counts
gdf["n_holes"] = shapely.get_num_interior_rings(geoms)

# force 2-D, dropping any Z
gdf["geometry"] = shapely.force_2d(geoms)

# snap coordinates to a grid, vectorised
gdf["geometry"] = shapely.set_precision(geoms, grid_size=0.001)

# a specific vertex from every geometry
first_points = shapely.get_point(shapely.get_exterior_ring(geoms), 0)

# pairwise between two arrays of equal length
gdf["overlap_m2"] = shapely.area(shapely.intersection(geoms, other.geometry.values))

print(f"{gdf['n_vertices'].sum():,} vertices, "
      f"{(gdf['n_parts'] > 1).sum():,} multipart, "
      f"{(gdf['n_holes'] > 0).sum():,} with holes")

473,520,912 vertices, 41,204 multipart, 8,412 with holes

Timed against the loop equivalent:

timed("shapely.get_num_coordinates", lambda: shapely.get_num_coordinates(geoms))
timed("apply(len(g.exterior.coords))",
      lambda: gdf.geometry.apply(lambda g: len(g.exterior.coords)))

shapely.get_num_coordinates          0.081 s
apply(len(g.exterior.coords))       18.442 s

Note that shapely.get_num_coordinates handles multipart geometries and holes correctly, while len(g.exterior.coords) raises on a MultiPolygon. The vectorised function is faster and more correct — a common pattern, because the array API was written to handle the general case.

Example 3: parallelising what remains #

When an operation genuinely cannot be vectorised — a per-feature call to an external tool, a complex custom algorithm — parallelise across processes:

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

def _worker(wkb_chunk):
    """Runs in a separate process. Takes and returns WKB, not Shapely objects."""
    import shapely
    geoms = shapely.from_wkb(wkb_chunk)
    results = [genuinely_irreducible(g) for g in geoms]
    return results

def parallel_apply(gdf, workers=4, chunks=None):
    chunks = chunks or workers * 4
    wkb = shapely.to_wkb(gdf.geometry.values)
    parts = np.array_split(wkb, chunks)

    with ProcessPoolExecutor(max_workers=workers) as pool:
        results = list(pool.map(_worker, parts))
    return np.concatenate([np.asarray(r) for r in results])

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

Two things matter here.

Serialise to WKB, not Shapely objects. Shapely geometries pickle, but WKB is smaller and faster to move between processes, and it avoids version-mismatch problems between parent and child interpreters.

Chunk generously. More chunks than workers gives the pool something to rebalance with when some geometries are far more expensive than others — which is normal, since cost scales with vertex count.

Parallelism is the last resort, and the ordering matters: vectorise first, because a 300× vectorisation beats a 6× parallelisation and does not multiply memory by the worker count. The memory reasoning is in parallel batch processing.

Explanation #

Vectorisation is about where the loop runs, not about whether there is a loop. There is always a loop over four million geometries; the question is whether Python executes it or C does.

gdf.geometry.area calls into Shapely 2.0, which iterates the GeometryArray in compiled code and calls GEOS per element without returning to the interpreter. One Python-level call, four million C-level operations.

gdf.geometry.apply(lambda g: g.area) makes four million Python function calls. Each has real overhead — frame creation, argument binding, attribute lookup, return — typically a microsecond or so, which is invisible once and four seconds at scale.

gdf.apply(lambda r: r.geometry.area, axis=1) is worse again, and the reason is the part people miss: it constructs a pandas Series for every row. Every column's value for that row is copied into a new object with an index. On a forty-column frame that is forty copies plus index construction, per row, before your function even runs. That is why it is 360× slower than the vectorised call while a list comprehension is only 3× slower — the gap is not "loops are slow", it is "building a Series per row is slow".

iterrows() has the same problem for the same reason, which is why itertuples() — returning a lightweight namedtuple with no index and no copying — is roughly thirty times faster.

Conditional logic feels unvectorisable and is not, because the vectorised form changes the strategy rather than the result. A loop with if/elif/else short-circuits: it evaluates the second condition only when the first fails. np.select evaluates every condition over the whole array and then selects. That is strictly more arithmetic, and still hundreds of times faster, because arithmetic in C is enormously cheaper than branching in Python.

Lookup loops are a different case, and the answer is a join. A loop searching another layer per feature is O(n×m) in Python; sjoin is an indexed algorithm in compiled code. The improvement there — 2,841 s to 1.9 s in the measurement above — is two effects multiplying: the complexity class changed and the constant factor changed. This is why "vectorise the loop" and "use a spatial index" are usually the same refactor, and why it is the highest-value change available.

Finally, the limit. Vectorisation trades memory for speed: gdf.geometry.buffer(50) materialises four million new geometries at once, where a loop would hold one. On data near the memory limit that trade can be the wrong one, and chunked vectorisation — vectorise within a chunk, loop over chunks — gets most of the speed at bounded memory. That is the same structure as processing a large GeoPackage in chunks, applied to computation rather than to reading.

Edge cases or notes #

• apply(axis=1) builds a Series per row. Its cost scales with column count; dropping unused columns helps before you even vectorise.
• itertuples beats iterrows by ~30× and is the right escape hatch for an irreducible loop.
• Iterating gdf.geometry is only ~3× slower than vectorised. Iterating rows is 300×.
• np.select evaluates all conditions, so avoid it when a branch would raise on inputs the earlier branch excludes.
• .map(dict) is a C hash lookup; .map(function) is a Python call per element. Prefer the dict.
• sjoin can return more rows than it received. Deduplicate on the index if one-to-one is required.
• sjoin_nearest can return several rows on exact distance ties. Same fix.
• Assigning a Series aligns on the index. After a join that reordered rows, reindex first or the values attach to the wrong rows.
• Vectorised operations materialise all results at once, so peak memory is higher than a loop's. Chunk if memory is tight.
• gdf.geometry.values is a GeometryArray that every shapely.* function accepts directly.

Internal links #

• Why GeoPandas is slow: the four real bottlenecks — where this fits among them
• GeoPandas .apply() takes hours: how to fix it — the specific symptom
• How to use the spatial index directly in GeoPandas (sindex) — for the loops sjoin does not cover
• How to perform a spatial join in Python (GeoPandas) — the replacement for search loops
• How to profile a slow Python GIS script — confirming the loop is the bottleneck
• How to speed up GeoPandas: tips for large datasets — the wider checklist
• Spatial join cardinality explained — why sjoin can add rows
• How to speed up batch GIS jobs with parallel processing — the last resort

FAQ #

Is apply vectorised? #

No. gdf.geometry.apply(f) makes one Python call per geometry, and gdf.apply(f, axis=1) additionally builds a Series per row. Neither is vectorised; the GeoSeries methods are.

What is the difference between iterrows and itertuples? #

iterrows yields a Series per row, copying every column. itertuples yields a lightweight namedtuple. That makes itertuples roughly 30× faster for identical work.

How do I vectorise an if/else? #

np.where for two branches, np.select for several, pd.cut for ranges of one value. All evaluate over the whole array, which is more work and far less time.

My loop searches another layer for each feature. What replaces it? #

gpd.sjoin for containment or intersection, gpd.sjoin_nearest for nearest. Both use a spatial index, so the improvement is algorithmic as well as vectorised.

What if the operation genuinely cannot be vectorised? #

Loop over gdf.geometry or itertuples, build a list, and assign once. Then consider ProcessPoolExecutor, passing WKB rather than Shapely objects between processes.

Does vectorising use more memory? #

Yes — it materialises every result at once, where a loop holds one. On data near the memory limit, vectorise within chunks and loop over the chunks.

Are there operations with no GeoSeries method? #

Several, and shapely.* functions accept gdf.geometry.values directly — get_num_coordinates, force_2d, set_precision and others. They are vectorised and often handle multipart geometry more correctly than a hand-written loop.