From e74c5708d7f9e99632ad508160f31be4027f7cdc Mon Sep 17 00:00:00 2001 From: Carles Fernandez Date: Wed, 30 Apr 2025 10:36:46 +0200 Subject: [PATCH] Add a python script that generates a skyplot from a RINEX navigation file --- docs/CHANGELOG.md | 6 + utils/skyplot/.gitignore | 46 ++++ utils/skyplot/README.md | 78 +++++++ utils/skyplot/skyplot.py | 475 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 605 insertions(+) create mode 100644 utils/skyplot/.gitignore create mode 100644 utils/skyplot/README.md create mode 100755 utils/skyplot/skyplot.py diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index c8d6a477f..8f8f86685 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -25,6 +25,12 @@ All notable changes to GNSS-SDR will be documented in this file. - Fix building option `-DENABLE_ION=ON` when using CMake >= 4.0. +### Improvements in Usability: + +- Added a GNSS skyplot visualization utility at `utils/skyplot/skyplot.py`, + which generates a skyplot from a RINEX navigation file and saves the image in + PDF format. It requires `numpy` and `matplotlib`. + See the definitions of concepts and metrics at https://gnss-sdr.org/design-forces/ diff --git a/utils/skyplot/.gitignore b/utils/skyplot/.gitignore new file mode 100644 index 000000000..120b9eecc --- /dev/null +++ b/utils/skyplot/.gitignore @@ -0,0 +1,46 @@ +# SPDX-License-Identifier: GPL-3.0-or-later +# SPDX-FileCopyrightText: 2025 Carles Fernandez-Prades + +# Temporary/backup files +*~ +*.bak +*.tmp +.DS_Store + +# RINEX Navigation files +*.??n +*.??p +*.??P +*.??g +*.??G +*.??h +*.??H +*.??q +*.??Q +*.??b +*.??B +*.??f +*.??F +*.??l +*.??L +*.rnx +*.RNX + +# RINEX v3+ navigation files +*_*.*n +*_*.*N +*_*.*p +*_*.*P + +# Compressed variants +*.gz +*.Z +*.zip + +# Output files +*.pdf +*.png +*.jpg +*.svg +*.eps +*.ps \ No newline at end of file diff --git a/utils/skyplot/README.md b/utils/skyplot/README.md new file mode 100644 index 000000000..3e97b2351 --- /dev/null +++ b/utils/skyplot/README.md @@ -0,0 +1,78 @@ +# GNSS Skyplot utility + + +[comment]: # ( +SPDX-License-Identifier: GPL-3.0-or-later +) + +[comment]: # ( +SPDX-FileCopyrightText: 2025 Carles Fernandez-Prades +) + + +A Python script that generates polar skyplots from RINEX navigation files, +showing satellite visibility over time. + +## Features + +- Processes RINEX navigation files. +- Calculates satellite positions using broadcast ephemeris. +- Plots satellite tracks in azimuth-elevation coordinates. +- Color-codes satellites by constellation (GPS, Galileo, GLONASS, BeiDou). +- Customizable observer location. +- Outputs high-quality image in PDF format. +- Non-interative mode for CI jobs (with `--no-show` flag). + +## Requirements + +- Python 3.6+ +- Required packages: + - `numpy` + - `matplotlib` + +## Usage + +### Basic Command + +``` +./skyplot.py [LATITUDE] [LONGITUDE] [ALTITUDE] [--no-show] +``` + +### Arguments + +| Argument | Type | Units | Description | Default | +| ------------ | -------- | ----------- | ---------------------- | -------- | +| `RINEX_FILE` | Required | - | RINEX nav file path | - | +| `LATITUDE` | Optional | degrees (°) | North/South position | 41.275°N | +| `LONGITUDE` | Optional | degrees (°) | East/West position | 1.9876°E | +| `ALTITUDE` | Optional | meters (m) | Height above sea level | 80.0 m | +| `--no-show` | Optional | - | Do not show plot | - | + +### Examples + +- Skyplot from default location (Castelldefels, Spain): + ``` + ./skyplot.py brdc0010.22n + ``` +- Skyplot from custom location (New York City, USA): + ``` + ./skyplot.py brdc0010.22n 40.7128 -74.0060 10.0 + ``` +- Skyplot from custom location (Santiago, Chile): + ``` + ./skyplot.py brdc0010.22n -33.4592 -70.6453 520.0 + ``` +- Non-interactive mode (for CI jobs): + ``` + ./skyplot.py brdc0010.22n -33.4592 -70.6453 520.0 --no-show + ``` + +## Output + +The script generates a PDF file named `skyplot_.pdf` (with dots in +`` replaced by `_`) with: + +- Satellite trajectories over all epochs in the file. +- Color-coded by constellation. +- Observer location in title. +- Time range in footer. diff --git a/utils/skyplot/skyplot.py b/utils/skyplot/skyplot.py new file mode 100755 index 000000000..8ff14dd5d --- /dev/null +++ b/utils/skyplot/skyplot.py @@ -0,0 +1,475 @@ +#!/usr/bin/env python +""" + skyplot.py + + Reads a RINEX navigation file and generates a skyplot + + Usage: python skyplot.py [observer_lat] [observer_lon] [observer_alt] + + ----------------------------------------------------------------------------- + + GNSS-SDR is a Global Navigation Satellite System software-defined receiver. + This file is part of GNSS-SDR. + + SPDX-FileCopyrightText: 2025 Carles Fernandez-Prades cfernandez(at)cttc.es + SPDX-License-Identifier: GPL-3.0-or-later + + ----------------------------------------------------------------------------- +""" + +import argparse +import re +import sys +from math import sin, cos, sqrt, atan2 +from datetime import datetime, timedelta +import numpy as np +import matplotlib.pyplot as plt + + +def parse_rinex_float(s): + """Parse RINEX formatted float string which may contain D or E exponent and compact spacing""" + # Handle empty string + if not s.strip(): + return 0.0 + + # Replace D exponent with E (some RINEX files use D instead of E) + s = s.replace('D', 'E').replace('d', 'e') + + # Handle cases where exponent lacks E (e.g., "12345-3") + if re.match(r'[+-]?\d+[+-]\d+', s.strip()): + s = s.replace('+', 'E+').replace('-', 'E-') + + try: + return float(s) + except ValueError: + # Handle cases where the number runs into the next field + # Try to split at the exponent if present + if 'E' in s: + base, exp = s.split('E')[:2] + # Take first character of exponent if needed + if exp and exp[0] in '+-' and len(exp) > 1: + return float(base + 'E' + exp[0] + exp[1:].split()[0]) + return 0.0 # Default if parsing fails + + +def read_rinex_nav(filename): + """Read RINEX v3.0 navigation file""" + satellites = {} + with open(filename, 'r', encoding='utf-8') as f: + # Skip header + while True: + line = f.readline() + if not line: + return satellites # Empty file + if "END OF HEADER" in line: + break + + # Read ephemeris data + current_line = f.readline() + while current_line: + if len(current_line) < 23: + current_line = f.readline() + continue + + prn = current_line[:3].strip() + + try: + # Parse epoch fields with careful position handling + year = int(current_line[4:8]) + month = int(current_line[9:11]) + day = int(current_line[12:14]) + hour = int(current_line[15:17]) + minute = int(current_line[18:20]) + second = int(float(current_line[21:23])) + + year += 2000 if year < 80 else 0 + epoch = datetime(year, month, day, hour, minute, second) + + # Read the next 7 lines + lines = [current_line] + for _ in range(7): + next_line = f.readline() + if not next_line: + break + lines.append(next_line) + + if len(lines) < 8: + current_line = f.readline() + continue + + # Parse all ephemeris parameters with robust float handling + ephemeris = { + 'prn': prn, + 'epoch': epoch, + 'sv_clock_bias': parse_rinex_float(lines[0][23:42]), + 'sv_clock_drift': parse_rinex_float(lines[0][42:61]), + 'sv_clock_drift_rate': parse_rinex_float(lines[0][61:80]), + 'iode': parse_rinex_float(lines[1][4:23]), + 'crs': parse_rinex_float(lines[1][23:42]), + 'delta_n': parse_rinex_float(lines[1][42:61]), + 'm0': parse_rinex_float(lines[1][61:80]), + 'cuc': parse_rinex_float(lines[2][4:23]), + 'ecc': parse_rinex_float(lines[2][23:42]), + 'cus': parse_rinex_float(lines[2][42:61]), + 'sqrt_a': parse_rinex_float(lines[2][61:80]), + 'toe': parse_rinex_float(lines[3][4:23]), + 'cic': parse_rinex_float(lines[3][23:42]), + 'omega0': parse_rinex_float(lines[3][42:61]), + 'cis': parse_rinex_float(lines[3][61:80]), + 'i0': parse_rinex_float(lines[4][4:23]), + 'crc': parse_rinex_float(lines[4][23:42]), + 'omega': parse_rinex_float(lines[4][42:61]), + 'omega_dot': parse_rinex_float(lines[4][61:80]), + 'idot': parse_rinex_float(lines[5][4:23]), + 'codes_l2': parse_rinex_float(lines[5][23:42]), + 'gps_week': parse_rinex_float(lines[5][42:61]), + 'l2p_flag': parse_rinex_float(lines[5][61:80]), + 'sv_accuracy': parse_rinex_float(lines[6][4:23]), + 'sv_health': parse_rinex_float(lines[6][23:42]), + 'tgd': parse_rinex_float(lines[6][42:61]), + 'iodc': parse_rinex_float(lines[6][61:80]), + 'transmission_time': parse_rinex_float(lines[7][4:23]), + 'fit_interval': ( + parse_rinex_float(lines[7][23:42])) if len(lines[7]) > 23 else 0.0 + } + + if prn not in satellites: + satellites[prn] = [] + satellites[prn].append(ephemeris) + + except (ValueError, IndexError) as e: + print(f"Error parsing PRN {prn} at {epoch}: {e}") + # Skip to next block by reading until next PRN + while current_line and not current_line.startswith(prn[0]): + current_line = f.readline() + continue + + current_line = f.readline() + + return satellites + + +def calculate_satellite_position(ephemeris, transmit_time): + """Calculate satellite position in ECEF coordinates at given transmission time""" + # Constants + mu = 3.986005e14 # Earth's gravitational constant (m^3/s^2) + omega_e_dot = 7.2921151467e-5 # Earth rotation rate (rad/s) + + # Semi-major axis + a = ephemeris['sqrt_a'] ** 2 + + # Time difference from ephemeris reference time + tk = transmit_time - ephemeris['toe'] + + # Corrected mean motion + n0 = sqrt(mu / (a ** 3)) + n = n0 + ephemeris['delta_n'] + + # Mean anomaly + mk = ephemeris['m0'] + n * tk + + # Solve Kepler's equation for eccentric anomaly (Ek) + ek = mk + for _ in range(10): + ek_old = ek + ek = mk + ephemeris['ecc'] * sin(ek) + if abs(ek - ek_old) < 1e-12: + break + + # True anomaly + nu_k = atan2(sqrt(1 - ephemeris['ecc']**2) * sin(ek), cos(ek) - ephemeris['ecc']) + + # Argument of latitude + phi_k = nu_k + ephemeris['omega'] + + # Second harmonic perturbations + delta_uk = ephemeris['cus'] * sin(2 * phi_k) + ephemeris['cuc'] * cos(2 * phi_k) + delta_rk = ephemeris['crs'] * sin(2 * phi_k) + ephemeris['crc'] * cos(2 * phi_k) + delta_ik = ephemeris['cis'] * sin(2 * phi_k) + ephemeris['cic'] * cos(2 * phi_k) + + # Corrected argument of latitude, radius and inclination + uk = phi_k + delta_uk + rk = a * (1 - ephemeris['ecc'] * cos(ek)) + delta_rk + ik = ephemeris['i0'] + delta_ik + ephemeris['idot'] * tk + + # Positions in orbital plane + xk_prime = rk * cos(uk) + yk_prime = rk * sin(uk) + + # Corrected longitude of ascending node + omega_k = ( + ephemeris['omega0'] + + (ephemeris['omega_dot'] - omega_e_dot) * tk + - omega_e_dot * ephemeris['toe'] + ) + + # Earth-fixed coordinates + xk = xk_prime * cos(omega_k) - yk_prime * cos(ik) * sin(omega_k) + yk = xk_prime * sin(omega_k) + yk_prime * cos(ik) * cos(omega_k) + zk = yk_prime * sin(ik) + + return xk, yk, zk + + +def calculate_satellite_positions(ephemeris, start_time, end_time, step_min=15): + """Generate multiple positions over time for a single satellite + between start_time and end_time. + """ + positions = [] + current_time = start_time + + while current_time <= end_time: + transmit_time = (current_time - ephemeris['epoch']).total_seconds() + + # Only calculate positions within ephemeris validity (typically 4 hours) + if abs(transmit_time) <= 4 * 3600: + x, y, z = calculate_satellite_position(ephemeris, transmit_time) + positions.append((current_time, x, y, z)) + + current_time += timedelta(minutes=step_min) + + return positions + + +def ecef_to_az_el(x, y, z, obs_lat, obs_lon, obs_alt): + """Convert ECEF coordinates to azimuth and elevation""" + # WGS-84 parameters + a = 6378137.0 # semi-major axis + e_sq = 6.69437999014e-3 # first eccentricity squared + + # Convert geodetic coordinates to ECEF + n = a / sqrt(1 - e_sq * sin(obs_lat)**2) + obs_x = (n + obs_alt) * cos(obs_lat) * cos(obs_lon) + obs_y = (n + obs_alt) * cos(obs_lat) * sin(obs_lon) + obs_z = (n * (1 - e_sq) + obs_alt) * sin(obs_lat) + + # Vector from observer to satellite + dx = x - obs_x + dy = y - obs_y + dz = z - obs_z + + # Convert to local ENU (East, North, Up) coordinates + enu_x = -sin(obs_lon) * dx + cos(obs_lon) * dy + enu_y = -sin(obs_lat) * cos(obs_lon) * dx - sin(obs_lat) * sin(obs_lon) * dy + cos(obs_lat) * dz + enu_z = cos(obs_lat) * cos(obs_lon) * dx + cos(obs_lat) * sin(obs_lon) * dy + sin(obs_lat) * dz + + # Calculate azimuth and elevation + azimuth = atan2(enu_x, enu_y) + elevation = atan2(enu_z, sqrt(enu_x**2 + enu_y**2)) + + # Convert to degrees and adjust azimuth to 0-360 + azimuth = np.degrees(azimuth) % 360 + elevation = np.degrees(elevation) + + return azimuth, elevation + + +def plot_satellite_tracks(satellites, obs_lat, obs_lon, obs_alt, + footer_text=None, filename=None, + show_plot=True): + """Plot trajectories for all visible satellites""" + plt.rcParams["font.family"] = "Times New Roman" + fig = plt.figure(figsize=(8, 8)) + ax = fig.add_subplot(111, projection='polar') + ax.tick_params(labelsize=16, pad=7) + + # Polar plot setup + ax.set_theta_zero_location('N') + ax.set_theta_direction(-1) + ax.set_ylim(0, 90) + + # Elevation ticks + ax.set_yticks(range(0, 91, 15)) + ax.set_yticklabels(['90°', '', '60°', '', '30°', '', '0°'], fontsize=14) + + # Color scheme by constellation + system_colors = { + 'G': 'blue', # GPS + 'E': 'green', # Galileo + 'R': 'red', # GLONASS + 'C': 'orange' # BeiDou + } + + for prn, ephemeris_list in satellites.items(): + color = system_colors.get(prn[0], 'purple') + + # Get the most recent ephemeris + if not ephemeris_list: + continue + ephemeris = max(ephemeris_list, key=lambda x: x['epoch']) + + # Calculate trajectory + all_epochs = sorted({e['epoch'] for prn_data in satellites.values() for e in prn_data}) + start_time = min(all_epochs) + end_time = max(all_epochs) + positions = calculate_satellite_positions(ephemeris, start_time, end_time) + + az = [] + el = [] + for _, x, y, z in positions: + azimuth, elevation = ecef_to_az_el(x, y, z, obs_lat, obs_lon, obs_alt) + if elevation > 0: # Above horizon + az.append(azimuth) + el.append(elevation) + + if len(az) > 1: + # Convert to polar coordinates + theta = np.radians(az) + r = 90 - np.array(el) + + # Plot trajectory + ax.plot(theta, r, '-', color=color, alpha=0.7, linewidth=2.5) + + # Label at midpoint + mid_idx = len(theta)//2 + ax.text(theta[mid_idx], r[mid_idx], prn, + fontsize=12, ha='center', va='center', + bbox={"facecolor": "white", "alpha": 0.8, "pad": 2}) + + # Add legend and metadata + legend_elements = [plt.Line2D([0], [0], marker='o', color='w', + label=f'{sys} ({name})', markerfacecolor=color, markersize=10) + for sys, (name, color) in [ + ('G', ('GPS', 'blue')), + ('E', ('Galileo', 'green')), + ('R', ('GLONASS', 'red')), + ('C', ('BeiDou', 'orange')) + ]] + + ax.legend(handles=legend_elements, loc='upper right', + bbox_to_anchor=(1.3, 1.1), fontsize=14) + + lat_deg = np.degrees(obs_lat) + lon_deg = np.degrees(obs_lon) + lat_hemisphere = 'N' if lat_deg >= 0 else 'S' + lon_hemisphere = 'E' if lon_deg >= 0 else 'W' + + plt.title( + f"GNSS skyplot from {abs(lat_deg):.2f}° {lat_hemisphere}, " + f"{abs(lon_deg):.2f}° {lon_hemisphere}", + pad=25, + fontsize=20 + ) + + if footer_text: + fig.text(0.42, 0.05, footer_text, ha='center', va='center', fontsize=16) + + plt.tight_layout() + + if filename: + filename_no_dots = filename.replace('.', '_') + output_name = f"skyplot_{filename_no_dots}.pdf" + else: + output_name = "skyplot.pdf" + + plt.savefig(output_name, format='pdf', bbox_inches='tight') + print(f"Image saved as {output_name}") + if show_plot: + plt.show() + else: + plt.close() + + +def main(): + """Generate the skyplot""" + +# Set up argument parser + parser = argparse.ArgumentParser( + description='Generate GNSS skyplot from RINEX navigation file', + add_help=False # We'll handle help manually to preserve your current style + ) + + # Add only the no-show flag + parser.add_argument( + '--no-show', + action='store_true', + help='Run without displaying plot window' + ) + + # Parse known args (this ignores other positional args) + args, remaining_args = parser.parse_known_args() + + # Handle help manually + if '-h' in remaining_args or '--help' in remaining_args: + print(""" +Usage: python skyplot.py [LATITUDE] [LONGITUDE] [ALTITUDE] [--no-show] + +Example: + python skyplot.py brdc0010.22n 41.275 1.9876 80.0 --no-show +""") + sys.exit(0) + + if len(remaining_args) < 1: + print("Error: RINEX file required") + sys.exit(1) + + filename = remaining_args[0] + + # Default observer location (Castelldefels, Barcelona, Spain) + obs_lat = np.radians(41.2750) + obs_lon = np.radians(1.9876) + obs_alt = 80.0 + + # Override with command line arguments if provided + if len(remaining_args) >= 4: + try: + obs_lat = np.radians(float(remaining_args[1])) + obs_lon = np.radians(float(remaining_args[2])) + if len(remaining_args) >= 5: + obs_alt = float(remaining_args[3]) + except ValueError: + print("Invalid observer coordinates. Using defaults.") + + # Read RINEX file + print(f"Reading {filename}...") + try: + satellites = read_rinex_nav(filename) + except FileNotFoundError: + print(f"Error: File '{filename}' not found.") + return + + if not satellites: + print("No satellite data found in the file.") + return + + # Print summary information + all_epochs = sorted(list(set( + e['epoch'] for prn, ephemerides in satellites.items() for e in ephemerides + ))) + print("\nFile contains:") + print(f"- {len(satellites)} unique satellites") + print(f"- {len(all_epochs)} unique epochs") + print(f"- From {all_epochs[0]} to {all_epochs[-1]}") + + # Calculate and print satellite counts by system + system_counts = {} + for prn in satellites: + system = prn[0] + system_counts[system] = system_counts.get(system, 0) + 1 + + print("\nSatellite systems found:") + for system, count in sorted(system_counts.items()): + system_name = { + 'G': 'GPS', + 'R': 'GLONASS', + 'E': 'Galileo', + 'C': 'BeiDou' + }.get(system, 'Unknown') + print(f"- {system_name} ({system}): {count} satellites") + + # Generate the combined skyplot + print("\nGenerating skyplot...") + footer = f"From {all_epochs[0]} to {all_epochs[-1]} UTC" + plot_satellite_tracks( + satellites, + obs_lat, + obs_lon, + obs_alt, + footer_text=footer, + filename=filename, + show_plot=not args.no_show + ) + + +if __name__ == "__main__": + main()