GeoSight MCP Server

""" Geospatial utility functions. Provides coordinate transformations, bounding box calculations, and area computations for satellite imagery analysis. """ import math from dataclasses import dataclass from typing import Tuple, List, Optional import numpy as np @dataclass class BoundingBox: """Geographic bounding box.""" west: float # Min longitude south: float # Min latitude east: float # Max longitude north: float # Max latitude @property def center(self) -> Tuple[float, float]: """Get center point (lat, lon).""" return ( (self.south + self.north) / 2, (self.west + self.east) / 2, ) @property def width_deg(self) -> float: """Width in degrees.""" return self.east - self.west @property def height_deg(self) -> float: """Height in degrees.""" return self.north - self.south def to_list(self) -> List[float]: """Convert to [west, south, east, north] list.""" return [self.west, self.south, self.east, self.north] def to_tuple(self) -> Tuple[float, float, float, float]: """Convert to (west, south, east, north) tuple.""" return (self.west, self.south, self.east, self.north) def to_polygon_coords(self) -> List[Tuple[float, float]]: """Convert to polygon coordinates [(lon, lat), ...].""" return [ (self.west, self.south), (self.east, self.south), (self.east, self.north), (self.west, self.north), (self.west, self.south), ] def buffer(self, degrees: float) -> "BoundingBox": """Create a new bbox buffered by given degrees.""" return BoundingBox( west=self.west - degrees, south=max(-90, self.south - degrees), east=self.east + degrees, north=min(90, self.north + degrees), ) # Earth's radius in kilometers EARTH_RADIUS_KM = 6371.0 def create_bbox_from_point( latitude: float, longitude: float, radius_km: float, ) -> BoundingBox: """ Create a bounding box around a center point. Args: latitude: Center latitude (-90 to 90) longitude: Center longitude (-180 to 180) radius_km: Radius in kilometers Returns: BoundingBox object """ # Validate inputs if not -90 <= latitude <= 90: raise ValueError(f"Latitude must be between -90 and 90, got {latitude}") if not -180 <= longitude <= 180: raise ValueError(f"Longitude must be between -180 and 180, got {longitude}") if radius_km <= 0: raise ValueError(f"Radius must be positive, got {radius_km}") # Convert to radians lat_rad = math.radians(latitude) # Calculate degree offsets # Latitude: 1 degree ≈ 111 km (constant) lat_offset = radius_km / 111.0 # Longitude: varies with latitude (cos(lat) factor) # At equator: 1 degree ≈ 111 km # At poles: approaches 0 lon_offset = radius_km / (111.0 * math.cos(lat_rad)) # Handle edge cases near poles if abs(latitude) > 89: lon_offset = 180.0 # Cover all longitudes near poles return BoundingBox( west=max(-180, longitude - lon_offset), south=max(-90, latitude - lat_offset), east=min(180, longitude + lon_offset), north=min(90, latitude + lat_offset), ) def calculate_area(bbox: BoundingBox) -> float: """ Calculate approximate area of a bounding box in square kilometers. Uses the spherical Earth approximation. Args: bbox: Bounding box Returns: Area in square kilometers """ # Convert to radians lat1 = math.radians(bbox.south) lat2 = math.radians(bbox.north) lon1 = math.radians(bbox.west) lon2 = math.radians(bbox.east) # Calculate area using spherical geometry # Area = R^2 * |sin(lat1) - sin(lat2)| * |lon1 - lon2| area = (EARTH_RADIUS_KM ** 2) * abs(math.sin(lat1) - math.sin(lat2)) * abs(lon2 - lon1) return area def haversine_distance( lat1: float, lon1: float, lat2: float, lon2: float, ) -> float: """ Calculate the great circle distance between two points. Args: lat1, lon1: First point coordinates lat2, lon2: Second point coordinates Returns: Distance in kilometers """ # Convert to radians lat1, lon1 = math.radians(lat1), math.radians(lon1) lat2, lon2 = math.radians(lat2), math.radians(lon2) # Haversine formula dlat = lat2 - lat1 dlon = lon2 - lon1 a = math.sin(dlat / 2) ** 2 + math.cos(lat1) * math.cos(lat2) * math.sin(dlon / 2) ** 2 c = 2 * math.asin(math.sqrt(a)) return EARTH_RADIUS_KM * c def pixel_to_coordinates( pixel_x: int, pixel_y: int, bbox: BoundingBox, image_width: int, image_height: int, ) -> Tuple[float, float]: """ Convert pixel coordinates to geographic coordinates. Args: pixel_x, pixel_y: Pixel coordinates bbox: Image bounding box image_width, image_height: Image dimensions Returns: Tuple of (latitude, longitude) """ # Calculate geographic coordinates lon = bbox.west + (pixel_x / image_width) * bbox.width_deg lat = bbox.north - (pixel_y / image_height) * bbox.height_deg return lat, lon def coordinates_to_pixel( latitude: float, longitude: float, bbox: BoundingBox, image_width: int, image_height: int, ) -> Tuple[int, int]: """ Convert geographic coordinates to pixel coordinates. Args: latitude, longitude: Geographic coordinates bbox: Image bounding box image_width, image_height: Image dimensions Returns: Tuple of (pixel_x, pixel_y) """ pixel_x = int((longitude - bbox.west) / bbox.width_deg * image_width) pixel_y = int((bbox.north - latitude) / bbox.height_deg * image_height) # Clamp to valid range pixel_x = max(0, min(image_width - 1, pixel_x)) pixel_y = max(0, min(image_height - 1, pixel_y)) return pixel_x, pixel_y def meters_per_pixel(latitude: float, zoom_level: int) -> float: """ Calculate meters per pixel at a given latitude and zoom level. Based on Web Mercator projection (EPSG:3857). Args: latitude: Latitude in degrees zoom_level: Map zoom level (0-22) Returns: Meters per pixel """ # At zoom 0, world is 256 pixels # Resolution halves with each zoom level meters_per_pixel_equator = 156543.03392 / (2 ** zoom_level) # Adjust for latitude return meters_per_pixel_equator * math.cos(math.radians(latitude)) def calculate_optimal_resolution(bbox: BoundingBox, max_pixels: int = 2048) -> int: """ Calculate optimal resolution for a given bounding box. Args: bbox: Bounding box max_pixels: Maximum pixels on longest side Returns: Resolution in meters """ # Calculate actual dimensions center_lat = (bbox.south + bbox.north) / 2 width_km = haversine_distance(center_lat, bbox.west, center_lat, bbox.east) height_km = haversine_distance(bbox.south, bbox.west, bbox.north, bbox.west) max_dim_km = max(width_km, height_km) max_dim_m = max_dim_km * 1000 # Calculate resolution resolution = max_dim_m / max_pixels # Round to standard resolutions standard_resolutions = [10, 20, 30, 60, 100, 250, 500, 1000] for std_res in standard_resolutions: if resolution <= std_res: return std_res return standard_resolutions[-1] def reproject_bbox( bbox: BoundingBox, from_crs: str = "EPSG:4326", to_crs: str = "EPSG:3857", ) -> BoundingBox: """ Reproject a bounding box to a different coordinate reference system. Note: This is a simplified implementation. For production use, consider using pyproj for accurate transformations. Args: bbox: Input bounding box from_crs: Source CRS (default: WGS84) to_crs: Target CRS (default: Web Mercator) Returns: Reprojected bounding box """ if from_crs == to_crs: return bbox if from_crs == "EPSG:4326" and to_crs == "EPSG:3857": # WGS84 to Web Mercator def to_mercator(lon: float, lat: float) -> Tuple[float, float]: x = lon * 20037508.34 / 180 y = math.log(math.tan((90 + lat) * math.pi / 360)) / (math.pi / 180) y = y * 20037508.34 / 180 return x, y west, south = to_mercator(bbox.west, bbox.south) east, north = to_mercator(bbox.east, bbox.north) return BoundingBox(west=west, south=south, east=east, north=north) raise NotImplementedError(f"Reprojection from {from_crs} to {to_crs} not implemented")