VALLONGOL/SXXXXXXX_FlightMonitor
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
0092abd909
SXXXXXXX_FlightMonitor/flightmonitor/map/map_utils.py
VALLONGOL 0092abd909 add map visualization
2025-05-16 12:23:48 +02:00

1249 lines
62 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, Dict # 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
# MODIFIED: Added new utility function to calculate the combined geographic bounds from a list of tile infos.
# WHY: Needed for calculating the overall geographic extent of a group of DEM tiles for map display aspect ratio.
# HOW: Iterate through the list, get bounds for each tile, and find the min/max lat/lon.
def get_combined_geographic_bounds_from_tile_info_list(tile_info_list: List[Dict]) -> Optional[Tuple[float, float, float, float]]:
"""
Calculates the minimum bounding box that encompasses all tiles in the provided list.
Args:
tile_info_list: A list of tile information dictionaries (e.g., from ElevationManager.get_area_tile_info).
Each dict must contain 'latitude_coord' and 'longitude_coord'.
Returns:
A tuple (west_lon, south_lat, east_lon, north_lat) encompassing all tiles,
or None if the list is empty or invalid info is found.
"""
if not tile_info_list:
logger.warning("Tile info list is empty, cannot calculate combined geographic bounds.")
return None
min_lon_combined, min_lat_combined, max_lon_combined, max_lat_combined = 180.0, 90.0, -180.0, -90.0
initialized = False
for tile_info in tile_info_list:
lat_coord = tile_info.get("latitude_coord")
# MODIFIED: Corrected typo in key name.
# WHY: To correctly access the longitude coordinate from the dictionary.
# HOW: Changed "longitude_longitude" to "longitude_coord".
lon_coord = tile_info.get("longitude_coord")
if lat_coord is None or lon_coord is None:
logger.warning(f"Skipping tile info entry due to missing coordinates: {tile_info}")
continue # Skip this entry if coordinates are missing
try:
# Get the precise geographic bounds for this HGT tile
tile_bounds = get_hgt_tile_geographic_bounds(lat_coord, lon_coord)
if tile_bounds: # Ensure bounds were calculated successfully
if not initialized:
min_lon_combined, min_lat_combined, max_lon_combined, max_lat_combined = tile_bounds
initialized = True
else:
min_lon_combined = min(min_lon_combined, tile_bounds[0])
min_lat_combined = min(min_lat_combined, tile_bounds[1])
max_lon_combined = max(max_lon_combined, tile_bounds[2])
max_lat_combined = max(max_lat_combined, tile_bounds[3])
else:
logger.warning(f"Could not get geographic bounds for tile ({lat_coord},{lon_coord}), skipping.")
except Exception as e_get_tile_bounds:
logger.warning(f"Error getting geographic bounds for tile ({lat_coord},{lon_coord}): {e_get_tile_bounds}. Skipping this tile.")
continue # Skip tile with invalid bounds
if not initialized:
logger.warning("No valid tile coordinates found in the list to calculate combined bounds.")
return None # No valid tiles processed
# Final validation of combined bounds (e.g., if it spans the whole globe)
# The bounds from HGT tiles are 1x1 degree, so combining shouldn't create invalid wrap-around issues easily,
# but defensive check.
# A more robust check might be needed for edge cases near antimeridian, but for 1x1 tiles this is usually fine.
# For simplicity, assume min_lat < max_lat and min_lon < max_lon unless it crosses the antimeridian.
# The get_hgt_tile_geographic_bounds clamps, so min/max comparison should work generally.
# If the combined area crosses the antimeridian, min_lon_combined will be > max_lon_combined.
# We might need to represent this differently if we need a bbox spanning the antimeridian,
# but for simple min/max extent, the current logic is okay for non-global areas.
# Let's just check for egregious errors like south > north.
if min_lat_combined > max_lat_combined:
logger.warning(f"Calculated invalid combined latitude range: S={min_lat_combined}, N={max_lat_combined}. Returning None.")
return None
# Longitude range check is tricky with antimeridian. Assuming for typical areas the min/max works.
logger.debug(f"Calculated combined geographic bounds: ({min_lon_combined:.6f}, {min_lat_combined:.6f}, {max_lon_combined:.6f}, {max_lat_combined:.6f})")
return (min_lon_combined, min_lat_combined, max_lon_combined, max_lat_combined)
# MODIFIED: Added new utility function to calculate a geographic bounding box
# given a center point, desired pixel dimensions, and zoom level.
# WHY: Necessary for interactive zoom to determine the map area to fetch and display
# while maintaining a relatively constant pixel size for the map window.
# HOW: Uses the inverse of the meters per pixel calculation to find the geographic
# width and height corresponding to the desired pixel dimensions at the given
# latitude and zoom. Then uses pyproj to find the geographic coordinates
# of the corners of a bounding box centered on the input point with these
# calculated geographic width/height.
def calculate_geographic_bbox_from_pixel_size_and_zoom(
center_latitude_deg: float,
center_longitude_deg: float,
target_pixel_width: int,
target_pixel_height: int,
zoom_level: int,
tile_pixel_size: int = 256
) -> Optional[Tuple[float, float, float, float]]: # (west_lon, south_lat, east_lon, north_lat)
"""
Calculates the geographic bounding box centered on a point that corresponds
to a specific pixel size at a given zoom level and latitude.
Args:
center_latitude_deg: Latitude of the center point in degrees.
center_longitude_deg: Longitude of the center point in degrees.
target_pixel_width: Desired pixel width for the bounding box.
target_pixel_height: Desired pixel height for the bounding box.
zoom_level: Map zoom level.
tile_pixel_size: The size of one side of a map tile in pixels.
Returns:
A tuple (west_lon, south_lat, east_lon, north_lat) in degrees,
or None if calculation fails or dependencies are unavailable.
"""
if not PYPROJ_AVAILABLE:
logger.error(
"'pyproj' library is required for bbox calculation but is not found."
)
return None
if not MERCANTILE_AVAILABLE_UTILS:
logger.error(
"'mercantile' library is required for bbox calculation but is not found."
)
return None
if not (-90.0 <= center_latitude_deg <= 90.0):
logger.error(f"Invalid center_latitude_deg for bbox calc: {center_latitude_deg}")
return None
if not (-180.0 <= center_longitude_deg <= 180.0):
logger.error(f"Invalid center_longitude_deg for bbox calc: {center_longitude_deg}")
return None
if not (isinstance(target_pixel_width, int) and target_pixel_width > 0 and
isinstance(target_pixel_height, int) and target_pixel_height > 0):
logger.error(f"Invalid target pixel dimensions ({target_pixel_width}x{target_pixel_height}) for bbox calc.")
return None
if not (isinstance(zoom_level, int) and 0 <= zoom_level <= 25): # Clamp to practical max zoom
logger.warning(f"Invalid zoom level for bbox calc: {zoom_level}. Clamping to [0, 20].")
zoom_level = max(0, min(zoom_level, 20)) # Clamp
if not (isinstance(tile_pixel_size, int) and tile_pixel_size > 0):
logger.error(f"Invalid tile_pixel_size for bbox calc: {tile_pixel_size}")
return None
logger.debug(
f"Calculating geographic bbox for center ({center_latitude_deg:.6f}, {center_longitude_deg:.6f}) "
f"at zoom {zoom_level} for {target_pixel_width}x{target_pixel_height}px."
)
try:
# Calculate the ground resolution (meters per pixel) at the center latitude for the given zoom
resolution_m_px = calculate_meters_per_pixel(center_latitude_deg, zoom_level, tile_pixel_size)
if resolution_m_px is None or resolution_m_px <= 0 or not math.isfinite(resolution_m_px):
logger.error("Could not calculate meters per pixel for bbox calculation.")
return None
# Calculate the geographic width and height in meters corresponding to the target pixel dimensions
geographic_width_meters = target_pixel_width * resolution_m_px
geographic_height_meters = target_pixel_height * resolution_m_px
# Now, use pyproj to find the geographic coordinates of the corners of a bounding box
# centered on the input point with these geographic dimensions.
geodetic_calculator = pyproj.Geod(ellps="WGS84") # type: ignore
# Calculate points by projecting from the center along cardinal directions
half_width_meters = geographic_width_meters / 2.0
half_height_meters = geographic_height_meters / 2.0
# 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_height_meters
)
_, south_boundary_lat, _ = geodetic_calculator.fwd(
lons=center_longitude_deg, lats=center_latitude_deg, az=180.0, dist=half_height_meters
)
east_boundary_lon, _, _ = geodetic_calculator.fwd(
lons=center_longitude_deg, lats=center_latitude_deg, az=90.0, dist=half_width_meters
)
west_boundary_lon, _, _ = geodetic_calculator.fwd(
lons=center_longitude_deg, lats=center_latitude_deg, az=270.0, dist=half_width_meters
)
# Clamp latitude boundaries to WGS84 range
north_boundary_lat = min(90.0, max(-90.0, north_boundary_lat))
south_boundary_lat = min(90.0, max(-90.0, south_boundary_lat))
# Longitude wrapping should be handled by geod.fwd/inv, but ensure the order is correct for a bbox
# If west_lon > east_lon after projection, it implies crossing the antimeridian.
# The min/max approach might not correctly represent a bbox spanning the antimeridian.
# For simplicity, if west > east, we assume it's a normal bbox and swap,
# which is okay for areas not crossing the antimeridian widely.
# A more robust solution for antimeridian spanning areas would involve mercantile's bbox handling.
# Let's stick to the simpler approach assuming areas don't span the antimeridian widely.
logger.debug(
f"Calculated BBox for {target_pixel_width}x{target_pixel_height}px at zoom {zoom_level}: "
f"W={west_boundary_lon:.6f}, S={south_boundary_lat:.6f}, "
f"E={east_boundary_lon:.6f}, N={north_boundary_lat:.6f}"
)
# Check for invalid bbox (e.g. zero size)
if west_boundary_lon == east_boundary_lon or south_boundary_lat == north_boundary_lat:
logger.warning("Calculated zero-size geographic bbox. Returning None.")
return None
return (west_boundary_lon, south_boundary_lat, east_boundary_lon, north_boundary_lat)
except Exception as e_bbox_calc:
logger.exception(f"Error calculating geographic bbox from pixel size and zoom: {e_bbox_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
# MODIFIED: Added new utility function to calculate the combined geographic bounds from a list of tile infos.
# WHY: Needed for calculating the overall geographic extent of a group of DEM tiles for map display aspect ratio.
# HOW: Iterate through the list, get bounds for each tile, and find the min/max lat/lon.
def get_combined_geographic_bounds_from_tile_info_list(tile_info_list: List[Dict]) -> Optional[Tuple[float, float, float, float]]:
"""
Calculates the minimum bounding box that encompasses all tiles in the provided list.
Args:
tile_info_list: A list of tile information dictionaries (e.g., from ElevationManager.get_area_tile_info).
Each dict must contain 'latitude_coord' and 'longitude_coord'.
Returns:
A tuple (west_lon, south_lat, east_lon, north_lat) encompassing all tiles,
or None if the list is empty or invalid info is found.
"""
if not tile_info_list:
logger.warning("Tile info list is empty, cannot calculate combined geographic bounds.")
return None
min_lon_combined, min_lat_combined, max_lon_combined, max_lat_combined = 180.0, 90.0, -180.0, -90.0
initialized = False
for tile_info in tile_info_list:
lat_coord = tile_info.get("latitude_coord")
# MODIFIED: Corrected typo in key name.
# WHY: To correctly access the longitude coordinate from the dictionary.
# HOW: Changed "longitude_longitude" to "longitude_coord".
lon_coord = tile_info.get("longitude_coord")
if lat_coord is None or lon_coord is None:
logger.warning(f"Skipping tile info entry due to missing coordinates: {tile_info}")
continue # Skip this entry if coordinates are missing
try:
# Get the precise geographic bounds for this HGT tile
tile_bounds = get_hgt_tile_geographic_bounds(lat_coord, lon_coord)
if tile_bounds: # Ensure bounds were calculated successfully
if not initialized:
min_lon_combined, min_lat_combined, max_lon_combined, max_lat_combined = tile_bounds
initialized = True
else:
min_lon_combined = min(min_lon_combined, tile_bounds[0])
min_lat_combined = min(min_lat_combined, tile_bounds[1])
max_lon_combined = max(max_lon_combined, tile_bounds[2])
max_lat_combined = max(max_lat_combined, tile_bounds[3])
except Exception as e_get_tile_bounds:
logger.warning(f"Error getting geographic bounds for tile ({lat_coord},{lon_coord}): {e_get_tile_bounds}. Skipping this tile.")
continue # Skip tile with invalid bounds
if not initialized:
logger.warning("No valid tile coordinates found in the list to calculate combined bounds.")
return None # No valid tiles processed
# Final validation of combined bounds (e.g., if it spans the whole globe)
# The bounds from HGT tiles are 1x1 degree, so combining shouldn't create invalid wrap-around issues easily,
# but defensive check.
# A more robust check might be needed for edge cases near antimeridian, but for 1x1 tiles this is usually okay.
# For simplicity, assume min_lat < max_lat and min_lon < max_lon unless it crosses the antimeridian.
# The get_hgt_tile_geographic_bounds clamps, so min/max comparison should work generally.
# If the combined area crosses the antimeridian, min_lon_combined will be > max_lon_combined.
# We might need to represent this differently if we need a bbox spanning the antimeridian,
# but for simple min/max extent, the current logic is okay for non-global areas.
# Let's just check for egregious errors like south > north.
if min_lat_combined > max_lat_combined:
logger.warning(f"Calculated invalid combined latitude range: S={min_lat_combined}, N={max_lat_combined}. Returning None.")
return None
# Longitude range check is tricky with antimeridian. Assuming for typical areas the min/max works.
logger.debug(f"Calculated combined geographic bounds: ({min_lon_combined:.6f}, {min_lat_combined:.6f}, {max_lon_combined:.6f}, {max_lat_combined:.6f})")
return (min_lon_combined, min_lat_combined, max_lon_combined, max_lat_combined)
# MODIFIED: Added new utility function to calculate a geographic bounding box
# given a center point, desired pixel dimensions, and zoom level.
# WHY: Necessary for interactive zoom to determine the map area to fetch and display
# while maintaining a relatively constant pixel size for the map window.
# HOW: Uses the inverse of the meters per pixel calculation to find the geographic
# width and height corresponding to the desired pixel dimensions at the given
# latitude and zoom. Then uses pyproj to find the geographic coordinates
# of the corners of a bounding box centered on the input point with these
# calculated geographic width/height.
def calculate_geographic_bbox_from_pixel_size_and_zoom(
center_latitude_deg: float,
center_longitude_deg: float,
target_pixel_width: int,
target_pixel_height: int,
zoom_level: int,
tile_pixel_size: int = 256
) -> Optional[Tuple[float, float, float, float]]: # (west_lon, south_lat, east_lon, north_lat)
"""
Calculates the geographic bounding box centered on a point that corresponds
to a specific pixel size at a given zoom level and latitude.
Args:
center_latitude_deg: Latitude of the center point in degrees.
center_longitude_deg: Longitude of the center point in degrees.
target_pixel_width: Desired pixel width for the bounding box.
target_pixel_height: Desired pixel height for the bounding box.
zoom_level: Map zoom level.
tile_pixel_size: The size of one side of a map tile in pixels.
Returns:
A tuple (west_lon, south_lat, east_lon, north_lat) in degrees,
or None if calculation fails or dependencies are unavailable.
"""
if not PYPROJ_AVAILABLE:
logger.error(
"'pyproj' library is required for bbox calculation but is not found."
)
return None
if not MERCANTILE_AVAILABLE_UTILS:
logger.error(
"'mercantile' library is required for bbox calculation but is not found."
)
return None
if not (-90.0 <= center_latitude_deg <= 90.0):
logger.error(f"Invalid center_latitude_deg for bbox calc: {center_latitude_deg}")
return None
if not (-180.0 <= center_longitude_deg <= 180.0):
logger.error(f"Invalid center_longitude_deg for bbox calc: {center_longitude_deg}")
return None
if not (isinstance(target_pixel_width, int) and target_pixel_width > 0 and
isinstance(target_pixel_height, int) and target_pixel_height > 0):
logger.error(f"Invalid target pixel dimensions ({target_pixel_width}x{target_pixel_height}) for bbox calc.")
return None
if not (isinstance(zoom_level, int) and 0 <= zoom_level <= 25): # Clamp to practical max zoom
logger.warning(f"Invalid zoom level for bbox calc: {zoom_level}. Clamping to [0, 20].")
zoom_level = max(0, min(zoom_level, 20)) # Clamp
if not (isinstance(tile_pixel_size, int) and tile_pixel_size > 0):
logger.error(f"Invalid tile_pixel_size for bbox calc: {tile_pixel_size}")
return None
logger.debug(
f"Calculating geographic bbox for center ({center_latitude_deg:.6f}, {center_longitude_deg:.6f}) "
f"at zoom {zoom_level} for {target_pixel_width}x{target_pixel_height}px."
)
try:
# Calculate the ground resolution (meters per pixel) at the center latitude for the given zoom
resolution_m_px = calculate_meters_per_pixel(center_latitude_deg, zoom_level, tile_pixel_size)
if resolution_m_px is None or resolution_m_px <= 0 or not math.isfinite(resolution_m_px):
logger.error("Could not calculate meters per pixel for bbox calculation.")
return None
# Calculate the geographic width and height in meters corresponding to the target pixel dimensions
geographic_width_meters = target_pixel_width * resolution_m_px
geographic_height_meters = target_pixel_height * resolution_m_px
# Now, use pyproj to find the geographic coordinates of the corners of a bounding box
# centered on the input point with these geographic dimensions.
geodetic_calculator = pyproj.Geod(ellps="WGS84") # type: ignore
# Calculate points by projecting from the center along cardinal directions
half_width_meters = geographic_width_meters / 2.0
half_height_meters = geographic_height_meters / 2.0
# 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_height_meters
)
_, south_boundary_lat, _ = geodetic_calculator.fwd(
lons=center_longitude_deg, lats=center_latitude_deg, az=180.0, dist=half_height_meters
)
east_boundary_lon, _, _ = geodetic_calculator.fwd(
lons=center_longitude_deg, lats=center_latitude_deg, az=90.0, dist=half_width_meters
)
west_boundary_lon, _, _ = geodetic_calculator.fwd(
lons=center_longitude_deg, lats=center_latitude_deg, az=270.0, dist=half_width_meters
)
# Clamp latitude boundaries to WGS84 range
north_boundary_lat = min(90.0, max(-90.0, north_boundary_lat))
south_boundary_lat = min(90.0, max(-90.0, south_boundary_lat))
# Longitude wrapping should be handled by geod.fwd/inv, but ensure the order is correct for a bbox
# If west_lon > east_lon after projection, it implies crossing the antimeridian.
# The min/max approach might not correctly represent a bbox spanning the antimeridian.
# For simplicity, if west > east, we assume it's a normal bbox and swap,
# which is okay for areas not crossing the antimeridian widely.
# A more robust solution for antimeridian spanning areas would involve mercantile's bbox handling.
# Let's stick to the simpler approach assuming areas don't span the antimeridian widely.
logger.debug(
f"Calculated BBox for {target_pixel_width}x{target_pixel_height}px at zoom {zoom_level}: "
f"W={west_boundary_lon:.6f}, S={south_boundary_lat:.6f}, "
f"E={east_boundary_lon:.6f}, N={north_boundary_lat:.6f}"
)
# Check for invalid bbox (e.g. zero size)
if west_boundary_lon == east_boundary_lon or south_boundary_lat == north_boundary_lat:
logger.warning("Calculated zero-size geographic bbox. Returning None.")
return None
return (west_boundary_lon, south_boundary_lat, east_boundary_lon, north_boundary_lat) # Return as (west, south, east, north)
except Exception as e_bbox_calc:
logger.exception(f"Error calculating geographic bbox from pixel size and zoom: {e_bbox_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
# MODIFIED: Added new utility function to calculate the combined geographic bounds from a list of tile infos.
# WHY: Needed for calculating the overall geographic extent of a group of DEM tiles for map display aspect ratio.
# HOW: Iterate through the list, get bounds for each tile, and find the min/max lat/lon.
def get_combined_geographic_bounds_from_tile_info_list(tile_info_list: List[Dict]) -> Optional[Tuple[float, float, float, float]]:
"""
Calculates the minimum bounding box that encompasses all tiles in the provided list.
Args:
tile_info_list: A list of tile information dictionaries (e.g., from ElevationManager.get_area_tile_info).
Each dict must contain 'latitude_coord' and 'longitude_coord'.
Returns:
A tuple (west_lon, south_lat, east_lon, north_lat) encompassing all tiles,
or None if the list is empty or invalid info is found.
"""
if not tile_info_list:
logger.warning("Tile info list is empty, cannot calculate combined geographic bounds.")
return None
min_lon_combined, min_lat_combined, max_lon_combined, max_lat_combined = 180.0, 90.0, -180.0, -90.0
initialized = False
for tile_info in tile_info_list:
lat_coord = tile_info.get("latitude_coord")
# MODIFIED: Corrected typo in key name.
# WHY: To correctly access the longitude coordinate from the dictionary.
# HOW: Changed "longitude_longitude" to "longitude_coord".
lon_coord = tile_info.get("longitude_coord")
if lat_coord is None or lon_coord is None:
logger.warning(f"Skipping tile info entry due to missing coordinates: {tile_info}")
continue # Skip this entry if coordinates are missing
try:
# Get the precise geographic bounds for this HGT tile
tile_bounds = get_hgt_tile_geographic_bounds(lat_coord, lon_coord)
if tile_bounds: # Ensure bounds were calculated successfully
if not initialized:
min_lon_combined, min_lat_combined, max_lon_combined, max_lat_combined = tile_bounds
initialized = True
else:
min_lon_combined = min(min_lon_combined, tile_bounds[0])
min_lat_combined = min(min_lat_combined, tile_bounds[1])
max_lon_combined = max(max_lon_combined, tile_bounds[2])
max_lat_combined = max(max_lat_combined, tile_bounds[3])
except Exception as e_get_tile_bounds:
logger.warning(f"Error getting geographic bounds for tile ({lat_coord},{lon_coord}): {e_get_tile_bounds}. Skipping this tile.")
continue # Skip tile with invalid bounds
if not initialized:
logger.warning("No valid tile coordinates found in the list to calculate combined bounds.")
return None # No valid tiles processed
# Final validation of combined bounds (e.g., if it spans the whole globe)
# The bounds from HGT tiles are 1x1 degree, so combining shouldn't create invalid wrap-around issues easily,
# but defensive check.
# A more robust check might be needed for edge cases near antimeridian, but for 1x1 tiles this is usually okay.
# For simplicity, assume min_lat < max_lat and min_lon < max_lon unless it crosses the antimeridian.
# The get_hgt_tile_geographic_bounds clamps, so min/max comparison should work generally.
# If the combined area crosses the antimeridian, min_lon_combined will be > max_lon_combined.
# We might need to represent this differently if we need a bbox spanning the antimeridian,
# but for simple min/max extent, the current logic is okay for non-global areas.
# Let's just check for egregious errors like south > north.
if min_lat_combined > max_lat_combined:
logger.warning(f"Calculated invalid combined latitude range: S={min_lat_combined}, N={max_lat_combined}. Returning None.")
return None
# Longitude range check is tricky with antimeridian. Assuming for typical areas the min/max works.
logger.debug(f"Calculated combined geographic bounds: ({min_lon_combined:.6f}, {min_lat_combined:.6f}, {max_lon_combined:.6f}, {max_lat_combined:.6f})")
return (min_lon_combined, min_lat_combined, max_lon_combined, max_lat_combined)
Reference in New Issue View Git Blame Copy Permalink