How to Replace a Row Loop with Vectorised GeoPandas Operations

Problem statement

On 200,000 rows this takes 2.6 s; on four million, about 52 s:

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

This does the same thing on the 200,000 rows in 0.004 s — 600 times faster:

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

The timings on this page come from tools/bench/vectorise.py in the site's repository, run on an AMD Ryzen 9 7900 with pandas 3.0.6, GeoPandas 1.1.4 and Shapely 2.1.2. The pipeline examples use real data: Delaware's 20,198 census blocks as parcels, its 262 tracts as zones and the TIGER 2024 primary and secondary roads.

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.

loop:    46.5 s    (20,198 points in 262 tracts; timed on 300 and scaled)
sjoin:    0.016 s

Geometry methods, np.where, map and spatial joins for loops

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

ZONE_NAMES = dict(zip(zones["GEOID"], zones["NAMELSAD"]))   # "10003014502" → "Census Tract 145.02"

# ❌
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:

lookup_df = zones.assign(county=zones["GEOID"].str[2:5], zone_code=zones["GEOID"],
                         land_km2=zones["ALAND"] / 1e6)[["county", "zone_code", "land_km2"]]
gdf["county"] = gdf["GEOID20"].str[2:5]
gdf = gdf.merge(lookup_df, on=["county", "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:

import time

# ❌ n × m comparisons, in Python — timed on 300 points and scaled, because the full run takes minutes
sample = points.head(300).copy()
start = time.perf_counter()
for idx, row in sample.iterrows():
    for _, zone in zones.iterrows():
        if zone.geometry.contains(row.geometry):
            sample.loc[idx, "zone"] = zone["zone_name"]
            break
loop = (time.perf_counter() - start) * len(points) / len(sample)

# ✅ indexed and vectorised
start = time.perf_counter()
points = gpd.sjoin(points, zones[["zone_name", "geometry"]],
                   how="left", predicate="within").drop(columns="index_right")
print(f"loop:  {loop:6.1f} s    ({len(points):,} points in {len(zones)} tracts; timed on {len(sample)} and scaled)")
print(f"sjoin: {time.perf_counter() - start:8.3f} s")

loop:    46.5 s    (20,198 points in 262 tracts; timed on 300 and scaled)
sjoin:    0.016 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["tract"], 0)
    totals[row["tract"]] += row.geometry.area

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

# ✅ and when the geometry should be merged too
merged = gdf.dissolve(by="tract", aggfunc={"POP20": "sum", "land_km2": "first"})
print(f"{len(merged)} tracts, {merged['POP20'].sum():,} people")

262 tracts, 989,948 people

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.

timed("❌ iterrows", lambda: [r.geometry.area for _, r in gdf.iterrows()])
timed("⚠ apply(axis=1)", lambda: gdf.apply(lambda r: r.geometry.area, axis=1))
timed("⚠ itertuples", lambda: [row.geometry.area for row in gdf.itertuples()])
timed("⚠ geometry.apply(f)", lambda: gdf.geometry.apply(lambda g: g.area))
timed("✅ iterate the GeoSeries", lambda: [g.area for g in gdf.geometry])
timed("✅ vectorised", lambda: gdf.geometry.area)

❌ iterrows                           0.2787 s
⚠ apply(axis=1)                      0.1114 s
⚠ itertuples                         0.0973 s
⚠ geometry.apply(f)                  0.0322 s
✅ iterate the GeoSeries              0.0395 s
✅ vectorised                         0.0027 s

Three tiers, and they held across runs. iterrows is slowest by far, with apply(axis=1) and itertuples well behind it but close to each other: all three build a row object. Iterating the geometries alone — a comprehension or geometry.apply — is two to three times cheaper again. The vectorised call is another order of magnitude below that. itertuples beats iterrows because it builds a plain namedtuple rather than a Series, but it still packs every column into each tuple, so select the columns you need first.

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

def expensive(geom):
    """A stand-in for a per-geometry call with no vectorised form."""
    return shapely.minimum_rotated_rectangle(geom).area

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

A Delaware pipeline before and after vectorising

The data used on this page

Delaware's 20,198 2020 census blocks stand in for parcels, its 262 census tracts for zones, and the TIGER 2024 primary and secondary roads of its three counties for roads — all from the US Census Bureau, all downloaded here. Every snippet on the page runs against these.

import time
import numpy as np
import pandas as pd
import geopandas as gpd
import shapely

TIGER = "https://www2.census.gov/geo/tiger"
parcels = gpd.read_file(f"{TIGER}/TIGER2020/TABBLOCK20/tl_2020_10_tabblock20.zip",
                        columns=["GEOID20", "POP20"]).to_crs(5070)
zones = gpd.read_file(f"{TIGER}/TIGER2020/TRACT/tl_2020_10_tract.zip",
                      columns=["GEOID", "NAMELSAD", "ALAND"]).to_crs(5070)
zones["zone_name"] = zones["GEOID"]
roads = pd.concat([gpd.read_file(f"{TIGER}/TIGER2024/ROADS/tl_2024_{county}_roads.zip")
                   for county in ("10001", "10003", "10005")])
roads = roads[roads["MTFCC"].isin(["S1100", "S1200"])].to_crs(5070).reset_index(drop=True)
roads["road_name"] = roads["FULLNAME"].fillna("unnamed")

gdf = parcels.copy()
gdf["zone_code"] = gdf["tract"] = gdf["GEOID20"].str[:11]     # a block's tract is its first 11 digits
points = parcels.set_geometry(parcels.geometry.representative_point())[["GEOID20", "geometry"]]
centre = city_centre = target = parcels.geometry.union_all().centroid


def timed(label, fn, repeat=3):
    """Print the best of `repeat` runs — the one least disturbed by other work."""
    times = []
    for _ in range(repeat):
        t0 = time.perf_counter()
        fn()
        times.append(time.perf_counter() - t0)
    print(f"{label:<36} {min(times):.4f} s")


print(f"{len(parcels):,} blocks, {len(zones)} tracts, {len(roads):,} road segments, EPSG:{parcels.crs.to_epsg()}")

20,198 blocks, 262 tracts, 1,170 road segments, EPSG:5070

Example 1: a real pipeline, before and after

import time
import numpy as np
import geopandas as gpd

# ── before: 4 loops, ~13 minutes on 20,198 census blocks ────────────────────
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, ~28 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

sample = parcels.sample(40, random_state=0)     # the loops would take minutes on all of them
t0 = time.perf_counter()
slow = process_slow(sample, zones, roads, centre)
t_slow = (time.perf_counter() - t0) * len(parcels) / len(sample)
t0 = time.perf_counter()
fast = process_fast(parcels, zones, roads, centre)
t_fast = time.perf_counter() - t0
print(f"loops        {t_slow:>9.1f} s   (timed on {len(sample)} blocks and scaled)")
print(f"vectorised   {t_fast:>9.1f} s")

loops            749.9 s   (timed on 40 blocks and scaled)
vectorised        28.8 s

26× — less than anywhere else on this page, and the reason is worth knowing: nearly all of the vectorised time is sjoin_nearest, which took 26.8 s on its own in tools/bench/vectorise.py. 7,584 of the blocks touch a primary or secondary road, so their nearest distance is zero and several roads tie; the join returned 34,148 extra rows before the deduplication below removed them. On the benchmark's 40-block sample every area, compactness and distance matched the loop. The road name differed for 18 blocks, all at identical distances: the loop keeps the first road it meets, sjoin_nearest another.

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 — defined for Polygons, so count on each block's
# first part (every block is stored as a MultiPolygon; only one has two parts)
gdf["n_holes"] = shapely.get_num_interior_rings(shapely.get_geometry(geoms, 0))

# 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: each block against its own tract
tract_of_block = zones.set_index("GEOID").geometry.loc[gdf["zone_code"]].values
gdf["overlap_m2"] = shapely.area(shapely.intersection(geoms, tract_of_block))

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

1,083,368 vertices, 1 multipart, 645 with holes

Timed against the loop equivalent:

parts = shapely.get_geometry(geoms, 0)            # the first part of each geometry
timed("shapely.get_num_coordinates", lambda: shapely.get_num_coordinates(parts))
timed("apply(len(g.exterior.coords))",
      lambda: gpd.GeoSeries(parts).apply(lambda g: len(g.exterior.coords)))

shapely.get_num_coordinates          0.0003 s
apply(len(g.exterior.coords))        0.1188 s

That is on the census blocks, and the comparison had to be made on their first parts: every block is stored as a MultiPolygon, even the 20,197 with one part, so gdf.geometry.apply(lambda g: len(g.exterior.coords)) raised AttributeError: 'MultiPolygon' object has no attribute 'exterior' on the first row. 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. The worker goes in its own module, because the child processes have to import it by name: on Windows, on macOS, and since Python 3.14 on Linux too, they start fresh rather than forking, so a function defined in a notebook or a __main__ script is not there for them to find.

from pathlib import Path

Path("irreducible.py").write_text('''
import shapely

def genuinely_irreducible(geom):
    """A stand-in for a per-geometry call with no vectorised form."""
    return shapely.minimum_rotated_rectangle(geom).area

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

import numpy as np
import shapely
from concurrent.futures import ProcessPoolExecutor
from irreducible import worker


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)
print(f"{len(gdf):,} results, same as a plain loop: "
      f"{np.allclose(gdf['result'], [shapely.minimum_rotated_rectangle(g).area for g in gdf.geometry])}")

20,198 results, same as a plain loop: True

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 600× vectorisation beats a 5× parallelisation and does not multiply memory by the worker count. The memory reasoning is in parallel batch processing.

Where the loop runs, not whether it exists

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. That is an allocation, a copy per column and an index, per row, before your function even runs. That is why, on the 20,198 census blocks, it took 0.111 s against 0.040 s for a list comprehension over the geometries — and iterrows, which builds a Series per row too, 0.279 s. The vectorised call took 0.003 s: next to it, every Python loop is slow, and the ones that build a row object are slowest.

iterrows() has the same problem for the same reason, which is why itertuples() — returning a lightweight namedtuple with no index and no copying — was about three times faster on the census blocks. The gap grows with the number of columns, because those are what iterrows copies.

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 — 46.5 s to 0.016 s in the measurement above, about 2,900× — 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.

Per-row Series cost, iterrows and np.select

• apply(axis=1) builds a Series per row. Its cost scales with column count; dropping unused columns helps before you even vectorise.
• itertuples beat iterrows by about 3× on the census blocks and is the right escape hatch for an irreducible loop.
• Iterating gdf.geometry was over ten times slower than vectorised; iterating rows with iterrows, over a hundred times.
• 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 made itertuples about three times faster for identical work: 0.097 s against 0.279 s on Delaware's 20,198 census blocks.

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.