357 lines
12 KiB
Python
357 lines
12 KiB
Python
# FlightMonitor/map/map_utils.py
|
|
|
|
import logging
|
|
import math
|
|
from typing import Tuple, Optional, List, Set, Dict, Any
|
|
|
|
try:
|
|
import pyproj
|
|
|
|
PYPROJ_MODULE_LOCALLY_AVAILABLE = True
|
|
except ImportError:
|
|
pyproj = None
|
|
PYPROJ_MODULE_LOCALLY_AVAILABLE = False
|
|
logging.warning(
|
|
"MapUtils: 'pyproj' not found. Geographic calculations will be impaired."
|
|
)
|
|
|
|
try:
|
|
import mercantile
|
|
|
|
MERCANTILE_MODULE_LOCALLY_AVAILABLE = True
|
|
except ImportError:
|
|
mercantile = None
|
|
MERCANTILE_MODULE_LOCALLY_AVAILABLE = False
|
|
logging.warning("MapUtils: 'mercantile' not found. Tile operations will fail.")
|
|
|
|
try:
|
|
from flightmonitor.utils.logger import get_logger
|
|
|
|
logger = get_logger(__name__)
|
|
except ImportError:
|
|
logger = logging.getLogger(__name__)
|
|
if not logger.hasHandlers():
|
|
logging.basicConfig(level=logging.DEBUG)
|
|
logger.warning("MapUtils using fallback standard Python logger.")
|
|
|
|
try:
|
|
from flightmonitor.map import map_constants
|
|
except ImportError:
|
|
logger.error("MapUtils: map_constants not found. Using mock constants.")
|
|
|
|
class MockMapConstants:
|
|
DMS_DEGREE_SYMBOL, DMS_MINUTE_SYMBOL, DMS_SECOND_SYMBOL = "°", "'", "''"
|
|
MIN_ZOOM_LEVEL, DEFAULT_MAX_ZOOM_FALLBACK = 0, 20
|
|
GRATICULE_MIN_PIXEL_SPACING = 80
|
|
|
|
map_constants = MockMapConstants()
|
|
|
|
|
|
class MapCalculationError(Exception):
|
|
"""Custom exception for map calculations."""
|
|
|
|
pass
|
|
|
|
|
|
def _is_valid_bbox_dict(bbox_dict: Dict[str, Any]) -> bool:
|
|
if not isinstance(bbox_dict, dict):
|
|
return False
|
|
required_keys = ["lat_min", "lon_min", "lat_max", "lon_max"]
|
|
if not all(key in bbox_dict for key in required_keys):
|
|
return False
|
|
try:
|
|
lat_min, lon_min = float(bbox_dict["lat_min"]), float(bbox_dict["lon_min"])
|
|
lat_max, lon_max = float(bbox_dict["lat_max"]), float(bbox_dict["lon_max"])
|
|
if not (
|
|
-90.0 <= lat_min <= 90.0
|
|
and -90.0 <= lat_max <= 90.0
|
|
and -180.0 <= lon_min <= 180.0
|
|
and -180.0 <= lon_max <= 180.0
|
|
):
|
|
return False
|
|
if lat_min >= lat_max or lon_min >= lon_max:
|
|
return False
|
|
return True
|
|
except (ValueError, TypeError):
|
|
return False
|
|
|
|
|
|
def calculate_graticule_interval(
|
|
span_degrees: float, image_pixel_span: int, min_pixel_spacing: Optional[int] = None
|
|
) -> float:
|
|
"""
|
|
Calculates a reasonable interval for graticule lines based on the map's
|
|
current scale, ensuring a minimum visual separation in pixels.
|
|
"""
|
|
if min_pixel_spacing is None:
|
|
min_pixel_spacing = map_constants.GRATICULE_MIN_PIXEL_SPACING
|
|
|
|
if span_degrees <= 0 or image_pixel_span <= 0:
|
|
return 5.0
|
|
|
|
degrees_per_pixel = span_degrees / image_pixel_span
|
|
min_degree_interval = degrees_per_pixel * min_pixel_spacing
|
|
|
|
possible_intervals = [
|
|
90,
|
|
45,
|
|
30,
|
|
20,
|
|
15,
|
|
10,
|
|
5,
|
|
2,
|
|
1,
|
|
0.5,
|
|
0.25,
|
|
0.1,
|
|
0.05,
|
|
0.02,
|
|
0.01,
|
|
0.005,
|
|
0.002,
|
|
0.001,
|
|
0.0005,
|
|
0.0001,
|
|
]
|
|
|
|
for interval in reversed(possible_intervals):
|
|
if interval >= min_degree_interval:
|
|
return interval
|
|
|
|
return possible_intervals[0]
|
|
|
|
|
|
def get_bounding_box_from_center_size(
|
|
center_latitude_deg: float, center_longitude_deg: float, area_size_km: float
|
|
) -> Optional[Tuple[float, float, float, float]]:
|
|
if not PYPROJ_MODULE_LOCALLY_AVAILABLE or not pyproj:
|
|
return None
|
|
try:
|
|
geod = pyproj.Geod(ellps="WGS84")
|
|
half_side_m = (area_size_km / 2.0) * 1000.0
|
|
_, north_lat, _ = geod.fwd(
|
|
center_longitude_deg, center_latitude_deg, 0.0, half_side_m
|
|
)
|
|
_, south_lat, _ = geod.fwd(
|
|
center_longitude_deg, center_latitude_deg, 180.0, half_side_m
|
|
)
|
|
east_lon, _, _ = geod.fwd(
|
|
center_longitude_deg, center_latitude_deg, 90.0, half_side_m
|
|
)
|
|
west_lon, _, _ = geod.fwd(
|
|
center_longitude_deg, center_latitude_deg, 270.0, half_side_m
|
|
)
|
|
return (west_lon, south_lat, east_lon, north_lat)
|
|
except Exception as e:
|
|
logger.exception(f"Error calculating BBox from center/size: {e}")
|
|
return None
|
|
|
|
|
|
def get_tile_ranges_for_bbox(
|
|
bounding_box_deg: Tuple[float, float, float, float], zoom_level: int
|
|
) -> Optional[Tuple[Tuple[int, int], Tuple[int, int]]]:
|
|
if not MERCANTILE_MODULE_LOCALLY_AVAILABLE or not mercantile:
|
|
return None
|
|
try:
|
|
west, south, east, north = bounding_box_deg
|
|
tiles = list(mercantile.tiles(west, south, east, north, zooms=[zoom_level]))
|
|
if not tiles:
|
|
return None
|
|
x_coords, y_coords = [t.x for t in tiles], [t.y for t in tiles]
|
|
return ((min(x_coords), max(x_coords)), (min(y_coords), max(y_coords)))
|
|
except Exception as e:
|
|
logger.exception(f"Error calculating tile ranges: {e}")
|
|
return None
|
|
|
|
|
|
def calculate_meters_per_pixel(
|
|
latitude_degrees: float, zoom_level: int, tile_pixel_size: int = 256
|
|
) -> Optional[float]:
|
|
try:
|
|
EARTH_CIRCUMFERENCE_METERS = 40075016.686
|
|
clamped_lat = max(-85.05112878, min(85.05112878, latitude_degrees))
|
|
resolution = (
|
|
EARTH_CIRCUMFERENCE_METERS * math.cos(math.radians(clamped_lat))
|
|
) / (tile_pixel_size * (2**zoom_level))
|
|
return resolution
|
|
except Exception as e:
|
|
logger.exception(f"Error calculating meters per pixel: {e}")
|
|
return None
|
|
|
|
|
|
def calculate_geographic_bbox_size_km(
|
|
bounding_box_deg: Tuple[float, float, float, float],
|
|
) -> Optional[Tuple[float, float]]:
|
|
if not PYPROJ_MODULE_LOCALLY_AVAILABLE or not pyproj:
|
|
return None
|
|
try:
|
|
west, south, east, north = bounding_box_deg
|
|
geod = pyproj.Geod(ellps="WGS84")
|
|
_, _, width_m = geod.inv(west, (south + north) / 2, east, (south + north) / 2)
|
|
_, _, height_m = geod.inv((west + east) / 2, south, (west + east) / 2, north)
|
|
return (abs(width_m) / 1000.0, abs(height_m) / 1000.0)
|
|
except Exception as e:
|
|
logger.exception(f"Error calculating BBox size: {e}")
|
|
return None
|
|
|
|
|
|
def calculate_zoom_level_for_geographic_size(
|
|
latitude_degrees: float,
|
|
geographic_span_meters: float,
|
|
target_pixel_span: int,
|
|
tile_pixel_size: int = 256,
|
|
) -> Optional[int]:
|
|
try:
|
|
if target_pixel_span <= 0:
|
|
return None
|
|
required_res = geographic_span_meters / target_pixel_span
|
|
EARTH_CIRCUMFERENCE_METERS = 40075016.686
|
|
clamped_lat = max(-85.05112878, min(85.05112878, latitude_degrees))
|
|
term = (EARTH_CIRCUMFERENCE_METERS * math.cos(math.radians(clamped_lat))) / (
|
|
tile_pixel_size * required_res
|
|
)
|
|
if term <= 0:
|
|
return None
|
|
precise_zoom = math.log2(term)
|
|
return max(
|
|
map_constants.MIN_ZOOM_LEVEL,
|
|
min(int(round(precise_zoom)), map_constants.DEFAULT_MAX_ZOOM_FALLBACK),
|
|
)
|
|
except Exception as e:
|
|
logger.exception(f"Error calculating zoom level: {e}")
|
|
return None
|
|
|
|
|
|
def deg_to_dms_string(degree_value: float, coord_type: str) -> str:
|
|
if not isinstance(degree_value, (int, float)) or not math.isfinite(degree_value):
|
|
return "N/A"
|
|
abs_val = abs(degree_value)
|
|
degrees = int(abs_val)
|
|
minutes = int((abs_val - degrees) * 60)
|
|
seconds = (((abs_val - degrees) * 60) - minutes) * 60
|
|
if coord_type.lower() == "lat":
|
|
suffix = "N" if degree_value >= 0 else "S"
|
|
elif coord_type.lower() == "lon":
|
|
suffix = "E" if degree_value >= 0 else "W"
|
|
else:
|
|
suffix = ""
|
|
return (
|
|
f"{degrees}{map_constants.DMS_DEGREE_SYMBOL} "
|
|
f"{minutes:02d}{map_constants.DMS_MINUTE_SYMBOL} "
|
|
f"{seconds:05.2f}{map_constants.DMS_SECOND_SYMBOL} {suffix}"
|
|
)
|
|
|
|
|
|
def calculate_geographic_bbox_from_pixel_size_and_zoom(
|
|
center_latitude_deg: float,
|
|
center_longitude_deg: float,
|
|
target_pixel_width: int,
|
|
target_pixel_height: int,
|
|
zoom_level: int,
|
|
tile_pixel_size: int = 256,
|
|
) -> Optional[Tuple[float, float, float, float]]:
|
|
if not (MERCANTILE_MODULE_LOCALLY_AVAILABLE and mercantile):
|
|
return None
|
|
try:
|
|
res = calculate_meters_per_pixel(
|
|
center_latitude_deg, zoom_level, tile_pixel_size
|
|
)
|
|
if not res:
|
|
return None
|
|
center_merc_x, center_merc_y = mercantile.xy(
|
|
center_longitude_deg, center_latitude_deg
|
|
)
|
|
half_w_merc, half_h_merc = (target_pixel_width / 2) * res, (
|
|
target_pixel_height / 2
|
|
) * res
|
|
ul_x, ul_y = center_merc_x - half_w_merc, center_merc_y + half_h_merc
|
|
lr_x, lr_y = center_merc_x + half_w_merc, center_merc_y - half_h_merc
|
|
west, north = mercantile.lnglat(ul_x, ul_y)
|
|
east, south = mercantile.lnglat(lr_x, lr_y)
|
|
return (west, south, east, north)
|
|
except Exception as e:
|
|
logger.exception(f"Error calculating geo BBox from pixel size/zoom: {e}")
|
|
return None
|
|
|
|
|
|
def geo_to_pixel_on_unscaled_map(
|
|
latitude_deg: float,
|
|
longitude_deg: float,
|
|
current_map_geo_bounds: Optional[Tuple[float, float, float, float]],
|
|
current_stitched_map_pixel_shape: Optional[Tuple[int, int]],
|
|
) -> Optional[Tuple[int, int]]:
|
|
if not (MERCANTILE_MODULE_LOCALLY_AVAILABLE and mercantile):
|
|
return None
|
|
if not (current_map_geo_bounds and current_stitched_map_pixel_shape):
|
|
return None
|
|
|
|
img_w, img_h = current_stitched_map_pixel_shape
|
|
if img_w <= 0 or img_h <= 0:
|
|
return None
|
|
|
|
map_west, map_south, map_east, map_north = current_map_geo_bounds
|
|
|
|
try:
|
|
map_ul_merc_x, map_ul_merc_y = mercantile.xy(map_west, map_north)
|
|
map_lr_merc_x, map_lr_merc_y = mercantile.xy(map_east, map_south)
|
|
merc_w, merc_h = abs(map_lr_merc_x - map_ul_merc_x), abs(
|
|
map_ul_merc_y - map_lr_merc_y
|
|
)
|
|
if merc_w < 1e-9 or merc_h < 1e-9:
|
|
return None
|
|
|
|
target_merc_x, target_merc_y = mercantile.xy(longitude_deg, latitude_deg)
|
|
|
|
rel_x = (target_merc_x - map_ul_merc_x) / merc_w
|
|
rel_y = (map_ul_merc_y - target_merc_y) / merc_h
|
|
|
|
return (int(round(rel_x * img_w)), int(round(rel_y * img_h)))
|
|
except Exception as e:
|
|
logger.error(f"Error converting geo to pixel: {e}", exc_info=False)
|
|
return None
|
|
|
|
|
|
def pixel_to_geo(
|
|
canvas_x: int,
|
|
canvas_y: int,
|
|
current_map_geo_bounds: Optional[Tuple[float, float, float, float]],
|
|
current_stitched_map_pixel_shape: Optional[Tuple[int, int]],
|
|
) -> Tuple[Optional[float], Optional[float]]:
|
|
if not (MERCANTILE_MODULE_LOCALLY_AVAILABLE and mercantile):
|
|
return None, None
|
|
if not (current_map_geo_bounds and current_stitched_map_pixel_shape):
|
|
return None, None
|
|
|
|
img_w, img_h = current_stitched_map_pixel_shape
|
|
if img_w <= 0 or img_h <= 0:
|
|
return None, None
|
|
|
|
map_west, map_south, map_east, map_north = current_map_geo_bounds
|
|
|
|
try:
|
|
map_ul_merc_x, map_ul_merc_y = mercantile.xy(
|
|
map_west, map_north, truncate=False
|
|
)
|
|
map_lr_merc_x, map_lr_merc_y = mercantile.xy(
|
|
map_east, map_south, truncate=False
|
|
)
|
|
|
|
total_merc_w, total_merc_h = abs(map_lr_merc_x - map_ul_merc_x), abs(
|
|
map_ul_merc_y - map_lr_merc_y
|
|
)
|
|
if total_merc_w < 1e-9 or total_merc_h < 1e-9:
|
|
return None, None
|
|
|
|
rel_x, rel_y = canvas_x / img_w, canvas_y / img_h
|
|
|
|
target_merc_x = map_ul_merc_x + (rel_x * total_merc_w)
|
|
target_merc_y = map_ul_merc_y - (rel_y * total_merc_h)
|
|
|
|
lon, lat = mercantile.lnglat(target_merc_x, target_merc_y, truncate=False)
|
|
|
|
return lon, lat
|
|
except Exception as e:
|
|
logger.error(f"Error converting pixel to geo: {e}", exc_info=True)
|
|
return None, None
|