feat: stage5 — Satellite tiles (F04) and Coordinates (F13)

src/gps_denied/utils/mercator.py

@@ -0,0 +1,48 @@
+"""Web Mercator utility functions (Component H06)."""
+
+import math
+
+from gps_denied.schemas import GPSPoint
+from gps_denied.schemas.satellite import TileBounds, TileCoords
+
+
+def latlon_to_tile(lat: float, lon: float, zoom: int) -> TileCoords:
+    """Convert GPS coordinates to Web Mercator tile coordinates."""
+    n = 2.0 ** zoom
+    x = int((lon + 180.0) / 360.0 * n)
+    lat_rad = math.radians(lat)
+    y = int((1.0 - math.log(math.tan(lat_rad) + (1.0 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
+    return TileCoords(x=x, y=y, zoom=zoom)
+
+
+def tile_to_latlon(x: float, y: float, zoom: int) -> GPSPoint:
+    """Convert tile coordinates (can be fractional) back to GPS WGS84."""
+    n = 2.0 ** zoom
+    lon_deg = x / n * 360.0 - 180.0
+    lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * y / n)))
+    lat_deg = math.degrees(lat_rad)
+    return GPSPoint(lat=lat_deg, lon=lon_deg)
+
+
+def compute_tile_bounds(coords: TileCoords) -> TileBounds:
+    """Compute the GPS bounds and GSD for a given tile."""
+    nw = tile_to_latlon(coords.x, coords.y, coords.zoom)
+    se = tile_to_latlon(coords.x + 1, coords.y + 1, coords.zoom)
+    center = tile_to_latlon(coords.x + 0.5, coords.y + 0.5, coords.zoom)
+    
+    ne = GPSPoint(lat=nw.lat, lon=se.lon)
+    sw = GPSPoint(lat=se.lat, lon=nw.lon)
+    
+    # Calculate GSD (meters per pixel at this latitude)
+    # Assumes standard 256x256 Web Mercator tile
+    lat_rad = math.radians(center.lat)
+    gsd = 156543.03392 * math.cos(lat_rad) / (2 ** coords.zoom)
+    
+    return TileBounds(
+        nw=nw,
+        ne=ne,
+        sw=sw,
+        se=se,
+        center=center,
+        gsd=gsd
+    )