233 lines
8.8 KiB
Python
233 lines
8.8 KiB
Python
# map_utils.py
|
|
"""
|
|
Provides utility functions for map-related calculations, such as determining
|
|
bounding boxes from geographic centers and sizes, and finding the necessary
|
|
map tiles to cover an area using the mercantile library.
|
|
"""
|
|
|
|
# Standard library imports
|
|
import logging
|
|
import math
|
|
from typing import Tuple, Optional, List, Set # <<< AGGIUNTO QUESTO IMPORT
|
|
|
|
# Third-party imports
|
|
try:
|
|
import pyproj # For geodetic calculations
|
|
except ImportError:
|
|
pyproj = None
|
|
try:
|
|
import mercantile # For tile calculations
|
|
except ImportError:
|
|
mercantile = None
|
|
|
|
# Local application imports
|
|
# (No direct imports from other app modules needed here usually)
|
|
|
|
|
|
class MapCalculationError(Exception):
|
|
"""Custom exception for errors during map calculations."""
|
|
|
|
pass
|
|
|
|
|
|
def get_bounding_box_from_center_size(
|
|
center_lat: float, center_lon: float, size_km: float
|
|
) -> Optional[
|
|
Tuple[float, float, float, float]
|
|
]: # <<< 'Optional' e 'Tuple' ora sono definiti
|
|
"""
|
|
Calculates the geographic bounding box (WGS84) given a center point and size.
|
|
|
|
Args:
|
|
center_lat (float): Latitude of the center point (degrees).
|
|
center_lon (float): Longitude of the center point (degrees).
|
|
size_km (float): The width and height of the desired square area in kilometers.
|
|
|
|
Returns:
|
|
Optional[Tuple[float, float, float, float]]: A tuple containing
|
|
(west_lon, south_lat, east_lon, north_lat) in degrees, or None on error.
|
|
|
|
Raises:
|
|
ImportError: If pyproj is not installed.
|
|
"""
|
|
log_prefix = "[MapUtils BBox]"
|
|
if pyproj is None:
|
|
logging.critical(
|
|
f"{log_prefix} 'pyproj' library is required for bounding box calculation but not found."
|
|
)
|
|
raise ImportError("'pyproj' library not found.")
|
|
if not (isinstance(size_km, (int, float)) and size_km > 0):
|
|
logging.error(
|
|
f"{log_prefix} Invalid size_km provided: {size_km}. Must be a positive number."
|
|
)
|
|
return None
|
|
|
|
logging.debug(
|
|
f"{log_prefix} Calculating bounding box for center ({center_lat:.6f}, {center_lon:.6f}), size {size_km} km."
|
|
)
|
|
|
|
try:
|
|
# Define the geodetic calculator (WGS84 ellipsoid)
|
|
geod = pyproj.Geod(ellps="WGS84")
|
|
|
|
# Calculate half-size in meters
|
|
half_size_meters = (size_km / 2.0) * 1000.0
|
|
|
|
# Calculate corner points by moving from the center along cardinals + diagonals
|
|
# Azimuth angles: 0=N, 90=E, 180=S, 270=W
|
|
# We need the points North, South, East, West from center to define the box easily.
|
|
lon_n, lat_n, _ = geod.fwd(center_lon, center_lat, 0, half_size_meters)
|
|
lon_s, lat_s, _ = geod.fwd(center_lon, center_lat, 180, half_size_meters)
|
|
lon_e, lat_e, _ = geod.fwd(center_lon, center_lat, 90, half_size_meters)
|
|
lon_w, lat_w, _ = geod.fwd(center_lon, center_lat, 270, half_size_meters)
|
|
|
|
# The bounding box is defined by the min/max latitudes and longitudes calculated
|
|
north_lat = lat_n
|
|
south_lat = lat_s
|
|
east_lon = lon_e
|
|
west_lon = lon_w
|
|
|
|
# Handle longitude wrapping around the 180 meridian (simple check)
|
|
if (east_lon - west_lon) > 180:
|
|
logging.warning(
|
|
f"{log_prefix} Potential longitude wrap detected (East={east_lon}, West={west_lon}). BBox might be inaccurate near dateline."
|
|
)
|
|
# This scenario is complex. A simple bbox might not work well.
|
|
# For moderate sizes not crossing the dateline, this shouldn't happen.
|
|
# If it does, we might need a more sophisticated approach or limit the area size.
|
|
# For now, we proceed but log the warning.
|
|
|
|
# Handle latitude crossing the pole (less likely for typical SAR sizes)
|
|
if abs(north_lat) > 90 or abs(south_lat) > 90:
|
|
logging.warning(
|
|
f"{log_prefix} Calculated latitude exceeds +/- 90 degrees. Clamping."
|
|
)
|
|
north_lat = min(north_lat, 90.0)
|
|
south_lat = max(south_lat, -90.0)
|
|
|
|
logging.debug(
|
|
f"{log_prefix} Calculated BBox: W={west_lon:.6f}, S={south_lat:.6f}, E={east_lon:.6f}, N={north_lat:.6f}"
|
|
)
|
|
return (west_lon, south_lat, east_lon, north_lat)
|
|
|
|
except Exception as e:
|
|
logging.exception(f"{log_prefix} Error calculating bounding box:")
|
|
return None
|
|
|
|
|
|
def get_tile_ranges_for_bbox(
|
|
bbox: Tuple[float, float, float, float], zoom: int
|
|
) -> Optional[
|
|
Tuple[Tuple[int, int], Tuple[int, int]]
|
|
]: # <<< 'Optional' e 'Tuple' ora sono definiti
|
|
"""
|
|
Calculates the required X and Y tile ranges for a given bounding box and zoom level.
|
|
|
|
Args:
|
|
bbox (Tuple[float, float, float, float]): Bounding box tuple
|
|
(west_lon, south_lat, east_lon, north_lat) in degrees.
|
|
zoom (int): The desired map zoom level.
|
|
|
|
Returns:
|
|
Optional[Tuple[Tuple[int, int], Tuple[int, int]]]: A tuple containing
|
|
((min_x, max_x), (min_y, max_y)) tile coordinates, or None on error.
|
|
|
|
Raises:
|
|
ImportError: If mercantile is not installed.
|
|
"""
|
|
log_prefix = "[MapUtils Tiles]"
|
|
if mercantile is None:
|
|
logging.critical(
|
|
f"{log_prefix} 'mercantile' library is required for tile calculation but not found."
|
|
)
|
|
raise ImportError("'mercantile' library not found.")
|
|
|
|
west, south, east, north = bbox
|
|
logging.debug(
|
|
f"{log_prefix} Calculating tile ranges for zoom {zoom} and BBox W={west:.4f}, S={south:.4f}, E={east:.4f}, N={north:.4f}"
|
|
)
|
|
|
|
try:
|
|
# Use mercantile.tiles to get all tiles intersecting the bounding box
|
|
# Note: mercantile expects (west, south, east, north)
|
|
tiles_generator = mercantile.tiles(
|
|
west, south, east, north, zooms=zoom
|
|
) # Pass zoom as integer
|
|
|
|
# Collect tiles as a list to find min/max easily
|
|
tiles = list(tiles_generator)
|
|
|
|
if not tiles:
|
|
logging.warning(
|
|
f"{log_prefix} No tiles found for the given bounding box and zoom level {zoom}."
|
|
)
|
|
# Attempt to find the single tile containing the center point as a fallback
|
|
center_lon = (west + east) / 2.0
|
|
center_lat = (south + north) / 2.0
|
|
center_tile = mercantile.tile(center_lon, center_lat, zoom)
|
|
logging.warning(
|
|
f"{log_prefix} Using fallback: single tile at center ({center_lon:.4f},{center_lat:.4f}) -> {center_tile}"
|
|
)
|
|
min_x = max_x = center_tile.x
|
|
min_y = max_y = center_tile.y
|
|
tile_count = 1
|
|
# return None # Original behavior: return None if mercantile.tiles is empty
|
|
|
|
else:
|
|
# Extract x and y coordinates and find min/max
|
|
x_coords = [tile.x for tile in tiles]
|
|
y_coords = [tile.y for tile in tiles]
|
|
|
|
min_x = min(x_coords)
|
|
max_x = max(x_coords)
|
|
min_y = min(y_coords)
|
|
max_y = max(y_coords)
|
|
tile_count = len(tiles)
|
|
|
|
logging.debug(
|
|
f"{log_prefix} Calculated tile ranges for zoom {zoom}: X=[{min_x}, {max_x}], Y=[{min_y}, {max_y}] ({tile_count} tiles)"
|
|
)
|
|
return ((min_x, max_x), (min_y, max_y))
|
|
|
|
except Exception as e:
|
|
logging.exception(f"{log_prefix} Error calculating tile ranges:")
|
|
return None
|
|
|
|
|
|
def calculate_meters_per_pixel(latitude_deg: float, zoom: int) -> Optional[float]:
|
|
"""
|
|
Calculates the approximate ground resolution (meters per pixel) at a given
|
|
latitude and zoom level for the Web Mercator projection (used by OSM/Google).
|
|
|
|
Args:
|
|
latitude_deg (float): Latitude in degrees.
|
|
zoom (int): Map zoom level.
|
|
|
|
Returns:
|
|
Optional[float]: Meters per pixel, or None if calculation fails.
|
|
"""
|
|
log_prefix = "[MapUtils Scale]"
|
|
try:
|
|
# Validate inputs
|
|
if not (-90 <= latitude_deg <= 90):
|
|
logging.warning(f"{log_prefix} Invalid latitude: {latitude_deg}")
|
|
return None
|
|
if not (0 <= zoom <= 22): # Practical zoom range limit
|
|
logging.warning(f"{log_prefix} Invalid zoom level: {zoom}")
|
|
return None
|
|
|
|
# Formula based on Earth circumference and tile size (usually 256px)
|
|
# Meters per pixel = (Earth Circumference * cos(latitude)) / (Tile Size * 2^zoom)
|
|
# Earth Circumference approx 40075016.686 meters at equator
|
|
C = 40075016.686 # meters
|
|
TILE_SIZE = 256 # pixels
|
|
latitude_rad = math.radians(latitude_deg)
|
|
meters_per_pixel = (C * math.cos(latitude_rad)) / (TILE_SIZE * (2**zoom))
|
|
logging.debug(
|
|
f"{log_prefix} Calculated meters/pixel at lat {latitude_deg:.4f}, zoom {zoom}: {meters_per_pixel:.4f}"
|
|
)
|
|
return meters_per_pixel
|
|
except Exception as e:
|
|
logging.exception(f"{log_prefix} Error calculating meters per pixel:")
|
|
return None
|