Shapely Integration

pip install GeoAlchemy2[shapely]
# or
pip install GeoAlchemy2 shapely

def to_shape(element):
    """
    Convert a GeoAlchemy2 spatial element to a Shapely geometry.
    
    Parameters:
    - element: Union[WKBElement, WKTElement], spatial element to convert
    
    Returns:
    Shapely geometry object (Point, LineString, Polygon, etc.)
    
    Raises:
    ImportError: if Shapely is not installed
    """

from geoalchemy2.shape import to_shape
from geoalchemy2 import WKTElement

# Convert WKT element to Shapely
wkt_elem = WKTElement('POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))', srid=4326)
shapely_geom = to_shape(wkt_elem)

# Query database and convert result
lake = session.query(Lake).first()
shapely_lake = to_shape(lake.geom)  # lake.geom is WKBElement

# Use Shapely operations
print(f"Area: {shapely_lake.area}")
print(f"Centroid: {shapely_lake.centroid}")
print(f"Bounds: {shapely_lake.bounds}")

# Advanced Shapely operations
simplified = shapely_lake.simplify(0.1)
buffered = shapely_lake.buffer(100)
convex_hull = shapely_lake.convex_hull

# Check geometric properties
print(f"Is valid: {shapely_lake.is_valid}")
print(f"Is simple: {shapely_lake.is_simple}")
print(f"Geometry type: {shapely_lake.geom_type}")

def from_shape(shape, srid: int = -1, extended: bool = None):
    """
    Convert a Shapely geometry to a GeoAlchemy2 spatial element.
    
    Parameters:
    - shape: Shapely geometry object
    - srid: int, spatial reference system identifier (default: -1)
    - extended: bool, use extended WKB format (default: None)
    
    Returns:
    WKBElement representing the Shapely geometry
    
    Raises:
    ImportError: if Shapely is not installed
    """

from geoalchemy2.shape import from_shape
from shapely.geometry import Point, Polygon, LineString
from shapely import wkt

# Create Shapely geometries
point = Point(1, 2)
polygon = Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])
line = LineString([(0, 0), (1, 1), (2, 0)])

# Convert to GeoAlchemy2 elements
point_elem = from_shape(point, srid=4326)
polygon_elem = from_shape(polygon, srid=4326)
line_elem = from_shape(line, srid=4326)

# Store in database
new_lake = Lake(name='Shapely Lake', geom=polygon_elem)
session.add(new_lake)
session.commit()

# Load from WKT and convert
shapely_from_wkt = wkt.loads('MULTIPOINT((0 0), (1 1), (2 2))')
multipoint_elem = from_shape(shapely_from_wkt, srid=4326)

HAS_SHAPELY: bool
"""Boolean indicating whether Shapely is available."""

def check_shapely():
    """
    Context manager to verify Shapely availability.
    
    Raises:
    ImportError: if Shapely is not installed with helpful message
    """

from geoalchemy2.shape import HAS_SHAPELY, check_shapely

# Check availability before using
if HAS_SHAPELY:
    from geoalchemy2.shape import to_shape, from_shape
    # Use Shapely functions
else:
    print("Shapely not available")

# Use context manager for safe operations
try:
    with check_shapely():
        shapely_geom = to_shape(spatial_element)
        # Perform Shapely operations
except ImportError as e:
    print(f"Shapely required: {e}")

from geoalchemy2.shape import to_shape
from shapely.ops import unary_union, cascaded_union
from shapely.geometry import Point

# Query spatial data from database
buildings = session.query(Building).filter(
    Building.district_id == district_id
).all()

# Convert all geometries to Shapely
shapely_buildings = [to_shape(b.geom) for b in buildings]

# Perform complex Shapely operations
total_footprint = unary_union(shapely_buildings)
building_centroids = [geom.centroid for geom in shapely_buildings]

# Analysis with Shapely
total_area = total_footprint.area
avg_building_area = sum(g.area for g in shapely_buildings) / len(shapely_buildings)

# Find buildings near a point
query_point = Point(x, y)
nearby_buildings = [
    building for building, geom in zip(buildings, shapely_buildings)
    if geom.distance(query_point) < 100
]

from geoalchemy2.shape import from_shape
from shapely.geometry import Polygon
from shapely.ops import transform
import pyproj

# Create complex geometry with Shapely
exterior = [(0, 0), (10, 0), (10, 10), (5, 15), (0, 10)]
holes = [[(2, 2), (8, 2), (8, 8), (2, 8)]]
complex_polygon = Polygon(exterior, holes)

# Perform geometric operations
simplified = complex_polygon.simplify(0.5)
buffered = simplified.buffer(1.0)

# Transform coordinate system using pyproj
transformer = pyproj.Transformer.from_crs('EPSG:4326', 'EPSG:3857', always_xy=True)
transformed = transform(transformer.transform, buffered)

# Convert back to GeoAlchemy2 and store
geom_elem = from_shape(transformed, srid=3857)
new_feature = SpatialFeature(name='Processed Area', geom=geom_elem)
session.add(new_feature)
session.commit()

from geoalchemy2.shape import to_shape, from_shape
from geoalchemy2.functions import ST_Transform, ST_Buffer
from shapely.geometry import Point

# Start with database operations for efficiency
query_point = Point(lon, lat)
query_elem = from_shape(query_point, srid=4326)

# Use database for initial filtering (more efficient)
candidate_areas = session.query(ProtectedArea).filter(
    ST_DWithin(
        ST_Transform(ProtectedArea.geom, 4326),
        query_elem,
        0.01  # ~1km in degrees
    )
).all()

# Convert to Shapely for complex analysis
shapely_areas = [to_shape(area.geom) for area in candidate_areas]
query_shapely = to_shape(query_elem)

# Perform complex Shapely operations
overlapping_areas = []
for area, shapely_area in zip(candidate_areas, shapely_areas):
    if shapely_area.contains(query_shapely):
        # Calculate overlap with 500m buffer
        buffer_zone = query_shapely.buffer(0.005)  # ~500m in degrees
        overlap = shapely_area.intersection(buffer_zone)
        if overlap.area > 0:
            overlapping_areas.append((area, overlap.area))

# Sort by overlap area
overlapping_areas.sort(key=lambda x: x[1], reverse=True)

from geoalchemy2.shape import to_shape, from_shape, check_shapely
from shapely.errors import ShapelyError

try:
    with check_shapely():
        # Convert database geometry
        shapely_geom = to_shape(db_element)
        
        # Validate geometry
        if not shapely_geom.is_valid:
            # Fix invalid geometry
            shapely_geom = shapely_geom.buffer(0)
        
        # Perform operations
        processed = shapely_geom.simplify(tolerance)
        
        # Convert back
        result_elem = from_shape(processed, srid=target_srid)
        
except ImportError:
    # Shapely not available - use database functions instead
    result_elem = db_element.ST_Simplify(tolerance)
    
except ShapelyError as e:
    # Handle Shapely-specific errors
    print(f"Geometric operation failed: {e}")
    result_elem = db_element  # Use original

# Efficient: Filter in database first, then use Shapely
candidates = session.query(Feature).filter(
    func.ST_DWithin(Feature.geom, query_point, 1000)
).all()
shapely_results = [to_shape(f.geom) for f in candidates]

# Inefficient: Load all data then filter with Shapely
all_features = session.query(Feature).all()
shapely_all = [to_shape(f.geom) for f in all_features]
filtered = [g for g in shapely_all if g.distance(query_shapely) < 1000]