SXXXXXXX_FlightMonitor/flightmonitor/map/map_utils.py

# 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