PyPIDataupdated 06 Aug 2026

shapely

Item: shapely
Rating: 4.3
Author: Keyoor

Shapely does computational geometry on flat, two-dimensional shapes: points, lines, polygons and collections of them. It answers questions like does this polygon contain that point, what is the overlapping area of these two regions, what does this road look like buffered by fifty metres, and which of these ten thousand shapes touch each other. Underneath it wraps GEOS, the same C++ library that powers PostGIS, so the answers match what your spatial database gives you. Version 2 added a second way to call everything: alongside the object-oriented interface where you write polygon.contains(point), there are module-level NumPy ufuncs where you write shapely.contains(polygons_array, points_array) and get a boolean array back, with the loop running in C and the GIL released. That vectorized layer is what makes it practical to work with millions of geometries. Two things Shapely deliberately does not do: it has no idea what a coordinate reference system is, and it does not read or write data files.

shapely · v2.1.2 · 18.4M/wk · ★ 4,485 · pushed 05 Aug 2026 · shapely/shapely · docs · reviewed by Keyoor

Verdict

Shapely is the correct answer for planar geometry in Python and there is no real competition; GEOS behind it means the results agree with PostGIS. Just remember it knows nothing about the curvature of the earth or about coordinate systems, so pair it with pyproj or geopandas before you trust a distance.

API stability4/5The 2.0 release in 2022 was a deliberate, large break: geometries became immutable and hashable, iterating over a MultiPolygon stopped working in favour of .geoms, and the array interface, ctypes access and cascaded_union were removed. Since then 2.1 has been additive, with M coordinate support and new constructive functions, and the one signature change to shapely.transform is being phased in across releases with a deprecation warning first.

Docs4/5shapely.readthedocs.io documents every function with its GEOS version requirement, includes a proper 1.8 to 2.0 migration guide, and the manual explains the difference between the scalar Geometry interface and the vectorized ufuncs. What is under-signposted is the biggest trap: that all measurements are in input units and there is no CRS anywhere, which is a footnote rather than a banner.

Maintenance4/5The repo was pushed on 5 August 2026 with 2.2 in active development on main, and 200 open issues (236 counting PRs) on a project running since 2011. Release cadence is slow though: 2.1.2 came out on 24 September 2025 and is still the current release, so fixes sit on main for a long time before you can pip install them.

Ecosystem5/5Around 18M weekly downloads and it is the geometry type that the entire Python GIS stack agrees on: geopandas, pyogrio, rasterio, movingpandas and the rest all take and return Shapely objects. Wheels with GEOS bundled removed the historical install pain, and the __geo_interface__ protocol makes it interoperate with tools that do not depend on it directly.

Use it if

You need correct geometric predicates and set operations: intersects, contains, touches, intersection, difference, union, buffer, with GEOS behind them rather than your own point-in-polygon attempt
You have arrays of geometries and want operations to run in C across all of them at once, which is what the module-level ufuncs in Shapely 2 exist for
You need a spatial index: STRtree builds one over a geometry array and query returns integer indices, which turns an O(n*m) nested loop over candidate pairs into something that finishes
You are building on top of it: geopandas, pyogrio and most Python GIS tooling speak Shapely geometry objects, so producing them is how you interoperate

Skip it if

Your coordinates are longitude and latitude and you want real distances or areas: Shapely is planar and unit-agnostic, so distance between two lon/lat points comes back in degrees and area comes back in square degrees, both of which are meaningless. Project to a suitable metric CRS with pyproj first, or use pyproj.Geod for geodesic measurements
You need to read or write spatial files: there is no shapefile, GeoPackage, or file-level GeoJSON support at all, only WKT, WKB and per-geometry GeoJSON strings, so you still need pyogrio, fiona or geopandas for IO
Your geometries have attributes and belong in a table: geopandas wraps Shapely and adds the CRS tracking, IO, spatial joins and groupby that you would otherwise write by hand, and using bare Shapely means rebuilding that layer
Your dataset does not fit in memory: every geometry is a Python object holding a GEOS pointer, so tens of millions of polygons is a memory problem, and a spatial database or DuckDB with its spatial extension handles that shape better
You are still on Shapely 1.8: the 2.0 upgrade removed iteration over multi-part geometries, made geometries immutable, dropped the array interface, ctypes access and cascaded_union, so this is a code migration rather than a version bump
You need spherical geometry rather than projected: there are no great-circle or geodesic operations in Shapely, so anything crossing the antimeridian or spanning a continent will be wrong in ways that look plausible

Setup reality

pip install shapely is genuinely easy now: the 2.1.2 release ships 56 wheels covering the usual Linux, macOS and Windows targets with GEOS statically bundled inside, so you do not install GEOS yourself and the version you get is whatever the wheel was built against. NumPy 1.21 or newer and Python 3.10 or newer. The complication only appears when you mix sources. If you install Shapely from a wheel and something else, such as an older GDAL or a conda build of another geospatial package, links its own GEOS, you can end up with two copies of the library in the process, which produces segfaults rather than errors; the fix is to install the whole geospatial stack from one channel, conda-forge being the usual choice. Check shapely.geos_version at runtime rather than guessing, because several functions need a minimum GEOS version and simply are not available below it. Building from source needs a GEOS 3.10 or newer development install plus a C compiler, and is worth avoiding. Nothing here validates your input: Shapely will happily construct a self-intersecting polygon, and the failure shows up much later as a GEOSException from an unrelated operation.

Patterns

Build the basic shapescreate-geometries

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

p = Point(2.0, 3.0)
line = LineString([(0, 0), (1, 1), (2, 0)])
square = box(0, 0, 10, 10)                      # minx, miny, maxx, maxy
ring_with_hole = Polygon(
    shell=[(0, 0), (0, 10), (10, 10), (10, 0)],
    holes=[[(2, 2), (2, 4), (4, 4), (4, 2)]],
)
print(square.area, ring_with_hole.area)

Geometries have been immutable since 2.0, so there is no way to move a point after construction; build a new one. Polygon closes the ring for you, and a shell with fewer than four distinct points raises. Coordinates are plain floats with no units attached, which is the root of most Shapely bugs.

Ask whether two shapes relatespatial-predicates

square.contains(p)      # strictly inside, boundary excluded
square.covers(p)        # inside or on the boundary
square.intersects(line)
square.touches(line)    # boundaries meet, interiors do not
square.disjoint(line)

# full DE-9IM when the named predicates are not precise enough
square.relate(line)                 # e.g. '1010F0212'
square.relate_pattern(line, 'T*F**F***')

contains and covers differ only on the boundary, and that difference is the source of a lot of off-by-one-point bugs: a point exactly on a polygon edge is covered but not contained. Predicates are exact-arithmetic-free floating point, so two shapes that should share an edge often do not, which is what set_precision and snap are for.

Run an operation across a whole array of geometriesvectorized-operations

import numpy as np
import shapely

points = shapely.points(np.random.rand(1_000_000, 2) * 100)
region = shapely.box(10, 10, 20, 20)

inside = shapely.contains(region, points)      # bool array, loop in C
areas = shapely.area(shapely.buffer(points[:1000], 0.5))
merged = shapely.union_all(shapely.buffer(points[:1000], 0.5))

This is the Shapely 2 addition worth the upgrade: the module-level functions are true NumPy ufuncs, so they broadcast, accept and return arrays, and release the GIL. A Python for loop calling point.within(region) a million times is roughly an order of magnitude slower for the same answer. Note that shapely.points() takes an (N, 2) coordinate array, not a list of Point objects.

Find candidate pairs without a nested loopspatial-index

from shapely import STRtree

tree = STRtree(parcels)                     # array of polygons

# which parcels intersect each building?
building_idx, parcel_idx = tree.query(buildings, predicate="intersects")

# nearest parcel to each building, with distance
idx, dist = tree.query_nearest(buildings, return_distance=True)

query returns integer indices, not geometries, so index back into your arrays. Without predicate= you get bounding box hits only, which is a superset you must filter yourself; passing predicate makes GEOS do the exact test. The tree is built once and is read-only, so mutating the source array afterwards leaves the index describing shapes that no longer exist.

Speed up repeated tests against one complex shapeprepare-geometry

import shapely

shapely.prepare(coastline)          # mutates in place, builds an internal index
hits = shapely.contains(coastline, points)
shapely.destroy_prepared(coastline)

Preparing builds a spatial index of the geometry's edges, which turns a many-vertex polygon from slow to fast for repeated predicate calls. In 2.x this happens in place on the existing object rather than returning a wrapper, so shapely.is_prepared() is how you check. It only helps the geometry on the left of the predicate, and only for predicates, not for set operations.

Detect and fix self-intersectionsrepair-invalid-geometry

from shapely import is_valid, make_valid, set_precision
from shapely.validation import explain_validity

bowtie = Polygon([(0, 0), (2, 2), (2, 0), (0, 2)])
print(is_valid(bowtie))          # False
print(explain_validity(bowtie))  # 'Self-intersection[1 1]'

fixed = make_valid(bowtie)       # returns a MultiPolygon here
snapped = set_precision(fixed, 0.001)

make_valid can change the geometry type: a broken Polygon can come back as a MultiPolygon or a GeometryCollection, so code that assumes a Polygon breaks downstream. The old buffer(0) trick also fixes some cases but silently deletes parts of the shape. set_precision rounds coordinates to a grid, which is the reliable fix for shapes that should share edges but differ in the last floating point bits.

Move geometries in and out as WKT, WKB and GeoJSONread-write-formats

import shapely

wkt = shapely.to_wkt(square, rounding_precision=6)
wkb = shapely.to_wkb(square, hex=True)
gj = shapely.to_geojson(square)

back = shapely.from_wkt(wkt)
from_db = shapely.from_wkb(row["geom"])
lenient = shapely.from_wkt(bad_text, on_invalid="warn")

These are the only serialization formats; there is no file reader. to_wkt applies rounding_precision=6 by default, which silently loses precision on a round trip, so pass -1 to keep everything. from_wkb is what you use for PostGIS results, and the hex form is what psycopg hands you. on_invalid accepts raise, warn or ignore, and ignore yields None for the bad entries rather than aborting a bulk load.

Get real distances by projecting firstreproject-coordinates

from pyproj import Transformer
from shapely.ops import transform

# WGS84 lon/lat to Web Mercator metres
to_metres = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True)
projected = transform(to_metres.transform, lonlat_geom)
print(projected.length)   # metres, not degrees

This is the step people skip. buffer(0.001) on lon/lat coordinates makes a buffer of one thousandth of a degree, which is roughly 111 metres at the equator and roughly zero near the poles. always_xy=True on the Transformer matters, because EPSG:4326 officially orders coordinates as latitude then longitude while Shapely always means x then y. Web Mercator distorts area badly away from the equator; use a local UTM zone or an equal-area projection when the number matters.

Measure, simplify and densifymeasure-and-simplify

import shapely

shapely.length(road)
shapely.area(parcel)
shapely.distance(parcel, well)

small = shapely.simplify(coastline, tolerance=50, preserve_topology=True)
dense = shapely.segmentize(road, max_segment_length=10)
hull = shapely.concave_hull(points_multi, ratio=0.3)   # GEOS 3.11+

tolerance is in coordinate units, so the right value depends entirely on your projection. preserve_topology=False uses the faster Douglas-Peucker variant and can produce self-intersecting output. segmentize is what you run before reprojecting a long straight line, because a two-point line stays a two-point line through the transform and ends up in the wrong place on a curved projection.

Get at the pieces of a geometryaccess-parts-and-coordinates

for poly in multi.geoms:            # NOT: for poly in multi
    print(poly.exterior.coords[:])
    for hole in poly.interiors:
        print(len(hole.coords))

import shapely
coords = shapely.get_coordinates(geoms_array)     # (N, 2) float array
parts = shapely.get_parts(geoms_array)            # explode multiparts

Iterating a multi-part geometry directly raises TypeError since 2.0; you must go through .geoms. That single change is the most common reason 1.8 code fails on 2.x. get_coordinates flattens every geometry in the array into one big coordinate array and loses the boundaries between them, so pair it with return_index=True if you need to know which coordinate came from which shape.

Combine and subtract shapesoverlay-operations

import shapely

union = shapely.union_all(parcels)                # one call, not a fold
clipped = shapely.intersection(parcel, flood_zone)
remaining = shapely.difference(parcel, easement)
ring = shapely.buffer(road, 25).difference(shapely.buffer(road, 15))
edge = shapely.buffer(site, -5, join_style="mitre")  # negative shrinks

union_all does a cascaded union internally and is dramatically faster than reducing with union in a Python loop; the old shapely.ops.cascaded_union was removed in 2.0 in its favour. A negative buffer can make a narrow polygon disappear entirely and return an empty geometry, so check is_empty before using the result. Overlay results are frequently GeometryCollections when inputs touch at a point or along a line.

Work with positions along a linelinear-referencing

from shapely.ops import nearest_points, substring, split

d = road.project(incident)              # distance along the line
snapped = road.interpolate(d)           # point at that distance
half = road.interpolate(0.5, normalized=True)

a, b = nearest_points(road, incident)
segment = substring(road, 100, 250)
pieces = split(road, crossing_line)

project and interpolate are inverses of each other and are how you snap a GPS fix to a route. The distance is in coordinate units unless you pass normalized=True, which switches to a 0 to 1 fraction. split requires the splitter to actually intersect the line; a splitter that misses by a floating point hair returns the original geometry unchanged with no warning, so snap them together first.

Alternatives

Package	Registry	Pick it when
geopandas	PyPI	Your geometries have attribute columns and a coordinate reference system, and you want file IO and spatial joins without writing them yourself
pyproj	PyPI	You need coordinate reference system transforms or true geodesic distances and areas on the ellipsoid, neither of which Shapely does
duckdb	PyPI	The dataset is too large to hold as Python objects and you would rather do spatial joins in SQL against Parquet
pyogrio	PyPI	The job is reading and writing shapefiles, GeoPackages or FlatGeobuf, which Shapely does not touch