VALLONGOL/SXXXXXXX_GeoElevation
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
aa7fb18626
SXXXXXXX_GeoElevation/geoelevation/map_viewer/map_utils.py
VALLONGOL 6925e8e323 ok view point on map
2025-05-13 16:01:44 +02:00

465 lines
22 KiB
Python
Raw Blame History

# 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] # Pass zoom as a list
)
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
# Clamp center_lat to avoid mercantile issues near poles if bbox extends beyond valid range
center_lat = max(-85.0, min(85.0, center_lat)) # Mercantile limits
center_lon = max(-180.0, min(180.0, center_lon)) # Mercantile limits
# 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))
# Avoid returning non-finite values if cos(latitude) is near zero at poles
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
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
# MODIFIED: Added function to calculate geographic size of a bounding box.
# WHY: Needed to report the displayed map area size in the GUI.
# HOW: Implemented logic using pyproj to calculate distances for width and height.
def calculate_geographic_bbox_size_km(
bounding_box_deg: Tuple[float, float, float, float] # (west, south, east, north)
) -> Optional[Tuple[float, float]]: # Returns (approx_width_km, approx_height_km)
"""
Calculates the approximate geographic width and height of a bounding box in kilometers.
Uses pyproj if available. Width is calculated along the center latitude, height along center longitude.
Args:
bounding_box_deg: A tuple (west_lon, south_lat, east_lon, north_lat) in degrees.
Returns:
A tuple (approx_width_km, approx_height_km), or None if calculation fails or pyproj is unavailable.
"""
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
# Basic validation
if not (-90.0 <= south_lat <= north_lat <= 90.0):
logger.warning(f"Invalid latitude range for size calculation: {south_lat}, {north_lat}")
return None
# For longitude, check range is within a reasonable bound, but allow west > east for dateline crossing
# The geographic width should be calculated carefully to handle wrapping around the globe.
# The height calculation is simpler as latitude is bounded.
try:
geodetic_calculator = pyproj.Geod(ellps="WGS84") # type: ignore
# Calculate approximate width along the center latitude
# The 'inv' method from pyproj.Geod is suitable for this, it handles the antimeridian.
center_lat = (south_lat + north_lat) / 2.0
# Clamp center lat to avoid issues near poles for geod.inv
center_lat = max(-89.9, min(89.9, center_lat))
# geod.inv returns (forward_azimuth, backward_azimuth, distance)
_, _, width_meters = geodetic_calculator.inv(west_lon, center_lat, east_lon, center_lat)
# Calculate approximate height along the center longitude
# This is simpler, distance between south_lat and north_lat at center_lon
# Need to handle potential longitude wrap around - use inv carefully
# The straight line distance calculation below should be generally fine for height.
center_lon = (west_lon + east_lon) / 2.0 # Use the average longitude for height calculation line
# Clamp center lon
center_lon = max(-179.9, min(179.9, center_lon))
_, _, height_meters = geodetic_calculator.inv(center_lon, south_lat, center_lon, north_lat)
approx_width_km = abs(width_meters) / 1000.0 # Ensure positive distance
approx_height_km = abs(height_meters) / 1000.0
# Add a sanity check: if width or height are zero, something is wrong (e.g., points are identical)
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
logger.debug(
f"Calculated BBox size for {bounding_box_deg}: "
f"Approx. {approx_width_km:.2f}km W x {approx_height_km:.2f}km H"
)
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
# MODIFIED: Added function to calculate the geographic bounds of an HGT tile from its integer coordinates.
# WHY: Needed to get the exact bounds of the DEM tile to determine the map fetch area and potentially draw the boundary.
# HOW: Based on HGT tile naming conventions (e.g., N45E007 covers 7E-8E, 45N-46N), the integer coordinates are the southwest corner.
def get_hgt_tile_geographic_bounds(lat_coord: int, lon_coord: int) -> Tuple[float, float, float, float]:
"""
Calculates the precise geographic bounding box (W, S, E, N) for an HGT tile
based on its integer latitude and longitude coordinates.
Assumes standard HGT 1x1 degree tile coverage where lat_coord is the
southern boundary latitude and lon_coord is the western boundary longitude.
E.g., tile N45E007 (lat_coord=45, lon_coord=7) covers 7E-8E, 45N-46N.
Args:
lat_coord: The integer latitude coordinate of the tile (e.g., 45 for N45).
lon_coord: The integer longitude coordinate of the tile (e.g., 7 for E007).
Returns:
A tuple (west_lon, south_lat, east_lon, north_lat) in degrees.
"""
west_lon = float(lon_coord)
south_lat = float(lat_coord)
east_lon = float(lon_coord + 1)
north_lat = float(lat_coord + 1)
# Clamping to strict WGS84 bounds for sanity, though tile coords should respect this
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)
# MODIFIED: Added function to calculate required zoom level for a target geographic size to fit in a target pixel size.
# WHY: To determine the appropriate zoom level for displaying the 1x1 degree DEM tile area within a manageable pixel size.
# HOW: Used the inverse of the calculate_meters_per_pixel formula to solve for the zoom level.
def calculate_zoom_level_for_geographic_size(
latitude_degrees: float,
geographic_height_meters: float, # The height of the area you want to fit in pixels
target_pixel_height: int, # The desired pixel height for that geographic area
tile_pixel_size: int = 256
) -> Optional[int]:
"""
Calculates the approximate Web Mercator zoom level required for a given
geographic height (in meters at a specific latitude) to span a target pixel height.
This is useful for determining the zoom needed to fit a known geographic area
(like a DEM tile's height) into a certain number of pixels on a map composed of tiles.
Args:
latitude_degrees: The latitude at which the geographic_height_meters is measured.
geographic_height_meters: The actual height of the geographic area in meters at that latitude.
target_pixel_height: The desired height in pixels for that geographic area on the map.
tile_pixel_size: The size of one side of a map tile in pixels (usually 256).
Returns:
The approximate integer zoom level, or None if calculation fails or inputs are invalid.
"""
if not (-90.0 <= latitude_degrees <= 90.0):
logger.warning(f"Invalid latitude for zoom calculation: {latitude_degrees}")
return None
if not (isinstance(geographic_height_meters, (int, float)) and geographic_height_meters > 0):
logger.warning(f"Invalid geographic_height_meters for zoom calculation: {geographic_height_meters}. Must be positive.")
return None
if not (isinstance(target_pixel_height, int) and target_pixel_height > 0):
logger.warning(f"Invalid target_pixel_height for zoom calculation: {target_pixel_height}. Must be positive integer.")
return None
if not (isinstance(tile_pixel_size, int) and tile_pixel_size > 0):
logger.warning(f"Invalid tile_pixel_size for zoom calculation: {tile_pixel_size}. Must be positive integer.")
return None
try:
# Earth's equatorial circumference in meters (WGS84)
EARTH_CIRCUMFERENCE_METERS = 40075016.686
# Calculate the required meters per pixel (resolution) to fit the geographic height into the target pixel height
required_resolution_m_px = geographic_height_meters / target_pixel_height
# Avoid division by zero or non-finite results
if required_resolution_m_px <= 0 or not math.isfinite(required_resolution_m_px):
logger.warning(f"Calculated non-positive or non-finite required resolution ({required_resolution_m_px} m/px). Cannot calculate zoom.")
return None
# Use the inverse of the resolution formula: Resolution = (Circumference * cos(latitude_rad)) / (tile_size * 2^z)
# Rearranging for 2^z: 2^z = (Circumference * cos(latitude_rad)) / (tile_size * Resolution)
# Solving for z: z = log2( (Circumference * cos(latitude_rad)) / (tile_size * Resolution) )
latitude_radians = math.radians(latitude_degrees)
cos_lat = math.cos(latitude_radians)
# Avoid division by zero or log of zero/negative if cos_lat is near zero (at poles)
# The latitude validation should prevent hitting exactly 90/-90, but check for very small values.
if abs(cos_lat) < 1e-9: # Handle near poles
logger.warning(f"Latitude {latitude_degrees} is too close to a pole for reliable zoom calculation.")
# Return a very low zoom level as a fallback? Or None?
# Given the context (mapping a DEM tile), this likely won't happen as DEMs stop at 60 deg.
return None # Return None for latitudes very close to poles
term_for_log = (EARTH_CIRCUMFERENCE_METERS * cos_lat) / (tile_pixel_size * required_resolution_m_px)
# Ensure the argument for log2 is positive
if term_for_log <= 0:
logger.warning(f"Calculated non-positive term for log2 ({term_for_log}) during zoom calculation. Cannot calculate zoom.")
return None
# Calculate the precise zoom level (can be fractional)
precise_zoom = math.log2(term_for_log)
# We need an integer zoom level for tile fetching. Rounding to the nearest integer is common.
# Floor might get fewer tiles than needed, ceil might get more. Rounding is a good balance.
integer_zoom = int(round(precise_zoom))
# Clamp the calculated zoom level to a reasonable range (e.g., 0 to 20)
# A zoom level too high (e.g., > 22) might not be supported by map services.
# A level too low (negative) indicates an issue or a request for an impossibly large area.
# We should probably cap it to the map service's max zoom as well, but that info isn't available here.
# Let's clamp it to a general reasonable range.
clamped_zoom = max(0, min(integer_zoom, 20)) # Max zoom of 20 is usually safe for OSM
logger.debug(
f"Calculated zoom for {geographic_height_meters:.2f}m at Lat {latitude_degrees:.4f} "
f"to fit in {target_pixel_height}px: Precise Zoom {precise_zoom:.2f}, Clamped Integer Zoom {clamped_zoom}"
)
return clamped_zoom
except Exception as e_zoom_calc:
logger.exception(f"Error calculating zoom level: {e_zoom_calc}")
return None
Reference in New Issue View Git Blame Copy Permalink