172 lines
8.5 KiB
Python
172 lines
8.5 KiB
Python
"""Utility functions for map calculations (migrated from original project).
|
|
|
|
Contains bounding box, tile range and resolution helpers.
|
|
"""
|
|
import logging
|
|
import math
|
|
from typing import Tuple, Optional, List
|
|
|
|
try:
|
|
import pyproj
|
|
PYPROJ_AVAILABLE = True
|
|
except ImportError:
|
|
pyproj = None # type: ignore
|
|
PYPROJ_AVAILABLE = False
|
|
logging.warning("MapUtils: 'pyproj' library not found. Some geodetic calculations will fail.")
|
|
|
|
try:
|
|
import mercantile
|
|
MERCANTILE_AVAILABLE_UTILS = True
|
|
except ImportError:
|
|
mercantile = None # type: ignore
|
|
MERCANTILE_AVAILABLE_UTILS = False
|
|
logging.warning("MapUtils: 'mercantile' library not found. Tile calculations will fail.")
|
|
|
|
logger = logging.getLogger(__name__)
|
|
|
|
|
|
class MapCalculationError(Exception):
|
|
pass
|
|
|
|
|
|
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_AVAILABLE:
|
|
logger.error("'pyproj' library is required for bounding box calculation from center/size but is not found.")
|
|
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}. Must be a positive number.")
|
|
return None
|
|
if not (-90.0 <= center_latitude_deg <= 90.0):
|
|
logger.error(f"Invalid center_latitude_deg: {center_latitude_deg}. Must be in [-90, 90].")
|
|
return None
|
|
if not (-180.0 <= center_longitude_deg <= 180.0):
|
|
logger.error(f"Invalid center_longitude_deg: {center_longitude_deg}. Must be in [-180, 180].")
|
|
return None
|
|
|
|
try:
|
|
geodetic_calculator = pyproj.Geod(ellps="WGS84") # type: ignore
|
|
half_side_length_meters = (area_size_km / 2.0) * 1000.0
|
|
|
|
_, north_boundary_lat, _ = geodetic_calculator.fwd(lons=center_longitude_deg, lats=center_latitude_deg, az=0.0, dist=half_side_length_meters)
|
|
_, south_boundary_lat, _ = geodetic_calculator.fwd(lons=center_longitude_deg, lats=center_latitude_deg, az=180.0, dist=half_side_length_meters)
|
|
east_boundary_lon, _, _ = geodetic_calculator.fwd(lons=center_longitude_deg, lats=center_latitude_deg, az=90.0, dist=half_side_length_meters)
|
|
west_boundary_lon, _, _ = geodetic_calculator.fwd(lons=center_longitude_deg, lats=center_latitude_deg, az=270.0, dist=half_side_length_meters)
|
|
|
|
north_boundary_lat = min(90.0, max(-90.0, north_boundary_lat))
|
|
south_boundary_lat = min(90.0, max(-90.0, south_boundary_lat))
|
|
|
|
logger.debug(f"Calculated BBox: W={west_boundary_lon:.6f}, S={south_boundary_lat:.6f}, E={east_boundary_lon:.6f}, N={north_boundary_lat:.6f}")
|
|
return (west_boundary_lon, south_boundary_lat, east_boundary_lon, north_boundary_lat)
|
|
|
|
except Exception as e_bbox_calc:
|
|
logger.exception(f"Error calculating bounding box 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]]]:
|
|
if not MERCANTILE_AVAILABLE_UTILS:
|
|
logger.error("'mercantile' library is required for tile range calculation but is not found.")
|
|
return None
|
|
west_lon, south_lat, east_lon, north_lat = bounding_box_deg
|
|
try:
|
|
tiles_in_bbox_generator = mercantile.tiles(west_lon, south_lat, east_lon, north_lat, zooms=[zoom_level])
|
|
list_of_tiles = list(tiles_in_bbox_generator)
|
|
if not list_of_tiles:
|
|
logger.warning("No tiles found by mercantile.tiles for BBox at zoom {zoom_level}.")
|
|
clamped_west_lon = max(-180.0, min(180.0, west_lon))
|
|
clamped_east_lon = max(-180.0, min(180.0, east_lon))
|
|
clamped_south_lat = max(-90.0, min(90.0, south_lat))
|
|
clamped_north_lat = max(-90.0, min(90.0, north_lat))
|
|
center_lon = (clamped_west_lon + clamped_east_lon) / 2.0
|
|
center_lat = (clamped_south_lat + clamped_north_lat) / 2.0
|
|
center_lat = max(-85.0, min(85.0, center_lat))
|
|
if clamped_west_lon > clamped_east_lon:
|
|
center_lon = (clamped_west_lon + clamped_east_lon + 360) / 2.0
|
|
if center_lon > 180: center_lon -= 360
|
|
|
|
center_point_tile = mercantile.tile(center_lon, center_lat, zoom_level)
|
|
min_tile_x = center_point_tile.x
|
|
max_tile_x = center_point_tile.x
|
|
min_tile_y = center_point_tile.y
|
|
max_tile_y = center_point_tile.y
|
|
else:
|
|
x_coordinates = [tile.x for tile in list_of_tiles]
|
|
y_coordinates = [tile.y for tile in list_of_tiles]
|
|
min_tile_x = min(x_coordinates)
|
|
max_tile_x = max(x_coordinates)
|
|
min_tile_y = min(y_coordinates)
|
|
max_tile_y = max(y_coordinates)
|
|
|
|
return ((min_tile_x, max_tile_x), (min_tile_y, max_tile_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]:
|
|
try:
|
|
if not (-90.0 <= latitude_degrees <= 90.0):
|
|
logger.warning(f"Invalid latitude for m/px calc: {latitude_degrees}")
|
|
return None
|
|
if not (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
|
|
latitude_radians = math.radians(latitude_degrees)
|
|
resolution_m_px = (EARTH_CIRCUMFERENCE_METERS * math.cos(latitude_radians)) / (tile_pixel_size * (2**zoom_level))
|
|
if not math.isfinite(resolution_m_px) or resolution_m_px <= 0:
|
|
logger.warning(f"Calculated non-finite or non-positive m/px ({resolution_m_px}) at Lat {latitude_degrees}. Returning None.")
|
|
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]]:
|
|
if not PYPROJ_AVAILABLE:
|
|
logger.error("'pyproj' library is required for geographic size calculation but is not found.")
|
|
return None
|
|
west_lon, south_lat, east_lon, north_lat = bounding_box_deg
|
|
if not (-90.0 <= south_lat <= north_lat <= 90.0):
|
|
south_lat = max(-90.0, south_lat)
|
|
north_lat = min(90.0, north_lat)
|
|
if south_lat >= north_lat:
|
|
logger.error(f"Invalid latitude range after clamping: {south_lat}, {north_lat}. Cannot calculate size.")
|
|
return None
|
|
try:
|
|
geodetic_calculator = pyproj.Geod(ellps="WGS84")
|
|
center_lat = (south_lat + north_lat) / 2.0
|
|
center_lat = max(-89.9, min(89.9, center_lat))
|
|
_, _, width_meters = geodetic_calculator.inv(west_lon, center_lat, east_lon, center_lat)
|
|
center_lon_for_height = (west_lon + east_lon) / 2.0
|
|
center_lon_for_height = max(-180.0, min(180.0, center_lon_for_height))
|
|
_, _, height_meters = geodetic_calculator.inv(center_lon_for_height, south_lat, center_lon_for_height, north_lat)
|
|
approx_width_km = abs(width_meters) / 1000.0
|
|
approx_height_km = abs(height_meters) / 1000.0
|
|
if approx_width_km <= 0 or approx_height_km <= 0:
|
|
logger.warning(f"Calculated non-positive width or height for BBox {bounding_box_deg}. Result: ({approx_width_km:.2f}, {approx_height_km:.2f}). Returning None.")
|
|
return None
|
|
return (approx_width_km, approx_height_km)
|
|
except Exception as e_size_calc:
|
|
logger.exception(f"Error calculating geographic bounding box size: {e_size_calc}")
|
|
return None
|
|
|
|
|
|
def get_hgt_tile_geographic_bounds(lat_coord: int, lon_coord: int) -> Tuple[float, float, float, float]:
|
|
west_lon = float(lon_coord)
|
|
south_lat = float(lat_coord)
|
|
east_lon = float(lon_coord + 1)
|
|
north_lat = float(lat_coord + 1)
|
|
west_lon = max(-180.0, min(180.0, west_lon))
|
|
south_lat = max(-90.0, min(90.0, south_lat))
|
|
east_lon = max(-180.0, min(180.0, east_lon))
|
|
north_lat = max(-90.0, min(90.0, north_lat))
|
|
logger.debug(f"Calculated HGT tile bounds for ({lat_coord},{lon_coord}): ({west_lon:.6f}, {south_lat:.6f}, {east_lon:.6f}, {north_lat:.6f})")
|
|
return (west_lon, south_lat, east_lon, north_lat)
|