shapely

# Shapely - Planar Geometry Shapely is the engine behind GeoPandas and many other GIS tools. It focuses on the geometry itself: calculating intersections, unions, distances, and checking spatial relationships (like "is this point inside this polygon?"). ## When to Use - Precise manipulation of 2D geometric shapes. - Performing set-theoretic operations (Intersection, Union, Difference). - Checking spatial predicates (Contains, Within, Intersects, Touches). - Cleaning and validating "dirty" geometry (fixing self-intersections). - Calculating geometric properties (Area, Length, Centroid, Bounds). - Generating buffers or simplifying complex lines. - Linear referencing (finding points along a line). ## Reference Documentation **Official docs**: https://shapely.readthedocs.io/ **GEOS (Engine)**: https://libgeos.org/ **Search patterns**: `shapely.geometry`, `shapely.ops.unary_union`, `shapely.validation.make_valid` ## Core Principles ### Geometric Objects Objects are immutable. Once created, you don't change them; you perform an operation that returns a new object. - **Points**: 0-dimensional. - **LineStrings**: 1-dimensional curves. - **Polygons**: 2-dimensional surfaces with optional holes. ### Cartesian Geometry Shapely operates in a Cartesian plane. It does not know about Earth's curvature, latitudes, or longitudes. Distance is sqrt(dx² + dy²). ### Vectorization (Shapely 2.0+) Modern Shapely supports vectorized operations on NumPy arrays of geometry objects, making it significantly faster than older versions. ## Quick Reference ### Installation ```bash pip install shapely numpy ``` ### Standard Imports ```python import numpy as np from shapely import Point, LineString, Polygon, MultiPoint, MultiPolygon from shapely import ops, wkt, wkb import shapely ``` ### Basic Pattern - Creation and Analysis ```python from shapely.geometry import Point, Polygon # 1. Create objects p = Point(0, 0) poly = Polygon([(0, 0), (2, 0), (2, 2), (0, 2)]) # 2. Check relationships is_inside = p.within(poly) # True is_on_border = p.touches(poly) # False (interior counts as within) # 3. Calculate print(f"Area: {poly.area}") print(f"Distance: {p.distance(Point(10, 10))}") ``` ## Critical Rules ### ✅ DO - **Check Validity** - Use `.is_valid` before complex operations. Invalid geometry (like a self-intersecting polygon) will cause errors. - **Use unary_union** - When merging many polygons, `ops.unary_union([list])` is orders of magnitude faster than a loop of `p1.union(p2)`. - **Prefer Vectorized Functions** - Use `shapely.intersects(array_a, array_b)` instead of loops for performance. - **Use prepare()** - If you are checking many points against the same polygon, use `shapely.prepare(poly)` to speed up subsequent queries. - **Simplify for Analysis** - Use `.simplify(tolerance)` for complex boundaries to improve performance if high precision isn't required. ### ❌ DON'T - **Mix Cartesian and Spherical** - Don't calculate distance on Lat/Lon points; the result will be in meaningless "degrees". - **Assume Polygon Orientation** - While Shapely handles it, remember that exterior rings should ideally be counter-clockwise and holes clockwise. - **Use for 3D Math** - Shapely supports Z-coordinates for storage, but most operations (like area, intersection) ignore Z and project onto the XY plane. - **Loop over large collections** - Use NumPy-style vectorization provided in Shapely 2.0+. ## Anti-Patterns (NEVER) ```python from shapely.geometry import Point, Polygon from shapely.ops import unary_union # ❌ BAD: Merging geometries in a loop (O(n²) complexity) result = geometries[0] for g in geometries[1:]: result = result.union(g) # ✅ GOOD: Use unary_union (O(n log n) complexity) result = unary_union(geometries) # ❌ BAD: Checking many points without preparing for p in many_points: if complex_poly.contains(p): # Slow for complex shapes pass # ✅ GOOD: Prepare the geometry (Builds a spatial index) from shapely import prepare prepare(complex_poly) for p in many_points: if complex_poly.contains(p): # Much faster pass ``` ## Geometry Types and Creation ### Standard Primitives ```python # Point (x, y, z) pt = Point(1.0, 2.0) # LineString (Ordered sequence of points) line = LineString([(0, 0), (1, 1), (2, 0)]) # Polygon (Shell, [Holes]) shell = [(0, 0), (10, 0), (10, 10), (0, 10)] hole = [(2, 2), (2, 4), (4, 4), (4, 2)] poly = Polygon(shell, [hole]) # Multi-Geometries (Collections) points = MultiPoint([(0,0), (1,1)]) ``` ## Spatial Predicates (Relationships) ### Checking how objects relate ```python a = Point(1, 1).buffer(1.5) # A circle b = Polygon([(0,0), (2,0), (2,2), (0,2)]) # A square print(a.intersects(b)) # Shared space? print(a.contains(b)) # B entirely inside A? print(a.disjoint(b)) # No shared space? print(a.overlaps(b)) # Same dimension, shared space, but not within? print(a.touches(b)) # Only boundaries share space? print(a.crosses(b)) # Line crossing a polygon? ``` ## Set-Theoretic Operations ### Creating new geometries from old ones ```python poly1 = Point(0, 0).buffer(1) poly2 = Point(1, 0).buffer(1) # Intersection (Shared area) inter = poly1.intersection(poly2) # Union (Combined area) union = poly1.union(poly2) # Difference (Area in poly1 NOT in poly2) diff = poly1.difference(poly2) # Symmetric Difference (Area in either but NOT both) sdiff = poly1.symmetric_difference(poly2) ``` ## Constructive Methods ### Buffering, Splicing, and Simplifying ```python # Buffer: Expand/shrink geometry # cap_style: 1=Round, 2=Flat, 3=Square line_thick = line.buffer(0.5, cap_style=2) # Centroid: Geometric center center = poly.centroid # Representative Point: Guaranteed to be INSIDE the geometry # Useful for label placement in U-shaped polygons label_pt = poly.representative_point() # Simplify: Reduce number of vertices simple_line = complex_line.simplify(tolerance=0.1, preserve_topology=True) # Convex Hull: Smallest convex box containing all points hull = MultiPoint(points).convex_hull ``` ## Linear Referencing ### Working with positions along a LineString ```python line = LineString([(0, 0), (0, 10), (10, 10)]) # Find distance along line to the point nearest to (5, 5) dist = line.project(Point(5, 5)) # returns 5.0 (it's at (0, 5)) # Find the actual point at a specific distance along the line pt = line.interpolate(15.0) # returns Point(5, 10) ``` ## I/O: WKT, WKB, and NumPy ### Serialization and Data Exchange ```python # WKT (Well-Known Text) - Human readable text = "POINT (10 20)" p = wkt.loads(text) print(p.wkt) # WKB (Well-Known Binary) - Fast and compact binary = wkb.dumps(p) p_new = wkb.loads(binary) # NumPy Integration (Shapely 2.0) points_array = np.array([Point(0,0), Point(1,1), Point(2,2)]) areas = shapely.area(points_array) # Returns array of zeros dist_matrix = shapely.distance(points_array[:, np.newaxis], points_array) ``` ## Practical Workflows ### 1. Cleaning Invalid Geometries ```python from shapely.validation import make_valid def safe_area(geom): """Calculates area even for invalid/self-intersecting polygons.""" if not geom.is_valid: geom = make_valid(geom) # After make_valid, a Polygon might become a MultiPolygon or GeometryCollection return geom.area ``` ### 2. Point-in-Polygon Search (Optimized) ```python from shapely import prepare def find_points_in_poly(points, poly): """Efficiently filters points inside a complex polygon.""" prepare(poly) # Builds internal STRtree or spatial index # Using vectorized intersection (much faster) mask = shapely.contains(poly, points) return points[mask] ``` ### 3. Splitting a Polygon by a Line ```python from shapely.ops import split def divide_land(polygon, line): """Splits a polygon into multiple parts using a LineString.""" result = split(polygon, line) # Returns a GeometryCollection of the resulting parts return list(result.geoms) ``` ## Performance Optimization ### STRtree for Nearest Neighbors If you have thousands of geometries and need to find which ones are near a point, use STRtree. ```python from shapely import STRtree tree = STRtree(geometries) # Find indices of geometries whose bounding boxes intersect the point's buffer indices = tree.query(Point(0,0).buffer(10)) # Find the single nearest geometry index nearest_idx = tree.nearest(Point(0,0)) ``` ## Common Pitfalls and Solutions ### The "Sliver Polygon" problem Calculations like intersection can sometimes produce tiny, almost invisible polygons due to floating-point errors. ```python # ✅ Solution: Filter by area intersection = p1.intersection(p2) if intersection.area < 1e-9: intersection = None ``` ### Latitude/Longitude Confusion Points are (x, y). In GIS, this usually means (Longitude, Latitude). ```python # ❌ Error: Point(Latitude, Longitude) # This will plot your maps sideways! # ✅ Solution: Always use (Lon, Lat) to match (X, Y) nyc = Point(-74.006, 40.7128) ``` ### GeometryCollections Operations like split or intersection can return GeometryCollection. This is a container for mixed types. ```python # ❌ Problem: Calling .area on a collection with Lines and Polygons # ✅ Solution: Filter for the type you want polys = [g for g in collection.geoms if g.geom_type == 'Polygon'] ``` ## Best Practices 1. Always validate geometry with `.is_valid` before complex operations 2. Use `unary_union` instead of looping over unions 3. Prepare geometries with `prepare()` when checking many points against the same shape 4. Use vectorized operations in Shapely 2.0+ for performance 5. Remember Shapely is Cartesian - don't use lat/lon directly for distance calculations 6. Filter sliver polygons by area threshold after geometric operations 7. Use STRtree for spatial indexing when working with many geometries 8. Simplify complex geometries when high precision isn't required 9. Handle GeometryCollections properly - filter by geometry type when needed 10. Use representative_point() for guaranteed interior points in complex polygons Shapely is a specialized, sharp tool. It doesn't care about your coordinate system or your file format — it only cares about the pure, mathematical relationship between shapes. Mastering it is the key to building advanced spatial algorithms.