""" Geographic Utility Functions Coordinate transformations, distance calculations, and tile math """ import math from typing import Tuple, Optional import pyproj from shapely.geometry import Point, Polygon, shape from shapely.ops import transform # Web Mercator (EPSG:3857) and WGS84 (EPSG:4326) transformers WGS84 = pyproj.CRS("EPSG:4326") WEB_MERCATOR = pyproj.CRS("EPSG:3857") ETRS89_LAEA = pyproj.CRS("EPSG:3035") # Equal area for Europe def lat_lon_to_tile(lat: float, lon: float, zoom: int) -> Tuple[int, int]: """ Convert latitude/longitude to tile coordinates (XYZ scheme). Args: lat: Latitude in degrees (-85.05 to 85.05) lon: Longitude in degrees (-180 to 180) zoom: Zoom level (0-22) Returns: Tuple of (x, y) tile coordinates """ n = 2 ** zoom x = int((lon + 180.0) / 360.0 * n) lat_rad = math.radians(lat) y = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n) # Clamp to valid range x = max(0, min(n - 1, x)) y = max(0, min(n - 1, y)) return x, y def tile_to_bounds(z: int, x: int, y: int) -> Tuple[float, float, float, float]: """ Convert tile coordinates to bounding box. Args: z: Zoom level x: Tile X coordinate y: Tile Y coordinate Returns: Tuple of (west, south, east, north) in degrees """ n = 2 ** z west = x / n * 360.0 - 180.0 east = (x + 1) / n * 360.0 - 180.0 north_rad = math.atan(math.sinh(math.pi * (1 - 2 * y / n))) south_rad = math.atan(math.sinh(math.pi * (1 - 2 * (y + 1) / n))) north = math.degrees(north_rad) south = math.degrees(south_rad) return west, south, east, north def tile_to_center(z: int, x: int, y: int) -> Tuple[float, float]: """ Get the center point of a tile. Returns: Tuple of (longitude, latitude) in degrees """ west, south, east, north = tile_to_bounds(z, x, y) return (west + east) / 2, (south + north) / 2 def calculate_distance( lat1: float, lon1: float, lat2: float, lon2: float ) -> float: """ Calculate the distance between two points using the Haversine formula. Args: lat1, lon1: First point coordinates in degrees lat2, lon2: Second point coordinates in degrees Returns: Distance in meters """ R = 6371000 # Earth's radius in meters lat1_rad = math.radians(lat1) lat2_rad = math.radians(lat2) delta_lat = math.radians(lat2 - lat1) delta_lon = math.radians(lon2 - lon1) a = ( math.sin(delta_lat / 2) ** 2 + math.cos(lat1_rad) * math.cos(lat2_rad) * math.sin(delta_lon / 2) ** 2 ) c = 2 * math.atan2(math.sqrt(a), math.sqrt(1 - a)) return R * c def transform_coordinates( geometry, from_crs: str = "EPSG:4326", to_crs: str = "EPSG:3857", ): """ Transform a Shapely geometry between coordinate reference systems. Args: geometry: Shapely geometry object from_crs: Source CRS (default WGS84) to_crs: Target CRS (default Web Mercator) Returns: Transformed geometry """ transformer = pyproj.Transformer.from_crs( from_crs, to_crs, always_xy=True, ) return transform(transformer.transform, geometry) def calculate_area_km2(geojson: dict) -> float: """ Calculate the area of a GeoJSON polygon in square kilometers. Uses ETRS89-LAEA projection for accurate area calculation in Europe. Args: geojson: GeoJSON geometry dict Returns: Area in square kilometers """ geom = shape(geojson) geom_projected = transform_coordinates(geom, "EPSG:4326", "EPSG:3035") return geom_projected.area / 1_000_000 def is_within_bounds( point: Tuple[float, float], bounds: Tuple[float, float, float, float], ) -> bool: """ Check if a point is within a bounding box. Args: point: (longitude, latitude) tuple bounds: (west, south, east, north) tuple Returns: True if point is within bounds """ lon, lat = point west, south, east, north = bounds return west <= lon <= east and south <= lat <= north def get_germany_bounds() -> Tuple[float, float, float, float]: """Get the bounding box of Germany.""" return (5.87, 47.27, 15.04, 55.06) def meters_per_pixel(lat: float, zoom: int) -> float: """ Calculate the ground resolution at a given latitude and zoom level. Args: lat: Latitude in degrees zoom: Zoom level Returns: Meters per pixel at that location and zoom """ # Earth's circumference at equator in meters C = 40075016.686 # Resolution at equator for this zoom level resolution_equator = C / (256 * (2 ** zoom)) # Adjust for latitude (Mercator projection) return resolution_equator * math.cos(math.radians(lat)) def simplify_polygon(geojson: dict, tolerance: float = 0.0001) -> dict: """ Simplify a polygon geometry to reduce complexity. Args: geojson: GeoJSON geometry dict tolerance: Simplification tolerance in degrees Returns: Simplified GeoJSON geometry """ from shapely.geometry import mapping geom = shape(geojson) simplified = geom.simplify(tolerance, preserve_topology=True) return mapping(simplified) def buffer_polygon(geojson: dict, distance_meters: float) -> dict: """ Buffer a polygon by a distance in meters. Args: geojson: GeoJSON geometry dict distance_meters: Buffer distance in meters Returns: Buffered GeoJSON geometry """ from shapely.geometry import mapping geom = shape(geojson) # Transform to metric CRS, buffer, transform back geom_metric = transform_coordinates(geom, "EPSG:4326", "EPSG:3035") buffered = geom_metric.buffer(distance_meters) geom_wgs84 = transform_coordinates(buffered, "EPSG:3035", "EPSG:4326") return mapping(geom_wgs84) def get_tiles_for_bounds( bounds: Tuple[float, float, float, float], zoom: int, ) -> list[Tuple[int, int]]: """ Get all tile coordinates that cover a bounding box. Args: bounds: (west, south, east, north) in degrees zoom: Zoom level Returns: List of (x, y) tile coordinates """ west, south, east, north = bounds # Get corner tiles x_min, y_max = lat_lon_to_tile(south, west, zoom) x_max, y_min = lat_lon_to_tile(north, east, zoom) # Generate all tiles in range tiles = [] for x in range(x_min, x_max + 1): for y in range(y_min, y_max + 1): tiles.append((x, y)) return tiles