Shapely 1.x vs Shapely 2.0 Vectorization

Shapely 2.0 was a ground-up rewrite around a NumPy-style vectorized API, and it changes both how fast geometry operations run and how you should write them. Code that loops over geometries the 1.x way still works, but leaves most of the speedup on the table. This guide compares the two idioms and shows how to convert per-object loops into array calls, and it is for anyone whose Shapely-heavy code feels slow on large feature sets. It sits under Shapely Geometry Operations in Mastering Core Geospatial Python Libraries, and it pairs with Shapely vs GeoPandas: When to Use Each for deciding whether to reach for the array API at all or let a GeoDataFrame do it for you.

Why This Approach / What Goes Wrong

In Shapely 1.x, every operation went through a bound Python method on a single geometry object, so processing a million parcel polygons meant a million Python-level GEOS round trips. Each call paid the cost of attribute lookup, argument marshalling, and a trip across the C boundary — overhead that dwarfs the actual GEOS computation for cheap operations like .area or .centroid. At scale, your program spends most of its time in the interpreter, not in the geometry engine.

Shapely 2.0 was rebuilt on pygeos, which stores geometries as opaque pointers inside NumPy object arrays and exposes module-level functions — shapely.area, shapely.buffer, shapely.intersects, shapely.distance — that accept those arrays and execute the loop in C. One Python call now covers the whole array, so the per-object overhead is amortized away. On bulk work the speedup is typically 5–100x, and it grows with the array size because the fixed dispatch cost is paid once instead of per feature.

Two things go wrong during migration. First, a handful of 1.x behaviors changed: iterating a multipart geometry no longer yields its parts (use shapely.get_parts), the old .ctypes/array interface is gone, and top-level names such as cascaded_union were renamed (shapely.union_all). Old code can break subtly rather than loudly. Second — and far more common — people port working 1.x code to 2.0 and keep the loop, calling .apply() or a Python for over object methods. That runs correctly but captures none of the gain, because the loop is still in Python. The whole point is to hand the array to a shapely.* function and let C do the iterating.

The concrete API changes, roughly in the order they bite during a port:

• shapely.ops.cascaded_union(list) → shapely.union_all(array). shapely.ops.unary_union still works and now forwards to the same implementation, so the rename is cosmetic — but only the array form avoids building an intermediate Python list first.
• Multipart geometries are no longer iterable. for part in footprint: raises TypeError in 2.0. Use footprint.geoms for one geometry, or shapely.get_parts(array) to explode a whole array at once.
• geom.type → geom.geom_type. The short form was deprecated in 2.0 and removed in 2.1, so code that only emits a warning today breaks on the next upgrade.
• == changed meaning, silently. In 1.x geom_a == geom_b delegated to .equals(), a topological test that calls two polygons equal even when their vertices are ordered differently. In 2.0 == is exact structural equality: same coordinates, same order. Deduplication built on ==, set(), or dict keys therefore returns different results after the upgrade without raising anything. Use shapely.equals() for the topological test and shapely.equals_exact(a, b, tolerance=...) when you need a tolerance.
• shapely.speedups.enable() does nothing. The module survives as a no-op stub because the fast path is now unconditional. Delete the call rather than guarding it.
• geom.ctypes, geom.array_interface(), and the asShape/asPoint adapters are gone. shapely.get_coordinates() replaces all of them and returns a plain (n, 2) array.
• Geometries are immutable. Assigning to geom.coords raises instead of mutating in place; construct a new geometry, or edit the coordinate array and rebuild.
• shapely.vectorized.contains(geom, x, y) was removed. Its replacement is the ordinary predicate: shapely.contains(geom, shapely.points(x, y)).

If you already work through a GeoDataFrame, GeoPandas routes its column operations through the Shapely 2.0 array functions for you, so parcels.geometry.area is already vectorized — the manual array API below matters most when you hold raw geometries outside a DataFrame or need a function GeoPandas does not surface as a method.

Prerequisites

• shapely>=2.0
• numpy>=1.26
• geopandas>=0.14 (optional, for the DataFrame comparison)

conda install -c conda-forge "shapely=2.0.*" "numpy=1.26.*" "geopandas=0.14.*"

Step-by-Step Implementation

The running example is a set of environmental sensor locations that we buffer to a fixed catchment radius and then test against a study-area boundary. All coordinates are in metres (ETRS89 / UTM 32N, EPSG:25832), so a buffer distance is a real distance — buffering in geographic degrees (EPSG:4326) would make the radius meaningless, the same axis-and-units discipline covered in Coordinate Systems with PyProj.

1. The 1.x idiom — a Python loop over objects (still valid, but slow at scale). Every .buffer and .area is a separate GEOS round trip.

from shapely.geometry import Point

# 250k sensor readings, easting/northing in metres (EPSG:25832)
sensor_xy = [(500_000 + i, 5_600_000 + i) for i in range(250_000)]
sensors = [Point(x, y) for x, y in sensor_xy]

# One GEOS round trip per object — the interpreter is the bottleneck
catchments_loop = [s.buffer(150) for s in sensors]          # 150 m catchments
catchment_areas_loop = [c.area for c in catchments_loop]     # m² per sensor

2. The 2.0 idiom — build the array once, operate on it in one call. shapely.points takes the x and y NumPy arrays directly and returns a geometry array; every subsequent call loops in C.

import numpy as np
import shapely

eastings = 500_000 + np.arange(250_000, dtype=float)
northings = 5_600_000 + np.arange(250_000, dtype=float)

sensors = shapely.points(eastings, northings)      # geometry array, no Python loop
catchments = shapely.buffer(sensors, 150)          # whole array, C-level loop
catchment_areas = shapely.area(catchments)         # returns a plain NumPy float array

3. Vectorized predicates replace per-pair loops. A spatial test that would be a nested for in 1.x becomes a single boolean array. Here we flag which sensor catchments fall inside a rectangular study area.

import numpy as np
import shapely

study_area = shapely.box(500_000, 5_600_000, 560_000, 5_660_000)  # metres
inside = shapely.contains(study_area, catchments)   # boolean NumPy array, one call
print("Catchments inside study area:", int(inside.sum()))

# Predicates broadcast, so element-wise pair tests are just as cheap:
gauges = shapely.points(eastings[::10], northings[::10])
near_gauge = shapely.dwithin(sensors[::10], gauges, 50)   # within 50 m, pairwise

4. Reductions collapse an array to one geometry without a Python accumulator. The 1.x cascaded_union over a list is now shapely.union_all over the array — a single call instead of a folding loop.

import shapely

# Merge overlapping catchments into a single coverage footprint
coverage = shapely.union_all(catchments)
print("Coverage type:", coverage.geom_type, "| area m²:", round(shapely.area(coverage)))

5. Inside a GeoDataFrame it is automatic — column operations already route through the 2.0 array functions, so you rarely touch shapely.* by hand once the data is tabular. Keep the CRS projected before any metric column.

import geopandas as gpd

parcels = gpd.read_file("parcels.gpkg").to_crs(epsg=25832)  # metres, not degrees
parcels["area_m2"] = parcels.geometry.area                  # Shapely 2.0 under the hood
parcels["catchment"] = parcels.geometry.buffer(150)         # vectorized, no .apply()

6. Coordinate math moves from a per-geometry callback to one array. The 1.x shapely.ops.transform(func, geom) invokes your function once per geometry; the 2.0 shapely.transform(geom, func) hands the function a single (n, 2) block covering the whole input and rebuilds the geometry structure around the result.

import numpy as np
import shapely
from shapely.ops import transform as transform_1x

# 1.x idiom: one call per geometry, callback receives loose x/y sequences
def shift_1x(x, y):
    return x + 100.0, y - 50.0

shifted_loop = [transform_1x(shift_1x, geom) for geom in catchments]

# 2.0 idiom: one call for the array, callback receives the coordinate block
shifted = shapely.transform(catchments, lambda xy: xy + [100.0, -50.0])

assert shapely.equals(shifted_loop[0], shifted[0])

This is the pattern that matters most when reprojecting geometry by hand outside a DataFrame: pyproj's Transformer.transform(xs, ys) is itself vectorized, so feeding it one coordinate block costs a single PROJ entry instead of one per feature. Build the transformer with always_xy=True so a geographic source CRS hands back (x, y) rather than (lat, lon) — the units-and-axis discipline in Coordinate Systems with PyProj applies unchanged, only the loop moves.

Verification

Confirm the vectorized path gives numerically identical results and is materially faster. The assert guards correctness; the timings prove the speedup is real rather than assumed.

import numpy as np, shapely, time
from shapely.geometry import Point

n = 200_000
eastings = 500_000 + np.arange(n, dtype=float)
northings = 5_600_000 + np.arange(n, dtype=float)

# 1.x-style loop: one GEOS round trip per sensor
t0 = time.perf_counter()
loop_areas = np.array([Point(x, y).buffer(5).area
                       for x, y in zip(eastings, northings)])
t_loop = time.perf_counter() - t0

# 2.0-style: build once, buffer once, measure once — all in C
t0 = time.perf_counter()
vec_areas = shapely.area(shapely.buffer(shapely.points(eastings, northings), 5))
t_vec = time.perf_counter() - t0

assert np.allclose(loop_areas, vec_areas), "Vectorized result must match the loop"
print(f"shapely {shapely.__version__}")
print(f"loop: {t_loop:.2f}s   vectorized: {t_vec:.2f}s   speedup: {t_loop/t_vec:.0f}x")
# shapely 2.0.6
# loop: 1.93s   vectorized: 0.07s   speedup: 28x

If shapely.__version__ prints 1.x the module functions will not exist and the import path above fails — that is itself the confirmation you are still on the old release.

Do not take the 28× from that run as a constant. The ratio tracks how much work GEOS does per geometry, not how many geometries there are. For cheap operations — area, bounds, centroid, is_valid, distance between points — almost all of the 1.x cost was Python dispatch, so removing it removes almost all of the runtime and 20–100× is normal. For expensive ones — buffering a 5,000-vertex administrative boundary, unioning complex polygons, make_valid on damaged rings — GEOS was already the bottleneck and the dispatch overhead was noise, so the same rewrite buys 1.1–1.5× and stops. Two practical consequences: benchmark on your own geometry rather than on points, and if you vectorized a buffer-heavy pipeline and saw nothing, the array API was never the constraint. Reduce the vertex count with simplify(tolerance, preserve_topology=True) or prune candidates with an index instead — the buffer-specific version of that argument is in Optimizing Buffer Operations for Large Datasets.

Memory moves the opposite way from speed. The 1.x loop discards each intermediate as soon as it has been appended or consumed, while a vectorized call materialises the entire output array before it returns: shapely.buffer(sensors, 150) over ten million points holds ten million polygons at once, each one a GEOS object with its own coordinate buffer. On wide geometry that is the difference between a job that finishes and an OOMKilled. Chunk the array — process a million at a time and reduce as you go — or push the whole workload into partitions with Scaling with Dask-GeoPandas. Converting an existing Python list with np.asarray(list_of_geoms) is cheap by comparison: it wraps the same GEOS pointers rather than copying geometry.

Edge Cases & Debugging

• AttributeError after upgrading. Renamed or removed 1.x helpers — cascaded_union → shapely.union_all, shapely.ops.unary_union still works but the array form is faster. Update the import rather than pinning back to 1.x.
• Iterating a multipolygon returns nothing new. In 1.x for part in multipolygon yielded the parts; in 2.0 a geometry is no longer iterable that way. Use shapely.get_parts(multipolygon) (or shapely.get_geometry(geom, i)) to reach components.
• Still slow on 2.0. You are looping with .apply(lambda g: g.area) or a Python for; pass the whole array to a shapely.* function so the loop runs in C. This is the single most common cause of "we upgraded and nothing got faster."
• Mixing object methods and module functions. catchment.area (method) and shapely.area(array) (function) both work, but only the array form vectorizes. Be consistent inside hot loops so you do not accidentally drop back to per-object dispatch.
• None or empty geometries in the array. Vectorized functions propagate missing values instead of raising: filter with shapely.is_missing / shapely.is_empty before a reduction like union_all, or the result can collapse to an empty geometry.
• Wrong units in buffers and areas. A buffer(150) on lon/lat coordinates buffers by 150 degrees, not metres. Reproject to a metric CRS such as EPSG:25832 first — never measure in EPSG:4326 or EPSG:3857 (Web Mercator distorts scale).
• GeoPandas resolved an older Shapely. A mixed conda/pip environment can leave shapely<2 installed even with recent GeoPandas; verify with shapely.__version__ and reinstall from conda-forge so GEOS stays ABI-compatible.
• Deduplication started keeping duplicates. == is structural equality in 2.0, not .equals(). Two rings with the same shape but a different start vertex now compare unequal. Switch the comparison to shapely.equals() or normalise with shapely.normalize() before hashing.
• TypeError: 'MultiPolygon' object is not iterable. 1.x part iteration was removed. geom.geoms for one geometry, shapely.get_parts(array) for a whole column.

Frequently Asked Questions

Do I have to rewrite working 1.x-style code after upgrading? Only where it is hot. Object-method code runs correctly on 2.0, so a cold path that touches a few hundred geometries is not worth touching. What does force a rewrite is the removed API — part iteration, .ctypes, mutation, shapely.vectorized — and the deprecations that 2.1 finished removing, such as geom.type. Convert the hot loops for speed and the removed calls for survival; leave the rest.

Can I pin shapely<2 and skip the migration entirely? Not for long. GeoPandas 1.0 requires shapely>=2.0, and 1.8 does not build against current NumPy or recent GEOS releases, so pinning back also pins your GeoPandas, your NumPy, and eventually your Python version. Treat the pin as a scheduling tool for a week, not an architecture decision.

How do I find which loops are worth converting? Profile with cProfile and sort by cumulative time, then look for frames inside shapely whose call count equals your feature count — that is a per-object method being dispatched once per geometry, and it is exactly what a shapely.* array call collapses. If the same profile shows the time landing inside GEOS rather than in dispatch, converting the loop will not help and the answer is fewer or simpler geometries.

Is .apply() on a GeoSeries ever faster than the array API? Never for pure geometry math — it reintroduces exactly the per-feature Python dispatch that 2.0 removed. It is the right tool only when the logic genuinely branches per shape, such as repairing invalid geometry conditionally; that trade-off is drawn in Shapely vs GeoPandas: When to Use Each.

My vectorized version is slower than the loop on a small array. Why? Below roughly a thousand geometries the fixed cost of allocating the output array and entering the ufunc machinery is comparable to the dispatch it saves, so the two are within noise of each other. The array API is a scaling win, not a micro-optimisation — if the array is small enough for the difference to be invisible, pick whichever form reads more clearly.