584 lines
28 KiB
Python
584 lines
28 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] # 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."
|
|
)
|
|
# Ensure coordinates are within valid WGS84 bounds before calculating center
|
|
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
|
|
|
|
# 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 - approx. latitude of tile row 0 or max
|
|
# Ensure center_lon wraps correctly if bbox spans the antimeridian
|
|
if clamped_west_lon > clamped_east_lon: # Bbox crosses antimeridian
|
|
center_lon = (clamped_west_lon + clamped_east_lon + 360) / 2.0
|
|
if center_lon > 180: center_lon -= 360
|
|
|
|
|
|
# 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}")
|
|
# Try clamping and continue, or return None? Let's clamp for robustness.
|
|
south_lat = max(-90.0, south_lat)
|
|
north_lat = min(90.0, north_lat)
|
|
if south_lat >= north_lat: # After clamping, check if still invalid
|
|
logger.error(f"Invalid latitude range after clamping: {south_lat}, {north_lat}. Cannot calculate size.")
|
|
return None
|
|
|
|
|
|
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.
|
|
# Using the average longitude for height calculation line might not be strictly necessary,
|
|
# distance between two points at the same longitude is simple geodetic distance.
|
|
# However, using inv is more general and handles the ellipsoid correctly.
|
|
center_lon_for_height = (west_lon + east_lon) / 2.0 # Use the average longitude for height calculation line
|
|
# Handle antimeridian crossing for height calculation by adjusting center_lon_for_height if needed.
|
|
# If the width calculation crossed the antimeridian (west_lon > east_lon originally), the average might be misleading.
|
|
# A simpler approach for height is just calculating distance between (center_lon, south) and (center_lon, north).
|
|
# Let's reuse the original logic using the average longitude, ensuring it's within range.
|
|
center_lon_for_height = max(-180.0, min(180.0, center_lon_for_height)) # Ensure within standard range
|
|
|
|
|
|
_, _, 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 # 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))
|
|
|
|
# Ensure west < east and south < north for valid bbox representation, handling wrap-around is complex
|
|
# For 1x1 degree tiles not crossing antimeridian or poles widely, this is fine.
|
|
# A more robust check might be needed for edge cases near antimeridian, but for 1x1 tiles this is usually okay.
|
|
if west_lon > east_lon:
|
|
# This case ideally shouldn't happen for 1x1 degree tiles unless lon_coord is 179 and it wraps,
|
|
# but let's handle defensively if coords are unusual. Could swap or adjust.
|
|
# For HGT tile boundaries, it's usually [lon, lon+1] so this shouldn't be an issue.
|
|
logger.warning(f"Calculated west > east for HGT tile bounds ({lat_coord},{lon_coord}): ({west_lon}, {east_lon}).")
|
|
# Assuming it's a simple swap needed
|
|
# west_lon, east_lon = east_lon, west_lon # This might not be correct if it actually wraps the globe
|
|
pass # Let's stick to the direct calculation based on HGT convention
|
|
|
|
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
|
|
|
|
# MODIFIED: Added utility function to convert decimal degrees to DMS string format.
|
|
# WHY: Needed for displaying coordinates in a user-friendly, copyable format in the GUI.
|
|
# HOW: Implemented standard conversion logic.
|
|
def deg_to_dms_string(degree_value: float, coord_type: str) -> str:
|
|
"""
|
|
Converts a decimal degree coordinate to a Degrees, Minutes, Seconds (DMS) string.
|
|
|
|
Args:
|
|
degree_value: The coordinate value in decimal degrees (float).
|
|
coord_type: 'lat' for latitude (determines N/S suffix), 'lon' for longitude (E/W suffix).
|
|
|
|
Returns:
|
|
A formatted DMS string (e.g., "45° 30' 15.23'' N"). Returns "N/A" for non-finite or invalid inputs.
|
|
"""
|
|
if not isinstance(degree_value, (int, float)) or not math.isfinite(degree_value):
|
|
# Handle None, NaN, Inf, etc.
|
|
return "N/A"
|
|
|
|
# Clamp to valid ranges for sanity, although conversion works outside.
|
|
# Note: DMS representation is strictly defined for lat [-90, 90] and lon [-180, 180]
|
|
if coord_type.lower() == 'lat':
|
|
if not -90.0 <= degree_value <= 90.0:
|
|
logger.warning(f"Latitude {degree_value} is outside valid range [-90, 90] for standard DMS.")
|
|
# Still convert, but maybe add a note or handle separately if needed.
|
|
# For now, just convert the clamped value.
|
|
degree_value = max(-90.0, min(90.0, degree_value))
|
|
elif coord_type.lower() == 'lon':
|
|
# Longitude wrapping - DMS usually represents within [-180, 180].
|
|
# Python's % operator handles negative numbers differently than some other languages,
|
|
# (a % b) has the same sign as b. So ((value + 180) % 360 - 180) ensures the result
|
|
# is in (-180, 180]. If the input is exactly 180, it becomes 180.
|
|
degree_value = ((degree_value + 180) % 360) - 180
|
|
# Check if exactly -180.0 and adjust to 180.0 if needed for conventional representation
|
|
# if degree_value == -180.0:
|
|
# degree_value = 180.0
|
|
|
|
else:
|
|
logger.warning(f"Unknown coordinate type '{coord_type}' for DMS conversion.")
|
|
return "N/A (Invalid Type)"
|
|
|
|
|
|
is_negative = degree_value < 0
|
|
abs_degree = abs(degree_value)
|
|
|
|
degrees = int(abs_degree)
|
|
minutes_decimal = (abs_degree - degrees) * 60
|
|
minutes = int(minutes_decimal)
|
|
seconds = (minutes_decimal - minutes) * 60
|
|
|
|
# Determine suffix (North/South, East/West)
|
|
suffix = ""
|
|
if coord_type.lower() == 'lat':
|
|
suffix = "N" if not is_negative else "S"
|
|
elif coord_type.lower() == 'lon':
|
|
# Longitude can be 0 or 180, technically not negative/positive in terms of E/W.
|
|
# Let's use the sign logic which is correct for the suffix E/W convention.
|
|
suffix = "E" if not is_negative else "W"
|
|
# Special case for exactly 0 or 180, maybe no suffix or specific suffix?
|
|
# Standard practice is usually to use E for 0 and 180 if positive/negative sign is ignored.
|
|
if degree_value == 0.0: suffix = "" # No suffix for 0? Or E? Let's use ""
|
|
elif abs(degree_value) == 180.0: suffix = "" # No suffix for 180? Or W? Let's use ""
|
|
|
|
# Format the string
|
|
# Use a consistent number of decimal places for seconds, e.g., 2
|
|
# Handle potential edge case where seconds round up to 60
|
|
if seconds >= 59.995: # Check if seconds are very close to 60 due to float precision
|
|
seconds = 0.0
|
|
minutes += 1
|
|
if minutes >= 60:
|
|
minutes = 0
|
|
degrees += 1
|
|
# Degree could potentially roll over if it was like 89 deg 59 min 59.99 sec
|
|
# This simple logic assumes it won't cross 90 or 180 significantly with typical inputs
|
|
# A more robust implementation might need to handle full degree rollover, but unlikely for standard inputs.
|
|
# For now, just increment degree if minutes rolled over from 59 to 60.
|
|
|
|
|
|
dms_string = f"{degrees}° {minutes}' {seconds:.2f}''"
|
|
|
|
# Add suffix only if it's meaningful (not empty)
|
|
if suffix:
|
|
dms_string += f" {suffix}"
|
|
|
|
|
|
return dms_string |