From acb1b541b5efed6b0bda0c7b98037b42996b85fd Mon Sep 17 00:00:00 2001 From: Carles Fernandez Date: Mon, 25 Aug 2025 09:34:37 +0200 Subject: [PATCH] skyplot: add usage of RINEX obs file approx position if found --- utils/skyplot/README.md | 7 ++++ utils/skyplot/skyplot.py | 82 ++++++++++++++++++++++++++++++++++++---- 2 files changed, 81 insertions(+), 8 deletions(-) diff --git a/utils/skyplot/README.md b/utils/skyplot/README.md index 6adf7a284..d6578ccbf 100644 --- a/utils/skyplot/README.md +++ b/utils/skyplot/README.md @@ -30,6 +30,13 @@ showing satellite visibility over time. to the receiver observation time via the`--use-obs` optional argument. - If this argument is set, the tool looks for a matching file following standard RINEX naming conventions, and uses it if found. + - If this argument is set and the observation file is found, the position + selection logic is: + 1. If user provides latitude, longitude, and altitude as positional + arguments, that position is always used. + 2. Otherwise, if the provided observation file contains a valid + `APPROX POSITION XYZ` field in its header, that position is used. + 3. Otherwise, the default position is used. ## Requirements diff --git a/utils/skyplot/skyplot.py b/utils/skyplot/skyplot.py index d0721a1db..90b19f6c5 100755 --- a/utils/skyplot/skyplot.py +++ b/utils/skyplot/skyplot.py @@ -43,6 +43,11 @@ except ImportError: __version__ = "1.0.0" +# Default position: Castelldefels, Barcelona +DEFAULT_LAT = 41.275 +DEFAULT_LON = 1.9876 +DEFAULT_ALT = 80.0 + def read_obs_time_bounds(obs_path: str) -> Tuple[Optional[datetime], Optional[datetime]]: """ Return (start_time, end_time) from a RINEX observation file (v2/3/4) @@ -188,6 +193,52 @@ def find_obs_for_nav(nav_file: str) -> Optional[str]: return None +def ecef_to_geodetic(x: float, y: float, z: float) -> Tuple[float, float, float]: + """Convert ECEF (X, Y, Z) to (lat, lon, h).""" + # WGS84 constants + a = 6378137.0 # semi-major axis + f = 1 / 298.257223563 + e2 = f * (2 - f) + + lon = atan2(y, x) + r = sqrt(x*x + y*y) + lat = atan2(z, r * (1 - e2)) # initial guess + + # Iterative improvement + for _ in range(5): + N = a / sqrt(1 - e2 * sin(lat)**2) + h = r / cos(lat) - N + lat = atan2(z, r * (1 - e2 * (N / (N + h)))) + + N = a / sqrt(1 - e2 * sin(lat)**2) + h = r / cos(lat) - N + + return np.degrees(lat), np.degrees(lon), h + + +def get_approx_position_from_obs(obs_file: str) -> Optional[Tuple[float, float, float]]: + """Read APPROX POSITION XYZ from an OBS RINEX file header. + Returns (lat, lon, h) in degrees/meters or None if not valid.""" + try: + with open(obs_file, "r") as f: + for line in f: + if "APPROX POSITION XYZ" in line: + parts = line.split() + if len(parts) >= 3: + try: + x, y, z = map(float, parts[0:3]) + except ValueError: + return None + # Check if position is valid (not all zeros) + if abs(x) < 1e-6 and abs(y) < 1e-6 and abs(z) < 1e-6: + return None + # Convert ECEF → geodetic + return ecef_to_geodetic(x, y, z) + return None + except OSError: + return None + + def parse_rinex_float(s: str) -> float: """Parse RINEX formatted float string which may contain D or E exponent and compact spacing""" # Handle empty string @@ -398,7 +449,7 @@ def read_rinex_nav(filename): if ephemeris: satellites.setdefault(prn, []).append(ephemeris) - except (ValueError, IndexError) as e: + except (ValueError, IndexError): # Skip malformed block; advance current_line = f.readline() line_number += 1 @@ -857,18 +908,18 @@ def main(): ) parser.add_argument( - 'lat', nargs='?', type=float, default=41.2750, - help='Observer latitude in degrees (default: 41.275° N)' + 'lat', nargs='?', type=float, default=DEFAULT_LAT, + help=f'Observer latitude in degrees (default: {DEFAULT_LAT}° N)' ) parser.add_argument( - 'lon', nargs='?', type=float, default=1.9876, - help='Observer longitude in degrees (default: 1.9876° E)' + 'lon', nargs='?', type=float, default=DEFAULT_LON, + help=f'Observer longitude in degrees (default: {DEFAULT_LON}° E)' ) parser.add_argument( - 'alt', nargs='?', type=float, default=80.0, - help='Observer altitude in meters (default: 80.0 m)' + 'alt', nargs='?', type=float, default=DEFAULT_ALT, + help=f'Observer altitude in meters (default: {DEFAULT_ALT} m)' ) # Optional arguments @@ -915,11 +966,18 @@ def main(): # Parse all arguments with full validation args = parser.parse_args() + filename = args.filename + user_provided_position = ( + (args.lat != DEFAULT_LAT) or + (args.lon != DEFAULT_LON) or + (args.alt != DEFAULT_ALT) + ) + # Convert coordinates to radians obs_lat = np.radians(args.lat) obs_lon = np.radians(args.lon) obs_alt = args.alt - filename = args.filename + # Read RINEX file print(f"Reading {filename} ...") @@ -989,6 +1047,14 @@ def main(): print(f"\nObservation window detected in {obs_path}: from {use_start} to {use_end}") else: print(f"\nWarning: Could not read valid times from {obs_path}. Using NAV span instead.") + if not user_provided_position: + approx_pos = get_approx_position_from_obs(obs_path) + if approx_pos is not None: + lat, lon, h = approx_pos + obs_lat = np.radians(lat) + obs_lon = np.radians(lon) + obs_alt = h + print(f"\nObserver position found in {obs_path}: lat {lat:.2f}°, lon {lon:.2f}°, height {h:.2f} m.") # Ensure at least two samples with the default 5-minute step if (use_end - use_start) < timedelta(minutes=5):