Try Live
Add Docs
Rankings
Pricing
Docs
Install
Theme
Install
Docs
Pricing
More...
More...
Try Live
Rankings
Enterprise
Create API Key
Add Docs
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
Context
Skills
Chat
Benchmark
94.5
Suggestions
Latest
Show doc for...
Code
Info
Show Results
Context Summary (auto-generated)
Raw
Copy
Link
# 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.