Shapely vs GeoPandas: When to Use Each for Spatial Tasks

Choosing between Shapely and GeoPandas comes down to one question — are you holding one geometry or a whole table of them — and getting it wrong means either hand-looping work that should be vectorized or spinning up a DataFrame to validate a single polygon. This guide is for anyone deciding how to structure a spatial pipeline: it sits under Shapely Geometry Operations in Mastering Core Geospatial Python Libraries, and pairs with Shapely 1.x vs Shapely 2.0 Vectorization, which explains why the array path is fast.

Why This Approach / What Goes Wrong

The two libraries are not competitors — GeoPandas is built on top of Shapely. Every geometry in a GeoDataFrame's active geometry column is a Shapely object, and GeoPandas' .area, .buffer(), .intersects() accessors dispatch to Shapely's vectorized array functions across the whole column in a single C call. So "Shapely vs GeoPandas" is really a question of altitude: do you operate on one geometry, or on a table of them plus their attributes?

The naive failures come from working at the wrong altitude. Reaching for GeoPandas to validate or transform a single geometry means paying for a DataFrame, an index, and a CRS registry you never use — noise, not error, but noise that clutters a codebase. The more expensive mistake runs the other way: staying in Shapely for tabular work. Writing a Python for loop over ten million geometries, or calling gdf.geometry.apply(lambda g: g.area), executes one Python-level GEOS call per feature. The vectorized gdf.geometry.area runs the same math as a single batched call and is typically 10–50× faster on large datasets. That gap is the whole reason GeoPandas exists — bypassing it defeats the point.

The correct pattern for real pipelines is hybrid. Let GeoPandas own what it is good at — file I/O, attribute joins, projection management, and vectorized column math — and drop to Shapely only for the genuinely per-shape logic that cannot be expressed as one array operation: conditional geometry repair, custom intersection branching, bespoke construction. The bridge between the two is GeoSeries.apply(), which hands each Shapely geometry to your function while preserving the DataFrame's schema and index.

Prerequisites

• geopandas>=1.0 — ships with pyogrio as the default I/O engine and Shapely 2.0 vectorization for column operations
• shapely>=2.0 — the geometry engine, exposing make_valid and array-level functions
• pyproj>=3.4 — CRS objects and transformations behind to_crs()
• pyogrio>=0.7 — the vectorized read/write engine GeoPandas delegates file I/O to

conda install -c conda-forge "geopandas=1.0.*" "shapely=2.0.*" "pyproj=3.6.*" "pyogrio=0.7.*"

Install from conda-forge rather than mixing pip wheels — GeoPandas, Shapely, and pyproj each bind the GEOS and PROJ C libraries, and letting conda resolve one consistent build avoids the ABI and PROJ-data mismatches that produce silent CRS failures on import.

The version boundary matters more than usual, because GeoPandas 1.0 removed most of what this decision used to hinge on:

• The PyGEOS backend is gone. Versions 0.10–0.13 shipped an optional array backend toggled with geopandas.options.use_pygeos; 1.0 removed the option because Shapely 2.0 is the backend. Advice written against that flag no longer applies.
• unary_union became union_all(). The old property is deprecated; the method form is what dissolves and overlays call internally.
• pyogrio replaced fiona as the default I/O engine, which is why file reads got several times faster on 1.0 without any code change. engine="fiona" still exists for drivers pyogrio does not cover.
• sjoin(op=...) is now sjoin(predicate=...). The old keyword was deprecated across the 0.x line and removed.
• Several operations gained vectorized accessors — make_valid(), segmentize(), concave_hull(), sample_points() among them. Each of those is a .apply() you no longer need to write, so check the current GeoSeries API before reaching for a lambda; a surprising amount of "I have to drop to Shapely for this" advice predates 1.0.

Step-by-Step Implementation

The canonical pipeline shows both altitudes in one flow: GeoPandas handles the file and the projection, Shapely handles the per-feature repair, and .apply() stitches the result back into the table.

1. Load the table and inspect its CRS — GeoPandas owns I/O.

import geopandas as gpd

# GeoPandas reads the file and carries the CRS off disk automatically
urban_zones = gpd.read_file("urban_zones.gpkg")
print("Source CRS:", urban_zones.crs)          # e.g. EPSG:4326 (degrees!)

2. Reproject to a metric CRS before any area or distance work.

Shapely computes everything in the plane using raw coordinate numbers — it has no concept of a datum. If those numbers are longitude/latitude degrees, .area returns square degrees, which are meaningless. Project to an appropriate UTM zone first; never measure in geographic degrees or in Web Mercator (EPSG:3857), whose area distortion grows severely away from the equator. estimate_utm_crs() picks the right zone from the data's extent — the mechanics are covered in choosing a UTM zone automatically in Python, and the broader transformation model in Coordinate Systems with PyProj.

# Reproject to the UTM zone that fits the data — metres, not degrees
metric_crs = urban_zones.estimate_utm_crs()
urban_zones = urban_zones.to_crs(metric_crs)
print("Working CRS:", urban_zones.crs.axis_info[0].unit_name)   # 'metre'

3. Do the per-feature work in Shapely, bridged by .apply().

Area repair needs branching — skip nulls and empties, fix invalid rings, then measure — which cannot be written as a single vectorized call. This is exactly where Shapely earns the .apply() bridge. make_valid calls GEOS to repair self-intersections and other topology errors covered in topology validation and repair.

from shapely.validation import make_valid

def valid_area_m2(geom) -> float:
    """Per-geometry repair-then-measure — genuine branching, so .apply() is correct."""
    if geom is None or geom.is_empty:
        return 0.0
    if not geom.is_valid:
        geom = make_valid(geom)          # GEOS repairs self-intersections
    return geom.area

urban_zones["area_m2"] = urban_zones.geometry.apply(valid_area_m2)

4. Prefer the vectorized accessor when no branching is needed.

If the geometries are already clean, do not .apply() at all — the vectorized column accessor runs the identical math in one batched C call and is far faster. Use .apply() only for the conditional path above.

# When repair is not required, this is the fast path — one C call, whole column
urban_zones["area_fast_m2"] = urban_zones.geometry.area
print(urban_zones[["zone_id", "area_m2"]].head())

5. Drop to shapely.ops for the operations a column accessor cannot express.

There is a clean rule for what GeoPandas will never surface: a GeoSeries accessor is row-aligned, so it can only produce one geometry per input row from parameters that are constant or broadcastable. Anything that changes the row count, or that needs a per-row parameter drawn from another column, lives in shapely.ops and needs an explicit bridge.

import geopandas as gpd
from shapely.ops import substring, polygonize

pipes = gpd.read_file("pipes.gpkg").to_crs(epsg=25832)   # metres

# Per-row parameters: each pipe is cut around its own recorded chainage
def inspection_reach(row):
    start = max(row.chainage_m - 50.0, 0.0)
    end = min(row.chainage_m + 50.0, row.geometry.length)
    return substring(row.geometry, start, end)

pipes["reach"] = gpd.GeoSeries(
    pipes.apply(inspection_reach, axis=1), crs=pipes.crs
)

# Cardinality change: 3,400 noded street segments become 812 enclosed blocks
streets = gpd.read_file("streets.gpkg").to_crs(pipes.crs)
blocks = gpd.GeoDataFrame(
    geometry=list(polygonize(streets.geometry.values)), crs=streets.crs
)
print(len(streets), "->", len(blocks))       # 3400 -> 812

split, substring, polygonize, and nearest_points all fall on that side of the line, and the gpd.GeoSeries(..., crs=...) wrapper in the first example is not optional: apply() returns a plain pandas Series of Shapely objects with no CRS at all, and assigning that straight into a GeoDataFrame produces a column that later silently refuses to reproject.

6. Recognise when the answer is neither library.

Both Shapely and GeoPandas are in-memory. The decision to move somewhere else is a memory-budget question, and it is worth measuring rather than guessing, because deserialised geometry costs several times its on-disk size:

import geopandas as gpd

parcels = gpd.read_parquet("parcels.parquet")
resident_gb = parcels.memory_usage(deep=True).sum() / 1e9
print(f"{len(parcels):,} rows, {resident_gb:.2f} GB resident")
# 4,182,996 rows, 6.41 GB resident      (the GeoParquet file on disk: 1.9 GB)

Expect roughly three to four times the file size once compact WKB has been inflated into individual GEOS objects with their own coordinate buffers, plus a pandas index. An overlay or join then needs headroom for the result on top of that. The practical threshold: once the load alone consumes a third of available memory, stop reaching for a bigger machine. Push the filter and the join into a query engine that streams — DuckDB Spatial Analytics for analytical scans over Parquet, PostGIS Integration with Python when the data already lives in a database and an index can do the work, or Scaling with Dask-GeoPandas when the workload really is a GeoPandas workload that simply does not fit in one process. All three return small results to GeoPandas at the end, which is where the familiar API earns its place again.

Verification

Confirm the pipeline produced projected, positive, finite areas — and that the vectorized and .apply() paths agree on already-valid geometry, proving the only difference between them is speed, not result.

import numpy as np

# The working CRS must be projected, or every area is nonsense
assert urban_zones.crs.is_projected, "Reproject to a metric CRS before measuring"

# Areas must be finite and non-negative
areas = urban_zones["area_m2"].to_numpy()
assert np.isfinite(areas).all() and (areas >= 0).all()

# On valid geometry, the fast accessor equals the repaired path
valid = urban_zones[urban_zones.geometry.is_valid]
assert np.allclose(valid.geometry.area, valid["area_m2"].loc[valid.index])
print("OK:", len(urban_zones), "zones,", int((areas == 0).sum()), "empty/null")
# OK: 1284 zones, 3 empty/null

Add one more assertion wherever a .apply() bridge produced geometry rather than a number, because that is the failure the numeric checks above cannot see: a geometry column that came back from Shapely without a CRS still plots, still measures, and still joins — right up until something calls to_crs() on it.

import geopandas as gpd

# The `reach` column from step 5 was rebuilt from Shapely objects — re-check it
assert isinstance(pipes["reach"], gpd.GeoSeries), "apply() returned a bare Series"
assert pipes["reach"].crs == pipes.crs, "CRS lost crossing the bridge"
assert pipes["reach"].is_valid.all(), "substring did not produce valid geometry"
assert (pipes["reach"].length <= pipes.geometry.length + 1e-6).all(), "reach overran its pipe"

Run those together rather than separately: the type check catches a missing gpd.GeoSeries() wrapper, the CRS check catches a wrapper applied without crs=, the validity check confirms the per-row branch actually ran, and the length check catches a chainage column recorded in the wrong units — the one error that produces perfectly valid geometry in entirely the wrong place.

Edge Cases & Debugging

• .apply() used for pure math. gdf.geometry.apply(lambda g: g.area) is the classic anti-pattern — it forces per-object Python. Replace with gdf.geometry.area; the difference is typically 10–50× on large tables. Reserve .apply() for logic with real branching.
• CRS mismatch before a join or distance. sjoin and distance ops silently misbehave when layers disagree on CRS. Check gdf.crs on both sides and to_crs() one to match; use estimate_utm_crs() when the source is geographic or unknown.
• TopologyException on buffer or intersection. A self-intersecting polygon crashes GEOS mid-operation. Pre-clean the column with gdf.geometry.make_valid() (Shapely 2.0 exposes it as a vectorized accessor), or guard the single-geometry call with make_valid().
• AttributeError: 'NoneType' object has no attribute 'area'. A row has a None geometry — common after a failed join or a bad read. Guard with if geom is None or geom.is_empty inside any .apply() function, as in step 3.
• Empty geometry returns nan, not an error. An empty (but non-null) geometry yields nan from .area, which then poisons sums. Filter with ~gdf.geometry.is_empty or coerce to 0.0 explicitly before aggregating.
• A hand-written lambda that GeoPandas already vectorizes. make_valid, segmentize, concave_hull, force_2d and others became GeoSeries accessors in 1.0. Check dir(gpd.GeoSeries) before writing the .apply(); the accessor is the same GEOS call without the per-row dispatch.
• The CRS vanished after a Shapely round trip. Shapely objects carry no coordinate system, so any column rebuilt from them comes back with crs=None and to_crs() then raises instead of reprojecting. Wrap the result: gpd.GeoSeries(values, crs=gdf.crs).
• shapely.* array functions on a GeoSeries return a bare ndarray. shapely.buffer(zones.geometry.values, 50) is valid and fast, but the result is a plain object array — assigning it directly to a column loses the CRS for the same reason as above.
• .apply() on a mixed-type column branches inconsistently. A GeometryCollection left behind by make_valid hits the polygon branch of your function and returns something the next step cannot use. Normalise types before the bridge, not inside it.

Frequently Asked Questions

If GeoPandas is built on Shapely, is there ever a reason to import Shapely directly? Three, reliably. Constructing geometry from raw coordinates before any table exists; calling anything in shapely.ops, which has no column accessor; and operating on geometry arrays that have no attributes attached, where a DataFrame would add an index and a schema you never read. Everything else is better expressed as a column operation.

How many features before plain Shapely lists stop being the right choice? The count matters less than whether the geometries have attributes. If each shape carries fields you will filter, group, or join on, a GeoDataFrame is correct at a hundred rows. If there are no attributes at all, raw arrays stay comfortable into the millions — the array API is the same engine GeoPandas calls.

Does converting a GeoDataFrame to a Shapely array make anything faster? No, and it is not really a conversion. gdf.geometry.values is a GeometryArray already wrapping the object array of Shapely geometries, so np.asarray() hands back the same pointers rather than copying geometry. You gain nothing on speed; you only lose the CRS and the index. Convert when you want to leave the table behind, not to optimise.

Which one belongs in a web service that handles one geometry per request? Shapely, with pyproj for projection. Importing GeoPandas pulls pandas, pyarrow, and pyogrio into every worker process, which costs both start-up time and resident memory per worker for a validate-buffer-serialise endpoint that never builds a table. The rule holds even if the service is written by people who use GeoPandas for everything else.

When does the choice stop being between these two at all? When the working set stops fitting in memory, or when the same query runs often enough to want an index maintained for it. Both point away from a process-local library and toward a query engine, as step 6 lays out. The signal is not the row count on its own — it is a load step that already consumes a large share of RAM before any analysis has run.