This repository has been archived on 2026-02-15. You can view files and clone it. You cannot open issues or pull requests or push a commit.
Files
breakpilot-pwa/geo-service/utils/geo_utils.py
Benjamin Admin bfdaf63ba9 fix: Restore all files lost during destructive rebase
A previous `git pull --rebase origin main` dropped 177 local commits,
losing 3400+ files across admin-v2, backend, studio-v2, website,
klausur-service, and many other services. The partial restore attempt
(660295e2) only recovered some files.

This commit restores all missing files from pre-rebase ref 98933f5e
while preserving post-rebase additions (night-scheduler, night-mode UI,
NightModeWidget dashboard integration).

Restored features include:
- AI Module Sidebar (FAB), OCR Labeling, OCR Compare
- GPU Dashboard, RAG Pipeline, Magic Help
- Klausur-Korrektur (8 files), Abitur-Archiv (5+ files)
- Companion, Zeugnisse-Crawler, Screen Flow
- Full backend, studio-v2, website, klausur-service
- All compliance SDKs, agent-core, voice-service
- CI/CD configs, documentation, scripts

Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com>
2026-02-09 09:51:32 +01:00

263 lines
6.5 KiB
Python

"""
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