Improve Python observables plots

This commit is contained in:
Carles Fernandez
2026-06-26 10:36:12 +02:00
parent 60a1c43f2b
commit d332fdbb21
2 changed files with 158 additions and 81 deletions
+150 -59
View File
@@ -27,6 +27,7 @@ from lib.gnss_sdr_conf import (
ConfigError,
add_conf_argument,
load_gnss_sdr_conf,
signal_pretty_name,
)
from lib.plot_format import add_output_format_argument, apply_publication_style
from lib.read_hybrid_observables_dump import read_hybrid_observables_dump
@@ -104,6 +105,20 @@ def apply_conf_defaults(args):
else:
args.channels = DEFAULT_CHANNELS
args.channel_signals = channel_signal_labels(conf, args.channels)
def channel_signal_labels(conf, channels):
labels = [""] * channels
if conf is None:
return labels
for signal in conf.enabled_signals:
for channel in signal.channels:
if 0 <= channel < channels:
labels[channel] = signal.signal
return labels
def first_valid_observable(gnss_observables, channels):
min_tow_idx = None
@@ -122,6 +137,100 @@ def first_valid_observable(gnss_observables, channels):
return min_tow_idx, obs_idx
def valid_channel_data(gnss_observables, channel, observable):
valid_mask = np.array(gnss_observables["valid"][channel]) > 0
rx_time = np.array(gnss_observables["RX_time"][channel])
values = np.array(gnss_observables[observable][channel])
return rx_time[valid_mask], values[valid_mask]
def valid_channel_samples(gnss_observables, channel, observable):
valid_mask = np.array(gnss_observables["valid"][channel]) > 0
rx_time = np.array(gnss_observables["RX_time"][channel])
values = np.array(gnss_observables[observable][channel])
prns = np.array(gnss_observables["PRN"][channel])
return rx_time[valid_mask], values[valid_mask], prns[valid_mask]
def valid_time_bounds(gnss_observables, channels):
valid_times = []
for channel in range(channels):
rx_time, _ = valid_channel_data(gnss_observables, channel, "valid")
if len(rx_time) > 0:
valid_times.append(rx_time)
if not valid_times:
raise ValueError("No valid observables found in the dump.")
valid_times = np.concatenate(valid_times)
return valid_times.min() - 100, valid_times.max() + 100
def series_by_tracked_prn(gnss_observables, channels, observable, channel_signals):
parts = {}
for channel in range(channels):
rx_time, values, prns = valid_channel_samples(
gnss_observables, channel, observable
)
signal = channel_signals[channel] if channel < len(channel_signals) else ""
for prn in np.unique(prns):
if int(prn) == 0:
continue
selected = prns == prn
parts.setdefault((signal, int(prn)), []).append(
(rx_time[selected], values[selected])
)
series = {}
for series_key, chunks in parts.items():
rx_time = np.concatenate([chunk[0] for chunk in chunks])
values = np.concatenate([chunk[1] for chunk in chunks])
order = np.argsort(rx_time, kind="stable")
series[series_key] = (rx_time[order], values[order])
return series
def tracked_prn_label(series_key):
signal, prn = series_key
if signal:
return f"{signal_pretty_name(signal)} PRN {prn}"
return f"PRN {prn}"
def plot_observable_by_tracked_prn(
gnss_observables,
args,
*,
observable,
title,
ylabel,
window_title,
output_name,
time_start,
time_end,
):
plt.figure()
plt.title(title)
series = series_by_tracked_prn(
gnss_observables, args.channels, observable, args.channel_signals
)
for series_key in sorted(series):
rx_time, values = series[series_key]
plt.scatter(
rx_time,
values,
s=1,
label=tracked_prn_label(series_key),
)
plt.xlim(time_start, time_end)
plt.xlabel("TOW [s]")
plt.ylabel(ylabel)
plt.grid(True)
plt.legend()
plt.gcf().canvas.manager.set_window_title(window_title)
save_figure(args.fig_path, output_name, args.show, args.output_format)
def save_figure(fig_path, name, show, output_format):
plt.tight_layout()
plt.savefig(fig_path / f"{name}.{output_format}")
@@ -140,74 +249,56 @@ def main():
directory, base = resolve_dump_prefix(args.file_prefix, args.input_path)
observables_file = directory / f"{base}.dat"
gnss_observables = read_hybrid_observables_dump(args.channels, observables_file)
min_tow_idx, obs_idx = first_valid_observable(gnss_observables, args.channels)
time_start = gnss_observables["RX_time"][obs_idx][min_tow_idx] - 100
time_end = gnss_observables["RX_time"][obs_idx][-1] + 100
first_valid_observable(gnss_observables, args.channels)
time_start, time_end = valid_time_bounds(gnss_observables, args.channels)
plt.figure()
plt.title("Pseudorange")
for channel in range(args.channels):
plt.scatter(
gnss_observables["RX_time"][channel][min_tow_idx:],
gnss_observables["Pseudorange_m"][channel][min_tow_idx:],
s=1,
label=f"Channel {channel}",
)
plt.xlim(time_start, time_end)
plt.grid(True)
plt.xlabel("TOW [s]")
plt.ylabel("Pseudorange [m]")
plt.legend()
plt.gcf().canvas.manager.set_window_title("Pseudorange.png")
save_figure(args.fig_path, "Pseudorange", args.show, args.output_format)
plot_observable_by_tracked_prn(
gnss_observables,
args,
observable="Pseudorange_m",
title="Pseudorange",
ylabel="Pseudorange [m]",
window_title="Pseudorange.png",
output_name="Pseudorange",
time_start=time_start,
time_end=time_end,
)
plt.figure()
plt.title("Carrier Phase")
for channel in range(args.channels):
plt.scatter(
gnss_observables["RX_time"][channel][min_tow_idx:],
gnss_observables["Carrier_phase_hz"][channel][min_tow_idx:],
s=1,
label=f"Channel {channel}",
)
plt.xlim(time_start, time_end)
plt.xlabel("TOW [s]")
plt.ylabel("Accumulated Carrier Phase [cycles]")
plt.grid(True)
plt.legend()
plt.gcf().canvas.manager.set_window_title("AccumulatedCarrierPhase.png")
save_figure(args.fig_path, "AccumulatedCarrierPhase", args.show, args.output_format)
plot_observable_by_tracked_prn(
gnss_observables,
args,
observable="Carrier_phase_hz",
title="Carrier Phase",
ylabel="Accumulated Carrier Phase [cycles]",
window_title="AccumulatedCarrierPhase.png",
output_name="AccumulatedCarrierPhase",
time_start=time_start,
time_end=time_end,
)
plt.figure()
plt.title("Doppler Effect")
for channel in range(args.channels):
plt.scatter(
gnss_observables["RX_time"][channel][min_tow_idx:],
gnss_observables["Carrier_Doppler_hz"][channel][min_tow_idx:],
s=1,
label=f"Channel {channel}",
)
plt.xlim(time_start, time_end)
plt.xlabel("TOW [s]")
plt.ylabel("Doppler Frequency [Hz]")
plt.grid(True)
plt.legend()
plt.gcf().canvas.manager.set_window_title("DopplerFrequency.png")
save_figure(args.fig_path, "DopplerFrequency", args.show, args.output_format)
plot_observable_by_tracked_prn(
gnss_observables,
args,
observable="Carrier_Doppler_hz",
title="Doppler Effect",
ylabel="Doppler Frequency [Hz]",
window_title="DopplerFrequency.png",
output_name="DopplerFrequency",
time_start=time_start,
time_end=time_end,
)
plt.figure()
plt.title("GNSS Channels captured")
for channel in range(args.channels):
label = "unknown"
for prn in gnss_observables["PRN"][channel][min_tow_idx:]:
if int(prn) != 0:
label = str(int(prn))
break
rx_time, prns = valid_channel_data(gnss_observables, channel, "PRN")
if len(rx_time) == 0:
continue
plt.scatter(
gnss_observables["RX_time"][channel][min_tow_idx:],
gnss_observables["PRN"][channel][min_tow_idx:],
rx_time,
prns,
s=1,
label=f"PRN {channel} = {label}",
label=f"Channel {channel}",
)
plt.xlim(time_start, time_end)
plt.xlabel("TOW [s]")
@@ -29,13 +29,13 @@ def read_hybrid_observables_dump(channels, filename):
double_size_bytes = 8
bytes_shift = 0
RX_time = [[] for _ in range(channels+1)]
d_TOW_at_current_symbol = [[] for _ in range(channels+1)]
Carrier_Doppler_hz = [[] for _ in range(channels+1)]
Carrier_phase_hz = [[] for _ in range(channels+1)]
Pseudorange_m = [[] for _ in range(channels+1)]
PRN = [[] for _ in range(channels+1)]
valid = [[] for _ in range(channels+1)]
RX_time = [[] for _ in range(channels)]
d_TOW_at_current_symbol = [[] for _ in range(channels)]
Carrier_Doppler_hz = [[] for _ in range(channels)]
Carrier_phase_hz = [[] for _ in range(channels)]
Pseudorange_m = [[] for _ in range(channels)]
PRN = [[] for _ in range(channels)]
valid = [[] for _ in range(channels)]
f = open(filename, 'rb')
if f is None:
@@ -43,8 +43,7 @@ def read_hybrid_observables_dump(channels, filename):
else:
while True:
try:
# There is an empty channel at the end (Channel-6)
for N in range(0, channels+1):
for N in range(0, channels):
f.seek(bytes_shift, 0)
RX_time[N].append(struct.unpack(
@@ -86,19 +85,6 @@ def read_hybrid_observables_dump(channels, filename):
# Reached a partial record at end of file: stop reading.
break
# Delete the trailing empty channel written after the configured ones.
RX_time = [row for i, row in enumerate(RX_time) if i != channels]
d_TOW_at_current_symbol = [row for i, row in enumerate(
d_TOW_at_current_symbol)if i != channels]
Carrier_Doppler_hz = [row for i, row in enumerate(
Carrier_Doppler_hz) if i != channels]
Carrier_phase_hz = [row for i, row in enumerate(
Carrier_phase_hz) if i != channels]
Pseudorange_m = [
row for i, row in enumerate(Pseudorange_m) if i != channels]
PRN = [row for i, row in enumerate(PRN) if i != channels]
valid = [row for i, row in enumerate(valid) if i != channels]
observables = {
'RX_time': RX_time,
'd_TOW_at_current_symbol': d_TOW_at_current_symbol,