Shapely: geometric objects and operations

Introduction

Shapely is an open-source Python library designed for creating, analyzing, and manipulating planar geometric objects. It wraps the robust C-based GEOS computational geometry library, which also powers QGIS, PostGIS, and GeoPandas for spatial topological operations.

Every geometric record inside a GeoPandas GeoDataFrame is natively a Shapely object. Understanding Shapely clarifies how spatial buffering, intersections, or point-in-polygon queries are computed at the algorithmic level. The package is installed with pip install shapely and installs automatically alongside GeoPandas.

Geometry types

Shapely implements foundational OGC Simple Features specifications: Point, LineString, and Polygon. Their multi-part counterparts—MultiPoint, MultiLineString, and MultiPolygon—model complex geometries consisting of several disconnected components, such as island archipelagos or fragmented districts.

from shapely import Point, LineString, Polygon

quduq = Point(420350, 4575120)
quvur = LineString([(420000, 4575000), (420600, 4575300)])
uchastka = Polygon([(420200, 4575000), (420500, 4575000),
                    (420500, 4575250), (420200, 4575250)])

print(uchastka.area, quvur.length)   # 75000.0  670.82...

A polygon exterior ring must close properly; interior voids (such as a pond inside a land parcel) are supplied as explicit interior ring sequences.

Geometries can be formatted and exchanged via standard Well-Known Text (WKT) strings: print(quduq.wkt) outputs POINT (420350 4575120). This exact notation is recognized by PostGIS SQL queries and QGIS expression builders, making WKT a versatile format for inspecting geometries across spatial tools.

Basic Shapely geometry types (Point, LineString, Polygon) and their multi-part representations

Spatial predicates

A spatial predicate is a boolean evaluation testing topological relationships between two geometries:

Predicate	Evaluated condition
`intersects`	Do features share at least one point in common?
`contains` / `within`	Is one feature completely enclosed inside another?
`touches`	Do features contact only along common boundaries?
`disjoint`	Do features share zero points across space?

print(uchastka.contains(quduq))     # True
print(quvur.intersects(uchastka))   # True

Predicates serve as primary query filters: for example, identifying land parcels intersecting a designated riparian protection buffer.

Spatial relationships between two geometries: intersects, contains, touches, and disjoint

Geometric operations

While predicates test conditions, geometric operations construct new vector geometries:

intersection — extracts the shared overlapping footprint (intersection);
union — fuses geometries into a combined feature (union);
difference — subtracts the second geometry from the first (difference);
buffer(d) — expands an equidistant boundary at distance d (buffer).

from shapely import box

zona = quduq.buffer(100)                  # 100 m sanitary buffer
tegishli = uchastka.intersection(zona)    # portion of parcel inside buffer
print(round(tegishli.area))               # 31365

qoshni = box(420500, 4575000, 420700, 4575250)   # adjacent parcel
print(uchastka.touches(qoshni))                  # True
print(uchastka.union(qoshni).area)               # 125000.0
print(round(uchastka.difference(zona).area))     # 43635

A buffer is an approximating polygon rather than an analytical curve: by default, each quadrant is subdivided into 16 line segments. Consequently, computed area is 31,365 m² rather than theoretical π·100² = 31,416 m²; specifying quduq.buffer(100, quad_segs=64) yields greater smoothness.

Outputs should be inspected immediately: empty shapes (is_empty) or unexpected types (lines instead of polygons) highlight spatial calculation flaws.

Geometric operations on polygons: intersection, union, difference, and buffer around a point

Working with multiple features

Production datasets process thousands of features. Since Shapely 2.0 (released December 2022), functions are vectorized at the C level: calling shapely.contains(uchastkalar, quduq) evaluates arrays in compiled loops much faster than native Python iterations.

For large datasets, pairwise cross-comparisons remain computationally expensive. Analysts employ an STRtree spatial index, which organizes bounding boxes into a tree hierarchy, returning only plausible overlapping candidates during initial passes:

import shapely
from shapely import STRtree

uchastkalar = [uchastka, qoshni]
print(shapely.contains(uchastkalar, quduq))    # [ True False]

daraxt = STRtree(uchastkalar)
idx = daraxt.query(zona, predicate="intersects")
print(idx)                                     # [0]
tanlangan = [uchastkalar[i] for i in idx]

query returns integer indices referencing input lists rather than cloning geometries, allowing fast attribute table lookups.

This exact indexing engine powers GeoPandas sjoin. Grasping Shapely indexing mechanics helps engineers evaluate computational bottlenecks before deploying spatial pipelines.

CRS and invalid geometry

Shapely is CRS-agnostic: it evaluates purely planar numeric coordinates. If geometries are defined in decimal degrees, buffer(100) builds a 100-degree circle enclosing continents, calculating areas in square degrees. Therefore, project coordinates into metric projections—such as UTM zone 42N (EPSG:32642) for central Uzbekistan—using pyproj before buffering:

from pyproj import Transformer
from shapely.ops import transform

t = Transformer.from_crs("EPSG:4326", "EPSG:32642", always_xy=True)

nuqta_grad = Point(69.2401, 41.2995)            # Tashkent in degrees
xato = nuqta_grad.buffer(500)                   # radius of 500 degrees — invalid
print(xato.bounds[0])                           # -430.7599 (longitude less than -180)

nuqta_m = transform(t.transform, nuqta_grad)    # projected to meters
print(round(nuqta_m.buffer(500).area))          # 784137 m²

A second critical hazard is invalid geometry: self-intersecting boundary rings or duplicated vertices produce topological failures. Geometries are verified using is_valid, while make_valid from Shapely 2 documentation resolves topological defects:

from shapely import is_valid_reason, make_valid

kapalak = Polygon([(0, 0), (2, 2), (2, 0), (0, 2)])   # self-intersecting bowtie polygon
print(kapalak.is_valid)              # False
print(is_valid_reason(kapalak))      # Self-intersection[1 1]

tuzatilgan = make_valid(kapalak)
print(tuzatilgan.geom_type, tuzatilgan.area)   # MultiPolygon 2.0

The uncorrected bowtie polygon evaluates to zero area because counter-rotating triangles cancel out mathematically. Applying make_valid decomposes the figure into a valid MultiPolygon with correct total area.

Conclusion

Shapely provides the computational engine for geometric analysis in Python: instantiating vector primitives, testing topological predicates, and generating derived overlays.

Reliable spatial modeling requires two practices: computing distance and area metrics within metric projections, and validating geometric validity prior to running overlay operations.