VALLONGOL/SXXXXXXX_FlightMonitor
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
029b3dc5e2
SXXXXXXX_FlightMonitor/flightmonitor/map/map_utils.py

726 lines
36 KiB
Python
Raw Blame History

# 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
import logging
logging.warning("MapUtils: 'pyproj' not found.")
try:
import mercantile
MERCANTILE_MODULE_LOCALLY_AVAILABLE = True
except ImportError:
mercantile = None
MERCANTILE_MODULE_LOCALLY_AVAILABLE = False
import logging
logging.warning("MapUtils: 'mercantile' not found.")
try:
from ..utils.logger import get_logger
from . import map_constants
except ImportError:
logger = logging.getLogger(__name__)
logger.warning("MapUtils: App logger or map_constants not available.")
class MockMapConstants:
COORDINATE_DECIMAL_PLACES = 4
DMS_DEGREE_SYMBOL = " deg"
DMS_MINUTE_SYMBOL = " min"
DMS_SECOND_SYMBOL = " sec"
map_constants = MockMapConstants()
logger = get_logger(__name__)
class MapCalculationError(Exception):
"""Custom exception for map calculations."""
pass
def _is_valid_bbox_dict(bbox_dict: Dict[str, Any]) -> bool:
"""Validate bounding box dict structure and ranges."""
if not isinstance(bbox_dict, dict):
logger.debug("BBox validation failed: not 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):
logger.debug(f"BBox validation failed: missing keys {required_keys}.")
return False
try:
lat_min = float(bbox_dict["lat_min"])
lon_min = float(bbox_dict["lon_min"])
lat_max = float(bbox_dict["lat_max"])
lon_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):
logger.debug(f"BBox validation failed: coordinates out of range {bbox_dict}.")
return False
if lat_min >= lat_max:
logger.debug(f"BBox validation failed: lat_min ({lat_min}) >= lat_max ({lat_max}).")
return False
# Note: For services not supporting antimeridian crossing, lon_min must be < lon_max.
# Current validation keeps this check for OpenSky compatibility.
if lon_min >= lon_max:
logger.debug(f"BBox validation failed: lon_min ({lon_min}) >= lon_max ({lon_max}).")
return False
if not all(math.isfinite(coord) for coord in [lat_min, lon_min, lat_max, lon_max]):
logger.debug(f"BBox validation failed: non-finite values {bbox_dict}.")
return False
return True
except (ValueError, TypeError) as e:
logger.debug(f"BBox validation failed: invalid types/values {bbox_dict}. Error: {e}")
return False
except Exception as e:
logger.error(f"Unexpected error validating BBox: {bbox_dict}. Error: {e}", exc_info=True)
return False
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]
]:
"""Calc WGS84 BBox from center and square size (km). Requires pyproj."""
if not PYPROJ_MODULE_LOCALLY_AVAILABLE or pyproj is None:
logger.error("'pyproj' required for BBox from center/size.")
return None
if not (isinstance(area_size_km, (int, float)) and area_size_km > 0):
logger.error(f"Invalid area_size_km: {area_size_km}.")
return None
if not (isinstance(center_latitude_deg, (int, float)) and -90.0 <= center_latitude_deg <= 90.0):
logger.error(f"Invalid center_latitude_deg: {center_latitude_deg}.")
return None
if not (isinstance(center_longitude_deg, (int, float)) and -180.0 <= center_longitude_deg <= 180.0):
logger.error(f"Invalid center_longitude_deg: {center_longitude_deg}.")
return None
try:
geodetic_calculator = pyproj.Geod(ellps="WGS84")
half_side_m = (area_size_km / 2.0) * 1000.0
_, north_lat, _ = geodetic_calculator.fwd(center_longitude_deg, center_latitude_deg, 0.0, half_side_m)
_, south_lat, _ = geodetic_calculator.fwd(center_longitude_deg, center_latitude_deg, 180.0, half_side_m)
east_lon, _, _ = geodetic_calculator.fwd(center_longitude_deg, center_latitude_deg, 90.0, half_side_m)
west_lon, _, _ = geodetic_calculator.fwd(center_longitude_deg, center_latitude_deg, 270.0, half_side_m)
north_lat = min(90.0, max(-90.0, north_lat))
south_lat = min(90.0, max(-90.0, south_lat))
if not all(math.isfinite(coord) for coord in [west_lon, south_lat, east_lon, north_lat]):
logger.error("Calculated BBox contains non-finite values.")
return None
return (west_lon, south_lat, east_lon, north_lat)
except Exception as e_bbox_calc:
logger.exception(f"Error calculating BBox from center/size: {e_bbox_calc}")
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]]
]:
"""Calc X/Y tile ranges covering BBox at zoom. Requires mercantile."""
if not MERCANTILE_MODULE_LOCALLY_AVAILABLE or mercantile is None:
logger.error("'mercantile' required for tile ranges.")
return None
west_lon, south_lat, east_lon, north_lat = bounding_box_deg
clamped_south_lat = max(-85.05112877, min(85.05112877, south_lat))
clamped_north_lat = max(-85.05112877, min(85.05112877, north_lat))
if clamped_south_lat >= clamped_north_lat:
logger.warning(f"Clamped lat range invalid: S={clamped_south_lat:.4f}, N={clamped_north_lat:.4f}. Cannot get tile ranges.")
return None
if not (isinstance(zoom_level, int) and 0 <= zoom_level <= 25):
logger.error(f"Invalid zoom level: {zoom_level}.")
return None
try:
# mercantile.tiles handles longitude wrapping automatically
tiles_in_bbox_generator = mercantile.tiles(west_lon, clamped_south_lat, east_lon, clamped_north_lat, zooms=[zoom_level])
list_of_tiles = list(tiles_in_bbox_generator)
if not list_of_tiles:
logger.warning(f"No tiles found by mercantile for BBox at zoom {zoom_level}. Using fallback: tile at BBox center.")
center_lon = (west_lon + east_lon) / 2.0
# Handle antimeridian crossing for center longitude
if west_lon > east_lon:
center_lon = (west_lon + east_lon + 360) / 2.0
if center_lon > 180: center_lon -= 360
center_lat = (south_lat + north_lat) / 2.0
center_lat = max(-85.05112877, min(85.05112877, center_lat))
try:
center_tile = mercantile.tile(center_lon, center_lat, zoom_level)
min_x, max_x = center_tile.x, center_tile.x
min_y, max_y = center_tile.y, center_tile.y
logger.debug(f"Fallback: Using tile at BBox center ({center_lon:.4f}, {center_lat:.4f}): ({min_x}, {min_y}, {zoom_level}).")
except Exception as e_center_tile:
logger.error(f"Fallback failed: Error getting tile for BBox center: {e_center_tile}", exc_info=True)
return None
else:
x_coordinates = [tile.x for tile in list_of_tiles]
y_coordinates = [tile.y for tile in list_of_tiles]
min_x, max_x = min(x_coordinates), max(x_coordinates)
min_y, max_y = min(y_coordinates), max(y_coordinates)
return ((min_x, max_x), (min_y, max_y))
except Exception as e_tile_range_calc:
logger.exception(f"Error calculating tile ranges: {e_tile_range_calc}")
return None
def calculate_meters_per_pixel(
latitude_degrees: float, zoom_level: int, tile_pixel_size: int = 256
) -> Optional[float]:
"""Calc approx ground resolution (m/px) at lat/zoom (Web Mercator)."""
try:
if not (isinstance(latitude_degrees, (int, float)) and -90.0 <= latitude_degrees <= 90.0):
logger.warning(f"Invalid latitude for m/px calc: {latitude_degrees}")
return None
if not (isinstance(zoom_level, int) and 0 <= zoom_level <= 25):
logger.warning(f"Invalid zoom level for m/px calc: {zoom_level}")
return None
if not (isinstance(tile_pixel_size, int) and tile_pixel_size > 0):
logger.warning(f"Invalid tile_pixel_size for m/px calc: {tile_pixel_size}")
return None
EARTH_CIRCUMFERENCE_METERS = 40075016.686
clamped_lat = max(-85.05112878, min(85.05112878, latitude_degrees))
lat_radians = math.radians(clamped_lat)
denominator = tile_pixel_size * (2**zoom_level)
if denominator <= 1e-9:
logger.warning("Denominator for m/px calc is zero/near-zero.")
return None
resolution_m_px = (EARTH_CIRCUMFERENCE_METERS * math.cos(lat_radians)) / denominator
if not math.isfinite(resolution_m_px) or resolution_m_px <= 0:
logger.warning(f"Calculated non-finite/non-positive m/px ({resolution_m_px}) at Lat {latitude_degrees}.")
return None
return resolution_m_px
except Exception as e_mpp_calc:
logger.exception(f"Error calculating meters per pixel: {e_mpp_calc}")
return None
def calculate_geographic_bbox_size_km(
bounding_box_deg: Tuple[float, float, float, float],
) -> Optional[Tuple[float, float]]:
"""Calc approx W/H (km) of BBox. Requires pyproj."""
if not PYPROJ_MODULE_LOCALLY_AVAILABLE or pyproj is None:
logger.error("'pyproj' required for BBox size.")
return None
west_lon, south_lat, east_lon, north_lat = bounding_box_deg
if not (isinstance(west_lon, (int, float)) and isinstance(south_lat, (int, float)) and isinstance(east_lon, (int, float)) and isinstance(north_lat, (int, float))):
logger.warning(f"Invalid coordinate types for size calc: {bounding_box_deg}.")
return None
if not (-90.0 <= south_lat <= north_lat <= 90.0):
logger.warning(f"Invalid lat range for size calc: {south_lat}, {north_lat}. Clamping.")
south_lat = max(-90.0, south_lat); north_lat = min(90.0, north_lat)
if south_lat >= north_lat: logger.error(f"Invalid lat range after clamping: {south_lat}, {north_lat}."); return None
try:
geod = pyproj.Geod(ellps="WGS84")
center_lat = (south_lat + north_lat) / 2.0
center_lat = max(-89.9, min(89.9, center_lat))
# Calculate width along the center latitude. Handles antimeridian crossing correctly in pyproj.Geod.inv.
_, _, width_m = geod.inv(west_lon, center_lat, east_lon, center_lat)
width_m = abs(width_m)
# Calculate height along the center longitude (or 0/180 meridian if BBox crosses).
# For simplicity, using center longitude here.
center_lon_for_height = (west_lon + east_lon) / 2.0
if west_lon > east_lon:
center_lon_for_height = (west_lon + east_lon + 360) / 2.0
if center_lon_for_height > 180: center_lon_for_height -= 360
center_lon_for_height = max(-179.9, min(179.9, center_lon_for_height)) # Clamp for geod.inv
_, _, height_m = geod.inv(center_lon_for_height, south_lat, center_lon_for_height, north_lat)
height_m = abs(height_m)
approx_width_km = width_m / 1000.0
approx_height_km = height_m / 1000.0
if (approx_width_km <= 1e-9 or approx_height_km <= 1e-9):
logger.warning(f"Calc non-positive W/H for BBox {bounding_box_deg}: ({approx_width_km:.6f}, {approx_height_km:.6f}).")
return None
if not all(math.isfinite(val) for val in [approx_width_km, approx_height_km]):
logger.warning(f"Calc non-finite W/H: ({approx_width_km}, {approx_height_km}).")
return None
return (approx_width_km, approx_height_km)
except Exception as e_size_calc:
logger.exception(f"Error calculating BBox size: {e_size_calc}")
return None
def get_hgt_tile_geographic_bounds(
lat_coord: int, lon_coord: int
) -> Optional[Tuple[float, float, float, float]]:
"""Calc WGS84 BBox (W,S,E,N) for HGT tile integer coords (S, W)."""
if not (isinstance(lat_coord, int) and isinstance(lon_coord, int)):
logger.warning(f"Invalid input types for HGT bounds: {type(lat_coord)}, {type(lon_coord)}.")
return None
# Note: HGT tile coordinates are typically S/W corner integers, e.g., N007E008 is lat_coord=7, lon_coord=8.
# But ranges are [-90, 89] for lat, [-180, 179] for lon based on standard grid.
if not (-90 <= lat_coord <= 89): logger.warning(f"Lat coord {lat_coord} outside expected HGT range [-90, 89].")
if not (-180 <= lon_coord <= 179): logger.warning(f"Lon coord {lon_coord} outside expected HGT range [-180, 179].")
west_lon = float(lon_coord)
south_lat = float(lat_coord)
east_lon = float(lon_coord + 1)
north_lat = float(lat_coord + 1)
# Clamp to valid WGS84 ranges defensively
west_lon = max(-180.0, west_lon); south_lat = max(-90.0, south_lat); east_lon = min(180.0, east_lon); north_lat = min(90.0, north_lat)
if south_lat >= north_lat: logger.warning(f"Calc HGT bounds have invalid lat range: S={south_lat}, N={north_lat}."); return None
# Longitude wrapping for a single HGT tile shouldn't happen unless input coords are out of range initially.
# The check `west_lon > east_lon` is only meaningful if not crossing antimeridian.
# For HGT tiles, they represent 1x1 degree squares not crossing antimeridian.
# if west_lon > east_lon and abs(west_lon - east_lon) > 1e-9: logger.warning(f"Calc HGT bounds have W > E: W={west_lon}, E={east_lon}."); return None
return (west_lon, south_lat, east_lon, north_lat)
def calculate_zoom_level_for_geographic_size(
latitude_degrees: float,
geographic_height_meters: float,
target_pixel_height: int,
tile_pixel_size: int = 256,
) -> Optional[int]:
"""Calc Web Mercator zoom for geographic height to fit pixel height."""
if not (isinstance(latitude_degrees, (int, float)) and -90.0 <= latitude_degrees <= 90.0):
logger.warning(f"Invalid latitude for zoom calc: {latitude_degrees}")
return None
if not (isinstance(geographic_height_meters, (int, float)) and geographic_height_meters > 0):
logger.warning(f"Invalid geographic_height_meters: {geographic_height_meters}.")
return None
if not (isinstance(target_pixel_height, int) and target_pixel_height > 0):
logger.warning(f"Invalid target_pixel_height: {target_pixel_height}.")
return None
if not (isinstance(tile_pixel_size, int) and tile_pixel_size > 0):
logger.warning(f"Invalid tile_pixel_size: {tile_pixel_size}.")
return None
try:
req_res_m_px = geographic_height_meters / target_pixel_height
if req_res_m_px <= 1e-9 or not math.isfinite(req_res_m_px):
logger.warning(f"Calc non-positive/non-finite required resolution ({req_res_m_px}).")
return None
EARTH_CIRCUMFERENCE_METERS = 40075016.686
clamped_lat = max(-85.05112878, min(85.05112878, latitude_degrees))
lat_radians = math.radians(clamped_lat)
cos_lat = math.cos(lat_radians)
if abs(cos_lat) < 1e-9:
logger.warning(f"Latitude {latitude_degrees} too close to pole for reliable zoom calc.")
return None
if tile_pixel_size <= 0: logger.error("Invalid tile_pixel_size (<= 0)."); return None
# res_m_px = (EARTH_CIRCUMFERENCE_METERS * cos_lat) / (tile_pixel_size * 2^Z)
# 2^Z = (EARTH_CIRCUMFERENCE_METERS * cos_lat) / (tile_pixel_size * req_res_m_px)
# Z = log2(...)
term_for_log = (EARTH_CIRCUMFERENCE_METERS * cos_lat) / (tile_pixel_size * req_res_m_px)
if term_for_log <= 1e-9: logger.warning(f"Calc non-positive term for log2 ({term_for_log:.2e})."); return None
precise_zoom = math.log2(term_for_log)
integer_zoom = int(round(precise_zoom))
clamped_zoom = max(0, min(integer_zoom, 20)) # Clamp to 0-20 range (arbitrary upper cap, services vary)
return clamped_zoom
except Exception as e_zoom_calc:
logger.exception(f"Error calculating zoom level: {e_zoom_calc}")
return None
def deg_to_dms_string(degree_value: float, coord_type: str) -> str:
"""Convert decimal degrees to DMS string."""
if not isinstance(degree_value, (int, float)) or not math.isfinite(degree_value): return "N/A"
if coord_type.lower() not in ("lat", "lon"): logger.warning(f"Unknown coord type '{coord_type}'."); return "N/A (Invalid Type)"
val_for_suffix = degree_value
if coord_type.lower() == "lat":
if not -90.0 <= degree_value <= 90.0:
logger.warning(f"Lat {degree_value} outside range [-90, 90] for suffix. Clamping.")
val_for_suffix = max(-90.0, min(90.0, degree_value))
elif coord_type.lower() == "lon":
val_for_suffix = (degree_value + 180) % 360 - 180
if val_for_suffix == -180.0 and degree_value != -180.0: val_for_suffix = 180.0
abs_deg = abs(degree_value)
degrees = int(abs_deg)
minutes_decimal = (abs_deg - degrees) * 60
minutes = int(minutes_decimal)
seconds = (minutes_decimal - minutes) * 60
# Round seconds and carry over if necessary
seconds = round(seconds, 2)
if seconds >= 60.0: # Can happen after rounding, e.g., 59.996 rounds to 60.00
seconds = 0.0
minutes += 1
if minutes >= 60:
minutes = 0
degrees += 1
# Handle degrees reaching 90 or 180 due to rounding, for suffix logic
if coord_type.lower() == "lat" and degrees > 90: degrees = 90
elif coord_type.lower() == "lon" and degrees > 180: degrees = 180
suffix = ""
if coord_type.lower() == "lat": suffix = "N" if val_for_suffix >= 0 else "S"
elif coord_type.lower() == "lon":
suffix = "E" if val_for_suffix >= 0 else "W"
if val_for_suffix == 0.0: suffix = "" # 0 longitude doesn't need E/W
try: # Use constants for symbols
dms_string = f"{degrees}{map_constants.DMS_DEGREE_SYMBOL} {minutes}{map_constants.DMS_MINUTE_SYMBOL} {seconds:.2f}{map_constants.DMS_SECOND_SYMBOL}"
except AttributeError: logger.warning("map_constants missing DMS symbols, using plain text."); dms_string = f"{degrees} Deg {minutes} Min {seconds:.2f} Sec"
except Exception as e: logger.error(f"Error formatting DMS string: {e}. Plain text.", exc_info=False); dms_string = f"{degrees} Deg {minutes} Min {seconds:.2f} Sec"
if suffix: dms_string += f" {suffix}"
return dms_string
def get_combined_geographic_bounds_from_tile_info_list(
tile_info_list: List[Dict],
) -> Optional[Tuple[float, float, float, float]]:
"""Calc min BBox encompassing all tiles in list."""
if not tile_info_list:
logger.warning("Tile info list empty for combined bounds.")
return None
min_lon, min_lat, max_lon, max_lat = 180.0, 90.0, -180.0, -90.0
initialized = False
for tile_info in tile_info_list:
lat_coord = tile_info.get("latitude_coord")
lon_coord = tile_info.get("longitude_coord")
if lat_coord is None or lon_coord is None:
logger.warning(f"Skipping tile info: missing coords {tile_info}.")
continue
try:
tile_bounds = get_hgt_tile_geographic_bounds(int(lat_coord), int(lon_coord))
if tile_bounds:
if not initialized:
min_lon, min_lat, max_lon, max_lat = tile_bounds
initialized = True
else:
min_lon = min(min_lon, tile_bounds[0])
min_lat = min(min_lat, tile_bounds[1])
max_lon = max(max_lon, tile_bounds[2])
max_lat = max(max_lat, tile_bounds[3])
else: logger.warning(f"Could not get geo bounds for tile ({lat_coord},{lon_coord}), skipping.")
except (ValueError, TypeError) as e: logger.warning(f"Invalid coord type in tile info ({lat_coord},{lon_coord}): {e}. Skipping.", exc_info=False)
except Exception as e: logger.warning(f"Error getting bounds for tile ({lat_coord},{lon_coord}): {e}. Skipping.", exc_info=True)
if not initialized: logger.warning("No valid tiles in list for combined bounds."); return None
if min_lat >= max_lat: logger.warning(f"Calc invalid combined lat range: S={min_lat}, N={max_lat}."); return None
# Longitude wrapping in combined bounds from multiple tiles is possible, but harder to represent as simple (W,S,E,N) tuple.
# The current min/max calculation assumes no wrapping. If tiles cross antimeridian, this needs refinement.
# Given the expected use case (local map views), assuming non-wrapping longitude range is usually fine.
return (min_lon, min_lat, max_lon, max_lat)
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]
]:
"""Calc geo BBox centered at point corresponding to pixel size/zoom."""
if not PYPROJ_MODULE_LOCALLY_AVAILABLE or pyproj is None:
logger.error("'pyproj' required for BBox from pixel size/zoom.")
return None
if not MERCANTILE_MODULE_LOCALLY_AVAILABLE or mercantile is None:
logger.error("'mercantile' required for BBox from pixel size/zoom (for mercator conversion).")
return None
if not (isinstance(center_latitude_deg, (int, float)) and -90.0 <= center_latitude_deg <= 90.0): logger.error(f"Invalid center_latitude_deg: {center_latitude_deg}."); return None
if not (isinstance(center_longitude_deg, (int, float)) and -180.0 <= center_longitude_deg <= 180.0): logger.error(f"Invalid center_longitude_deg: {center_longitude_deg}."); return None
if not (isinstance(target_pixel_width, int) and target_pixel_width > 0 and isinstance(target_pixel_height, int) and target_pixel_height > 0): logger.error(f"Invalid pixel dims ({target_pixel_width}x{target_pixel_height})."); return None
if not (isinstance(zoom_level, int) and 0 <= zoom_level <= 25): logger.warning(f"Invalid zoom level: {zoom_level}. Clamping to."); zoom_level = max(0, min(zoom_level, 25)) # Clamp zoom level
if not (isinstance(tile_pixel_size, int) and tile_pixel_size > 0): logger.error(f"Invalid tile_pixel_size: {tile_pixel_size}."); return None
try:
# Convert center geo to mercator pixel at zoom 0
# Origin (0,0) in mercator pixel space at zoom 0 is the top-left corner of the earth map.
# X increases east, Y increases south. Total size is 256x256 pixels at zoom 0.
# Using the tile_pixel_size gives the size relative to that.
origin_mercator_x_at_z0, origin_mercator_y_at_z0 = mercantile.xy(-180.0, 85.05112878) # Top-left corner of web mercator projection
center_mercator_x, center_mercator_y = mercantile.xy(center_longitude_deg, center_latitude_deg)
# Calculate the pixel coordinates of the center point at the target zoom level
# Pixel coordinates at zoom Z are related to mercator coordinates and tile size.
# Total earth pixel size at zoom Z is tile_pixel_size * 2^Z
total_earth_pixels_z = tile_pixel_size * (2**zoom_level)
# Mercator X/Y range from (-PI * R, PI * R)
# Total width/height in mercator units is 2 * PI * R (~40075016.686 meters)
# Resolution at equator at zoom Z is (2 * PI * R) / (tile_pixel_size * 2^Z) meters/pixel
# A simpler way using mercantile x/y:
# X pixel = (mercator_x - mercator_min_x) / mercator_total_width * total_earth_pixels_z
# Y pixel = (mercator_max_y - mercator_y) / mercator_total_height * total_earth_pixels_z (Y is inverted)
# mercantile.xy returns coordinates in "Web Mercator" units (meters), relative to the origin of the projection.
# The origin of the projection is at (-180, +85.0511) geographic, corresponding to mercator (-20037508.34, 20037508.34).
# The full extent is 2 * 20037508.34 meters in both directions.
mercator_extent_meters = 20037508.34 * 2 # Total width/height of the Mercator projection in meters
center_pixel_x_z0 = (center_mercator_x - (-mercator_extent_meters/2)) / mercator_extent_meters * tile_pixel_size # Pixel relative to 256x256 tile at z0
center_pixel_y_z0 = ( (mercator_extent_meters/2) - center_mercator_y) / mercator_extent_meters * tile_pixel_size # Pixel relative to 256x256 tile at z0 (inverted Y)
# Pixel coordinates at zoom Z are pixel_at_z0 * 2^Z
center_pixel_x_z = center_pixel_x_z0 * (2**zoom_level)
center_pixel_y_z = center_pixel_y_z0 * (2**zoom_level)
# The target BBox pixel center is (target_pixel_width/2, target_pixel_height/2) relative to its top-left corner.
# The top-left pixel of the target BBox corresponds to a mercator coordinate.
# Let BBox_center_pixel_x, BBox_center_pixel_y be the center pixel of the *rendered* BBox (at the center of the canvas).
# The top-left pixel of the BBox is (BBox_center_pixel_x - target_pixel_width/2, BBox_center_pixel_y - target_pixel_height/2)
# The bottom-right pixel is (BBox_center_pixel_x + target_pixel_width/2, BBox_center_pixel_y + target_pixel_height/2)
# Assuming the stitched image is centered on the geo coordinate center_latitude_deg, center_longitude_deg,
# the center of the stitched image is at pixel (img_width/2, img_height/2).
# The pixel coordinates relative to the *top-left of the entire stitched map at zoom Z*
# that corresponds to the target BBox's top-left corner would be:
# ul_pixel_x_z = center_pixel_x_z - target_pixel_width / 2.0
# ul_pixel_y_z = center_pixel_y_z - target_pixel_height / 2.0
# lr_pixel_x_z = center_pixel_x_z + target_pixel_width / 2.0
# lr_pixel_y_z = center_pixel_y_z + target_pixel_height / 2.0
# Convert these pixel coordinates back to mercator, then to geo.
# mercator_x = mercator_min_x + (pixel_x_z / total_earth_pixels_z) * mercator_total_width
# mercator_y = mercator_max_y - (pixel_y_z / total_earth_pixels_z) * mercator_total_height # Inverted Y
ul_pixel_x_z = center_pixel_x_z - target_pixel_width / 2.0
ul_pixel_y_z = center_pixel_y_z - target_pixel_height / 2.0
lr_pixel_x_z = center_pixel_x_z + target_pixel_width / 2.0
lr_pixel_y_z = center_pixel_y_z + target_pixel_height / 2.0
ul_mercator_x = (-mercator_extent_meters/2) + (ul_pixel_x_z / total_earth_pixels_z) * mercator_extent_meters
ul_mercator_y = (mercator_extent_meters/2) - (ul_pixel_y_z / total_earth_pixels_z) * mercator_extent_meters # Inverted Y
lr_mercator_x = (-mercator_extent_meters/2) + (lr_pixel_x_z / total_earth_pixels_z) * mercator_extent_meters
lr_mercator_y = (mercator_extent_meters/2) - (lr_pixel_y_z / total_earth_pixels_z) * mercator_extent_meters # Inverted Y
# Convert mercator corners back to geo (lon, lat)
ul_lon, ul_lat = mercantile.lnglat(ul_mercator_x, ul_mercator_y)
lr_lon, lr_lat = mercantile.lnglat(lr_mercator_x, lr_mercator_y)
# The bounding box is (west_lon, south_lat, east_lon, north_lat)
west_lon = ul_lon
north_lat = ul_lat
east_lon = lr_lon
south_lat = lr_lat
# Ensure latitudes are within valid range for WGS84
north_lat = min(90.0, max(-90.0, north_lat))
south_lat = min(90.0, max(-90.0, south_lat))
# Handle longitude wrapping (e.g., if the canvas width spans the antimeridian)
# mercantile.lnglat handles the -180/180 wrap, but the calculated east_lon might be < west_lon.
# For consistency, ensure west_lon <= east_lon if not crossing antimeridian.
# If it crosses, the standard (W,S,E,N) tuple doesn't represent it well.
# Given our assumption for OpenSky, we expect non-wrapping BBoxes.
# If lr_lon < ul_lon, it might indicate wrapping, but for a single canvas view,
# we usually expect ul_lon < lr_lon unless the canvas is extremely wide and centered on antimeridian.
# For this calculation, we just return the calculated ul_lon and lr_lon as west and east.
if (abs(west_lon - east_lon) < 1e-9 or abs(south_lat - north_lat) < 1e-9):
logger.warning("Calc zero-size geo BBox from pixel size/zoom.")
return None
if not all(math.isfinite(coord) for coord in [west_lon, south_lat, east_lon, north_lat]):
logger.error("Calc geo BBox from pixel size/zoom contains non-finite values.")
return None
return (west_lon, south_lat, east_lon, north_lat)
except Exception as e_bbox_calc:
logger.exception(f"Error calculating geo BBox from pixel size/zoom: {e_bbox_calc}")
return None
def calculate_inner_geographic_bbox(
outer_bbox_deg: Tuple[float, float, float, float],
reduction_distance_km: float,
) -> Optional[
Tuple[float, float, float, float]
]:
"""Calc inner geo BBox by reducing outer BBox by distance (km). Requires pyproj."""
if not PYPROJ_MODULE_LOCALLY_AVAILABLE or pyproj is None:
logger.error("'pyproj' required for inner BBox.")
return None
west_lon, south_lat, east_lon, north_lat = outer_bbox_deg
reduct_dist_m = reduction_distance_km * 1000.0
if not (isinstance(reduction_distance_km, (int, float)) and reduction_distance_km >= 0): logger.error(f"Invalid reduct_dist_km: {reduction_distance_km}."); return None
if not (isinstance(west_lon, (int, float)) and isinstance(south_lat, (int, float)) and isinstance(east_lon, (int, float)) and isinstance(north_lat, (int, float))): logger.warning(f"Invalid coord types for outer_bbox: {outer_bbox_deg}."); return None
if not (-90.0 <= south_lat <= north_lat <= 90.0): logger.warning(f"Invalid lat range in outer_bbox: {south_lat}, {north_lat}."); return None
if not (-180.0 <= west_lon <= 180.0 and -180.0 <= east_lon <= 180.0): logger.warning(f"Invalid lon range in outer_bbox: {west_lon}, {east_lon}."); return None # Assuming non-wrapping input BBox
if reduction_distance_km == 0: return outer_bbox_deg
try:
geod = pyproj.Geod(ellps="WGS84")
outer_size_km = calculate_geographic_bbox_size_km(outer_bbox_deg)
if outer_size_km is None: logger.error("Could not calc outer BBox size for inner BBox."); return None
outer_width_km, outer_height_km = outer_size_km
inner_width_km = outer_width_km - (2 * reduction_distance_km)
inner_height_km = outer_height_km - (2 * reduction_distance_km)
if (inner_width_km <= 1e-9 or inner_height_km <= 1e-9):
logger.warning(f"Reduct dist {reduction_distance_km}km too large. Inner BBox size {inner_width_km:.2f}x{inner_height_km:.2f}km is zero or negative.")
return None
center_lat = (south_lat + north_lat) / 2.0
center_lon = (west_lon + east_lon) / 2.0
if west_lon > east_lon:
center_lon = (west_lon + east_lon + 360) / 2.0
if center_lon > 180: center_lon -= 360
center_lat_clamped = max(-89.9, min(89.9, center_lat))
center_lon_clamped = max(-179.9, min(179.9, center_lon)) # Clamp center_lon for fwd if needed
half_inner_width_m = inner_width_km * 500.0
half_inner_height_m = inner_height_km * 500.0
# Calculate inner BBox corners by moving from the center
_, inner_north_lat, _ = geod.fwd(center_lon_clamped, center_lat_clamped, 0.0, half_inner_height_m)
_, inner_south_lat, _ = geod.fwd(center_lon_clamped, center_lat_clamped, 180.0, half_inner_height_m)
inner_east_lon, _, _ = geod.fwd(center_lon_clamped, center_lat_clamped, 90.0, half_inner_width_m)
inner_west_lon, _, _ = geod.fwd(center_lon_clamped, center_lat_clamped, 270.0, half_inner_width_m)
inner_north_lat = min(90.0, max(-90.0, inner_north_lat)); inner_south_lat = min(90.0, max(-90.0, inner_south_lat))
if (abs(inner_west_lon - inner_east_lon) < 1e-9 or abs(inner_south_lat - inner_north_lat) < 1e-9):
logger.warning("Calc zero-size inner geo BBox.")
return None
if not all(math.isfinite(coord) for coord in [inner_west_lon, inner_south_lat, inner_east_lon, inner_north_lat]):
logger.error("Calc inner geo BBox contains non-finite values.")
return None
return (inner_west_lon, inner_south_lat, inner_east_lon, inner_north_lat)
except Exception as e_inner_bbox_calc:
logger.exception(f"Error calculating inner geo BBox: {e_inner_bbox_calc}")
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]], # (width, height)
) -> Tuple[Optional[float], Optional[float]]:
"""
Convert canvas pixel (x, y) to geo (lon, lat). Requires mercantile.
Uses current map geo bounds and pixel shape for context.
"""
if not MERCANTILE_MODULE_LOCALLY_AVAILABLE or mercantile is None:
logger.error("_pixel_to_geo: Mercantile missing.")
return None, None
if current_map_geo_bounds is None:
logger.warning("_pixel_to_geo: Map geo bounds missing.")
return None, None
if current_stitched_map_pixel_shape is None:
logger.warning("_pixel_to_geo: Map pixel shape missing.")
return None
img_w, img_h = current_stitched_map_pixel_shape
if img_w <= 0 or img_h <= 0:
logger.warning(f"_pixel_to_geo: Map pixel dims invalid ({img_w}x{img_h}).")
return None, None
map_west_lon, map_south_lat, map_east_lon, map_north_lat = current_map_geo_bounds
try:
map_ul_merc_x, map_ul_merc_y = mercantile.xy(map_west_lon, map_north_lat)
map_lr_merc_x, map_lr_merc_y = mercantile.xy(map_east_lon, map_south_lat)
total_map_width_merc = abs(map_lr_merc_x - map_ul_merc_x)
total_map_height_merc = abs(map_ul_merc_y - map_lr_merc_y)
if total_map_width_merc <= 1e-9 or total_map_height_merc <= 1e-9:
logger.warning(f"_pixel_to_geo: Map mercator dims zero/near-zero ({total_map_width_merc:.2e}x{total_map_height_merc:.2e}).")
return None, None
# Calculate the pixel coordinates relative to the top-left of the full mercator projection at the current zoom.
# Need the zoom level to know the total pixel size of the earth at that zoom.
# This function doesn't receive the zoom level directly, but it's used by MapCanvasManager which has it.
# However, the conversion from mercator to pixel relative to the *stitched image* only needs the stitched image geo bounds and pixel size.
# Convert pixel position (canvas_x, canvas_y) relative to stitched image to mercator coordinate.
# relative_x_on_map_pixels is 0 at the left edge of the stitched image, 1 at the right edge.
relative_x_on_map_pixels = canvas_x / img_w
# relative_y_on_map_pixels is 0 at the top edge, 1 at the bottom.
relative_y_on_map_pixels = canvas_y / img_h
# Convert relative pixel position to relative mercator position
# X: 0.0 at map_west_lon (ul_merc_x), 1.0 at map_east_lon (lr_merc_x)
# Y: 0.0 at map_north_lat (ul_merc_y), 1.0 at map_south_lat (lr_merc_y)
target_merc_x = map_ul_merc_x + (relative_x_on_map_pixels * total_map_width_merc)
target_merc_y = map_ul_merc_y - (relative_y_on_map_pixels * total_map_height_merc) # Y is inverted relative to map_ul_merc_y
# Convert mercator back to geo (lon, lat)
clicked_lon, clicked_lat = mercantile.lnglat(target_merc_x, target_merc_y)
# Clamp latitude to Mercator limits for safety
MAX_MERCATOR_LATITUDE = 85.05112878
clicked_lat = max(-MAX_MERCATOR_LATITUDE, min(MAX_MERCATOR_LATITUDE, clicked_lat))
return clicked_lon, clicked_lat
except Exception as e_pixel_to_geo:
logger.error(
f"Error converting pixel ({canvas_x}, {canvas_y}) to geo. Error: {e_pixel_to_geo}",
exc_info=True,
)
return None, None
Reference in New Issue View Git Blame Copy Permalink