Shapely Geometry Operations: Predicates, Set Operations & Buffers

Shapely is the geometry engine of the Python geospatial stack — a thin, GEOS-backed layer that constructs, validates, and manipulates planar geometries one object or one array at a time. This stage of Mastering Core Geospatial Python Libraries covers the operations practitioners reach for daily: building Point, LineString, and Polygon objects, testing spatial predicates like intersects and contains, and computing buffers, unions, and clips. It supplies the topological primitives that GeoPandas DataFrames delegates to across whole tables, and it operates strictly in Cartesian space — so every metric result depends on first projecting with Coordinate Systems with PyProj.

Shapely geometry type hierarchy The three base geometry types — Point, LineString and Polygon — each have a multipart counterpart, all unified under GeometryCollection. Three primitives, each with a multipart form Point x, y (z) LineString vertices Polygon ring + holes MultiPoint MultiLineString MultiPolygon
Every Shapely object is one of these types — knowing the hierarchy explains why operations return Multi* and GeometryCollection results.

Architecture & Data Structures

Every Shapely geometry is an immutable Python wrapper around a GEOS C++ object. The three primitives — Point, LineString, and Polygon — each have a multipart counterpart (MultiPoint, MultiLineString, MultiPolygon), and any mix of them can be held in a GeometryCollection. Understanding this hierarchy is what explains otherwise surprising results: a Polygon.intersection() that clips through a hole returns a MultiPolygon, and a line-on-line intersection can return a GeometryCollection of both points and segments.

A Polygon is not a flat ring of coordinates — it is one exterior ring plus zero or more interior rings (holes), each stored as a closed LinearRing. Coordinates are exposed as immutable sequences; to "edit" a geometry you build a new one. Since Shapely 2.0 the objects are hashable and thread-safe, and construction is cheap because coordinates are stored in a contiguous NumPy buffer rather than per-vertex Python objects.

from shapely import Point, LineString, Polygon, MultiPolygon
from shapely import GeometryCollection

# A parcel with a courtyard hole: one exterior ring + one interior ring
parcel = Polygon(
    shell=[(0, 0), (40, 0), (40, 30), (0, 30)],
    holes=[[(10, 10), (20, 10), (20, 20), (10, 20)]],
)

print(parcel.geom_type)          # 'Polygon'
print(len(parcel.interiors))     # 1  -> the courtyard
print(parcel.exterior.is_ccw)    # ring orientation matters for validity

# A collection of unrelated geometries stays addressable by index
network = GeometryCollection([Point(5, 5), LineString([(0, 0), (40, 30)])])
print([g.geom_type for g in network.geoms])   # ['Point', 'LineString']

Prefer the top-level shapely.Point / shapely.Polygon imports (Shapely 2.0+) over the legacy shapely.geometry.Point path; both still work, but the flat namespace is the documented API and aligns with the vectorized functions in the next sections.

Two structural details decide how the rest of the API behaves. The first is that coordinates are reachable as a flat (n, 2) NumPy array regardless of how many rings or parts a geometry has, and shapely.transform() rebuilds the same structure around a modified array — which is how you move, scale, or reproject geometry without unpacking rings by hand. The second is that "no geometry" has two distinct spellings, and they do not behave alike:

import numpy as np
import shapely
from shapely import Polygon

parcel = Polygon(
    shell=[(0, 0), (40, 0), (40, 30), (0, 30)],
    holes=[[(10, 10), (20, 10), (20, 20), (10, 20)]],
)

coords = shapely.get_coordinates(parcel)     # (n, 2) — exterior then interior rings
print(coords.shape)                          # (10, 2)

# Structure-preserving edit: shift the parcel 5 m east, hole and all
moved = shapely.transform(parcel, lambda xy: xy + [5.0, 0.0])
print(len(moved.interiors))                  # 1 — the courtyard survived

# Missing (None) is not the same thing as empty (a valid geometry with no points)
layer = np.array([parcel, None, Polygon()], dtype=object)
print(shapely.is_missing(layer))             # [False  True False]
print(shapely.is_empty(layer))               # [False False  True]
print(shapely.area(layer))                   # [1100.   nan    0.]

The nan in that last line is the detail that costs people an afternoon: a missing geometry propagates nan through every measurement instead of raising, so a single None left by a failed join turns a total area into nan with no traceback. Filter with shapely.is_missing before any reduction, and decide deliberately whether an empty geometry should contribute 0.0 or be dropped.

Environment Configuration & Dependency Resolution

Shapely's behaviour is tied to the GEOS build it links against, so version pinning matters. Shapely 2.0 was a hard break from 1.x: geometries became immutable, the array-based (shapely.*) functions replaced the old shapely.vectorized and shapely.ops.cascaded_union helpers, and per-object attribute mutation was removed. Code written for 1.x can silently misbehave under 2.x — see Shapely 1.x vs Shapely 2.0 Vectorization for the migration specifics.

Install the whole native stack from a single channel. Mixing a pip Shapely wheel (which bundles its own GEOS) with a conda GEOS, or with a pip pyproj/fiona that bundle their own PROJ and GDAL, is the leading cause of DLL and symbol-clash errors.

# One channel, pinned — reproducible across machines
conda create -n geo python=3.12 \
    "shapely>=2.0" "pyproj>=3.6" "fiona>=1.9" "geopandas>=1.0" "numpy>=1.26" \
    -c conda-forge
conda activate geo

Verify the link at runtime before trusting any topology result — a mismatched or ancient GEOS explains bugs that look like data problems:

import shapely
print(shapely.__version__)        # e.g. 2.0.6
print(shapely.geos_version)       # (3, 12, 1)  -> the C engine actually in use
assert shapely.geos_version >= (3, 10, 0), "make_valid and coverage ops need GEOS 3.10+"

The GEOS version — not the Shapely version — gates which functions exist at all, which is why the same pip install shapely behaves differently on two machines. Shapely exposes a Python-level name for an operation, but calling it on an older engine raises UnsupportedGEOSVersionError rather than silently degrading. The thresholds worth knowing:

Probe rather than assume when a library has to run on someone else's environment:

import shapely

def hull_of(observations):
    """Tight hull where the engine supports it, convex hull as the floor."""
    if shapely.geos_version >= (3, 11, 0):
        return shapely.concave_hull(observations, ratio=0.3)
    return shapely.convex_hull(observations)

One more source of confusion: a pip wheel bundles a private GEOS build, so shapely.geos_version can legitimately disagree with whatever geos-config --version reports on the system, and upgrading the system GEOS changes nothing. Conda-forge builds link the shared library instead, which is why a single-channel environment is also the only one where the two numbers agree.

Vectorized Operations & Core Workflow

The defining feature of Shapely 2.0 is that the top-level functions are NumPy ufuncs: pass an array of geometries and the loop runs in C, not Python. This is what makes Shapely fast enough to sit under GeoPandas DataFrames, which is essentially an array of Shapely geometries plus a CRS. Working directly with arrays is the right call when you need geometry math without the DataFrame overhead — for example, filtering a stream of sensor readings against a study boundary.

The end-to-end pattern below buffers a set of point observations, tests which ones fall inside a floodplain boundary, and measures distances — all vectorized:

import numpy as np
import shapely
from shapely import Point, Polygon

# 5 sensor locations (already in a projected, metric CRS — see next section)
sensors = shapely.points(
    np.array([[312, 118], [340, 205], [455, 260], [500, 90], [610, 300]])
)

floodplain_boundary = Polygon([(300, 80), (520, 80), (520, 280), (300, 280)])

# Vectorized predicate: which sensors sit inside the floodplain?
inside = shapely.contains(floodplain_boundary, sensors)   # boolean array
print(inside)                                             # [ True  True False  True False]

# Vectorized measurement: distance from each sensor to the boundary edge
edge = floodplain_boundary.exterior
dist_to_edge = shapely.distance(edge, sensors)           # float array, metres
print(np.round(dist_to_edge, 1))

# Vectorized buffering: a 25 m alert radius around every sensor at once
alert_zones = shapely.buffer(sensors, 25.0)
merged_alert = shapely.union_all(alert_zones)            # single dissolved geometry
print(merged_alert.geom_type)                            # 'Polygon' or 'MultiPolygon'

Two rules keep this fast. First, reach for shapely.union_all(array) rather than folding .union() in a Python loop — the former hands the whole set to GEOS's cascaded-union sweep-line at once. Second, when you test one fixed reference geometry against many candidates, wrap it in shapely.prepared.prep() to build a spatial index over its edges (covered below).

Geometry & Data Processing Details

Topology validation

Invalid geometry is the silent killer of spatial pipelines. Self-intersecting rings, bowtie polygons, duplicate consecutive vertices, and wrong ring orientation cause boolean operations to fail or return wrong answers with no exception raised. Gate every incoming layer on .is_valid and repair with make_valid() before analysis — the same discipline the dedicated Topology Validation & Repair stage enforces across a full ingestion pipeline.

from shapely import Polygon, make_valid
from shapely.validation import explain_validity

# A bowtie: the ring crosses itself, so area and overlays are unreliable
bowtie = Polygon([(0, 0), (20, 20), (0, 20), (20, 0)])

if not bowtie.is_valid:
    print(explain_validity(bowtie))     # 'Self-intersection[10 10]'
    repaired = make_valid(bowtie)       # -> MultiPolygon of two clean triangles
    print(repaired.geom_type, round(repaired.area, 1))

Spatial predicates and the DE-9IM model

Predicates — intersects, contains, within, crosses, touches, overlaps, disjoint — are all evaluations of the DE-9IM intersection matrix under the hood. Two gotchas trip up newcomers: contains is strict about the boundary (a point exactly on a polygon edge is not contained but does intersect), and .distance() returns 0.0 for any geometries that touch or overlap. For tolerance-based proximity, use shapely.dwithin(a, b, distance) (GEOS 3.10+) instead of a.buffer(tol).intersects(b) — dwithin skips building the buffer geometry entirely and is markedly faster.

Predicate truth table for a point inside, on, and outside a polygon On the left, a floodplain boundary polygon carries three marked points: A sits in the interior, B sits exactly on the boundary, and C sits outside. On the right, a four-column table gives the result of intersects, contains, touches and disjoint for each point. Point A is true for intersects and contains only. Point B is true for intersects and touches but false for contains. Point C is true only for disjoint. Boundary cases: which predicate is true for a point on the edge? floodplain_boundary A B C A inside · B exactly on the boundary · C outside predicate(floodplain_boundary, point) point intersects contains touches disjoint A inside true true false false B on edge true false true false C outside false false false true contains() excludes the boundary; intersects() includes it.
Point B is the row that catches people out: it is on the polygon, so intersects and touches are true while contains is false.
from shapely import LineString, Point, Polygon
import shapely

pipeline = LineString([(0, 0), (500, 500)])
inspection_pt = Point(300, 100)

print(shapely.dwithin(pipeline, inspection_pt, 150.0))   # within 150 m? True/False
print(round(pipeline.distance(inspection_pt), 1))        # exact metric distance

# nearest_points resolves the actual closest coordinate pair
from shapely.ops import nearest_points
on_line, on_point = nearest_points(pipeline, inspection_pt)
print(on_line.x, on_line.y)

Set operations and generalization

intersection, union, difference, and symmetric_difference derive analytical boundaries — clipping features to a study area, cutting exclusion zones, or computing overlap. buffer, convex_hull, and minimum_rotated_rectangle generalize features for zoning or footprint analysis. Two practical warnings: never chain .buffer() calls (each pass adds vertices and topological noise — compute once with an explicit join_style and mitre_limit), and use simplify(tolerance, preserve_topology=True) rather than raw Douglas–Peucker when you must keep rings valid.

from shapely import Polygon
from shapely.affinity import rotate, translate

base = Polygon([(0, 0), (100, 0), (100, 100), (0, 100)])
overlap = Polygon([(60, 60), (160, 60), (160, 160), (60, 160)])

study_area = base.intersection(overlap)        # the shared region
exclusion = base.difference(study_area)        # base minus the overlap

# Single, explicit buffer — never chained
setback = base.buffer(5.0, join_style="mitre", mitre_limit=2.0)
footprint = rotate(translate(setback, xoff=10, yoff=-10), 30, origin="centroid")

Precision, snapping, and the overlay engine

Most TopologyException failures are not bad data in the human sense — they are two boundaries that were digitised to coincide but differ in the eleventh decimal place. GEOS has to node every edge intersection exactly before it can build an overlay result, and floating-point coordinates that are almost-but-not-quite equal produce intersection points that cannot be represented consistently. The result is either an exception or a hairline sliver polygon along what should have been a shared edge.

Since GEOS 3.9 the fix is a fixed-precision overlay: pass grid_size and GEOS snaps both operands to that grid inside the operation, leaving the inputs untouched.

import shapely
from shapely import Polygon

# Two parcels that share an edge — except for 0.1 nanometres of float noise
west = Polygon([(0, 0), (100, 0), (100, 100), (0, 100)])
east = Polygon([(100.0000000001, 0), (200, 0), (200, 100), (100.0000000001, 100)])

noisy = shapely.union(west, east)
print(len(shapely.get_parts(noisy)))            # 1, but with a hairline interior ring

# Fixed precision: snap to a 1 mm grid during the overlay only
clean = shapely.union(west, east, grid_size=0.001)
print(round(clean.area, 3), len(clean.interiors))   # 20000.0 0

Choose grid_size from the survey accuracy of the data, not from how small you can make it: 1 mm for cadastral parcels, 0.1 m for digitised imagery, 1 m for consumer GPS traces. A grid finer than the noise you are trying to remove removes nothing. The alternative, shapely.set_precision(geom, grid_size), rewrites the coordinates permanently — useful as an ingestion-time normalisation, dangerous mid-analysis because collapsing a grid can shorten edges below the grid size and delete them.

Repair has two strategies as of Shapely 2.1. The default make_valid(geom, method="linework") — the only behaviour available in 2.0 — preserves every input edge, which means an invalid polygon can legitimately come back as a GeometryCollection containing both polygons and dangling lines. method="structure" rebuilds the areal interior instead and returns polygonal output only, and with keep_collapsed=False it discards zero-area shards. When the repaired geometry has to stay in a polygon layer, structure is the one that will not surprise a later to_file() call with a mixed-type collection.

Linear referencing and geometry surgery

Network and utility work is mostly positions along a line rather than positions in the plane, and Shapely handles that with a pair of inverse operations: project() maps a point to its distance along a line, interpolate() maps a distance back to a point. Everything else — cutting a maintenance reach, splitting a run at a junction, rebuilding a route from noded segments — is built on those two.

from shapely import LineString, Point
from shapely.ops import substring, split, snap

# A pipeline in EPSG:25832 metres: 400 m east, then 900 m north
pipeline = LineString([(0, 0), (400, 0), (400, 900)])
print(round(pipeline.length, 1))                 # 1300.0

leak = Point(400, 310)
chainage = pipeline.project(leak)                # 710.0 m along the run
print(chainage, pipeline.interpolate(chainage).wkt)   # POINT (400 310)

# The 50 m of pipe either side of the leak — the repair order's extent
repair_reach = substring(pipeline, chainage - 50, chainage + 50)
print(round(repair_reach.length, 1))             # 100.0

# Splitting needs the cut point to lie exactly ON the line
survey_pt = Point(400.0004, 310)                 # 0.4 mm off — split does nothing
on_line = snap(survey_pt, pipeline, tolerance=0.01)
upstream, downstream = split(pipeline, on_line).geoms
print(round(upstream.length, 1), round(downstream.length, 1))   # 710.0 590.0

The survey_pt line is the failure practitioners actually hit: split() cuts only where the splitter genuinely touches the line, so a point that misses by a fraction of a millimetre returns the original geometry as a single-part collection and the pipeline silently never gets divided. Snap it on first, or derive the cut point with interpolate() so it is exact by construction. project() also takes normalized=True to work in 0–1 fractions instead of map units, which is the safer form when the same code runs against lines of very different lengths. All of it is planar distance, so the line must already be in a metric CRS before a chainage means anything.

Prepared geometries for high-volume filtering

When one static reference geometry is tested against many candidates, repeated predicate calls scale O(n) with full edge comparisons each time. shapely.prepared.prep() builds a cached edge index once, accelerating repeated contains/intersects by 10–100×.

from shapely.prepared import prep

reference_zone = Polygon([(0, 0), (1000, 0), (1000, 1000), (0, 1000)])
fast_zone = prep(reference_zone)               # index built once

hits = [pt for pt in sensor_points if fast_zone.contains(pt)]

STRtree for many-against-many queries

prep() solves one shape against many. The opposite problem — many shapes against many shapes — needs an index over the whole candidate set, and that is shapely.STRtree: a packed Sorted-Tile-Recursive R-tree built once over an array of geometries. It is the same structure GeoPandas exposes as sindex, and reaching for it directly is the right move when you have geometry arrays but no attribute table.

import numpy as np
import shapely
from shapely import STRtree

# 40k parcel polygons and 5k sensor points, all in EPSG:25832 metres
tree = STRtree(parcels)                       # O(n log n), built once, then immutable

# Which parcel does each sensor fall in? Returns index pairs, not geometries.
sensor_idx, parcel_idx = tree.query(sensors, predicate="within")
print(sensor_idx[:3], parcel_idx[:3])         # [0 1 2] [118 118 4402]

matched_parcels = parcels[parcel_idx]         # index back into the source array

# Nearest parcel within 250 m, distance included
near_idx, distance = tree.query_nearest(
    sensors, max_distance=250.0, return_distance=True
)

Three behaviours are worth internalising. Without predicate=, query() returns envelope hits — a cheap candidate set that still needs an exact test, which is exactly the two-phase filter-and-refine that makes large overlays tractable. With predicate=, GEOS runs the exact relate on those candidates for you, so the returned pairs are final. And the result has an inner-join shape: a sensor that matches nothing simply does not appear in the output arrays, so reconstructing a per-sensor answer means seeding an output array with -1 and scattering the hits into it, rather than assuming the result lines up positionally with the input.

The tree is a snapshot. Mutating or replacing entries in the array you built it from leaves the index pointing at the old geometries, with no error and quietly wrong answers — rebuild after any repair pass. Build cost is real but amortises fast: below a few thousand comparisons a brute-force loop wins, and above that the tree wins by orders of magnitude. Once the same query needs to travel with attributes and joins, this is the machinery underneath GeoPandas sjoin, and the tabular form is documented in Spatial Joins & Merging.

Dimensionality: what silently drops the z coordinate

Shapely stores a third ordinate when you give it one, but GEOS is a 2D engine and most operations ignore or discard z. That is fine until an elevation-bearing dataset passes through a buffer and comes back flat.

import shapely
from shapely import LineString

# A sewer run with invert levels: 300 m horizontally, 2.6 m of fall
sewer_run = LineString([(0, 0, 12.4), (300, 0, 9.8)])
print(sewer_run.has_z, round(sewer_run.length, 1))    # True 300.0  <- planar length

corridor = shapely.buffer(sewer_run, 5.0)
print(corridor.has_z)                                  # False — buffer drops z

print(shapely.get_coordinates(sewer_run, include_z=True))
# [[  0.    0.   12.4]
#  [300.    0.    9.8]]

.length is the horizontal run, not the slope distance — for the true 3D length compute it from get_coordinates(..., include_z=True) yourself. Predicates are 2D as well, which means a road and the tunnel beneath it intersects as True no matter how far apart their z values are; grade-separated networks need an explicit level attribute, never a reliance on elevation. Buffer, overlay, and convex hull return 2D results, so put any z-dependent calculation before those steps. On the storage side, WKB and GeoParquet carry z faithfully, and shapely.force_2d() is the explicit way to strip it when a consumer cannot handle three ordinates.

CRS Alignment & Projection Pipeline

Shapely has no concept of a coordinate reference system — it treats every coordinate as a plain Cartesian number. Feed it EPSG:4326 longitude/latitude and .area returns square degrees, .distance() returns degrees, and a buffer(100) inflates by 100 degrees. These are mathematically valid and geographically meaningless. Always project to an appropriate metric CRS before any area, length, or buffer operation, using Coordinate Systems with PyProj.

Two PyProj gotchas recur constantly. Set always_xy=True so the transformer accepts and returns (longitude, latitude) / (easting, northing) order rather than the authority's declared axis order — otherwise your x and y silently swap. And never use Web Mercator (EPSG:3857) for measurement: it distorts area and distance badly away from the equator. Pick the local UTM zone (or a national grid) instead.

from pyproj import Transformer
from shapely import Point
from shapely.ops import transform

# WGS84 lon/lat -> UTM zone 33N (metres). always_xy keeps (lon, lat) order.
to_utm = Transformer.from_crs("EPSG:4326", "EPSG:32633", always_xy=True)

gauge_wgs84 = Point(13.404, 52.520)            # (lon, lat), Berlin
gauge_utm = transform(to_utm.transform, gauge_wgs84)

# NOW metric operations are meaningful
catchment = gauge_utm.buffer(500.0)            # a true 500 m radius
print(round(catchment.area))                   # ~785398 m^2 (pi * 500^2)

For arrays, transform the coordinate columns with the vectorized Transformer.transform(xs, ys) and rebuild geometries with shapely.points() / shapely.set_coordinates() rather than looping — the same principle behind GeoPandas' to_crs(). If a layer arrives without a declared CRS, you cannot guess it: confirm the source before projecting, because relabeling and reprojecting are different operations with very different results.

Production Export & Integration

Shapely produces geometry; getting it into a durable, interoperable form is a separate concern. Serialize with the built-in WKB/WKT and GeoJSON functions, and stream vector I/O through fiona so you never hold an entire dataset in memory.

import fiona
import shapely
from shapely.geometry import shape, mapping

# Round-trip a geometry through the interchange formats
wkb_blob = shapely.to_wkb(catchment)           # compact binary, ideal for PostGIS/Parquet
geojson_geom = mapping(catchment)              # dict, ready for web maps

# Stream features, filter, and write — constant memory
with fiona.open("gauges.gpkg", layer="stations") as src:
    profile = src.profile
    with fiona.open("catchments.gpkg", "w", **profile) as dst:
        for feat in src:
            geom = shape(feat["geometry"])
            if geom.is_valid:
                dst.write({"geometry": mapping(geom.buffer(500.0)),
                           "properties": feat["properties"]})

For analytical storage, hand geometries to a GeoDataFrame and write GeoParquet — it preserves the CRS in file metadata and is far faster than Shapefile. For a spatial database, load WKB into a PostGIS geometry column (planar, SRID-tagged) rather than geography unless you specifically need ellipsoidal math. When a Shapely-heavy workflow starts iterating over thousands of features to attach attributes or run joins, that is the signal to move up to the DataFrame layer — the trade-offs are laid out in Shapely vs GeoPandas: When to Use Each.

WKB has two flavours and picking the wrong one costs an SRID. to_wkb(flavor="iso") writes the standards-conformant form, where Z and M are encoded in the geometry type code and there is no room for a coordinate system. to_wkb(flavor="extended", srid=...) writes PostGIS EWKB, which sets high bits on the type word to flag Z, M, and an embedded SRID:

import shapely

# Interchange: ISO WKB, CRS lives in the container's metadata (GeoParquet, GeoPackage)
iso_blob = shapely.to_wkb(catchment, flavor="iso")

# PostGIS: extended WKB carries SRID 25832 with the geometry itself
ewkb_hex = shapely.to_wkb(catchment, flavor="extended", srid=25832, hex=True)

# Text output for logs and diffs — trim the digits nobody can act on
print(shapely.to_wkt(catchment, rounding_precision=3, trim=True)[:60])

Hand EWKB to a PostGIS geometry column and the SRID travels with the value; hand it ISO WKB and you must assert the SRID in SQL (ST_SetSRID) or the column silently accepts SRID 0. GeoParquet goes the other way — it expects ISO WKB with the CRS declared once in file metadata, so embedding an SRID per row is redundant and non-conformant. rounding_precision on to_wkt is a text-formatting control only; it does not change the geometry, which is what distinguishes it from set_precision.

Production checklist

Windows / Platform Edge Cases & Debugging

The stack behaves identically across platforms until the native GEOS/PROJ/GDAL libraries fail to resolve — and on Windows they fail more often. The recurring symptoms and their one-line fixes:

Decision flow for diagnosing a wrong Shapely result Check three things in order: is the geometry valid, is it in a metric CRS, and is the predicate or distance semantics correct. Each failing check points to its fix — make_valid, reproject with PyProj, or use intersects and dwithin — before the result can be trusted. Wrong result? Check three things, in order 1 is_valid == True? self-intersections, bowties 2 In a metric CRS? not raw EPSG:4326 degrees 3 Correct predicate? boundary & touch cases Result you can trust yes yes no no no make_valid() snap noise with set_precision() Reproject to a metric CRS PyProj · always_xy · UTM, not EPSG:3857 Fix the semantics on-edge → intersects · touching → dwithin()
Most "wrong" Shapely answers are one of three things: invalid input, a non-metric CRS, or a predicate/distance semantic that does not mean what you assumed.

Frequently Asked Questions

Does Shapely handle coordinate reference systems? No. Shapely operates purely in Cartesian space and ignores any CRS. Project with PyProj (or via GeoPandas' to_crs()) before measuring, then hand the projected geometry to Shapely.

Why does .distance() return 0 for shapes that clearly differ? Because they intersect or touch — planar distance between overlapping geometries is zero by definition. Use shapely.dwithin() for a tolerance test, or nearest_points() to inspect the closest pair.

When should I switch from Shapely to GeoPandas? Stay in Shapely for a handful of geometries or custom per-shape logic. Move to GeoPandas once you have hundreds-plus features with attributes, joins, or batch reprojections — see the dedicated comparison.

How do I fix a TopologyException during a union? An invalid geometry reached GEOS. Run make_valid() on the inputs first, and if precision noise persists, snap with shapely.set_precision() before retrying.

Is shapely.geometry.Point still correct in 2.0? It still works, but the documented API is the flat shapely.Point. Prefer the top-level imports, which pair naturally with the vectorized shapely.* array functions.

prep() or STRtree — which one do I actually need? Count the geometries on each side. One fixed reference against many candidates is prep(), which caches an edge index for that single shape. Many against many is STRtree, which indexes the whole candidate set so each query prunes to a handful of envelopes first. Using prep() inside a loop over both sides rebuilds the cache every iteration and is slower than doing nothing.

Can I use threads to speed up a large shapely.* call? No — a vectorized call is one single-threaded C loop, so it saturates exactly one core no matter how many you have. Geometries are immutable and safe to read from multiple threads, but throughput comes from splitting the array across processes or partitions instead; the pattern for doing that on tabular data is in Scaling with Dask-GeoPandas.

Why did union_all() return a GeometryCollection instead of a MultiPolygon? Because the input array was not all polygons. A union preserves whatever dimensions it was given, so one stray LineString or Point — often a degenerate polygon that make_valid reduced to a line — forces a mixed collection. Filter the array by shapely.get_type_id or repair with make_valid(method="structure") before the reduction if the output must stay polygonal.

Related