Add a python script that generates a skyplot from a RINEX navigation file

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/
 

utils/skyplot/.gitignore

@@ -0,0 +1,46 @@
+# SPDX-License-Identifier: GPL-3.0-or-later
+# SPDX-FileCopyrightText: 2025 Carles Fernandez-Prades <carles.fernandez@cttc.es>
+
+# 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

utils/skyplot/README.md

@@ -0,0 +1,78 @@
+# GNSS Skyplot utility
+
+<!-- prettier-ignore-start -->
+[comment]: # (
+SPDX-License-Identifier: GPL-3.0-or-later
+)
+
+[comment]: # (
+SPDX-FileCopyrightText: 2025 Carles Fernandez-Prades <carles.fernandez@cttc.es>
+)
+<!-- prettier-ignore-end -->
+
+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 <RINEX_FILE> [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_<RINEX_FILE>.pdf` (with dots in
+`<RINEX_FILE>` replaced by `_`) with:
+
+- Satellite trajectories over all epochs in the file.
+- Color-coded by constellation.
+- Observer location in title.
+- Time range in footer.

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 <RINEX_NAV_FILE> [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 <RINEX_FILE> [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()