diff --git a/utils/skyplot/skyplot.py b/utils/skyplot/skyplot.py index 01661bd87..547a429be 100755 --- a/utils/skyplot/skyplot.py +++ b/utils/skyplot/skyplot.py @@ -1,21 +1,21 @@ #!/usr/bin/env python """ -skyplot.py + skyplot.py -Reads a RINEX navigation file and generates a skyplot. Optionally, a RINEX observation file can -also be read to match the skyplot to the receiver processing time. + Reads a RINEX navigation file and generates a skyplot. Optionally, a RINEX observation file can + also be read to match the skyplot to the receiver processing time. -Usage: python skyplot.py [observer_lat] [observer_lon] [observer_alt] [--use-obs] + Usage: python skyplot.py [observer_lat] [observer_lon] [observer_alt] [--use-obs] ------------------------------------------------------------------------------ + ----------------------------------------------------------------------------- -GNSS-SDR is a Global Navigation Satellite System software-defined receiver. -This file is part of GNSS-SDR. + 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 + SPDX-FileCopyrightText: 2025 Carles Fernandez-Prades cfernandez(at)cttc.es + SPDX-License-Identifier: GPL-3.0-or-later ------------------------------------------------------------------------------ + ----------------------------------------------------------------------------- """ import argparse @@ -30,20 +30,19 @@ import numpy as np def _read_obs_time_bounds(obs_path): - """ - Return (start_time, end_time) from a RINEX 2/3 OBS file by scanning epoch lines. + """Return (start_time, end_time) from a RINEX 2/3 OBS file by scanning epoch lines. If parsing fails or the file is not OBS, return (None, None). """ try: - with open(obs_path, "r", encoding="utf-8", errors="ignore") as f: + with open(obs_path, 'r', encoding='utf-8', errors='ignore') as f: # Detect OBS in header and skip to END OF HEADER is_obs = False for line in f: if "RINEX VERSION / TYPE" in line: # Robust OBS detection: - # - In RINEX 2/3 the file-type letter at col 21 (0-based idx 20) + # In RINEX 2/3 the file-type letter at col 21 (0-based idx 20) tchar = line[20:21] - if tchar == "O" or "OBSERVATION DATA" in line: + if tchar == 'O' or 'OBSERVATION DATA' in line: is_obs = True if "END OF HEADER" in line: break @@ -56,27 +55,17 @@ def _read_obs_time_bounds(obs_path): for line in f: if not line.strip(): continue - if line.startswith(">"): # RINEX 3 epoch line - yyyy = int(line[2:6]) - mm = int(line[7:9]) - dd = int(line[10:12]) - hh = int(line[13:15]) - mi = int(line[16:18]) - ss = float(line[19:29]) - epoch = datetime(yyyy, mm, dd, hh, mi, int(ss), int((ss % 1) * 1e6)) + if line.startswith('>'): # RINEX 3 epoch line + yyyy = int(line[2:6]); mm = int(line[7:9]); dd = int(line[10:12]) + hh = int(line[13:15]); mi = int(line[16:18]); ss = float(line[19:29]) + epoch = datetime(yyyy, mm, dd, hh, mi, int(ss), int((ss % 1)*1e6)) else: # RINEX 2 epoch line try: - yy = int(line[1:3]) - mm = int(line[4:6]) - dd = int(line[7:9]) - hh = int(line[10:12]) - mi = int(line[13:15]) - ss = float(line[15:26]) + yy = int(line[1:3]); mm = int(line[4:6]); dd = int(line[7:9]) + hh = int(line[10:12]); mi = int(line[13:15]); ss = float(line[15:26]) yyyy = 1900 + yy if yy >= 80 else 2000 + yy - epoch = datetime( - yyyy, mm, dd, hh, mi, int(ss), int((ss % 1) * 1e6) - ) + epoch = datetime(yyyy, mm, dd, hh, mi, int(ss), int((ss % 1)*1e6)) except Exception: continue if start_time is None or epoch < start_time: @@ -95,22 +84,22 @@ def parse_rinex_float(s): return 0.0 # Replace D exponent with E (some RINEX files use D instead of E) - s = s.replace("D", "E").replace("d", "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-") + 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] + 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]) + 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 @@ -118,7 +107,7 @@ def read_rinex_nav(filename): """Read RINEX v3.0 navigation file""" satellites = {} line_number = 0 - with open(filename, "r", encoding="utf-8") as f: + with open(filename, 'r', encoding='utf-8') as f: # Skip header while True: line = f.readline() @@ -154,7 +143,7 @@ def read_rinex_nav(filename): # Read the next lines lines = [current_line] - line_count = 4 if system == "R" else 7 + line_count = 4 if system == 'R' else 7 for _ in range(line_count): next_line = f.readline() line_number += 1 @@ -167,64 +156,61 @@ def read_rinex_nav(filename): line_number += 1 continue - if system == "R": # GLONASS specific parsing + if system == 'R': # GLONASS specific parsing ephemeris = { - "prn": prn, - "epoch": epoch, - "sv_clock_bias": parse_rinex_float(lines[0][23:42]), - "sv_relative_freq_bias": parse_rinex_float(lines[0][42:61]), - "message_frame_time": parse_rinex_float(lines[0][61:80]), - "x": parse_rinex_float(lines[1][4:23]), # Position (km) - "x_vel": parse_rinex_float(lines[1][23:42]), # Velocity (km/s) - "x_acc": parse_rinex_float(lines[1][42:61]), - "health": parse_rinex_float(lines[1][61:80]), - "y": parse_rinex_float(lines[2][4:23]), - "y_vel": parse_rinex_float(lines[2][23:42]), - "y_acc": parse_rinex_float(lines[2][42:61]), - "freq_num": parse_rinex_float(lines[2][61:80]), - "z": parse_rinex_float(lines[3][4:23]), - "z_vel": parse_rinex_float(lines[3][23:42]), - "z_acc": parse_rinex_float(lines[3][42:61]), - "age": parse_rinex_float(lines[3][61:80]), + 'prn': prn, + 'epoch': epoch, + 'sv_clock_bias': parse_rinex_float(lines[0][23:42]), + 'sv_relative_freq_bias': parse_rinex_float(lines[0][42:61]), + 'message_frame_time': parse_rinex_float(lines[0][61:80]), + 'x': parse_rinex_float(lines[1][4:23]), # Position (km) + 'x_vel': parse_rinex_float(lines[1][23:42]), # Velocity (km/s) + 'x_acc': parse_rinex_float(lines[1][42:61]), + 'health': parse_rinex_float(lines[1][61:80]), + 'y': parse_rinex_float(lines[2][4:23]), + 'y_vel': parse_rinex_float(lines[2][23:42]), + 'y_acc': parse_rinex_float(lines[2][42:61]), + 'freq_num': parse_rinex_float(lines[2][61:80]), + 'z': parse_rinex_float(lines[3][4:23]), + 'z_vel': parse_rinex_float(lines[3][23:42]), + 'z_acc': parse_rinex_float(lines[3][42:61]), + 'age': parse_rinex_float(lines[3][61:80]) } else: # Parse all ephemeris parameters 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 - ), + '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: @@ -252,60 +238,52 @@ def read_rinex_nav(filename): def calculate_satellite_position(ephemeris, transmit_time): """Calculate satellite position in ECEF coordinates at given transmission time""" - system = ephemeris["prn"][0] + system = ephemeris['prn'][0] - if system == "R": # GLONASS - use position + velocity * time + if system == 'R': # GLONASS - use position + velocity * time dt = transmit_time # Convert km to meters - xk = ( - ephemeris["x"] + ephemeris["x_vel"] * dt + 0.5 * ephemeris["x_acc"] * dt**2 - ) * 1000 - yk = ( - ephemeris["y"] + ephemeris["y_vel"] * dt + 0.5 * ephemeris["y_acc"] * dt**2 - ) * 1000 - zk = ( - ephemeris["z"] + ephemeris["z_vel"] * dt + 0.5 * ephemeris["z_acc"] * dt**2 - ) * 1000 + xk = (ephemeris['x'] + ephemeris['x_vel'] * dt + 0.5 * ephemeris['x_acc'] * dt**2) * 1000 + yk = (ephemeris['y'] + ephemeris['y_vel'] * dt + 0.5 * ephemeris['y_acc'] * dt**2) * 1000 + zk = (ephemeris['z'] + ephemeris['z_vel'] * dt + 0.5 * ephemeris['z_acc'] * dt**2) * 1000 else: # 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 + a = ephemeris['sqrt_a'] ** 2 # Corrected mean motion - n0 = sqrt(mu / (a**3)) - n = n0 + ephemeris["delta_n"] + n0 = sqrt(mu / (a ** 3)) + n = n0 + ephemeris['delta_n'] # Mean anomaly - mk = ephemeris["m0"] + n * transmit_time + mk = ephemeris['m0'] + n * transmit_time # Solve Kepler's equation for eccentric anomaly (Ek) ek = mk for _ in range(10): ek_old = ek - ek = mk + ephemeris["ecc"] * sin(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"] - ) + nu_k = atan2(sqrt(1 - ephemeris['ecc']**2) * sin(ek), cos(ek) - ephemeris['ecc']) # Argument of latitude - phi_k = nu_k + ephemeris["omega"] + 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) + 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"] * transmit_time + rk = a * (1 - ephemeris['ecc'] * cos(ek)) + delta_rk + ik = ephemeris['i0'] + delta_ik + ephemeris['idot'] * transmit_time # Positions in orbital plane xk_prime = rk * cos(uk) @@ -313,9 +291,9 @@ def calculate_satellite_position(ephemeris, transmit_time): # Corrected longitude of ascending node omega_k = ( - ephemeris["omega0"] - + (ephemeris["omega_dot"] - omega_e_dot) * transmit_time - - omega_e_dot * ephemeris["toe"] + ephemeris['omega0'] + + (ephemeris['omega_dot'] - omega_e_dot) * transmit_time + - omega_e_dot * ephemeris['toe'] ) # Earth-fixed coordinates @@ -332,10 +310,10 @@ def calculate_satellite_positions(ephemeris, start_time, end_time, step_min=5): """ positions = [] current_time = start_time - system = ephemeris["prn"][0] - max_valid_time = 1800 if system == "R" else 14400 + system = ephemeris['prn'][0] + max_valid_time = 1800 if system == 'R' else 14400 while current_time <= end_time: - transmit_time = (current_time - ephemeris["epoch"]).total_seconds() + transmit_time = (current_time - ephemeris['epoch']).total_seconds() if abs(transmit_time) <= max_valid_time: x, y, z = calculate_satellite_position(ephemeris, transmit_time) @@ -353,7 +331,7 @@ def ecef_to_az_el(x, y, z, obs_lat, obs_lon, obs_alt): e_sq = 6.69437999014e-3 # first eccentricity squared # Convert geodetic coordinates to ECEF - n = a / sqrt(1 - e_sq * sin(obs_lat) ** 2) + 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) @@ -365,16 +343,8 @@ def ecef_to_az_el(x, y, z, obs_lat, obs_lon, obs_alt): # 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 - ) + 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) @@ -387,52 +357,45 @@ def ecef_to_az_el(x, y, z, obs_lat, obs_lon, obs_alt): return azimuth, elevation -def plot_satellite_tracks( - satellites, - obs_lat, - obs_lon, - obs_alt, - footer_text=None, - filename=None, - show_plot=True, - start_time=None, - end_time=None, -): +def plot_satellite_tracks(satellites, obs_lat, obs_lon, obs_alt, + footer_text=None, filename=None, + show_plot=True, start_time=None, + end_time=None): """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 = 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_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) + 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 - "J": "brown", # QZSS - "I": "pink", # IRNSS - "S": "gray", # SBAS + 'G': 'blue', # GPS + 'E': 'green', # Galileo + 'R': 'red', # GLONASS + 'C': 'orange', # BeiDou + 'J': 'brown', # QZSS + 'I': 'pink', # IRNSS + 'S': 'gray' # SBAS } # System names mapping system_names = { - "G": "GPS", - "E": "Galileo", - "R": "GLONASS", - "C": "BeiDou", - "J": "QZSS", - "I": "IRNSS", - "S": "SBAS", + 'G': 'GPS', + 'E': 'Galileo', + 'R': 'GLONASS', + 'C': 'BeiDou', + 'J': 'QZSS', + 'I': 'IRNSS', + 'S': 'SBAS' } # Find which systems are actually present @@ -440,32 +403,23 @@ def plot_satellite_tracks( # Plot each satellite for prn, ephemeris_list in satellites.items(): - color = system_colors.get( - prn[0], "purple" - ) # Default to purple for unknown systems + color = system_colors.get(prn[0], 'purple') # Default to purple for unknown systems # Get the most recent ephemeris if not ephemeris_list: continue mid_time = start_time + (end_time - start_time) / 2 - prev_eph = [ - e for e in ephemeris_list if e["epoch"] <= mid_time - ] # Previous ephemeris + prev_eph = [e for e in ephemeris_list if e['epoch'] <= mid_time] # Previous ephemeris if prev_eph: # Pick the most recent ephemeris before or at the midpoint - ephemeris = max(prev_eph, key=lambda e: e["epoch"]) - else: # pick the ephemeris whose epoch is closest in time to the midpoint - ephemeris = min( - ephemeris_list, - key=lambda e: abs((e["epoch"] - mid_time).total_seconds()), - ) + ephemeris = max(prev_eph, key=lambda e: e['epoch']) + else: # Pick the ephemeris whose epoch is closest in time to the midpoint + ephemeris = min(ephemeris_list, key=lambda e: abs((e['epoch'] - mid_time).total_seconds())) # Calculate trajectory if start_time is None or end_time is None: - all_epochs = sorted( - {e["epoch"] for prn_data in satellites.values() for e in prn_data} - ) + 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) @@ -485,7 +439,7 @@ def plot_satellite_tracks( r = 90 - np.array(el) # Plot trajectory - ax.plot(theta, r, "-", color=color, alpha=0.7, linewidth=2.5) + ax.plot(theta, r, '-', color=color, alpha=0.7, linewidth=2.5) # Add arrow at last point if len(theta) >= 2: # Need at least 2 points for direction @@ -496,43 +450,30 @@ def plot_satellite_tracks( arrow_length_factor = 1.3 extended_theta = theta[-2] + dx * arrow_length_factor extended_r = r[-2] + dy * arrow_length_factor - ax.annotate( - "", - xytext=(theta[-1], r[-1]), - xy=(extended_theta, extended_r), - arrowprops={ - "arrowstyle": "->", - "color": color, - "alpha": 0.9, - "linewidth": 1.5, - "shrinkA": 0, - "shrinkB": 0, - }, - ) + ax.annotate('', + xytext=(theta[-1], r[-1]), + xy=(extended_theta, extended_r), + arrowprops={ + 'arrowstyle': '->', + 'color': color, + 'alpha': 0.9, + 'linewidth': 1.5, + 'shrinkA': 0, + 'shrinkB': 0 + }) # 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}, - ) + 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}) # Create legend elements only for present systems legend_elements = [ - plt.Line2D( - [0], - [0], - marker="o", - color="w", - label=f"{system_names[sys]} ({sys})", - markerfacecolor=system_colors[sys], - markersize=10, - ) + plt.Line2D([0], [0], marker='o', color='w', + label=f'{system_names[sys]} ({sys})', + markerfacecolor=system_colors[sys], + markersize=10) for sys in present_systems ] @@ -540,36 +481,36 @@ def plot_satellite_tracks( if legend_elements: ax.legend( handles=legend_elements, - loc="upper right", + loc='upper right', bbox_to_anchor=(1.3, 1.1), - fontsize=14, + 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" + 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, + fontsize=20 ) if footer_text: - fig.text(0.42, 0.05, footer_text, ha="center", va="center", fontsize=16) + fig.text(0.42, 0.05, footer_text, ha='center', va='center', fontsize=16) plt.tight_layout() if filename: filename_no_path = Path(filename).name - filename_no_dots = filename_no_path.replace(".", "_") + filename_no_dots = filename_no_path.replace('.', '_') output_name = f"skyplot_{filename_no_dots}.pdf" else: output_name = "skyplot.pdf" - plt.savefig(output_name, format="pdf", bbox_inches="tight") + plt.savefig(output_name, format='pdf', bbox_inches='tight') print(f"Image saved as {output_name}") if show_plot: plt.show() @@ -581,33 +522,34 @@ def main(): """Generate the skyplot""" # Set up argument parser parser = argparse.ArgumentParser( - description="Generate GNSS skyplot from RINEX navigation file", add_help=False + description='Generate GNSS skyplot from RINEX navigation file', + add_help=False ) # Add the no-show flag parser.add_argument( - "--no-show", action="store_true", help="Run without displaying plot window" + '--no-show', + action='store_true', + help='Run without displaying plot window' ) # Add the observation filename parser.add_argument( - "--use-obs", - action="store_true", - help="Use corresponding RINEX observation file to bound the skyplot to the receiver time window", + '--use-obs', + action='store_true', + help='Use corresponding RINEX observation file to bound the skyplot to the receiver time 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( - """ + if '-h' in remaining_args or '--help' in remaining_args: + print(""" Usage: python skyplot.py [LATITUDE] [LONGITUDE] [ALTITUDE] [--use-obs] [--no-show] Example: python skyplot.py brdc0010.22n 41.275 1.9876 80.0 --use-obs --no-show -""" - ) +""") sys.exit(0) if len(remaining_args) < 1: @@ -644,15 +586,9 @@ Example: return # Print summary information - all_epochs = sorted( - list( - set( - e["epoch"] - for prn, ephemerides in satellites.items() - for e in ephemerides - ) - ) - ) + 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") @@ -666,9 +602,12 @@ Example: 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" - ) + system_name = { + 'G': 'GPS', + 'R': 'GLONASS', + 'E': 'Galileo', + 'C': 'BeiDou' + }.get(system, 'Unknown') print(f"- {system_name} ({system}): {count} satellites") # Generate the combined skyplot @@ -679,7 +618,7 @@ Example: obs_path = None stem = filename[:-1] - for s in ("O", "o"): # try uppercase then lowercase + for s in ('O', 'o'): # Try uppercase then lowercase candidate = stem + s tried.append(candidate) if Path(candidate).exists(): @@ -690,17 +629,11 @@ Example: obs_start, obs_end = _read_obs_time_bounds(obs_path) if obs_start and obs_end: use_start, use_end = obs_start, obs_end - print( - f"\nObservation window detected from {obs_path}: {use_start} → {use_end}" - ) + print(f"\nObservation window detected from {obs_path}: {use_start} → {use_end}") else: - print( - f"\nWarning: Could not read valid times from {obs_path}. Using NAV span instead." - ) + print(f"\nWarning: Could not read valid times from {obs_path}. Using NAV span instead.") else: - print( - f"\nOBS file not found. Tried: {', '.join(tried)}. Using NAV span instead." - ) + print(f"\nOBS file not found. Tried: {', '.join(tried)}. Using NAV span instead.") # Ensure at least two samples with the default 5-minute step if (use_end - use_start) < timedelta(minutes=5): @@ -719,7 +652,7 @@ Example: filename=filename, show_plot=not args.no_show, start_time=use_start, - end_time=use_end, + end_time=use_end )