251 lines
10 KiB
Python
251 lines
10 KiB
Python
# geoelevation/map_viewer/map_utils.py
|
|
"""
|
|
Provides utility functions for map-related calculations.
|
|
|
|
Includes functions for determining geographic bounding boxes from center points
|
|
and sizes, finding necessary map tile ranges to cover an area using the
|
|
'mercantile' library, and calculating ground resolution (meters per pixel).
|
|
"""
|
|
|
|
# Standard library imports
|
|
import logging
|
|
import math
|
|
from typing import Tuple, Optional, List, Set # Ensure all used types are imported
|
|
|
|
# Third-party imports
|
|
try:
|
|
import pyproj # For geodetic calculations (e.g., bounding box from center)
|
|
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 # For tile calculations (Web Mercator)
|
|
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.")
|
|
|
|
# Module-level logger
|
|
logger = logging.getLogger(__name__)
|
|
|
|
|
|
class MapCalculationError(Exception):
|
|
"""Custom exception for errors occurring during map-related calculations."""
|
|
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]]: # (west_lon, south_lat, east_lon, north_lat)
|
|
"""
|
|
Calculates a geographic bounding box (WGS84 decimal degrees) given a center
|
|
point and a size (width/height of a square area).
|
|
|
|
Args:
|
|
center_latitude_deg: Latitude of the center point in degrees.
|
|
center_longitude_deg: Longitude of the center point in degrees.
|
|
area_size_km: The side length of the desired square area in kilometers.
|
|
|
|
Returns:
|
|
A tuple (west_lon, south_lat, east_lon, north_lat) in degrees,
|
|
or None if an error occurs or pyproj is unavailable.
|
|
|
|
Raises:
|
|
ImportError: Propagated if pyproj is required but not installed.
|
|
(Note: function now returns None instead of raising directly for this)
|
|
"""
|
|
if not PYPROJ_AVAILABLE:
|
|
logger.error(
|
|
"'pyproj' library is required for bounding box calculation from center/size but is not found."
|
|
)
|
|
return None # Return None instead of raising ImportError here to allow graceful degradation
|
|
|
|
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
|
|
|
|
|
|
logger.debug(
|
|
f"Calculating bounding box for center ({center_latitude_deg:.6f}, {center_longitude_deg:.6f}), "
|
|
f"size {area_size_km} km."
|
|
)
|
|
|
|
try:
|
|
# Initialize a geodetic calculator using the WGS84 ellipsoid
|
|
geodetic_calculator = pyproj.Geod(ellps="WGS84") # type: ignore
|
|
|
|
half_side_length_meters = (area_size_km / 2.0) * 1000.0
|
|
|
|
# Calculate points by projecting from the center along cardinal directions
|
|
# Azimuths: 0=North, 90=East, 180=South, 270=West
|
|
# geod.fwd returns (end_longitude, end_latitude, back_azimuth)
|
|
_, 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
|
|
)
|
|
|
|
# Handle potential latitude clamping at poles
|
|
# pyproj should handle this correctly, but a sanity check can be useful
|
|
north_boundary_lat = min(90.0, max(-90.0, north_boundary_lat))
|
|
south_boundary_lat = min(90.0, max(-90.0, south_boundary_lat))
|
|
|
|
# Handle longitude wrapping if the area is extremely large (unlikely for typical use)
|
|
# pyproj fwd/inv handles dateline crossing correctly for points.
|
|
# If west_lon > east_lon after projection, it implies dateline crossing.
|
|
# For bounding box, we typically want west < east unless it spans more than 180 deg.
|
|
# This simple assignment should be fine for areas not spanning the dateline/poles widely.
|
|
|
|
logger.debug(
|
|
f"Calculated BBox: W={west_boundary_lon:.6f}, S={south_boundary_lat:.6f}, "
|
|
f"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], # (west, south, east, north)
|
|
zoom_level: int
|
|
) -> Optional[Tuple[Tuple[int, int], Tuple[int, int]]]: # ((min_x, max_x), (min_y, max_y))
|
|
"""
|
|
Calculates the X and Y tile ranges (min_x, max_x, min_y, max_y)
|
|
required to cover a given geographic bounding box at a specific zoom level.
|
|
|
|
Args:
|
|
bounding_box_deg: A tuple (west_lon, south_lat, east_lon, north_lat) in degrees.
|
|
zoom_level: The Web Mercator map zoom level.
|
|
|
|
Returns:
|
|
A tuple containing ((min_x, max_x), (min_y, max_y)) tile coordinates,
|
|
or None if an error occurs or mercantile is unavailable.
|
|
"""
|
|
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
|
|
logger.debug(
|
|
f"Calculating tile ranges for zoom {zoom_level} and BBox W={west_lon:.4f}, "
|
|
f"S={south_lat:.4f}, E={east_lon:.4f}, N={north_lat:.4f}"
|
|
)
|
|
|
|
try:
|
|
# mercantile.tiles expects (west, south, east, north) and zoom as integer
|
|
# It returns a generator of Tile(x, y, z) objects.
|
|
tiles_in_bbox_generator = mercantile.tiles( # type: ignore
|
|
west_lon, south_lat, east_lon, north_lat, zooms=zoom_level
|
|
)
|
|
|
|
list_of_tiles = list(tiles_in_bbox_generator)
|
|
|
|
if not list_of_tiles:
|
|
# If the bbox is very small or outside standard tile limits, mercantile.tiles might be empty.
|
|
# As a fallback, find the tile containing the center of the bbox.
|
|
logger.warning(
|
|
f"No tiles found by mercantile.tiles for BBox at zoom {zoom_level}. "
|
|
"Using fallback: tile at BBox center."
|
|
)
|
|
center_lon = (west_lon + east_lon) / 2.0
|
|
center_lat = (south_lat + north_lat) / 2.0
|
|
# mercantile.tile(lon, lat, zoom)
|
|
center_point_tile = mercantile.tile(center_lon, center_lat, zoom_level) # type: ignore
|
|
|
|
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
|
|
num_tiles_found = 1
|
|
else:
|
|
# Extract all x and y coordinates from the list of tiles
|
|
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)
|
|
num_tiles_found = len(list_of_tiles)
|
|
|
|
logger.debug(
|
|
f"Calculated tile ranges for zoom {zoom_level}: "
|
|
f"X=[{min_tile_x}, {max_tile_x}], Y=[{min_tile_y}, {max_tile_y}] "
|
|
f"({num_tiles_found} tiles)"
|
|
)
|
|
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]:
|
|
"""
|
|
Calculates the approximate ground resolution (meters per pixel) at a given
|
|
latitude and zoom level for the Web Mercator projection.
|
|
|
|
Args:
|
|
latitude_degrees: Latitude in degrees.
|
|
zoom_level: Map zoom level.
|
|
tile_pixel_size: The size of one side of a map tile in pixels (usually 256).
|
|
|
|
Returns:
|
|
Approximate meters per pixel, or None if calculation fails or inputs are invalid.
|
|
"""
|
|
try:
|
|
if not (-90.0 <= latitude_degrees <= 90.0):
|
|
logger.warning(f"Invalid latitude for m/px calc: {latitude_degrees}")
|
|
return None
|
|
# Zoom levels for Web Mercator are typically 0-22, some services go higher.
|
|
if not (0 <= zoom_level <= 25): # Practical upper limit
|
|
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's equatorial circumference in meters (WGS84)
|
|
EARTH_CIRCUMFERENCE_METERS = 40075016.686
|
|
|
|
latitude_radians = math.radians(latitude_degrees)
|
|
|
|
# Formula: Resolution = (Circumference * cos(latitude_rad)) / (tile_size * 2^zoom)
|
|
resolution_m_px = (EARTH_CIRCUMFERENCE_METERS * math.cos(latitude_radians)) / \
|
|
(tile_pixel_size * (2**zoom_level))
|
|
|
|
logger.debug(
|
|
f"Calculated meters/pixel at lat {latitude_degrees:.4f}, zoom {zoom_level}, "
|
|
f"tile_size {tile_pixel_size}px: {resolution_m_px:.4f} m/px"
|
|
)
|
|
return resolution_m_px
|
|
|
|
except Exception as e_mpp_calc:
|
|
logger.exception(f"Error calculating meters per pixel: {e_mpp_calc}")
|
|
return None |