"""Geo pipeline — processes raw GPS tracks into zone polygons. All geo math is done via Shapely (pure Python), so it works with both SQLite (dev) and PostgreSQL (prod). """ import json import math import uuid from datetime import datetime, timezone from rdp import rdp from shapely.geometry import Polygon, mapping from shapely import wkt from sqlmodel import Session from app.config import ( MIN_ZONE_AREA_M2, LOOP_CLOSURE_RADIUS_M, DOUGLAS_PEUCKER_EPSILON_M, ) from app.models.activity import Activity from app.models.zone import Zone from app.schemas.activity import GpsPoint, ActivityDetail def haversine_m(lat1: float, lon1: float, lat2: float, lon2: float) -> float: """Haversine distance in meters between two lat/lon points.""" R = 6_371_000 # Earth radius in meters phi1, phi2 = math.radians(lat1), math.radians(lat2) dphi = math.radians(lat2 - lat1) dlam = math.radians(lon2 - lon1) a = ( math.sin(dphi / 2) ** 2 + math.cos(phi1) * math.cos(phi2) * math.sin(dlam / 2) ** 2 ) return 2 * R * math.atan2(math.sqrt(a), math.sqrt(1 - a)) def compute_distance_m(points: list[GpsPoint]) -> float: """Total distance in meters along a GPS track.""" total = 0.0 for i in range(1, len(points)): total += haversine_m( points[i - 1].lat, points[i - 1].lon, points[i].lat, points[i].lon ) return total def simplify_track(points: list[GpsPoint], epsilon: float) -> list[GpsPoint]: """Douglas-Peucker simplification on lat/lon coordinates. epsilon is in approximate degrees; we convert from meters. ~0.00001 degree ≈ 1.11 m at the equator. """ coords = [(p.lat, p.lon) for p in points] eps_deg = epsilon / 111_000 # rough meters->degrees simplified = rdp(coords, epsilon=eps_deg) # Rebuild GpsPoint list (drop timestamps for simplified points) return [GpsPoint(lat=c[0], lon=c[1]) for c in simplified] def is_loop_closed(points: list[GpsPoint], radius_m: float) -> bool: """Check if start and end points are within radius_m meters.""" if len(points) < 4: return False return ( haversine_m(points[0].lat, points[0].lon, points[-1].lat, points[-1].lon) <= radius_m ) def build_polygon(points: list[GpsPoint]) -> Polygon | None: """Build a Shapely Polygon from GPS points. Uses (lon, lat) ordering for GeoJSON/WKT compatibility. Closes the ring if needed. """ coords = [(p.lon, p.lat) for p in points] if coords[0] != coords[-1]: coords.append(coords[0]) if len(coords) < 4: return None poly = Polygon(coords) if not poly.is_valid: poly = poly.buffer(0) # attempt to fix self-intersections return poly if poly.is_valid and not poly.is_empty else None def polygon_area_m2(poly: Polygon) -> float: """Approximate area in m² from a WGS84 polygon. Uses a simple cos(lat) correction. Accurate enough for city-scale zones. """ centroid = poly.centroid lat_rad = math.radians(centroid.y) # Degrees to meters conversion factors m_per_deg_lat = 111_320.0 m_per_deg_lon = 111_320.0 * math.cos(lat_rad) # Scale polygon to meters and compute area coords = list(poly.exterior.coords) scaled = [(x * m_per_deg_lon, y * m_per_deg_lat) for x, y in coords] return abs(Polygon(scaled).area) def wkt_to_geojson(wkt_str: str) -> dict: """Convert WKT polygon string to GeoJSON dict.""" geom = wkt.loads(wkt_str) return mapping(geom) def process_activity( user_id: uuid.UUID, activity_type: str, started_at: datetime, ended_at: datetime, gps_track: list[GpsPoint], session: Session, ) -> ActivityDetail: """Full geo pipeline: filter -> simplify -> validate loop -> build polygon -> persist. Returns an ActivityDetail with zone info if a zone was created. """ # 1. Create activity record distance = compute_distance_m(gps_track) activity = Activity( user_id=user_id, type=activity_type, started_at=started_at, ended_at=ended_at, distance_m=distance, raw_gpx=json.dumps([{"lat": p.lat, "lon": p.lon} for p in gps_track]), status="pending", ) session.add(activity) session.flush() # get activity.id zone_id = None area_m2 = None try: # 2. Filter outliers (skip points with hdop > 4 if present) filtered = [p for p in gps_track if p.hdop is None or p.hdop <= 4] if len(filtered) < 10: activity.status = "failed" session.commit() return _to_detail(activity, None, None) # 3. Simplify simplified = simplify_track(filtered, DOUGLAS_PEUCKER_EPSILON_M) # 4. Check loop closure if not is_loop_closed(simplified, LOOP_CLOSURE_RADIUS_M): activity.status = "failed" session.commit() return _to_detail(activity, None, None) # 5. Build polygon poly = build_polygon(simplified) if poly is None: activity.status = "failed" session.commit() return _to_detail(activity, None, None) # 6. Validate area area_m2 = polygon_area_m2(poly) if area_m2 < MIN_ZONE_AREA_M2: activity.status = "failed" session.commit() return _to_detail(activity, None, None) # 7. Persist zone zone = Zone( owner_id=user_id, activity_id=activity.id, polygon_wkt=poly.wkt, area_m2=area_m2, ) session.add(zone) activity.status = "completed" session.commit() session.refresh(zone) session.refresh(activity) zone_id = zone.id area_m2 = zone.area_m2 except Exception: activity.status = "failed" session.commit() return _to_detail(activity, zone_id, area_m2) def _to_detail( activity: Activity, zone_id: uuid.UUID | None, area_m2: float | None, ) -> ActivityDetail: return ActivityDetail( id=activity.id, user_id=activity.user_id, type=activity.type, started_at=activity.started_at, ended_at=activity.ended_at, distance_m=activity.distance_m, status=activity.status, created_at=activity.created_at, zone_id=zone_id, area_m2=area_m2, )