Shapely creates and analyzes geometric objects in Python. pip install shapely. Point: from shapely.geometry import Point; p = Point(1.0, 2.0). LineString: from shapely.geometry import LineString; line = LineString([(0,0),(1,1),(2,0)]). Polygon: from shapely.geometry import Polygon; poly = Polygon([(0,0),(1,0),(1,1),(0,1)]). With hole: Polygon(exterior, [hole_coords]). Multi: MultiPolygon([poly1, poly2]). WKT: from shapely import wkt; p = wkt.loads("POINT (1 2)"). GeoJSON: from shapely.geometry import shape; geom = shape(geojson_dict). to_wkt: poly.wkt. to_geojson: from shapely.geometry import mapping; mapping(poly). Contains: poly.contains(point). Intersects: a.intersects(b). Within: a.within(b). Crosses: a.crosses(b). Touches: a.touches(b). Disjoint: a.disjoint(b). Intersection: a.intersection(b). Union: a.union(b). Difference: a.difference(b). Sym diff: a.symmetric_difference(b). Buffer: p.buffer(1.0) — circle radius 1. Line buffer: line.buffer(0.5) — corridor. Simplify: poly.simplify(0.01, preserve_topology=True). Convex hull: poly.convex_hull. Centroid: poly.centroid → Point. Bounds: poly.bounds → (minx, miny, maxx, maxy). Area: poly.area. Length: line.length. Distance: a.distance(b). Nearest point: shapely.ops.nearest_points(a, b). Unary union: from shapely.ops import unary_union; merged = unary_union([a, b, c]). Affine: from shapely.affinity import translate, rotate, scale; rotate(poly, 45, origin="centroid"). STR-tree: from shapely.strtree import STRtree; tree = STRtree(geoms); tree.query(bbox). Claude Code generates Shapely spatial analyzers, GIS tools, collision detectors, and geometry preprocessors.
CLAUDE.md for Shapely
## Shapely Stack
- Version: shapely >= 2.0 | pip install shapely
- Create: Point(x,y) | LineString(coords) | Polygon(exterior, [holes])
- Load: wkt.loads("WKT") | shape(geojson_dict) | from_wkb(bytes)
- Predicates: contains | within | intersects | touches | crosses | disjoint
- Ops: a.intersection(b) | a.union(b) | a.buffer(r) | a.simplify(tol)
- Measure: .area | .length | .distance(b) | .bounds | .centroidShapely Geometry Pipeline
# app/geometry.py — Shapely shapes, predicates, set ops, spatial index, transforms, GeoJSON
from __future__ import annotations
import json
from typing import Any, Generator
from shapely import wkt, wkb
from shapely.affinity import (
affine_transform,
rotate,
scale,
translate,
)
from shapely.geometry import (
GeometryCollection,
LinearRing,
LineString,
MultiLineString,
MultiPoint,
MultiPolygon,
Point,
Polygon,
mapping,
shape,
)
from shapely.ops import (
nearest_points,
split,
unary_union,
polygonize,
transform,
)
from shapely.strtree import STRtree
# ─────────────────────────────────────────────────────────────────────────────
# 1. Construction helpers
# ─────────────────────────────────────────────────────────────────────────────
def make_point(x: float, y: float, z: float | None = None) -> Point:
"""
Create a 2D or 3D Point.
Example:
p2 = make_point(1.0, 2.0)
p3 = make_point(1.0, 2.0, 100.0)
"""
return Point(x, y) if z is None else Point(x, y, z)
def make_bbox(minx: float, miny: float, maxx: float, maxy: float) -> Polygon:
"""
Create a bounding-box Polygon.
Example:
bbox = make_bbox(-122.5, 37.7, -122.3, 37.9) # SF area
"""
return Polygon([
(minx, miny), (maxx, miny),
(maxx, maxy), (minx, maxy),
(minx, miny),
])
def from_geojson(geojson: dict | str) -> Any:
"""
Create a Shapely geometry from a GeoJSON dict or string.
Example:
geom = from_geojson({"type": "Point", "coordinates": [1.0, 2.0]})
geom = from_geojson(geojson_string)
"""
if isinstance(geojson, str):
geojson = json.loads(geojson)
if "geometry" in geojson: # Feature
geojson = geojson["geometry"]
return shape(geojson)
def to_geojson(geom) -> dict:
"""Convert Shapely geometry to GeoJSON dict."""
return mapping(geom)
def from_wkt_str(wkt_str: str):
"""Load geometry from WKT string."""
return wkt.loads(wkt_str)
def to_wkt_str(geom, rnd: int | None = 6) -> str:
"""Serialize geometry to WKT, rounding coordinates."""
if rnd is not None:
return geom.wkt # Shapely 2 handles precision via set_precision
return geom.wkt
# ─────────────────────────────────────────────────────────────────────────────
# 2. Spatial predicates
# ─────────────────────────────────────────────────────────────────────────────
def contains(container, item) -> bool:
"""Return True if container fully contains item (no boundary contact)."""
return container.contains(item)
def intersects_geom(a, b) -> bool:
"""Return True if a and b share any geometry."""
return a.intersects(b)
def within_distance(a, b, distance: float) -> bool:
"""Return True if the distance between a and b is <= distance."""
return a.distance(b) <= distance
def classify_relationship(a, b) -> str:
"""
Classify the DE-9IM spatial relationship between two geometries.
Returns one of: contains, within, intersects, touches, crosses, disjoint, equals
Example:
rel = classify_relationship(polygon, point) # "contains" or "disjoint"
"""
if a.equals(b): return "equals"
if a.contains(b): return "contains"
if a.within(b): return "within"
if a.crosses(b): return "crosses"
if a.touches(b): return "touches"
if a.intersects(b): return "intersects"
return "disjoint"
# ─────────────────────────────────────────────────────────────────────────────
# 3. Set operations
# ─────────────────────────────────────────────────────────────────────────────
def intersection(a, b):
"""Return the intersection geometry of a and b."""
return a.intersection(b)
def union(a, b):
"""Return the union of a and b."""
return a.union(b)
def merge_geometries(geoms: list) -> MultiPolygon | Polygon | GeometryCollection:
"""
Merge a list of geometries into a single geometry using unary_union.
Example:
block = merge_geometries([parcel1, parcel2, parcel3])
"""
return unary_union(geoms)
def clip(geom, bbox_or_mask):
"""
Clip a geometry to a bounding box or mask polygon.
Example:
clipped = clip(city_boundary, study_area_polygon)
"""
return geom.intersection(bbox_or_mask)
def subtract(base, other):
"""Remove other's area from base (difference)."""
return base.difference(other)
# ─────────────────────────────────────────────────────────────────────────────
# 4. Measurement
# ─────────────────────────────────────────────────────────────────────────────
def area(geom) -> float:
"""Return the area (in coordinate units squared)."""
return geom.area
def length(geom) -> float:
"""Return the length/perimeter (in coordinate units)."""
return geom.length
def distance_between(a, b) -> float:
"""Euclidean distance between two geometries."""
return a.distance(b)
def closest_points(a, b) -> tuple[Point, Point]:
"""
Return the pair of closest points between two geometries.
Example:
p1, p2 = closest_points(line, polygon)
dist = p1.distance(p2)
"""
return nearest_points(a, b)
def bounds(geom) -> tuple[float, float, float, float]:
"""Return (minx, miny, maxx, maxy) bounding box."""
return geom.bounds
def centroid(geom) -> Point:
"""Return the geometric centroid."""
return geom.centroid
def envelope(geom) -> Polygon:
"""Return the minimum bounding rectangle."""
return geom.envelope
# ─────────────────────────────────────────────────────────────────────────────
# 5. Transformations
# ─────────────────────────────────────────────────────────────────────────────
def buffer_geom(geom, distance: float, resolution: int = 16, cap_style: str = "round") -> Polygon:
"""
Buffer a geometry by distance.
cap_style: "round" | "flat" | "square"
resolution: segments per quarter circle (higher = smoother).
Example:
circle = buffer_geom(point, 100) # 100-unit circle
corridor = buffer_geom(road_line, 5) # 5-unit road buffer
expanded = buffer_geom(polygon, -10) # inset (negative buffer)
"""
cap_map = {"round": 1, "flat": 2, "square": 3}
return geom.buffer(distance, resolution=resolution, cap_style=cap_map.get(cap_style, 1))
def simplify_geom(geom, tolerance: float = 1.0, preserve_topology: bool = True):
"""
Simplify geometry using Douglas-Peucker algorithm.
Example:
simple = simplify_geom(complex_coastline, tolerance=0.001)
"""
return geom.simplify(tolerance, preserve_topology=preserve_topology)
def translate_geom(geom, x: float = 0, y: float = 0, z: float = 0):
"""Translate geometry by (x, y, z) offset."""
return translate(geom, xoff=x, yoff=y, zoff=z)
def rotate_geom(geom, angle_degrees: float, origin: str | tuple = "centroid"):
"""
Rotate geometry by angle_degrees around origin.
origin: "centroid" | "center" | (x, y)
Example:
rotated = rotate_geom(polygon, 45)
rotated_around = rotate_geom(polygon, 90, origin=(0, 0))
"""
return rotate(geom, angle_degrees, origin=origin)
def scale_geom(geom, x_factor: float = 1.0, y_factor: float = 1.0, origin: str | tuple = "centroid"):
"""Scale geometry by x_factor and y_factor."""
return scale(geom, xfact=x_factor, yfact=y_factor, origin=origin)
# ─────────────────────────────────────────────────────────────────────────────
# 6. Spatial indexing
# ─────────────────────────────────────────────────────────────────────────────
class SpatialIndex:
"""
Wrapper around Shapely STRtree for efficient spatial queries.
Example:
idx = SpatialIndex(parcels)
nearby = idx.query_intersects(search_polygon)
closest = idx.nearest(point, k=5)
"""
def __init__(self, geometries: list) -> None:
self._geoms = geometries
self._tree = STRtree(geometries)
def query_intersects(self, geom) -> list:
"""Return all geometries that intersect geom."""
indices = self._tree.query(geom, predicate="intersects")
return [self._geoms[i] for i in indices]
def query_within(self, geom) -> list:
"""Return geometries fully contained within geom."""
indices = self._tree.query(geom, predicate="within")
return [self._geoms[i] for i in indices]
def query_bbox(self, minx: float, miny: float, maxx: float, maxy: float) -> list:
"""Return geometries intersecting a bounding box."""
bbox = make_bbox(minx, miny, maxx, maxy)
return self.query_intersects(bbox)
def nearest(self, geom, k: int = 1) -> list:
"""Return k nearest geometries."""
indices = self._tree.nearest(geom) if k == 1 else self._tree.query_nearest(geom, max_distance=None, exclusive=False, all_matches=False)
if not hasattr(indices, "__iter__"):
indices = [indices]
return [self._geoms[i] for i in list(indices)[:k]]
def __len__(self) -> int:
return len(self._geoms)
# ─────────────────────────────────────────────────────────────────────────────
# 7. GeoJSON collection helpers
# ─────────────────────────────────────────────────────────────────────────────
def to_feature(geom, properties: dict | None = None) -> dict:
"""Wrap a geometry as a GeoJSON Feature."""
return {"type": "Feature", "geometry": to_geojson(geom), "properties": properties or {}}
def to_feature_collection(geoms: list, properties_list: list[dict] | None = None) -> dict:
"""
Convert a list of geometries to a GeoJSON FeatureCollection.
Example:
fc = to_feature_collection(polygons, [{"id": i} for i in range(len(polygons))])
"""
props = properties_list or [{} for _ in geoms]
features = [to_feature(g, p) for g, p in zip(geoms, props)]
return {"type": "FeatureCollection", "features": features}
# ─────────────────────────────────────────────────────────────────────────────
# Demo
# ─────────────────────────────────────────────────────────────────────────────
if __name__ == "__main__":
print("=== Construction ===")
pt = make_point(1.0, 2.0)
line = LineString([(0, 0), (4, 0), (4, 4)])
poly = Polygon([(0, 0), (3, 0), (3, 3), (0, 3)])
bbox = make_bbox(0, 0, 5, 5)
print(f" Point: {pt.wkt}")
print(f" Polygon: area={poly.area:.1f}, perimeter={poly.length:.1f}")
print(f" Bounds: {bbox.bounds}")
print("\n=== Predicates ===")
print(f" poly.contains(pt): {contains(poly, pt)}")
print(f" poly.intersects(line): {intersects_geom(poly, line)}")
print(f" classify: {classify_relationship(bbox, poly)}")
print("\n=== Set operations ===")
poly2 = Polygon([(2, 2), (5, 2), (5, 5), (2, 5)])
inter = intersection(poly, poly2)
unio = union(poly, poly2)
diff = subtract(poly, poly2)
print(f" Intersection area: {inter.area:.2f}")
print(f" Union area: {unio.area:.2f}")
print(f" Difference area: {diff.area:.2f}")
print("\n=== Measurements ===")
print(f" Distance pt to poly: {distance_between(make_point(5, 5), poly):.2f}")
p1, p2 = closest_points(make_point(5, 5), poly)
print(f" Closest points: {p1.wkt} — {p2.wkt}")
print("\n=== Transforms ===")
circle = buffer_geom(pt, 2.0)
print(f" Circle area (r=2): {circle.area:.2f}")
simple = simplify_geom(circle, tolerance=0.5)
print(f" Simplified coords: {len(list(simple.exterior.coords))}")
rotated = rotate_geom(poly, 45)
print(f" Rotated centroid: {rotated.centroid.wkt}")
print("\n=== Spatial index ===")
polys = [buffer_geom(make_point(i, i), 0.6) for i in range(10)]
idx = SpatialIndex(polys)
query_area = make_bbox(1.5, 1.5, 3.5, 3.5)
found = idx.query_intersects(query_area)
print(f" Indexed {len(idx)} circles, found {len(found)} in query bbox")
print("\n=== GeoJSON round-trip ===")
fc = to_feature_collection([poly, poly2], [{"id": 1}, {"id": 2}])
reloaded = [from_geojson(f["geometry"]) for f in fc["features"]]
print(f" FeatureCollection features: {len(fc['features'])}")
print(f" Reloaded area[0]: {reloaded[0].area:.1f}")For the GeoPandas alternative — GeoPandas builds on Shapely and Fiona to provide a DataFrame-style interface with .plot() visualization, CRS management via pyproj, and read/write support for Shapefiles, GeoPackages, and PostGIS; Shapely is the pure-geometry computational engine underneath — use Shapely directly when you need geometry operations in non-tabular pipelines (collision detection, API geometry validation, spatial joins in code), GeoPandas when you’re working with geographic datasets that need CRS transformations and tabular analysis. For the pyproj alternative — pyproj handles coordinate reference system (CRS) transformations between geographic and projected systems (WGS84, UTM, Web Mercator); Shapely is CRS-agnostic and works in pure coordinate space — combine Shapely with pyproj’s Transformer when real-world distances (meters) matter; use Shapely alone when your data is already in a consistent projected coordinate system. The Claude Skills 360 bundle includes Shapely skill sets covering make_point()/make_bbox()/from_geojson()/from_wkt_str() construction, contains()/intersects_geom()/within_distance()/classify_relationship() predicates, intersection()/union()/merge_geometries()/clip()/subtract() set ops, area()/length()/distance_between()/closest_points()/bounds()/centroid() measurement, buffer_geom()/simplify_geom()/translate_geom()/rotate_geom()/scale_geom() transforms, SpatialIndex with STRtree, and to_feature_collection() GeoJSON export. Start with the free tier to try spatial geometry and GIS analysis code generation.