Shapely (shapely/shapely)

Shapely

https://github.com/shapely/shapely
Admin
Manipulation and analysis of geometric objects

Tokens:33,173
Snippets:181
Trust Score:5.7
Update:5 months ago
Show doc for...
Context Summary (auto-generated)
Raw
# Shapely - Geometric Object Manipulation and Analysis

Shapely is a BSD-licensed Python package for manipulation and analysis of planar geometric objects in the Cartesian plane. It wraps the widely deployed GEOS (Geometry Engine Open Source) library, which is also the engine of PostGIS and a port of the Java Topology Suite (JTS). Shapely provides both a feature-rich Geometry interface for scalar geometries and high-performance NumPy ufuncs for vectorized operations on arrays of geometries.

The library focuses exclusively on geometric operations without handling data serialization formats, map projections, or file I/O, making it lightweight and easy to integrate with other spatial data tools. Shapely enables PostGIS-type geometry operations outside of a database, supports multithreading by releasing Python's GIL during GEOS operations, and follows the OGC Simple Features specification. It requires Python ≥3.11, GEOS ≥3.10, and NumPy ≥1.23.

## Point - Zero-dimensional geometry creation and operations

Create and manipulate point geometries representing single coordinate locations in 2D or 3D space.

```python
from shapely import Point
import shapely

# Create 2D and 3D points
point_2d = Point(0.0, 0.0)
point_3d = Point(1.5, 2.5, 3.5)

# Access coordinates
print(point_2d.x, point_2d.y)  # 0.0 0.0
print(point_3d.x, point_3d.y, point_3d.z)  # 1.5 2.5 3.5

# Check for Z coordinates
print(shapely.has_z(point_2d))  # False
print(shapely.has_z(point_3d))  # True

# Create buffer (circular area) around point
circle = point_2d.buffer(10.0)
print(circle.area)  # 313.6548490545941
print(circle.geom_type)  # 'Polygon'

# Vectorized point creation
import numpy as np
points_array = shapely.points([[0, 0], [1, 1], [2, 2], [3, 3]])
print(points_array)  # array([<POINT (0 0)>, <POINT (1 1)>, ...])
```

## LineString - One-dimensional curve geometry

Create and manipulate line geometries with measurement and interpolation capabilities.

```python
from shapely import LineString
import shapely

# Create linestring from coordinates
line = LineString([(0, 0), (1, 1), (2, 1), (3, 0)])

# Basic properties
print(line.length)  # 4.242640687119285
print(line.bounds)  # (0.0, 0.0, 3.0, 1.0)
print(line.is_simple)  # True (no self-intersections)
print(line.is_closed)  # False

# Linear referencing - interpolate point at distance
point_at_2 = shapely.line_interpolate_point(line, 2.0)
print(point_at_2)  # POINT (1.414... 1.0)

# Find distance along line to a point
point = shapely.Point(2, 1)
distance = shapely.line_locate_point(line, point)
print(distance)  # 2.414...

# Normalized (fractional) interpolation
midpoint = shapely.line_interpolate_point(line, 0.5, normalized=True)
print(midpoint)  # Point at 50% along the line

# Create parallel offset
offset = line.parallel_offset(0.5, 'left')
print(offset.geom_type)  # 'LineString'
```

## Polygon - Two-dimensional surface with holes

Create and manipulate polygon geometries with interior and exterior rings.

```python
from shapely import Polygon, Point
import shapely

# Create simple polygon
polygon = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)])
print(polygon.area)  # 4.0
print(polygon.length)  # 8.0 (perimeter)

# Create polygon with hole
exterior = [(0, 0), (4, 0), (4, 4), (0, 4)]
hole = [(1, 1), (3, 1), (3, 3), (1, 3)]
polygon_with_hole = Polygon(exterior, [hole])
print(polygon_with_hole.area)  # 12.0 (16 - 4)

# Access exterior and interior rings
print(polygon_with_hole.exterior.coords[:])  # [(0, 0), (4, 0), ...]
print(len(polygon_with_hole.interiors))  # 1

# Test point containment
point_inside = Point(2, 2)
point_in_hole = Point(2, 2)
point_outside = Point(5, 5)

print(shapely.contains(polygon, point_inside))  # True
print(shapely.contains(polygon_with_hole, point_in_hole))  # False (in hole)
print(shapely.contains(polygon, point_outside))  # False

# Centroid and bounds
print(polygon.centroid)  # POINT (1 1)
minx, miny, maxx, maxy = polygon.bounds
print(f"Bounds: ({minx}, {miny}) to ({maxx}, {maxy})")
```

## Box - Rectangular polygon creation

Create axis-aligned rectangular polygons efficiently.

```python
import shapely

# Create box from min/max coordinates
box = shapely.box(0, 0, 10, 5)
print(box.area)  # 50.0
print(box.bounds)  # (0.0, 0.0, 10.0, 5.0)

# Counter-clockwise ordering
box_ccw = shapely.box(0, 0, 5, 5, ccw=True)
print(shapely.is_ccw(box_ccw.exterior))  # True

# Vectorized box creation for spatial queries
import numpy as np
boxes = [shapely.box(i, i, i+1, i+1) for i in range(5)]
print(len(boxes))  # 5 boxes
```

## Spatial Predicates - Boolean geometry relationships

Test spatial relationships between geometries with vectorized operations.

```python
import shapely
from shapely import Point, Polygon, LineString
import numpy as np

# Create test geometries
polygon = shapely.box(0, 0, 5, 5)
point_inside = Point(2, 2)
point_outside = Point(7, 7)
point_on_edge = Point(0, 2)

# Binary predicates
print(shapely.contains(polygon, point_inside))  # True
print(shapely.contains(polygon, point_outside))  # False
print(shapely.intersects(polygon, point_on_edge))  # True
print(shapely.within(point_inside, polygon))  # True
print(shapely.touches(polygon, point_on_edge))  # True

# Vectorized predicate operations
points = shapely.points([[1, 1], [3, 3], [6, 6], [2.5, 2.5]])
contained = shapely.contains(polygon, points)
print(contained)  # array([ True,  True, False,  True])

# More predicates
line = LineString([(0, 0), (6, 6)])
print(shapely.crosses(line, polygon))  # True
print(shapely.overlaps(polygon, shapely.box(2, 2, 7, 7)))  # True
print(shapely.disjoint(polygon, shapely.box(10, 10, 15, 15)))  # True

# Unary predicates
complex_line = LineString([(0, 0), (1, 1), (1, 0), (0, 1)])
print(shapely.is_simple(complex_line))  # False (self-intersecting)
print(shapely.is_valid(polygon))  # True
print(shapely.is_empty(Point()))  # True
```

## Set Operations - Union, intersection, difference

Perform set-theoretic operations on geometries.

```python
import shapely
from shapely import Polygon

# Create overlapping polygons
poly1 = shapely.box(0, 0, 3, 3)
poly2 = shapely.box(2, 2, 5, 5)

# Basic set operations
union = shapely.union(poly1, poly2)
print(union.area)  # 21.0 (9 + 9 - 4 overlap)

intersection = shapely.intersection(poly1, poly2)
print(intersection.area)  # 1.0 (overlap area)

difference = shapely.difference(poly1, poly2)
print(difference.area)  # 8.0 (poly1 minus overlap)

sym_diff = shapely.symmetric_difference(poly1, poly2)
print(sym_diff.area)  # 17.0 (union minus intersection)

# Using operators (object-oriented API)
union_op = poly1 | poly2  # equivalent to union
intersection_op = poly1 & poly2  # equivalent to intersection
difference_op = poly1 - poly2  # equivalent to difference
sym_diff_op = poly1 ^ poly2  # equivalent to symmetric_difference

# Aggregating operations on multiple geometries
polygons = [shapely.box(i, i, i+2, i+2) for i in range(5)]
merged = shapely.union_all(polygons)
print(merged.geom_type)  # Could be 'Polygon' or 'MultiPolygon'
```

## Buffer - Create offset regions around geometries

Generate buffer zones with configurable styles and resolution.

```python
import shapely
from shapely import Point, LineString, Polygon

# Point buffer (circle)
point = Point(0, 0)
circle = shapely.buffer(point, 5.0)
print(circle.area)  # ~78.54 (π * 5²)

# Line buffer with different cap styles
line = LineString([(0, 0), (10, 0)])
buffer_round = shapely.buffer(line, 2.0, cap_style='round')
buffer_flat = shapely.buffer(line, 2.0, cap_style='flat')
buffer_square = shapely.buffer(line, 2.0, cap_style='square')

print(buffer_round.area)  # Larger (includes rounded ends)
print(buffer_flat.area)  # Smaller (flat ends)

# Join styles for corners
zigzag = LineString([(0, 0), (5, 5), (10, 0)])
buffer_round_join = shapely.buffer(zigzag, 1.0, join_style='round')
buffer_mitre_join = shapely.buffer(zigzag, 1.0, join_style='mitre')
buffer_bevel_join = shapely.buffer(zigzag, 1.0, join_style='bevel')

# Negative buffer (erosion)
polygon = shapely.box(0, 0, 10, 10)
eroded = shapely.buffer(polygon, -1.0)
print(eroded.bounds)  # (1.0, 1.0, 9.0, 9.0)

# High resolution buffer
high_res_circle = shapely.buffer(point, 5.0, quad_segs=32)
print(len(high_res_circle.exterior.coords))  # More points
```

## Simplify - Reduce geometry complexity

Simplify geometries using Douglas-Peucker algorithm.

```python
import shapely
from shapely import LineString, Polygon

# Create complex line
coords = [(i, i % 3) for i in range(100)]
complex_line = LineString(coords)
print(len(complex_line.coords))  # 100 points

# Simplify with tolerance
simple_line = shapely.simplify(complex_line, tolerance=1.0)
print(len(simple_line.coords))  # Much fewer points

# Topology preservation
polygon = Polygon([
    (0, 0), (0, 10), (5, 5), (10, 10), (10, 0)
])
simplified_preserve = shapely.simplify(polygon, tolerance=2.0, preserve_topology=True)
simplified_no_preserve = shapely.simplify(polygon, tolerance=2.0, preserve_topology=False)

print(shapely.is_valid(simplified_preserve))  # True
print(simplified_preserve.area <= polygon.area)  # True

# Vectorized simplification
import numpy as np
lines = [LineString([(i+j, j) for j in range(50)]) for i in range(5)]
simplified_lines = shapely.simplify(lines, tolerance=1.0)
print([len(line.coords) for line in simplified_lines])
```

## Convex Hull and Envelopes

Generate bounding shapes around geometries.

```python
import shapely
from shapely import MultiPoint, Point

# Convex hull from points
points = MultiPoint([
    Point(0, 0), Point(2, 2), Point(3, 1), Point(1, 3), Point(1, 1)
])
hull = shapely.convex_hull(points)
print(hull.geom_type)  # 'Polygon'
print(hull.area > points.bounds[2] * points.bounds[3] / 2)  # True

# Envelope (bounding box)
from shapely import LineString
line = LineString([(1, 1), (5, 3), (2, 7)])
envelope = shapely.envelope(line)
print(envelope.bounds)  # (1.0, 1.0, 5.0, 7.0)

# Minimum rotated rectangle
scattered_points = MultiPoint([Point(i, i**2) for i in range(10)])
min_rect = shapely.minimum_rotated_rectangle(scattered_points)
print(min_rect.area < shapely.envelope(scattered_points).area)  # Usually true

# Minimum bounding circle
circle_bounds = shapely.minimum_bounding_circle(scattered_points)
print(circle_bounds.geom_type)  # 'Polygon' (approximation of circle)
```

## Distance Calculations

Measure distances between geometries.

```python
import shapely
from shapely import Point, LineString, Polygon

# Point to point distance
p1 = Point(0, 0)
p2 = Point(3, 4)
print(shapely.distance(p1, p2))  # 5.0

# Point to line distance
line = LineString([(0, 0), (10, 0)])
point = Point(5, 5)
print(shapely.distance(line, point))  # 5.0

# Find nearest points
nearest = shapely.shortest_line(line, point)
print(nearest.coords[:])  # [(5.0, 0.0), (5.0, 5.0)]
print(nearest.length)  # 5.0

# Polygon to polygon distance
poly1 = shapely.box(0, 0, 1, 1)
poly2 = shapely.box(3, 3, 4, 4)
print(shapely.distance(poly1, poly2))  # ~2.828 (diagonal distance)

# Hausdorff distance (maximum distance)
line1 = LineString([(0, 0), (1, 0)])
line2 = LineString([(0, 1), (1, 1)])
print(shapely.hausdorff_distance(line1, line2))  # 1.0

# Within distance check
print(shapely.dwithin(poly1, poly2, 3.0))  # True (distance < 3.0)
print(shapely.dwithin(poly1, poly2, 2.0))  # False (distance > 2.0)

# Vectorized distance calculations
import numpy as np
points = shapely.points([[i, 0] for i in range(5)])
distances = shapely.distance(Point(0, 0), points)
print(distances)  # array([0., 1., 2., 3., 4.])
```

## Coordinate Operations

Extract, modify, and transform geometry coordinates.

```python
import shapely
from shapely import LineString, Polygon
import numpy as np

# Extract coordinates as numpy array
line = LineString([(0, 0), (1, 1), (2, 0)])
coords = shapely.get_coordinates(line)
print(coords)
# [[0. 0.]
#  [1. 1.]
#  [2. 0.]]

# Set new coordinates
new_coords = coords * 2  # Scale by 2
new_line = shapely.set_coordinates(line, new_coords)
print(new_line.coords[:])  # [(0.0, 0.0), (2.0, 2.0), (4.0, 0.0)]

# Count total coordinates
polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
print(shapely.count_coordinates(polygon))  # 5 (including closing point)

# Custom transformation function
def rotate_90(coords):
    """Rotate coordinates 90 degrees counter-clockwise."""
    return np.column_stack([-coords[:, 1], coords[:, 0]])

rotated_line = shapely.transform(line, rotate_90)
print(shapely.get_coordinates(rotated_line))
# [[ 0.  0.]
#  [-1.  1.]
#  [ 0.  2.]]

# Transform with external projection (example with scale)
def wgs84_to_mercator_approx(coords):
    """Simplified Web Mercator projection."""
    R = 6378137  # Earth radius in meters
    x = coords[:, 0] * R * np.pi / 180
    y = np.log(np.tan((90 + coords[:, 1]) * np.pi / 360)) * R
    return np.column_stack([x, y])

# Apply to geometry
lat_lon_point = shapely.Point(-122.4, 37.8)  # San Francisco
mercator_point = shapely.transform(lat_lon_point, wgs84_to_mercator_approx)
```

## Affine Transformations

Apply geometric transformations to geometries.

```python
from shapely import affinity, Point, Polygon
import shapely

# Create geometry
polygon = Polygon([(0, 0), (2, 0), (2, 1), (0, 1)])

# Translation
translated = affinity.translate(polygon, xoff=5, yoff=3)
print(translated.bounds)  # (5.0, 3.0, 7.0, 4.0)

# Rotation around origin
rotated_origin = affinity.rotate(polygon, 45, origin='center')
rotated_point = affinity.rotate(polygon, 90, origin=(0, 0))
rotated_radians = affinity.rotate(polygon, 1.5708, use_radians=True)

# Scaling
scaled = affinity.scale(polygon, xfact=2.0, yfact=3.0)
print(scaled.bounds)  # (0.0, 0.0, 4.0, 3.0)

# Scale from specific origin
scaled_centered = affinity.scale(polygon, xfact=0.5, yfact=0.5, origin='center')

# Skewing
skewed_x = affinity.skew(polygon, xs=30, ys=0)  # Shear along x-axis
skewed_y = affinity.skew(polygon, xs=0, ys=30)  # Shear along y-axis

# General affine transformation (matrix form)
# [x', y'] = [a b xoff] [x]
#            [d e yoff] [y]
#                       [1]
# Scale by 2, rotate, and translate
matrix = [2, 0, 0, 2, 10, 10]  # [a, b, d, e, xoff, yoff]
transformed = affinity.affine_transform(polygon, matrix)
print(transformed.bounds)  # (10.0, 10.0, 14.0, 12.0)
```

## STRtree - Spatial Indexing

Perform efficient spatial queries using R-tree indexing.

```python
import shapely
from shapely import STRtree, Point, box
import numpy as np

# Create spatial index from geometries
geometries = [box(i, i, i+1, i+1) for i in range(100)]
tree = STRtree(geometries)

# Bounding box query (fast)
query_box = box(5, 5, 10, 10)
indices = tree.query(query_box)
print(len(indices))  # Number of potentially intersecting boxes

# Precise predicate query
precise_indices = tree.query(query_box, predicate='intersects')
print(len(precise_indices))  # Exact intersections only

# Other predicates
contains_indices = tree.query(query_box, predicate='contains')
within_indices = tree.query(query_box, predicate='within')

# Nearest neighbor search
query_point = Point(50, 50)
nearest_idx = tree.query_nearest(query_point)
print(f"Nearest geometry at index: {nearest_idx}")

# Nearest within distance
nearest_within = tree.query_nearest(query_point, max_distance=5.0, return_distance=True)
indices, distances = nearest_within
print(f"Found {len(indices)} geometries within distance 5.0")

# Bulk query (query multiple geometries at once)
query_geoms = [box(i*10, i*10, i*10+5, i*10+5) for i in range(5)]
bulk_results = tree.query_bulk(query_geoms, predicate='intersects')
tree_indices, query_indices = bulk_results
print(f"Found {len(tree_indices)} intersections")

# Access indexed geometries
for idx in indices[:5]:
    geom = tree.geometries[idx]
    print(f"Geometry {idx}: {geom.bounds}")
```

## WKT/WKB Serialization

Convert between Shapely geometries and standard text/binary formats.

```python
import shapely
from shapely import Point, LineString, Polygon
from shapely.wkt import dumps as wkt_dumps, loads as wkt_loads
from shapely.wkb import dumps as wkb_dumps, loads as wkb_loads

# Well-Known Text (WKT)
point = Point(10.5, 20.3)
wkt_string = wkt_dumps(point)
print(wkt_string)  # 'POINT (10.5000000000000000 20.3000000000000000)'

# Control precision
wkt_short = shapely.to_wkt(point, rounding_precision=2)
print(wkt_short)  # 'POINT (10.5 20.3)'

# Parse WKT
geom_from_wkt = wkt_loads('POINT (0 0)')
print(geom_from_wkt)  # <POINT (0 0)>

# Complex geometries
polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
wkt_poly = wkt_dumps(polygon)
print(wkt_poly)  # 'POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))'

# Well-Known Binary (WKB)
wkb_bytes = wkb_dumps(point)
print(type(wkb_bytes))  # <class 'bytes'>
print(len(wkb_bytes))  # Binary size

# WKB hex encoding
wkb_hex = shapely.to_wkb(point, hex=True)
print(wkb_hex[:20])  # Hex string representation

# Parse WKB
geom_from_wkb = wkb_loads(wkb_bytes)
print(geom_from_wkb.equals(point))  # True

# Vectorized serialization
geoms = [Point(i, i) for i in range(5)]
wkt_array = shapely.to_wkt(geoms, rounding_precision=1)
print(wkt_array)  # array(['POINT (0 0)', 'POINT (1 1)', ...])

# Error handling
try:
    invalid = shapely.from_wkt('INVALID WKT', on_invalid='raise')
except Exception as e:
    print(f"Error: {e}")

# Ignore invalid geometries
result = shapely.from_wkt(['POINT (0 0)', 'INVALID'], on_invalid='ignore')
print(result)  # array([<POINT (0 0)>, None])
```

## GeoJSON Integration

Work with GeoJSON format using the __geo_interface__ protocol.

```python
import json
from shapely.geometry import shape, mapping, Point, LineString, Polygon
import shapely

# Shapely to GeoJSON dict
point = Point(10, 20)
geojson_dict = mapping(point)
print(geojson_dict)
# {'type': 'Point', 'coordinates': (10.0, 20.0)}

# Complex geometry to GeoJSON
polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
poly_geojson = mapping(polygon)
print(json.dumps(poly_geojson, indent=2))
# {
#   "type": "Polygon",
#   "coordinates": [[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0], [0.0, 0.0]]]
# }

# GeoJSON dict to Shapely
geojson_point = {'type': 'Point', 'coordinates': [5.0, 7.0]}
shapely_geom = shape(geojson_point)
print(shapely_geom)  # <POINT (5 7)>

# Parse GeoJSON string
geojson_str = '{"type": "LineString", "coordinates": [[0, 0], [1, 1], [2, 0]]}'
geojson_obj = json.loads(geojson_str)
line = shape(geojson_obj)
print(line.length)  # ~2.828

# Using __geo_interface__ directly
print(point.__geo_interface__)
# {'type': 'Point', 'coordinates': (10.0, 20.0)}

# Vectorized GeoJSON operations
geoms = [Point(i, i*2) for i in range(3)]
geojson_strings = shapely.to_geojson(geoms)
print(geojson_strings)
# array(['{"type":"Point","coordinates":[0,0]}', ...])

# Parse array of GeoJSON
geojson_array = [
    '{"type": "Point", "coordinates": [0, 0]}',
    '{"type": "Point", "coordinates": [1, 1]}'
]
geoms_from_json = shapely.from_geojson(geojson_array)
print(geoms_from_json)  # array([<POINT (0 0)>, <POINT (1 1)>])
```

## Multi-Geometries - Collections

Work with collections of geometries.

```python
from shapely import MultiPoint, MultiLineString, MultiPolygon, GeometryCollection
from shapely import Point, LineString, Polygon
import shapely

# MultiPoint
points = MultiPoint([Point(0, 0), Point(1, 1), Point(2, 2)])
print(len(points.geoms))  # 3
for pt in points.geoms:
    print(pt.coords[:])

# MultiLineString
lines = MultiLineString([
    [(0, 0), (1, 1)],
    [(2, 2), (3, 3)],
    [(0, 1), (1, 0)]
])
print(lines.length)  # Total length of all lines

# MultiPolygon
polys = MultiPolygon([
    Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
    Polygon([(2, 2), (3, 2), (3, 3), (2, 3)])
])
print(polys.area)  # Total area: 2.0

# GeometryCollection (mixed types)
collection = GeometryCollection([
    Point(0, 0),
    LineString([(0, 0), (1, 1)]),
    Polygon([(2, 2), (3, 2), (3, 3), (2, 3)])
])
print(len(collection.geoms))  # 3
print([g.geom_type for g in collection.geoms])
# ['Point', 'LineString', 'Polygon']

# Operations on collections
buffer_all = shapely.buffer(points, 0.5)
print(buffer_all.geom_type)  # 'MultiPolygon'

# Union of multi-geometry
merged = shapely.unary_union(polys)
print(merged.geom_type)  # 'MultiPolygon' or 'Polygon' if merged

# Access individual geometries
first_poly = polys.geoms[0]
print(first_poly.area)  # 1.0
```

## Polygonize and Line Merging

Construct polygons from lines and merge connected linestrings.

```python
import shapely
from shapely import LineString, MultiLineString

# Line merge - join connected linestrings
lines = [
    LineString([(0, 0), (1, 1)]),
    LineString([(1, 1), (2, 2)]),
    LineString([(2, 2), (3, 0)])
]
merged = shapely.line_merge(lines)
print(merged.geom_type)  # 'LineString' (all connected)
print(len(merged.coords))  # 4 points

# Disconnected lines remain separate
lines_disconnected = [
    LineString([(0, 0), (1, 1)]),
    LineString([(5, 5), (6, 6)])
]
result = shapely.line_merge(lines_disconnected)
print(result.geom_type)  # 'MultiLineString' (not connected)

# Polygonize - create polygons from dangling lines
boundary_lines = [
    LineString([(0, 0), (2, 0)]),
    LineString([(2, 0), (2, 2)]),
    LineString([(2, 2), (0, 2)]),
    LineString([(0, 2), (0, 0)])
]
polygons = shapely.polygonize(boundary_lines)
print(polygons.geom_type)  # 'GeometryCollection'
print(len(list(polygons.geoms)))  # 1 polygon created

# Complex polygonization (with holes)
complex_lines = [
    # Outer square
    LineString([(0, 0), (4, 0)]),
    LineString([(4, 0), (4, 4)]),
    LineString([(4, 4), (0, 4)]),
    LineString([(0, 4), (0, 0)]),
    # Inner square (hole)
    LineString([(1, 1), (3, 1)]),
    LineString([(3, 1), (3, 3)]),
    LineString([(3, 3), (1, 3)]),
    LineString([(1, 3), (1, 1)])
]
poly_result = shapely.polygonize(complex_lines)
print(f"Created {len(list(poly_result.geoms))} polygons")
```

## Validation and Repair

Check geometry validity and fix invalid geometries.

```python
import shapely
from shapely import Polygon, LineString

# Create invalid polygon (self-intersecting bowtie)
invalid_poly = Polygon([(0, 0), (2, 2), (2, 0), (0, 2), (0, 0)])
print(shapely.is_valid(invalid_poly))  # False

# Get reason for invalidity
reason = shapely.is_valid_reason(invalid_poly)
print(reason)  # 'Self-intersection[1 1]'

# Make valid (repair)
valid_poly = shapely.make_valid(invalid_poly)
print(shapely.is_valid(valid_poly))  # True
print(valid_poly.geom_type)  # Might be 'MultiPolygon' after repair

# Normalize geometry (canonical form)
polygon = Polygon([(0, 0), (0, 1), (1, 1), (1, 0), (0, 0)])
normalized = shapely.normalize(polygon)
print(normalized.exterior.coords[:])  # Standardized coordinate order

# Check various validity conditions
print(shapely.is_simple(LineString([(0, 0), (1, 1), (0, 1), (1, 0)])))  # False (crosses itself)
print(shapely.is_ring(LineString([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)])))  # True
print(shapely.is_closed(LineString([(0, 0), (1, 1)])))  # False

# Empty geometry checks
from shapely import Point
empty = Point()
print(shapely.is_empty(empty))  # True
print(shapely.is_valid(empty))  # True (empty is valid)
```

## Voronoi Diagrams and Delaunay Triangulation

Generate computational geometry constructs from point sets.

```python
import shapely
from shapely import MultiPoint, Point

# Create point set
points = MultiPoint([
    Point(0, 0), Point(1, 0), Point(0.5, 1), Point(2, 0.5)
])

# Voronoi diagram
voronoi = shapely.voronoi_polygons(points)
print(voronoi.geom_type)  # 'GeometryCollection'
print(len(list(voronoi.geoms)))  # One polygon per point

# Extended Voronoi (larger envelope)
voronoi_extended = shapely.voronoi_polygons(points, extend_to=shapely.box(-2, -2, 4, 4))

# Delaunay triangulation
triangles = shapely.delaunay_triangles(points)
print(triangles.geom_type)  # 'GeometryCollection'
print([g.geom_type for g in triangles.geoms])  # All 'Polygon' (triangles)

# Extract edges from triangulation
for triangle in triangles.geoms:
    print(f"Triangle area: {triangle.area}")
    print(f"Triangle coords: {list(triangle.exterior.coords)}")

# Single point - no triangulation possible
single = MultiPoint([Point(0, 0)])
result = shapely.delaunay_triangles(single)
print(len(list(result.geoms)))  # 0 triangles

# Collinear points
collinear = MultiPoint([Point(i, i) for i in range(5)])
tri = shapely.delaunay_triangles(collinear)
print(len(list(tri.geoms)))  # 0 (collinear, no area)
```

## Vectorized Array Operations

Perform high-performance operations on geometry arrays.

```python
import shapely
import numpy as np

# Create geometry arrays
points = shapely.points(np.random.rand(1000, 2) * 100)
polygons = [shapely.box(i*10, i*10, i*10+15, i*10+15) for i in range(100)]

# Vectorized contains check (broadcasting)
container = shapely.box(0, 0, 50, 50)
contained_mask = shapely.contains(container, points)
print(f"Points inside: {contained_mask.sum()}")

# Vectorized distance calculations
origin = shapely.Point(0, 0)
distances = shapely.distance(origin, points)
print(f"Average distance: {distances.mean():.2f}")
print(f"Max distance: {distances.max():.2f}")

# Filter geometries by predicate
large_polygons_mask = shapely.area(polygons) > 200
large_polygons = np.array(polygons)[large_polygons_mask]
print(f"Large polygons: {len(large_polygons)}")

# Vectorized buffer
buffered_points = shapely.buffer(points[:100], 5.0)
print(f"Buffered geometries: {len(buffered_points)}")

# Pairwise operations
poly_array1 = np.array([shapely.box(i, 0, i+2, 2) for i in range(5)])
poly_array2 = np.array([shapely.box(i+1, 0, i+3, 2) for i in range(5)])
intersections = shapely.intersection(poly_array1, poly_array2)
intersection_areas = shapely.area(intersections)
print(f"Intersection areas: {intersection_areas}")

# Aggregate operations
all_points = shapely.points([(i, i) for i in range(10)])
bounds = shapely.total_bounds(all_points)
print(f"Total bounds: {bounds}")  # [minx, miny, maxx, maxy]

# Boolean array indexing
valid_mask = shapely.is_valid(polygons)
valid_polygons = np.array(polygons)[valid_mask]
print(f"Valid polygons: {len(valid_polygons)}/{len(polygons)}")
```

---

## Use Cases and Integration Patterns

Shapely excels in computational geometry tasks across multiple domains including GIS applications, spatial analysis, geographic data processing, urban planning, environmental modeling, and robotics path planning. Common use cases include calculating area and perimeter statistics for land parcels, finding intersections between administrative boundaries, generating buffer zones around features like roads or hazard areas, performing point-in-polygon tests for geocoding, simplifying complex geometries for visualization, and building spatial indices for efficient queries on large datasets.

Integration with the Python geospatial ecosystem is straightforward and powerful. Shapely works seamlessly with GeoPandas for dataframe-based spatial operations, Fiona for reading/writing spatial file formats (Shapefiles, GeoJSON, etc.), Rasterio for raster-vector operations, and Matplotlib/Cartopy for visualization. The library's focus on pure geometric operations without handling projections or I/O makes it a perfect building block - use PyProj for coordinate transformations, Shapely for geometry manipulation, and specialized libraries for data persistence. The __geo_interface__ protocol ensures compatibility across the ecosystem, while the dual API (object-oriented for single geometries, vectorized ufuncs for arrays) provides flexibility to optimize performance. Shapely's thread-safe operations with GIL release enable true parallel processing, making it suitable for high-performance spatial analytics on large datasets.