diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 922b70f77..7aacd6660 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -105,6 +105,15 @@ All notable changes to GNSS-SDR will be documented in this file. Their generation can be disabled by setting `Observables.enable_E6=false`. This setting is now independent of `PVT.use_e6_for_pvt`, which keeps controlling whether E6 observables are used in the PVT solution. +- Reworked the Python plotting utilities under `utils/python` (acquisition, + tracking, telemetry, observables, and PVT diagnostics). Each script now + exposes a command-line interface (run with `--help`) and can be executed from + any directory with configurable input and output locations, instead of + requiring edits to the source to change file paths. The `--file-prefix` option + takes the value of the corresponding block's `dump_filename` configuration + parameter directly, reconstructing the dump file names the same way the + receiver does. A new `utils/python/README.md` documents all the utilities and + their options. See the definitions of concepts and metrics at https://gnss-sdr.org/design-forces/ diff --git a/utils/python/README.md b/utils/python/README.md new file mode 100644 index 000000000..06636f1ca --- /dev/null +++ b/utils/python/README.md @@ -0,0 +1,305 @@ + +[comment]: # ( +SPDX-License-Identifier: GPL-3.0-or-later +) + +[comment]: # ( +SPDX-FileCopyrightText: 2026 Carles Fernandez-Prades +) + + +# GNSS-SDR Python Plot Utilities + +These scripts read GNSS-SDR dump files and save diagnostic plots. They are +intended to be run from any directory; input and output locations are +configurable through command-line options. + +Every script accepts `-h`/`--help` (printing the full option list and exiting) +and saves PNG figures to a `--fig-path` directory, creating it if needed. Pass +`--show` to also display the figures interactively. + +## Common Options + +The following options are shared by every script. The per-script tables below +repeat them with the default value that applies to each tool. + +- `-h`, `--help`: print the option list and exit. +- `-i`, `--input-path`: directory containing the dump files (default: current + directory). +- `--file-prefix`: the value of the block's `dump_filename` configuration + parameter in GNSS-SDR. It may include a directory and an extension; each + script reconstructs the actual dump filenames from it exactly the way GNSS-SDR + does — the directory is split off, the basename extension is removed, and the + block-specific suffix is appended: + + - tracking / telemetry: `.dat` + - observables / PVT: `.dat` + - acquisition: + `___ch___sat_.mat` + + Any directory in `--file-prefix` is resolved relative to `--input-path` (an + absolute `--file-prefix` path overrides it). + +- `-o`, `--fig-path`: directory where plots are saved. Defaults to a per-script + subdirectory under `./plots/`. +- `--show`: display plots interactively after saving them. By default, scripts + save the PNG files and close the figures. When `--show` is given, all of a + run's figures are kept open and displayed together with a single blocking + `plt.show()` at the end; closing the windows lets the script exit. +- `--format`: saved figure format, one of `png` (default), `pdf`, `eps`, `svg`, + `jpg`. The vector formats (`pdf`, `eps`, `svg`) are written publication-ready + with embedded fonts — PDF and EPS embed TrueType fonts (font type 42) instead + of the default Type 3 fonts that many journals reject, so the output is + suitable for camera-ready figures. Raster formats are saved at 300 DPI. This + mirrors the `--format` option of `utils/skyplot/skyplot.py`. + +## Dependencies + +The scripts use Python 3 plus: + +- `numpy`, `matplotlib` — all scripts. +- `scipy`, `h5py` — acquisition grid (`plot_acq_grid.py`). +- `pyproj` — PVT UTM coordinate conversion. +- `folium` — PVT map generation (skipped with `--no-position`). + +## Acquisition Grid — `plot_acq_grid.py` + +Reads acquisition `.mat` dumps and plots the 2D/3D acquisition grid (search +space) of the test statistic. + +Input pattern (per dump): + +```text +/___ch___sat_.mat +``` + +Output PNGs (`_2D` and `_3D`): + +```text +/___ch___sat__2D.png +/___ch___sat__3D.png +``` + +| Option | Default | Description | +| ---------------------------------------- | ------------------------ | ----------------------------------------------------------------------------------------------------------------- | +| `-i`, `--input-path` | `.` | Directory containing the acquisition `.mat` dumps. | +| `--file-prefix` | `acquisition` | `Acquisition.dump_filename` value (may include a directory/extension). | +| `-o`, `--fig-path` | `plots/acquisition` | Output directory for PNGs. | +| `--output-basename`, `--output-prefix` | `--file-prefix` basename | Base name for the saved PNGs (directory and extension stripped). | +| `--all-files` / `--single-file` | `--all-files` | Plot every matching dump, or a single dump chosen by the selectors below. | +| `--positive-only` / `--include-negative` | `--positive-only` | Plot only positive acquisitions, or positive and negative. | +| `--lite-view` / `--full-view` | `--lite-view` | Interpolated light grid, or the raw acquisition grid. | +| `--sat` | `1` in single-file mode | Satellite PRN. Filters in all-files mode; selects the dump in single-file mode. | +| `--channel` | `0` in single-file mode | Acquisition channel number. Filter / selector. | +| `--execution` | `1` in single-file mode | Acquisition dump execution index. Filter / selector. | +| `--signal-type CODE` | `1C` (GPS L1 C/A) | GNSS-SDR signal code (see table below); case-insensitive. Filters in all-files mode; selects in single-file mode. | +| `--samples-per-chip` | `3` | Samples per chip used by the light-grid interpolation. | +| `--samples-per-code` | `25000` | Samples per code used for the 2D normalization. | +| `--input-power` | `100.0` | Input power used for the 2D normalization. | +| `--show` | off | Display figures after saving. | +| `--format` | `png` | Saved figure format: `png` (default), `pdf`, `eps`, `svg`, `jpg`. Vector formats embed fonts (publication-ready). | + +Signal codes accepted by `--signal-type`, using the same nomenclature as the +GNSS-SDR configuration and dump filenames (e.g. `--signal-type=E6`): + +| Code | System | Signal | Code length (chips) | +| ---- | ------- | ------ | ------------------- | +| `1C` | GPS | L1 C/A | 1023 | +| `2S` | GPS | L2C | 10230 | +| `L5` | GPS | L5 | 10230 | +| `1B` | Galileo | E1B | 4092 | +| `5X` | Galileo | E5a | 10230 | +| `7X` | Galileo | E5b | 10230 | +| `E6` | Galileo | E6 | 5115 | +| `1G` | GLONASS | L1 | 511 | +| `2G` | GLONASS | L2 | 511 | +| `B1` | BeiDou | B1I | 2046 | +| `B3` | BeiDou | B3I | 10230 | +| `J1` | QZSS | L1 C/A | 1023 | +| `J5` | QZSS | L5 | 10230 | + +In all-files mode `--sat`, `--channel`, `--execution`, and `--signal-type` +filter the matching dumps; if a `--file-prefix` matches no files, the available +dump prefixes found in the directory are listed. + +Examples: + +```bash +python3 plot_acq_grid.py +python3 plot_acq_grid.py --input-path ./out --fig-path ./plots/acquisition +python3 plot_acq_grid.py --file-prefix hackrf_l1_acq_dump --channel 0 +python3 plot_acq_grid.py --file-prefix hackrf_l1_acq_dump --sat 26 +python3 plot_acq_grid.py --signal-type E6 +python3 plot_acq_grid.py --single-file --signal-type J1 --sat 193 --channel 0 --execution 1 +python3 plot_acq_grid.py --single-file --sat 3 --channel 0 --execution 1 --full-view +python3 plot_acq_grid.py --include-negative --output-basename my_run +``` + +## DLL/PLL VEML Tracking — `dll_pll_veml_plot_sample.py` + +Reads per-channel tracking `.dat` dumps and plots VEML tracking diagnostics +(correlator outputs, discriminators, C/N0, ...), plus optional Doppler plots. + +Input pattern: `/.dat`. + +| Option | Default | Description | +| ---------------------- | ----------------------------- | ----------------------------------------------------------------------------------------------------------------- | +| `-i`, `--input-path` | `.` | Directory containing the tracking `.dat` dumps. | +| `--file-prefix` | `track_ch` | `Tracking.dump_filename` value (may include a directory/extension). | +| `-o`, `--fig-path` | `plots/dll-pll-veml-tracking` | Output directory for PNGs. | +| `--sampling-frequency` | `3000000.0` | Signal sampling frequency in Hz (sets the RX-time axis). | +| `--channels` | `5` | Number of channels to read. | +| `--first-channel` | `0` | First channel number in the dump filenames. | +| `--plot-last-outputs` | `0` | Only plot the last N outputs; `0` plots all of them. | +| `--no-doppler` | Doppler plots on | Skip the extra Doppler-only plots. | +| `--show` | off | Display figures after saving. | +| `--format` | `png` | Saved figure format: `png` (default), `pdf`, `eps`, `svg`, `jpg`. Vector formats embed fonts (publication-ready). | + +Examples: + +```bash +python3 dll_pll_veml_plot_sample.py --input-path ./out +python3 dll_pll_veml_plot_sample.py --file-prefix tracking_ch_ --channels 8 +python3 dll_pll_veml_plot_sample.py --plot-last-outputs 5000 --no-doppler +``` + +## GPS L1 C/A KF Tracking — `gps_l1_ca_kf_plot_sample.py` + +Reads per-channel Kalman-filter tracking `.dat` dumps produced by the +`GPS_L1_CA_KF_Tracking` block and plots tracking plus Kalman-filter variables. + +Input pattern: `/.dat`. + +| Option | Default | Description | +| ---------------------- | ------------------- | ----------------------------------------------------------------------------------------------------------------- | +| `-i`, `--input-path` | `.` | Directory containing the tracking `.dat` dumps. | +| `--file-prefix` | `track_ch` | `Tracking.dump_filename` value (may include a directory/extension). | +| `-o`, `--fig-path` | `plots/kf-tracking` | Output directory for PNGs. | +| `--sampling-frequency` | `4000000.0` | Signal sampling frequency in Hz. | +| `--channels` | `5` | Number of channels to read. | +| `--first-channel` | `0` | First channel number in the dump filenames. | +| `--code-period` | `0.001` | Code period in seconds. | +| `--show` | off | Display figures after saving. | +| `--format` | `png` | Saved figure format: `png` (default), `pdf`, `eps`, `svg`, `jpg`. Vector formats embed fonts (publication-ready). | + +> The dumps must come from the `GPS_L1_CA_KF_Tracking` block; standard DLL/PLL +> tracking dumps have a different record layout and will be misread. + +Examples: + +```bash +python3 gps_l1_ca_kf_plot_sample.py --input-path ./out +python3 gps_l1_ca_kf_plot_sample.py --file-prefix track_ch --channels 8 +python3 gps_l1_ca_kf_plot_sample.py --sampling-frequency 4000000 --code-period 0.001 +``` + +## Tracking Quality Indicators — `plot_tracking_quality_indicators.py` + +Plots the carrier-lock test and C/N0 indicators of all channels on two shared +figures. + +Input pattern: `/.dat`. + +| Option | Default | Description | +| -------------------- | ------------------------ | --------------------------------------------------------------------------------------------------------------------------------------------- | +| `-i`, `--input-path` | `.` | Directory containing the tracking `.dat` dumps. | +| `--file-prefix` | `track_ch` | `Tracking.dump_filename` value (may include a directory/extension). | +| `-o`, `--fig-path` | `plots/tracking-quality` | Output directory for PNGs. | +| `--channels` | `5` | Number of channels to read. | +| `--first-channel` | `0` | First channel number in the dump filenames. | +| `--signal-type` | `1C` | GNSS-SDR signal code used as the C/N0 legend label (case-insensitive). The tracking dump stores only the PRN, so the signal is supplied here. | +| `--show` | off | Display figures after saving. | +| `--format` | `png` | Saved figure format: `png` (default), `pdf`, `eps`, `svg`, `jpg`. Vector formats embed fonts (publication-ready). | + +The C/N0 legend labels each trace as ` PRN ` (e.g. +`1C PRN 9`); the carrier-lock legend labels traces as `SV `. + +Examples: + +```bash +python3 plot_tracking_quality_indicators.py --input-path ./out +python3 plot_tracking_quality_indicators.py --file-prefix hackrf_l1_tracking_ch_ --channels 8 +python3 plot_tracking_quality_indicators.py --file-prefix hackrf_l1_tracking_ch_ --channels 8 --signal-type 1C +``` + +## Telemetry — `gps_l1_ca_telemetry_plot_sample.py` + +Reads telemetry decoder `.dat` dumps and compares the TOW timing of two +channels. + +Input pattern: `/.dat`. + +| Option | Default | Description | +| -------------------- | ----------------- | ----------------------------------------------------------------------------------------------------------------- | +| `-i`, `--input-path` | `.` | Directory containing the telemetry `.dat` dumps. | +| `--file-prefix` | `telemetry` | `TelemetryDecoder.dump_filename` value (may include a directory/extension). | +| `-o`, `--fig-path` | `plots/telemetry` | Output directory for PNGs. | +| `--channels` | `18` | Number of channels to try reading; missing channels are skipped. | +| `--first-channel` | `0` | First channel number in the dump filenames. | +| `--channel-a` | `0` | First channel number to plot. | +| `--channel-b` | `5` | Second channel number to plot. | +| `--show` | off | Display figures after saving. | +| `--format` | `png` | Saved figure format: `png` (default), `pdf`, `eps`, `svg`, `jpg`. Vector formats embed fonts (publication-ready). | + +Examples: + +```bash +python3 gps_l1_ca_telemetry_plot_sample.py --input-path ./out +python3 gps_l1_ca_telemetry_plot_sample.py --file-prefix hackrf_l1_tlm_ch_ --channel-a 0 --channel-b 5 +``` + +## Hybrid Observables — `hybrid_observables_plot_sample.py` + +Reads a single observables dump and plots pseudorange, accumulated carrier +phase, Doppler frequency, and captured PRNs across channels. + +Input pattern: `/.dat` (single file). + +| Option | Default | Description | +| -------------------- | -------------------------- | ----------------------------------------------------------------------------------------------------------------- | +| `-i`, `--input-path` | `.` | Directory containing the observables dump. | +| `--file-prefix` | `observables.dat` | `Observables.dump_filename` value (may include a directory/extension). | +| `-o`, `--fig-path` | `plots/hybrid-observables` | Output directory for PNGs. | +| `--channels` | `5` | Number of observable channels in the dump. | +| `--show` | off | Display figures after saving. | +| `--format` | `png` | Saved figure format: `png` (default), `pdf`, `eps`, `svg`, `jpg`. Vector formats embed fonts (publication-ready). | + +Examples: + +```bash +python3 hybrid_observables_plot_sample.py +python3 hybrid_observables_plot_sample.py --input-path ./out --channels 8 +python3 hybrid_observables_plot_sample.py --file-prefix hackrf_l1_observables.dat +``` + +## PVT Raw Dump — `gps_l1_ca_pvt_raw_plot_sample.py` + +Reads a single PVT raw dump and plots navigation/position diagnostics. It also +generates folium HTML maps under `/maps` unless `--no-position` is +given. + +Input pattern: `/.dat` (single file). + +| Option | Default | Description | +| ----------------------------------- | ------------------ | ----------------------------------------------------------------------------------------------------------------- | +| `-i`, `--input-path` | `.` | Directory containing the PVT dump. | +| `--file-prefix` | `PVT.dat` | `PVT.dump_filename` value (may include a directory/extension). | +| `-o`, `--fig-path` | `plots/pvt` | Output directory for PNGs and maps. | +| `--nav-sol-period` | `10.0` | Navigation solution period in milliseconds. | +| `--true-position E_UTM N_UTM U_UTM` | unset (NaN) | Reference receiver position in UTM, used for error plots. | +| `--plot-skyplot` | off | Try to generate the sky-plot panel. | +| `--no-position` | position plot on | Skip the position/map diagnostic (avoids the `folium` dependency). | +| `--one-vs-time NAME` | `X_vel`, `Tot_Vel` | `navSolutions` variable to plot versus time. Repeatable. | +| `--show` | off | Display figures after saving. | +| `--format` | `png` | Saved figure format: `png` (default), `pdf`, `eps`, `svg`, `jpg`. Vector formats embed fonts (publication-ready). | +| `--open-maps` | off | Open the generated folium map HTML files in a browser. | + +Examples: + +```bash +python3 gps_l1_ca_pvt_raw_plot_sample.py +python3 gps_l1_ca_pvt_raw_plot_sample.py --input-path ./out --file-prefix PVT.dat +python3 gps_l1_ca_pvt_raw_plot_sample.py --true-position 431234.0 4587654.0 120.0 +python3 gps_l1_ca_pvt_raw_plot_sample.py --one-vs-time X_vel --one-vs-time Y_vel +python3 gps_l1_ca_pvt_raw_plot_sample.py --no-position +``` diff --git a/utils/python/dll_pll_veml_plot_sample.py b/utils/python/dll_pll_veml_plot_sample.py old mode 100644 new mode 100755 index c7eac6a52..78a660f0d --- a/utils/python/dll_pll_veml_plot_sample.py +++ b/utils/python/dll_pll_veml_plot_sample.py @@ -1,20 +1,12 @@ +#!/usr/bin/env python3 """ dll_pll_veml_plot_sample.py - Reads GNSS-SDR Tracking dump binary file using the provided function and - plots some internal variables + Reads GNSS-SDR DLL/PLL VEML tracking dump binary files and plots internal + tracking variables. - Irene Pérez Riega, 2023. iperrie@inta.es - - Modifiable in the file: - sampling_freq - Sampling frequency [Hz] - plot_last_outputs - If 0 -> process everything / number of items processed - channels - Number of channels - first_channel - Number of the first channel - doppler_opt - = 1 -> Plot // = 0 -> No plot - path - Path to folder which contains raw files - fig_path - Path where doppler plots will be save - 'trk_dump_ch' - Fixed part of the tracking dump files names + File format: + {input_path}/{file_prefix}{channel}.dat ----------------------------------------------------------------------------- @@ -27,84 +19,169 @@ ----------------------------------------------------------------------------- """ -import os -import numpy as np +import argparse +from pathlib import Path + import matplotlib.pyplot as plt +import numpy as np + from lib.dll_pll_veml_read_tracking_dump import dll_pll_veml_read_tracking_dump +from lib.dump_filename import resolve_dump_prefix +from lib.plot_format import add_output_format_argument, apply_publication_style from lib.plotVEMLTracking import plotVEMLTracking -trackResults = [] -settings = {} -GNSS_tracking = [] -# ---------- CHANGE HERE: -sampling_freq = 3000000 -plot_last_outputs = 0 -channels = 5 -first_channel = 0 -doppler_opt = 1 -settings['numberOfChannels'] = channels +def parse_args(): + parser = argparse.ArgumentParser( + description="Plot GNSS-SDR DLL/PLL VEML tracking dumps." + ) + parser.add_argument( + "-i", + "--input-path", + type=Path, + default=Path("."), + help="Directory containing tracking .dat dumps (default: .).", + ) + parser.add_argument( + "-o", + "--fig-path", + type=Path, + default=Path("plots/dll-pll-veml-tracking"), + help="Directory where plots are saved.", + ) + parser.add_argument( + "--file-prefix", + default="track_ch", + help="GNSS-SDR Tracking.dump_filename value (default: track_ch). May " + "include a directory and extension; the matching .dat " + "files are read, resolved against --input-path.", + ) + parser.add_argument( + "--sampling-frequency", + type=float, + default=3000000.0, + help="Signal sampling frequency in Hz.", + ) + parser.add_argument( + "--channels", + type=int, + default=5, + help="Number of channels to read.", + ) + parser.add_argument( + "--first-channel", + type=int, + default=0, + help="First channel number in the dump filenames.", + ) + parser.add_argument( + "--plot-last-outputs", + type=int, + default=0, + help="Only plot the last N outputs; 0 plots all outputs.", + ) + parser.add_argument( + "--no-doppler", + dest="plot_doppler", + action="store_false", + help="Do not generate the extra Doppler-only plots.", + ) + parser.add_argument( + "--show", + action="store_true", + help="Display figures interactively after saving them.", + ) + add_output_format_argument(parser) + parser.set_defaults(plot_doppler=True) + return parser.parse_args() -path = '/home/labnav/Desktop/TEST_IRENE/tracking' -fig_path = '/home/labnav/Desktop/TEST_IRENE/PLOTS/Doppler' -for N in range(1, channels+1): - tracking_log_path = os.path.join(path, - f'trk_dump_ch{N-1+first_channel}.dat') - GNSS_tracking.append(dll_pll_veml_read_tracking_dump(tracking_log_path)) +def read_tracking_dumps(args): + directory, base = resolve_dump_prefix(args.file_prefix, args.input_path) + dumps = [] + for channel in range(args.first_channel, args.first_channel + args.channels): + tracking_log_path = directory / f"{base}{channel}.dat" + dumps.append(dll_pll_veml_read_tracking_dump(tracking_log_path)) + return dumps -# GNSS-SDR format conversion to Python GPS receiver -for N in range (1, channels+1): - if 0 < plot_last_outputs < len(GNSS_tracking[N - 1].get("code_freq_hz")): - start_sample = (len(GNSS_tracking[N-1].get("code_freq_hz")) - - plot_last_outputs) - else: - start_sample = 0 - trackResult = { - 'status': 'T', # fake track - 'codeFreq': np.copy(GNSS_tracking[N-1]["code_freq_hz"][start_sample:]), - 'carrFreq': np.copy(GNSS_tracking[N-1]["carrier_doppler_hz"][start_sample:]), - 'dllDiscr': np.copy(GNSS_tracking[N-1]["code_error"][start_sample:]), - 'dllDiscrFilt': np.copy(GNSS_tracking[N-1]["code_nco"][start_sample:]), - 'pllDiscr': np.copy(GNSS_tracking[N-1]["carr_error"][start_sample:]), - 'pllDiscrFilt': np.copy(GNSS_tracking[N-1]["carr_nco"][start_sample:]), +def main(): + args = parse_args() + args.fig_path.mkdir(parents=True, exist_ok=True) - 'I_P': np.copy(GNSS_tracking[N-1]["P"][start_sample:]), - 'Q_P': np.zeros(len(GNSS_tracking[N-1]["P"][start_sample:])), + apply_publication_style() - 'I_VE': np.copy(GNSS_tracking[N-1]["VE"][start_sample:]), - 'I_E': np.copy(GNSS_tracking[N-1]["E"][start_sample:]), - 'I_L': np.copy(GNSS_tracking[N-1]["L"][start_sample:]), - 'I_VL': np.copy(GNSS_tracking[N-1]["VL"][start_sample:]), - 'Q_VE': np.zeros(len(GNSS_tracking[N-1]["VE"][start_sample:])), - 'Q_E': np.zeros(len(GNSS_tracking[N-1]["E"][start_sample:])), - 'Q_L': np.zeros(len(GNSS_tracking[N-1]["L"][start_sample:])), - 'Q_VL': np.zeros(len(GNSS_tracking[N-1]["VL"][start_sample:])), - 'data_I': np.copy(GNSS_tracking[N-1]["prompt_I"][start_sample:]), - 'data_Q': np.copy(GNSS_tracking[N-1]["prompt_Q"][start_sample:]), - 'PRN': np.copy(GNSS_tracking[N-1]["PRN"][start_sample:]), - 'CNo': np.copy(GNSS_tracking[N-1]["CN0_SNV_dB_Hz"][start_sample:]), - 'prn_start_time_s': np.copy(GNSS_tracking[N-1]["PRN_start_sample"] - [start_sample:]) / sampling_freq + gnss_tracking = read_tracking_dumps(args) + track_results = [] + settings = { + "numberOfChannels": args.channels, + "fig_path": args.fig_path, + "show": args.show, + "output_format": args.output_format, } - trackResults.append(trackResult) - # Plot results: - plotVEMLTracking(N,trackResults,settings) + for index, tracking in enumerate(gnss_tracking, start=1): + if 0 < args.plot_last_outputs < len(tracking["code_freq_hz"]): + start_sample = len(tracking["code_freq_hz"]) - args.plot_last_outputs + else: + start_sample = 0 - # Plot Doppler according to selected in doppler_opt variable: - if doppler_opt == 1: - if not os.path.exists(fig_path): - os.makedirs(fig_path) + track_result = { + "status": "T", + "codeFreq": np.copy(tracking["code_freq_hz"][start_sample:]), + "carrFreq": np.copy(tracking["carrier_doppler_hz"][start_sample:]), + "dllDiscr": np.copy(tracking["code_error"][start_sample:]), + "dllDiscrFilt": np.copy(tracking["code_nco"][start_sample:]), + "pllDiscr": np.copy(tracking["carr_error"][start_sample:]), + "pllDiscrFilt": np.copy(tracking["carr_nco"][start_sample:]), + "I_P": np.copy(tracking["P"][start_sample:]), + "Q_P": np.zeros(len(tracking["P"][start_sample:])), + "I_VE": np.copy(tracking["VE"][start_sample:]), + "I_E": np.copy(tracking["E"][start_sample:]), + "I_L": np.copy(tracking["L"][start_sample:]), + "I_VL": np.copy(tracking["VL"][start_sample:]), + "Q_VE": np.zeros(len(tracking["VE"][start_sample:])), + "Q_E": np.zeros(len(tracking["E"][start_sample:])), + "Q_L": np.zeros(len(tracking["L"][start_sample:])), + "Q_VL": np.zeros(len(tracking["VL"][start_sample:])), + "data_I": np.copy(tracking["prompt_I"][start_sample:]), + "data_Q": np.copy(tracking["prompt_Q"][start_sample:]), + "PRN": np.copy(tracking["PRN"][start_sample:]), + "CNo": np.copy(tracking["CN0_SNV_dB_Hz"][start_sample:]), + "prn_start_time_s": ( + np.copy(tracking["PRN_start_sample"][start_sample:]) + / args.sampling_frequency + ), + } + track_results.append(track_result) - plt.figure() - plt.plot(trackResults[N - 1]['prn_start_time_s'], - [x/1000 for x in GNSS_tracking[N - 1]['carrier_doppler_hz'] - [start_sample:]]) - plt.xlabel('Time(s)') - plt.ylabel('Doppler(KHz)') - plt.title('Doppler frequency channel ' + str(N)) + plotVEMLTracking(index, track_results, settings) - plt.savefig(os.path.join(fig_path, f'Doppler_freq_ch_{N}.png')) + if args.plot_doppler: + channel = args.first_channel + index - 1 + plt.figure() + plt.plot( + track_result["prn_start_time_s"], + [x / 1000 for x in tracking["carrier_doppler_hz"][start_sample:]], + ) + plt.xlabel("Time(s)") + plt.ylabel("Doppler(KHz)") + plt.title(f"Doppler frequency channel {channel}") + plt.savefig( + args.fig_path + / f"Doppler_freq_ch_{channel}.{args.output_format}" + ) + if not args.show: + plt.close() + + # Show all saved figures with a single plt.show() to avoid the repeated + # show()/close() cycle that can crash interactive backends on macOS. + if args.show: plt.show() + + +if __name__ == "__main__": + try: + main() + except OSError as exc: + raise SystemExit(f"Error: {exc}") diff --git a/utils/python/gps_l1_ca_kf_plot_sample.py b/utils/python/gps_l1_ca_kf_plot_sample.py old mode 100644 new mode 100755 index 9f310f1da..e755e4e55 --- a/utils/python/gps_l1_ca_kf_plot_sample.py +++ b/utils/python/gps_l1_ca_kf_plot_sample.py @@ -1,21 +1,12 @@ +#!/usr/bin/env python3 """ gps_l1_ca_kf_plot_sample.py - Reads a GPS L1 C/A Kalman-filter tracking dump binary file - (GPS_L1_CA_KF_Tracking) and plots some internal tracking and Kalman-filter - variables. + Reads GPS L1 C/A Kalman-filter tracking dump binary files + (GPS_L1_CA_KF_Tracking) and plots tracking and Kalman-filter variables. - Irene Pérez Riega, 2023. iperrie@inta.es - Minhaj Uddin Ahmad, 2026. mahmad12@crimson.ua.edu - - Modifiable in the file: - samplingFreq - Sampling frequency [Hz] - channels - Number of channels - first_channel - Number of the first channel - code_period - Code period [s] - path - Path to folder which contains the tracking dump files - fig_path - Path where the plots will be saved - file_prefix - Fixed part of the tracking dump file names + File format: + {input_path}/{file_prefix}{channel}.dat ----------------------------------------------------------------------------- @@ -28,80 +19,153 @@ ----------------------------------------------------------------------------- """ -import os +import argparse +from pathlib import Path + import numpy as np + +from lib.dump_filename import resolve_dump_prefix from lib.gps_l1_ca_kf_read_tracking_dump import gps_l1_ca_kf_read_tracking_dump -from lib.plotTracking import plotTracking +from lib.plot_format import add_output_format_argument, apply_publication_style from lib.plotKalman import plotKalman +from lib.plotTracking import plotTracking -GNSS_tracking = [] -trackResults = [] -kalmanResults = [] -# ---------- CHANGE HERE: -samplingFreq = 4000000 # must match SignalSource.sampling_frequency in the .conf -channels = 5 -first_channel = 0 -code_period = 0.001 -path = '../../../out/' -fig_path = '../../../PLOTS/KF_Tracking' -file_prefix = 'track_ch' +def parse_args(): + parser = argparse.ArgumentParser( + description="Plot GPS L1 C/A Kalman-filter tracking dumps." + ) + parser.add_argument( + "-i", + "--input-path", + type=Path, + default=Path("."), + help="Directory containing tracking .dat dumps (default: .).", + ) + parser.add_argument( + "-o", + "--fig-path", + type=Path, + default=Path("plots/kf-tracking"), + help="Directory where plots are saved.", + ) + parser.add_argument( + "--file-prefix", + default="track_ch", + help="GNSS-SDR Tracking.dump_filename value (default: track_ch). May " + "include a directory and extension; the matching .dat " + "files are read, resolved against --input-path.", + ) + parser.add_argument( + "--sampling-frequency", + type=float, + default=4000000.0, + help="Signal sampling frequency in Hz.", + ) + parser.add_argument( + "--channels", + type=int, + default=5, + help="Number of channels to read.", + ) + parser.add_argument( + "--first-channel", + type=int, + default=0, + help="First channel number in the dump filenames.", + ) + parser.add_argument( + "--code-period", + type=float, + default=0.001, + help="Code period in seconds.", + ) + parser.add_argument( + "--show", + action="store_true", + help="Display figures interactively after saving them.", + ) + add_output_format_argument(parser) + return parser.parse_args() -for N in range(1, channels + 1): - tracking_log_path = os.path.join(path, - f'{file_prefix}{N-1+first_channel}.dat') - GNSS_tracking.append(gps_l1_ca_kf_read_tracking_dump(tracking_log_path)) -for N in range(1, channels + 1): - trackResult = { - 'status': 'T', # fake track - 'codeFreq': np.copy(GNSS_tracking[N-1]["code_freq_hz"]), - 'carrFreq': np.copy(GNSS_tracking[N-1]["carrier_doppler_hz"]), - 'carrFreqRate': - np.copy(GNSS_tracking[N-1]["carrier_doppler_rate_hz2"]), - 'dllDiscr': np.copy(GNSS_tracking[N-1]["code_error"]), - 'dllDiscrFilt': np.copy(GNSS_tracking[N-1]["code_nco"]), - 'pllDiscr': np.copy(GNSS_tracking[N-1]["carr_error"]), - 'pllDiscrFilt': np.copy(GNSS_tracking[N-1]["carr_nco"]), +def read_tracking_dumps(args): + directory, base = resolve_dump_prefix(args.file_prefix, args.input_path) + dumps = [] + for channel in range(args.first_channel, args.first_channel + args.channels): + tracking_log_path = directory / f"{base}{channel}.dat" + dumps.append(gps_l1_ca_kf_read_tracking_dump(tracking_log_path)) + return dumps - 'I_P': np.copy(GNSS_tracking[N-1]["prompt_I"]), - 'Q_P': np.copy(GNSS_tracking[N-1]["prompt_Q"]), - 'I_E': np.copy(GNSS_tracking[N-1]["E"]), - 'I_L': np.copy(GNSS_tracking[N-1]["L"]), - 'Q_E': np.zeros(len(GNSS_tracking[N-1]["E"])), - 'Q_L': np.zeros(len(GNSS_tracking[N-1]["L"])), - 'PRN': np.copy(GNSS_tracking[N-1]["PRN"]), - 'CNo': np.copy(GNSS_tracking[N-1]["CN0_SNV_dB_Hz"]), - 'prn_start_time_s': - np.copy(GNSS_tracking[N-1]["PRN_start_sample"]) / samplingFreq - } +def main(): + args = parse_args() + args.fig_path.mkdir(parents=True, exist_ok=True) - # Kalman-filter internals. The KF dump stores the carrier discriminator - # (carr_error -> innovation) and the filter states, but no measurement - # noise covariance, so 'r_noise_cov' is intentionally left out (plotKalman - # skips that panel when it is absent). - kalmanResult = { - 'PRN': np.copy(GNSS_tracking[N-1]["PRN"]), - 'innovation': np.copy(GNSS_tracking[N-1]["carr_error"]), - 'state1': np.copy(GNSS_tracking[N-1]["carr_nco"]), - 'state2': np.copy(GNSS_tracking[N-1]["carrier_doppler_hz"]), - 'state3': np.copy(GNSS_tracking[N-1]["carrier_doppler_rate_hz2"]), - 'CNo': np.copy(GNSS_tracking[N-1]["CN0_SNV_dB_Hz"]) - } + apply_publication_style() - trackResults.append(trackResult) - kalmanResults.append(kalmanResult) + gnss_tracking = read_tracking_dumps(args) + track_results = [] + kalman_results = [] - settings = { - 'numberOfChannels': channels, - 'msToProcess': len(GNSS_tracking[N-1]['E']), - 'codePeriod': code_period, - 'timeStartInSeconds': 0, - 'fig_path': fig_path, - 'show': False # set True to display the figures interactively - } + for index, tracking in enumerate(gnss_tracking, start=1): + track_result = { + "status": "T", + "codeFreq": np.copy(tracking["code_freq_hz"]), + "carrFreq": np.copy(tracking["carrier_doppler_hz"]), + "carrFreqRate": np.copy(tracking["carrier_doppler_rate_hz2"]), + "dllDiscr": np.copy(tracking["code_error"]), + "dllDiscrFilt": np.copy(tracking["code_nco"]), + "pllDiscr": np.copy(tracking["carr_error"]), + "pllDiscrFilt": np.copy(tracking["carr_nco"]), + "I_P": np.copy(tracking["prompt_I"]), + "Q_P": np.copy(tracking["prompt_Q"]), + "I_E": np.copy(tracking["E"]), + "I_L": np.copy(tracking["L"]), + "Q_E": np.zeros(len(tracking["E"])), + "Q_L": np.zeros(len(tracking["L"])), + "PRN": np.copy(tracking["PRN"]), + "CNo": np.copy(tracking["CN0_SNV_dB_Hz"]), + "prn_start_time_s": ( + np.copy(tracking["PRN_start_sample"]) / args.sampling_frequency + ), + } - # Create and save graphics as PNG - plotTracking(N, trackResults, settings) - plotKalman(N, kalmanResults, settings) + kalman_result = { + "PRN": np.copy(tracking["PRN"]), + "innovation": np.copy(tracking["carr_error"]), + "state1": np.copy(tracking["carr_nco"]), + "state2": np.copy(tracking["carrier_doppler_hz"]), + "state3": np.copy(tracking["carrier_doppler_rate_hz2"]), + "CNo": np.copy(tracking["CN0_SNV_dB_Hz"]), + } + + track_results.append(track_result) + kalman_results.append(kalman_result) + + settings = { + "numberOfChannels": args.channels, + "msToProcess": len(tracking["E"]), + "codePeriod": args.code_period, + "timeStartInSeconds": 0, + "fig_path": args.fig_path, + "show": args.show, + "output_format": args.output_format, + } + + plotTracking(index, track_results, settings) + plotKalman(index, kalman_results, settings) + + # Show all saved figures with a single plt.show() to avoid the repeated + # show()/close() cycle that can crash interactive backends on macOS. + if args.show: + import matplotlib.pyplot as plt + + plt.show() + + +if __name__ == "__main__": + try: + main() + except OSError as exc: + raise SystemExit(f"Error: {exc}") diff --git a/utils/python/gps_l1_ca_pvt_raw_plot_sample.py b/utils/python/gps_l1_ca_pvt_raw_plot_sample.py old mode 100644 new mode 100755 index fb922b5e5..0158f3ba7 --- a/utils/python/gps_l1_ca_pvt_raw_plot_sample.py +++ b/utils/python/gps_l1_ca_pvt_raw_plot_sample.py @@ -1,23 +1,8 @@ +#!/usr/bin/env python3 """ gps_l1_ca_pvt_raw_plot_sample.py - Reads GNSS-SDR PVT raw dump binary file using the provided function and plots - some internal variables - - Irene Pérez Riega, 2023. iperrie@inta.es - - Modifiable in the file: - sampling_freq - Sampling frequency [Hz] - channels - Number of channels to check if they exist - path - Path to folder which contains raw file - pvt_raw_log_path - Completed path to PVT raw data file - nav_sol_period - Measurement period [ms] - plot_skyplot - = 1 -> Sky Plot (TO DO) // = 0 -> No Sky Plot - true_position - In settings, If not known enter all NaN's and mean - position will be used as a reference in UTM - coordinate system - plot_position - Optional function at the end - plot_oneVStime - Optional function at the end, select variable to plot + Reads a GNSS-SDR PVT raw dump and plots navigation/position diagnostics. ----------------------------------------------------------------------------- @@ -30,72 +15,189 @@ ----------------------------------------------------------------------------- """ -import utm +import argparse +from pathlib import Path + import numpy as np -from lib.gps_l1_ca_read_pvt_dump import gps_l1_ca_read_pvt_dump -from lib.plotNavigation import plotNavigation -import pyproj -from lib.plotPosition import plot_position, plot_oneVStime -settings = {} -utm_e = [] -utm_n = [] -E_UTM = [] -N_UTM = [] -utm_zone = [] - -# ---------- CHANGE HERE: -samplingFreq = 64e6 / 16 -channels = 5 -path = '/home/labnav/Desktop/TEST_IRENE/' -pvt_raw_log_path = path + 'PVT.dat' -nav_sol_period = 10 -plot_skyplot = 0 - -settings['true_position'] = {'E_UTM':np.nan,'N_UTM':np.nan,'U_UTM':np.nan} -settings['navSolPeriod'] = nav_sol_period - -navSolutions = gps_l1_ca_read_pvt_dump(pvt_raw_log_path) -X, Y, Z = navSolutions['X'], navSolutions['Y'], navSolutions['Z'] - -utm_coords = [] - -for i in range(len(navSolutions['longitude'])): - utm_coords.append(utm.from_latlon(navSolutions['latitude'][i], - navSolutions['longitude'][i])) +from lib.dump_filename import resolve_dump_prefix +from lib.plot_format import add_output_format_argument, apply_publication_style - -for i in range(len(utm_coords)): - utm_e.append(utm_coords[i][0]) - utm_n.append(utm_coords[i][1]) - utm_zone.append(utm_coords[i][2]) - -# Transform from Lat Long degrees to UTM coordinate system -# It's supposed utm_zone and letter will not change during tracking -input_projection = pyproj.CRS.from_string("+proj=longlat " - "+datum=WGS84 +no_defs") - -utm_e = [] -utm_n = [] -for i in range(len(navSolutions['longitude'])): - output_projection = pyproj.CRS (f"+proj=utm +zone={utm_zone[i]} " - f"+datum=WGS84 +units=m +no_defs") - transformer = pyproj.Transformer.from_crs(input_projection, - output_projection) - utm_e, utm_n = transformer.transform(navSolutions['longitude'][i], - navSolutions['latitude'][i]) - E_UTM.append(utm_e) - N_UTM.append(utm_n) +def parse_args(): + parser = argparse.ArgumentParser( + description="Plot GNSS-SDR GPS L1 C/A PVT raw dump data." + ) + parser.add_argument( + "-i", + "--input-path", + type=Path, + default=Path("."), + help="Directory containing the PVT dump (default: .).", + ) + parser.add_argument( + "--file-prefix", + default="PVT.dat", + help="GNSS-SDR PVT.dump_filename value (default: PVT.dat). May include " + "a directory and extension; the matching .dat file is read, " + "resolved against --input-path.", + ) + parser.add_argument( + "-o", + "--fig-path", + type=Path, + default=Path("plots/pvt"), + help="Directory where plots are saved.", + ) + parser.add_argument( + "--nav-sol-period", + type=float, + default=10.0, + help="Navigation solution period in milliseconds.", + ) + parser.add_argument( + "--true-position", + type=float, + nargs=3, + metavar=("E_UTM", "N_UTM", "U_UTM"), + default=[np.nan, np.nan, np.nan], + help="Reference receiver position in UTM coordinates.", + ) + parser.add_argument( + "--plot-skyplot", + action="store_true", + help="Placeholder: the sky plot panel is not available from a PVT " + "dump (no azimuth/elevation data). Use utils/skyplot/skyplot.py " + "instead.", + ) + parser.add_argument( + "--no-position", + dest="plot_position", + action="store_false", + help="Skip the position/map diagnostic plot.", + ) + parser.add_argument( + "--one-vs-time", + action="append", + default=None, + help="navSolutions variable to plot versus time. Can be repeated.", + ) + parser.add_argument( + "--show", + action="store_true", + help="Display figures interactively after saving them.", + ) + parser.add_argument( + "--open-maps", + action="store_true", + help="Open generated folium map HTML files in a browser.", + ) + add_output_format_argument(parser) + parser.set_defaults(plot_position=True) + return parser.parse_args() -navSolutions['E_UTM'] = E_UTM -navSolutions['N_UTM'] = N_UTM -navSolutions['U_UTM'] = navSolutions['Z'] +def add_utm_coordinates(nav_solutions): + try: + from pyproj import CRS, Transformer + from pyproj.aoi import AreaOfInterest + from pyproj.database import query_utm_crs_info + except ModuleNotFoundError as exc: + raise SystemExit( + "gps_l1_ca_pvt_raw_plot_sample.py requires pyproj for PVT " + "coordinate conversion." + ) from exc -plotNavigation(navSolutions,settings,plot_skyplot) + input_projection = CRS.from_epsg(4326) + e_utm = [] + n_utm = [] + for longitude, latitude in zip( + nav_solutions["longitude"], nav_solutions["latitude"] + ): + utm_crs_info = query_utm_crs_info( + datum_name="WGS 84", + area_of_interest=AreaOfInterest( + west_lon_degree=longitude, + south_lat_degree=latitude, + east_lon_degree=longitude, + north_lat_degree=latitude, + ), + ) + if utm_crs_info: + output_projection = CRS.from_epsg(utm_crs_info[0].code) + else: + zone = int((longitude + 180) // 6) + 1 + epsg = (32600 if latitude >= 0 else 32700) + zone + output_projection = CRS.from_epsg(epsg) -# OPTIONAL: Other plots -> -plot_position(navSolutions) -plot_oneVStime(navSolutions, 'X_vel') -plot_oneVStime(navSolutions, 'Tot_Vel') \ No newline at end of file + transformer = Transformer.from_crs( + input_projection, output_projection, always_xy=True) + east, north = transformer.transform(longitude, latitude) + e_utm.append(east) + n_utm.append(north) + + nav_solutions["E_UTM"] = e_utm + nav_solutions["N_UTM"] = n_utm + nav_solutions["U_UTM"] = nav_solutions["Z"] + + +def main(): + args = parse_args() + args.fig_path.mkdir(parents=True, exist_ok=True) + apply_publication_style() + + from lib.gps_l1_ca_read_pvt_dump import gps_l1_ca_read_pvt_dump + from lib.plotNavigation import plotNavigation + from lib.plotPosition import plot_oneVStime, plot_position + + settings = { + "true_position": { + "E_UTM": args.true_position[0], + "N_UTM": args.true_position[1], + "U_UTM": args.true_position[2], + }, + "navSolPeriod": args.nav_sol_period, + "fig_path": args.fig_path, + "show": args.show, + "output_format": args.output_format, + } + + directory, base = resolve_dump_prefix(args.file_prefix, args.input_path) + pvt_file = directory / f"{base}.dat" + nav_solutions = gps_l1_ca_read_pvt_dump(pvt_file) + add_utm_coordinates(nav_solutions) + + plotNavigation(nav_solutions, settings, int(args.plot_skyplot)) + + if args.plot_position: + plot_position( + nav_solutions, + fig_path=args.fig_path, + show=args.show, + open_maps=args.open_maps, + output_format=args.output_format, + ) + + one_vs_time = args.one_vs_time or ["X_vel", "Tot_Vel"] + for variable_name in one_vs_time: + plot_oneVStime( + nav_solutions, + variable_name, + fig_path=args.fig_path, + show=args.show, + output_format=args.output_format, + ) + + # Show all saved figures with a single plt.show() to avoid the repeated + # show()/close() cycle that can crash interactive backends on macOS. + if args.show: + import matplotlib.pyplot as plt + + plt.show() + + +if __name__ == "__main__": + try: + main() + except OSError as exc: + raise SystemExit(f"Error: {exc}") diff --git a/utils/python/gps_l1_ca_telemetry_plot_sample.py b/utils/python/gps_l1_ca_telemetry_plot_sample.py old mode 100644 new mode 100755 index 51e456d66..26aaefe7f --- a/utils/python/gps_l1_ca_telemetry_plot_sample.py +++ b/utils/python/gps_l1_ca_telemetry_plot_sample.py @@ -1,18 +1,12 @@ +#!/usr/bin/env python3 """ gps_l1_ca_telemetry_plot_sample.py - Reads GNSS-SDR Tracking dump binary file using the provided function and - plots some internal variables + Reads GNSS-SDR GPS L1 C/A telemetry dump binary files and compares telemetry + timing fields for two channels. - Irene Pérez Riega, 2023. iperrie@inta.es - - Modifiable in the file: - sampling_freq - Sampling frequency [Hz] - channels - Number of channels to check if they exist - doppler_opt - = 1 -> Plot // = 0 -> No plot - path - Path to folder which contains raw file - fig_path - Path where plots will be save - chn_num_a / b - Channel which will be plotted + File format: + {input_path}/{file_prefix}{channel}.dat ----------------------------------------------------------------------------- @@ -25,75 +19,157 @@ ----------------------------------------------------------------------------- """ -import os +import argparse +from pathlib import Path + import matplotlib.pyplot as plt + +from lib.dump_filename import resolve_dump_prefix from lib.gps_l1_ca_read_telemetry_dump import gps_l1_ca_read_telemetry_dump +from lib.plot_format import add_output_format_argument, apply_publication_style -GNSS_telemetry = [] -# ---------- CHANGE HERE: -sampling_freq = 2000000 -channels = list(range(18)) -path = '/home/labnav/Desktop/TEST_IRENE/' -fig_path = '/home/labnav/Desktop/TEST_IRENE/PLOTS/Telemetry' +def parse_args(): + parser = argparse.ArgumentParser( + description="Plot GNSS-SDR GPS L1 C/A telemetry dumps." + ) + parser.add_argument( + "-i", + "--input-path", + type=Path, + default=Path("."), + help="Directory containing telemetry .dat dumps (default: .).", + ) + parser.add_argument( + "-o", + "--fig-path", + type=Path, + default=Path("plots/telemetry"), + help="Directory where plots are saved.", + ) + parser.add_argument( + "--file-prefix", + default="telemetry", + help="GNSS-SDR TelemetryDecoder.dump_filename value (default: " + "telemetry). May include a directory and extension; the matching " + ".dat files are read, resolved against --input-path.", + ) + parser.add_argument( + "--channels", + type=int, + default=18, + help="Number of channels to try reading.", + ) + parser.add_argument( + "--first-channel", + type=int, + default=0, + help="First channel number in the dump filenames.", + ) + parser.add_argument( + "--channel-a", + type=int, + default=0, + help="First channel number to plot.", + ) + parser.add_argument( + "--channel-b", + type=int, + default=5, + help="Second channel number to plot.", + ) + parser.add_argument( + "--show", + action="store_true", + help="Display figures interactively after saving them.", + ) + add_output_format_argument(parser) + return parser.parse_args() -if not os.path.exists(fig_path): - os.makedirs(fig_path) -i = 0 -for N in channels: +def read_telemetry_dumps(args): + directory, base = resolve_dump_prefix(args.file_prefix, args.input_path) + telemetry = {} + for channel in range(args.first_channel, args.first_channel + args.channels): + telemetry_log_path = directory / f"{base}{channel}.dat" + try: + telemetry[channel] = gps_l1_ca_read_telemetry_dump(telemetry_log_path) + except OSError: + continue + return telemetry + + +def require_channel(telemetry, channel): + if channel not in telemetry: + available = ", ".join(str(ch) for ch in sorted(telemetry)) or "none" + raise ValueError(f"Channel {channel} was not read. Available: {available}") + return telemetry[channel] + + +def main(): + args = parse_args() + args.fig_path.mkdir(parents=True, exist_ok=True) + + apply_publication_style() + + telemetry = read_telemetry_dumps(args) + data_a = require_channel(telemetry, args.channel_a) + data_b = require_channel(telemetry, args.channel_b) + + plt.figure() + plt.gcf().canvas.manager.set_window_title( + f"Telem_Current_Symbol_TOW_{args.channel_a}_{args.channel_b}.png" + ) + plt.plot( + data_a["tracking_sample_counter"], + [x / 1000 for x in data_a["tow_current_symbol_ms"]], + "b", + ) + plt.plot( + data_b["tracking_sample_counter"], + [x / 1000 for x in data_b["tow_current_symbol_ms"]], + "r", + ) + plt.grid(True) + plt.xlabel("TRK Sampling Counter") + plt.ylabel("Current Symbol TOW [s]") + plt.legend([f"CHN-{args.channel_a}", f"CHN-{args.channel_b}"]) + plt.tight_layout() + plt.savefig( + args.fig_path + / f"Telem_Current_Symbol_TOW_{args.channel_a}_{args.channel_b}" + f".{args.output_format}" + ) + if not args.show: + plt.close() + + plt.figure() + plt.gcf().canvas.manager.set_window_title( + f"Telem_TRK_Sampling_Counter_{args.channel_a}_{args.channel_b}.png" + ) + plt.plot(data_a["tracking_sample_counter"], data_a["tow"], "b") + plt.plot(data_b["tracking_sample_counter"], data_b["tow"], "r") + plt.grid(True) + plt.xlabel("TRK Sampling Counter") + plt.ylabel("Decoded Nav TOW") + plt.legend([f"CHN-{args.channel_a}", f"CHN-{args.channel_b}"]) + plt.tight_layout() + plt.savefig( + args.fig_path + / f"Telem_TRK_Sampling_Counter_{args.channel_a}_{args.channel_b}" + f".{args.output_format}" + ) + if not args.show: + plt.close() + + # Show all saved figures with a single plt.show() to avoid the repeated + # show()/close() cycle that can crash interactive backends on macOS. + if args.show: + plt.show() + + +if __name__ == "__main__": try: - telemetry_log_path = os.path.join(path, f'telemetry{N}.dat') - telemetry_data = gps_l1_ca_read_telemetry_dump(telemetry_log_path) - GNSS_telemetry.append(telemetry_data) - i += 1 - except: - pass - -# ---------- CHANGE HERE: -chn_num_a = 0 -chn_num_b = 5 - -# Plotting values -if chn_num_a in range(i) and chn_num_b in range(i): - - # First Plot: - plt.figure() - plt.gcf().canvas.manager.set_window_title(f'Telem_Current_Simbol_TOW_' - f'{chn_num_a}_{chn_num_b}.png') - - plt.plot(GNSS_telemetry[chn_num_a]['tracking_sample_counter'], - [x / 1000 for x in GNSS_telemetry[chn_num_a] - ['tow_current_symbol_ms']], 'b') - plt.plot(GNSS_telemetry[chn_num_b]['tracking_sample_counter'], - GNSS_telemetry[chn_num_b]['tow_current_symbol_ms'], 'r') - - plt.grid(True) - plt.xlabel('TRK Sampling Counter') - plt.ylabel('Current Symbol TOW') - plt.legend([f'CHN-{chn_num_a-1}', f'CHN-{chn_num_b-1}']) - plt.tight_layout() - - plt.savefig(os.path.join(fig_path, f'Telem_Current_Simbol_TOW_{chn_num_a}' - f'_{chn_num_b}.png')) - plt.show() - - # Second Plot: - plt.figure() - plt.gcf().canvas.manager.set_window_title(f'Telem_TRK_Sampling_Counter_' - f'{chn_num_a}_{chn_num_b}.png') - - plt.plot(GNSS_telemetry[chn_num_a]['tracking_sample_counter'], - GNSS_telemetry[chn_num_a]['tow'], 'b') - plt.plot(GNSS_telemetry[chn_num_b]['tracking_sample_counter'], - GNSS_telemetry[chn_num_b]['tow'], 'r') - - plt.grid(True) - plt.xlabel('TRK Sampling Counter') - plt.ylabel('Decoded Nav TOW') - plt.legend([f'CHN-{chn_num_a-1}', f'CHN-{chn_num_b-1}']) - plt.tight_layout() - - plt.savefig(os.path.join(fig_path, f'Telem_TRK_Sampling_Counter_' - f'{chn_num_a}_{chn_num_b}.png')) - plt.show() + main() + except (OSError, ValueError) as exc: + raise SystemExit(f"Error: {exc}") diff --git a/utils/python/hybrid_observables_plot_sample.py b/utils/python/hybrid_observables_plot_sample.py old mode 100644 new mode 100755 index c4d932529..cb2a16f59 --- a/utils/python/hybrid_observables_plot_sample.py +++ b/utils/python/hybrid_observables_plot_sample.py @@ -1,17 +1,9 @@ +#!/usr/bin/env python3 """ hybrid_observables_plot_sample.py - Reads GNSS-SDR observables raw dump binary file using the provided function - and plots some internal variables - - Irene Pérez Riega, 2023. iperrie@inta.es - - Modifiable in the file: - sampling_freq - Sampling frequency [Hz] - channels - Number of channels to check if they exist - path - Path to folder which contains raw file - fig_path - Path where plots will be save - observables_log_path - Completed path to observables raw data file + Reads a GNSS-SDR hybrid observables raw dump and plots pseudorange, carrier + phase, Doppler, and PRN values. ----------------------------------------------------------------------------- @@ -24,116 +16,177 @@ ----------------------------------------------------------------------------- """ -import numpy as np +import argparse +from pathlib import Path + import matplotlib.pyplot as plt +import numpy as np + +from lib.dump_filename import resolve_dump_prefix +from lib.plot_format import add_output_format_argument, apply_publication_style from lib.read_hybrid_observables_dump import read_hybrid_observables_dump -import os -observables = {} -double_size_bytes = 8 -bytes_shift = 0 -# ---------- CHANGE HERE: -samplingFreq = 2000000 -channels = 5 -path = '/home/labnav/Desktop/TEST_IRENE/' -fig_path = '/home/labnav/Desktop/TEST_IRENE/PLOTS/HybridObservables/' -observables_log_path = path + 'observables.dat' +def parse_args(): + parser = argparse.ArgumentParser( + description="Plot GNSS-SDR hybrid observables dump data." + ) + parser.add_argument( + "-i", + "--input-path", + type=Path, + default=Path("."), + help="Directory containing the observables dump (default: .).", + ) + parser.add_argument( + "--file-prefix", + default="observables.dat", + help="GNSS-SDR Observables.dump_filename value (default: " + "observables.dat). May include a directory and extension; the " + "matching .dat file is read, resolved against --input-path.", + ) + parser.add_argument( + "-o", + "--fig-path", + type=Path, + default=Path("plots/hybrid-observables"), + help="Directory where plots are saved.", + ) + parser.add_argument( + "--channels", + type=int, + default=5, + help="Number of observable channels in the dump.", + ) + parser.add_argument( + "--show", + action="store_true", + help="Display figures interactively after saving them.", + ) + add_output_format_argument(parser) + return parser.parse_args() -if not os.path.exists(fig_path): - os.makedirs(fig_path) -GNSS_observables = read_hybrid_observables_dump(channels,observables_log_path) +def first_valid_observable(gnss_observables, channels): + min_tow_idx = None + obs_idx = 0 + for channel in range(channels): + valid_indices = np.where(np.array(gnss_observables["valid"][channel]) > 0)[0] + if len(valid_indices) == 0: + continue + idx = valid_indices[0] + if min_tow_idx is None or idx < min_tow_idx: + min_tow_idx = idx + obs_idx = channel -# Plot data -# --- optional: plot since it detect the first satellite connected -min_tow_idx = 1 -obs_idx = 1 + if min_tow_idx is None: + raise ValueError("No valid observables found in the dump.") + return min_tow_idx, obs_idx -for n in range(0, channels): - idx = np.where(np.array(GNSS_observables['valid'][n])>0)[0][0] - # Find the index from the first satellite connected - if n == 0: - min_tow_idx = idx - if min_tow_idx > idx: - min_tow_idx = idx - obs_idx = n +def save_figure(fig_path, name, show, output_format): + plt.tight_layout() + plt.savefig(fig_path / f"{name}.{output_format}") + # Close unless it will be shown; main() triggers a single plt.show() at the + # end. Avoids repeated show()/close() cycles, which can crash interactive + # matplotlib backends (e.g. macOS) on window close. + if not show: + plt.close() -# Plot observables from that index -# First plot -plt.figure() -plt.title('Pseudorange') -for i in range(channels): - plt.scatter(GNSS_observables['RX_time'][i][min_tow_idx:], - GNSS_observables['Pseudorange_m'][i][min_tow_idx:],s=1, - label=f'Channel {i}') -plt.xlim(GNSS_observables['RX_time'][obs_idx][min_tow_idx]-100, - GNSS_observables['RX_time'][obs_idx][-1]+100) -plt.grid(True) -plt.xlabel('TOW [s]') -plt.ylabel('Pseudorange [m]') -plt.legend() -plt.gcf().canvas.manager.set_window_title('Pseudorange.png') -plt.tight_layout() -plt.savefig(os.path.join(fig_path, 'Pseudorange.png')) -plt.show() -# Second plot -plt.figure() -plt.title('Carrier Phase') -for i in range(channels): - plt.scatter(GNSS_observables['RX_time'][i][min_tow_idx:], - GNSS_observables['Carrier_phase_hz'][i][min_tow_idx:],s=1, - label=f'Channel {i}') -plt.xlim(GNSS_observables['RX_time'][obs_idx][min_tow_idx]-100, - GNSS_observables['RX_time'][obs_idx][-1]+100) -plt.xlabel('TOW [s]') -plt.ylabel('Accumulated Carrier Phase [cycles]') -plt.grid(True) -plt.legend() -plt.gcf().canvas.manager.set_window_title('AccumulatedCarrierPhase.png') -plt.tight_layout() -plt.savefig(os.path.join(fig_path, 'AccumulatedCarrierPhase.png')) -plt.show() +def main(): + args = parse_args() + args.fig_path.mkdir(parents=True, exist_ok=True) + apply_publication_style() -# Third plot -plt.figure() -plt.title('Doppler Effect') -for i in range(channels): - plt.scatter(GNSS_observables['RX_time'][i][min_tow_idx:], - GNSS_observables['Carrier_Doppler_hz'][i][min_tow_idx:],s=1, - label=f'Channel {i}') -plt.xlim(GNSS_observables['RX_time'][obs_idx][min_tow_idx]-100, - GNSS_observables['RX_time'][obs_idx][-1]+100) -plt.xlabel('TOW [s]') -plt.ylabel('Doppler Frequency [Hz]') -plt.grid(True) -plt.legend() -plt.gcf().canvas.manager.set_window_title('DopplerFrequency.png') -plt.tight_layout() -plt.savefig(os.path.join(fig_path, 'DopplerFrequency.png')) -plt.show() + 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 -# Fourth plot -plt.figure() -plt.title('GNSS Channels captured') -for i in range(channels): - lab = 0 - a = 0 - while lab == 0: - lab = int(GNSS_observables["PRN"][i][min_tow_idx+a]) - a += 1 - plt.scatter(GNSS_observables['RX_time'][i][min_tow_idx:], - GNSS_observables['PRN'][i][min_tow_idx:], s=1, - label=f'PRN {i} = {lab}') -plt.xlim(GNSS_observables['RX_time'][obs_idx][min_tow_idx]-100, - GNSS_observables['RX_time'][obs_idx][-1]+100) -plt.xlabel('TOW [s]') -plt.ylabel('PRN') -plt.grid(True) -plt.legend() -plt.gcf().canvas.manager.set_window_title('PRNs.png') -plt.tight_layout() -plt.savefig(os.path.join(fig_path, 'PRNs.png')) -plt.show() \ No newline at end of file + 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) + + 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) + + 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) + + 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 + plt.scatter( + gnss_observables["RX_time"][channel][min_tow_idx:], + gnss_observables["PRN"][channel][min_tow_idx:], + s=1, + label=f"PRN {channel} = {label}", + ) + plt.xlim(time_start, time_end) + plt.xlabel("TOW [s]") + plt.ylabel("PRN") + plt.grid(True) + plt.legend() + plt.gcf().canvas.manager.set_window_title("PRNs.png") + save_figure(args.fig_path, "PRNs", args.show, args.output_format) + + # Show all saved figures with a single plt.show() to avoid the repeated + # show()/close() cycle that can crash interactive backends on macOS. + if args.show: + plt.show() + + +if __name__ == "__main__": + try: + main() + except (OSError, ValueError) as exc: + raise SystemExit(f"Error: {exc}") diff --git a/utils/python/lib/dump_filename.py b/utils/python/lib/dump_filename.py new file mode 100644 index 000000000..3464cbc6d --- /dev/null +++ b/utils/python/lib/dump_filename.py @@ -0,0 +1,46 @@ +""" + dump_filename.py + + Helper to interpret a GNSS-SDR ``dump_filename`` configuration value the same + way the C++ signal-processing blocks do, so the plotting utilities can locate + the dump files from the exact string the user wrote in the configuration. + + ----------------------------------------------------------------------------- + + GNSS-SDR is a Global Navigation Satellite System software-defined receiver. + This file is part of GNSS-SDR. + + Copyright (C) 2026 (see AUTHORS file for a list of contributors) + SPDX-License-Identifier: GPL-3.0-or-later + + ----------------------------------------------------------------------------- +""" + +from pathlib import Path + + +def resolve_dump_prefix(file_prefix, input_path="."): + """Split a GNSS-SDR ``dump_filename`` into ``(directory, basename)``. + + Mirrors the C++ dump-file logic shared by the tracking, telemetry, + observables, PVT, and acquisition blocks: the directory part (everything + before the last ``/``) is separated from the basename, and any extension + on the basename is removed. The directory is resolved relative to + ``input_path`` (an absolute ``file_prefix`` path overrides it). The + block-specific suffixes (channel number, ``.dat``, ``.mat``, ...) are + appended by the caller. + + Examples: + resolve_dump_prefix("track_ch") -> (Path("."), "track_ch") + resolve_dump_prefix("./observables.dat") -> (Path("."), "observables") + resolve_dump_prefix("run1/acq", "data") -> (Path("data/run1"), "acq") + """ + prefix_path = Path(file_prefix) + name = prefix_path.name + # Remove a trailing extension, ignoring a possible leading dot (matches the + # ``substr(1).find_last_of('.')`` logic used by the GNSS-SDR blocks). + dot = name[1:].rfind(".") + if dot != -1: + name = name[: dot + 1] + directory = Path(input_path) / prefix_path.parent + return directory, name diff --git a/utils/python/lib/plotKalman.py b/utils/python/lib/plotKalman.py index f2dffe245..a2a482d85 100644 --- a/utils/python/lib/plotKalman.py +++ b/utils/python/lib/plotKalman.py @@ -32,9 +32,8 @@ import os def plotKalman(channelNr, trackResults, settings): - # ---------- CHANGE HERE (or pass 'fig_path' in settings): - fig_path = settings.get('fig_path', - '../../../PLOTS/PlotKalman') + fig_path = settings.get('fig_path', 'plots/kalman') + output_format = settings.get('output_format', 'png') if not os.path.exists(fig_path): os.makedirs(fig_path) @@ -140,6 +139,9 @@ def plotKalman(channelNr, trackResults, settings): plt.savefig(os.path.join(fig_path, f'kalman_ch{channelNr}_PRN_' f'{trackResults[channelNr - 1]["PRN"][-1]}' - f'.png')) - if settings.get('show', True): - plt.show() + f'.{output_format}')) + # Close unless it will be shown; the caller triggers a single + # plt.show() at the end. Avoids repeated show()/close() cycles, which + # can crash interactive matplotlib backends (e.g. macOS) on close. + if not settings.get('show', True): + plt.close() diff --git a/utils/python/lib/plotNavigation.py b/utils/python/lib/plotNavigation.py index 661f38b29..c7df99d94 100644 --- a/utils/python/lib/plotNavigation.py +++ b/utils/python/lib/plotNavigation.py @@ -39,8 +39,8 @@ import os def plotNavigation(navSolutions, settings, plot_skyplot=0): - # ---------- CHANGE HERE: - fig_path = '/home/labnav/Desktop/TEST_IRENE/PLOTS/PlotNavigation' + fig_path = settings.get('fig_path', 'plots/pvt') + output_format = settings.get('output_format', 'png') if not os.path.exists(fig_path): os.makedirs(fig_path) @@ -67,9 +67,9 @@ def plotNavigation(navSolutions, settings, plot_skyplot=0): else: # Compute the mean error for static receiver ref_coord = { - 'E_UTM': settings.truePosition['E_UTM'], - 'N_UTM': settings.truePosition['N_UTM'], - 'U_UTM': settings.truePosition['U_UTM'] + 'E_UTM': settings['true_position']['E_UTM'], + 'N_UTM': settings['true_position']['N_UTM'], + 'U_UTM': settings['true_position']['U_UTM'] } mean_position = { @@ -96,9 +96,15 @@ def plotNavigation(navSolutions, settings, plot_skyplot=0): ax3 = plt.subplot(4, 2, (6, 8), projection='3d') # (ax1) Coordinate differences in UTM system from reference point - ax1.plot(np.vstack([navSolutions['E_UTM'] - ref_coord['E_UTM'], - navSolutions['N_UTM'] - ref_coord['N_UTM'], - navSolutions['U_UTM'] - ref_coord['U_UTM']]).T) + # Coerce to arrays so the offset works whether the reference is a + # float (from --true-position) or a NumPy mean of the solutions. + e_utm = np.asarray(navSolutions['E_UTM'], dtype=float) + n_utm = np.asarray(navSolutions['N_UTM'], dtype=float) + u_utm = np.asarray(navSolutions['U_UTM'], dtype=float) + d_east = e_utm - ref_coord['E_UTM'] + d_north = n_utm - ref_coord['N_UTM'] + d_up = u_utm - ref_coord['U_UTM'] + ax1.plot(np.vstack([d_east, d_north, d_up]).T) ax1.set_title('Coordinates variations in UTM system', fontweight='bold') ax1.legend(['E_UTM', 'N_UTM', 'U_UTM']) ax1.set_xlabel(f"Measurement period: {settings['navSolPeriod']} ms") @@ -106,19 +112,18 @@ def plotNavigation(navSolutions, settings, plot_skyplot=0): ax1.grid(True) ax1.axis('tight') - # (ax2) Satellite sky plot - if plot_skyplot: #todo posicion de los satelites - skyPlot(ax2, navSolutions['channel']['az'], - navSolutions['channel']['el'], - navSolutions['channel']['PRN'][:, 0]) - ax2.set_title(f'Sky plot (mean PDOP: ' - f'{np.nanmean(navSolutions["DOP"][1, :]):.1f})', - fontweight='bold') + # (ax2) Reserved for a satellite sky plot, which cannot be drawn from a + # PVT dump (it has no per-satellite azimuth/elevation). Hide the empty + # panel; use utils/skyplot/skyplot.py to plot a skyplot from a RINEX + # navigation file. + ax2.axis('off') + if plot_skyplot: + print("Warning: the sky plot panel is not available from a PVT " + "dump (no azimuth/elevation data). Use " + "utils/skyplot/skyplot.py instead.") # (ax3) Position plot in UTM system - ax3.scatter(navSolutions['E_UTM'] - ref_coord['E_UTM'], - navSolutions['N_UTM'] - ref_coord['N_UTM'], - navSolutions['U_UTM'] - ref_coord['U_UTM'], marker='+') + ax3.scatter(d_east, d_north, d_up, marker='+') ax3.scatter([0], [0], [0], color='r', marker='+', linewidth=1.5) ax3.view_init(0, 90) ax3.set_box_aspect([1, 1, 1]) @@ -130,5 +135,9 @@ def plotNavigation(navSolutions, settings, plot_skyplot=0): ax3.set_zlabel('Upping (m)') plt.tight_layout() - plt.savefig(os.path.join(fig_path, 'measures_UTM.png')) - plt.show() + plt.savefig(os.path.join(fig_path, f'measures_UTM.{output_format}')) + # Close unless it will be shown; the caller triggers a single + # plt.show() at the end. Avoids repeated show()/close() cycles, which + # can crash interactive matplotlib backends (e.g. macOS) on close. + if not settings.get('show', True): + plt.close() diff --git a/utils/python/lib/plotPosition.py b/utils/python/lib/plotPosition.py index 3284b16cb..a117a3da8 100644 --- a/utils/python/lib/plotPosition.py +++ b/utils/python/lib/plotPosition.py @@ -42,21 +42,27 @@ import math import os.path import webbrowser +from pathlib import Path import numpy as np import matplotlib.pyplot as plt -import folium -def plot_position(navSolutions): +def plot_position(navSolutions, fig_path='plots/pvt', show=True, + open_maps=False, output_format='png'): + try: + import folium + except ModuleNotFoundError as exc: + raise SystemExit( + "plot_position requires folium. Re-run with --no-position to skip " + "map generation." + ) from exc - # ---------- CHANGE HERE: - fig_path = '/home/labnav/Desktop/TEST_IRENE/PLOTS/PlotPosition/' - fig_path_maps = fig_path + 'maps/' + fig_path = os.fspath(fig_path) + fig_path_maps = os.path.join(fig_path, 'maps') filename_map = 'mapPlotPosition.html' - filename_map_t = 'mapTerrainPotPosition.html' if not os.path.exists(fig_path_maps): - os.mkdir(fig_path_maps) + os.makedirs(fig_path_maps) # Statics Positions: m_lat = sum(navSolutions['latitude']) / len(navSolutions['latitude']) @@ -95,8 +101,11 @@ def plot_position(navSolutions): icon=folium.Icon(color='red')).add_to(m) """ - m.save(fig_path_maps + filename_map) - webbrowser.open(fig_path_maps + filename_map) + map_path = os.path.join(fig_path_maps, filename_map) + m.save(map_path) + if open_maps: + # webbrowser needs a URL; pass an absolute file:// URI, not a bare path. + webbrowser.open(Path(map_path).resolve().as_uri()) # Optional: with terrain -> """ @@ -117,7 +126,7 @@ def plot_position(navSolutions): plt.figure(figsize=(1920 / 120, 1080 / 120)) plt.clf() - plt.suptitle(f'Plot file PVT process data results') + plt.suptitle('Plot file PVT process data results') # Latitude and Longitude plt.subplot(1, 2, 1) @@ -149,16 +158,20 @@ def plot_position(navSolutions): ax.set_title('Positions x-y-z') plt.tight_layout() - plt.savefig(os.path.join(fig_path, f'PVT_ProcessDataResults.png')) - plt.show() + plt.savefig(os.path.join(fig_path, f'PVT_ProcessDataResults.{output_format}')) + # Close unless it will be shown; the caller triggers a single plt.show() + # at the end. Avoids repeated show()/close() cycles, which can crash + # interactive matplotlib backends (e.g. macOS) on window close. + if not show: + plt.close() -def plot_oneVStime(navSolutions, name): +def plot_oneVStime(navSolutions, name, fig_path='plots/pvt', show=True, + output_format='png'): - # ---------- CHANGE HERE: - fig_path = '/home/labnav/Desktop/TEST_IRENE/PLOTS/PlotPosition/' + fig_path = os.fspath(fig_path) if not os.path.exists(fig_path): - os.mkdir(fig_path) + os.makedirs(fig_path) time = [] for i in range(len(navSolutions['TransmitTime'])): @@ -175,8 +188,12 @@ def plot_oneVStime(navSolutions, name): plt.ticklabel_format(style='plain', axis='both', useOffset=False) plt.tight_layout() - plt.savefig(os.path.join(fig_path, f'{name}VSTime.png')) - plt.show() + plt.savefig(os.path.join(fig_path, f'{name}VSTime.{output_format}')) + # Close unless it will be shown; the caller triggers a single plt.show() + # at the end. Avoids repeated show()/close() cycles, which can crash + # interactive matplotlib backends (e.g. macOS) on window close. + if not show: + plt.close() def calcularCEFP(percentil, navSolutions, m_lat, m_long): diff --git a/utils/python/lib/plotTracking.py b/utils/python/lib/plotTracking.py index 01f949b28..29966c4a8 100644 --- a/utils/python/lib/plotTracking.py +++ b/utils/python/lib/plotTracking.py @@ -32,9 +32,8 @@ import matplotlib.pyplot as plt def plotTracking(channelNr, trackResults, settings): - # ---------- CHANGE HERE (or pass 'fig_path' in settings): - fig_path = settings.get('fig_path', - '../../../PLOTS/PlotTracking') + fig_path = settings.get('fig_path', 'plots/tracking') + output_format = settings.get('output_format', 'png') if not os.path.exists(fig_path): os.makedirs(fig_path) @@ -187,6 +186,9 @@ def plotTracking(channelNr, trackResults, settings): plt.savefig(os.path.join(fig_path, f'trk_dump_ch{channelNr}_PRN_' f'{trackResults[channelNr - 1]["PRN"][-1]}' - f'.png')) - if settings.get('show', True): - plt.show() + f'.{output_format}')) + # Close unless it will be shown; the caller triggers a single + # plt.show() at the end. Avoids repeated show()/close() cycles, which + # can crash interactive matplotlib backends (e.g. macOS) on close. + if not settings.get('show', True): + plt.close() diff --git a/utils/python/lib/plotVEMLTracking.py b/utils/python/lib/plotVEMLTracking.py index 2af50e03a..01a4aba28 100644 --- a/utils/python/lib/plotVEMLTracking.py +++ b/utils/python/lib/plotVEMLTracking.py @@ -33,8 +33,8 @@ import os def plotVEMLTracking(channelNr, trackResults, settings): - # ---------- CHANGE HERE: - fig_path = '/home/labnav/Desktop/TEST_IRENE/PLOTS/VEMLTracking' + fig_path = settings.get('fig_path', 'plots/dll-pll-veml-tracking') + output_format = settings.get('output_format', 'png') if not os.path.exists(fig_path): os.makedirs(fig_path) @@ -43,7 +43,7 @@ def plotVEMLTracking(channelNr, trackResults, settings): plt.figure(figsize=(1920 / 120, 1080 / 120)) plt.clf() - plt.gcf().canvas.set_window_title( + plt.gcf().canvas.manager.set_window_title( f'Channel {channelNr} (PRN ' f'{trackResults[channelNr-1]["PRN"][0]}) results') plt.subplots_adjust(left=0.1, right=0.9, top=0.9, bottom=0.1, @@ -167,5 +167,9 @@ def plotVEMLTracking(channelNr, trackResults, settings): plt.savefig(os.path.join(fig_path, f'Ch{channelNr}_PRN' f'{trackResults[channelNr-1]["PRN"][0]}' - f'_results')) - plt.show() + f'_results.{output_format}')) + # Close the figure unless it will be shown; the caller triggers a single + # plt.show() at the end. Avoids repeated show()/close() cycles, which can + # crash interactive matplotlib backends (e.g. macOS) on window close. + if not settings.get('show', True): + plt.close() diff --git a/utils/python/lib/plot_format.py b/utils/python/lib/plot_format.py new file mode 100644 index 000000000..008d4001e --- /dev/null +++ b/utils/python/lib/plot_format.py @@ -0,0 +1,55 @@ +""" + plot_format.py + + Shared output-format helpers for the GNSS-SDR Python plotting utilities. + + Provides a common --format command-line argument (same set of formats as + utils/skyplot/skyplot.py) and a matplotlib style that makes the vector + outputs (PDF/EPS) publication-ready with embedded fonts. + + ----------------------------------------------------------------------------- + + GNSS-SDR is a Global Navigation Satellite System software-defined receiver. + This file is part of GNSS-SDR. + + Copyright (C) 2026 (see AUTHORS file for a list of contributors) + SPDX-License-Identifier: GPL-3.0-or-later + + ----------------------------------------------------------------------------- +""" + +# Output formats offered by every plotting utility, matching the --format +# option of utils/skyplot/skyplot.py. PNG is the default here. +OUTPUT_FORMATS = ["png", "pdf", "eps", "svg", "jpg"] +DEFAULT_OUTPUT_FORMAT = "png" + + +def add_output_format_argument(parser): + """Add the shared --format option to an argparse parser.""" + parser.add_argument( + "--format", + dest="output_format", + type=str.lower, + default=DEFAULT_OUTPUT_FORMAT, + choices=OUTPUT_FORMATS, + help="Saved figure format (default: png). Vector formats (pdf, eps, " + "svg) are written publication-ready with embedded fonts.", + ) + + +def apply_publication_style(): + """Configure matplotlib so saved figures are publication-ready. + + The key settings embed the fonts in PDF and EPS output (font type 42 = + TrueType) instead of leaving them as the default Type 3 fonts, which many + journals reject. SVG text is kept editable, and raster output uses a high + DPI. The active font family is left untouched, so plots keep their usual + look; whatever font is in use gets embedded. + """ + import matplotlib.pyplot as plt + + plt.rcParams["pdf.fonttype"] = 42 # embed TrueType fonts in PDF + plt.rcParams["ps.fonttype"] = 42 # embed TrueType fonts in EPS/PS + plt.rcParams["svg.fonttype"] = "none" # keep SVG text editable + plt.rcParams["savefig.dpi"] = 300 # high-resolution png/jpg + plt.rcParams["savefig.bbox"] = "tight" # trim surrounding whitespace diff --git a/utils/python/lib/read_hybrid_observables_dump.py b/utils/python/lib/read_hybrid_observables_dump.py index dc7cf09d2..9b37df886 100644 --- a/utils/python/lib/read_hybrid_observables_dump.py +++ b/utils/python/lib/read_hybrid_observables_dump.py @@ -82,20 +82,22 @@ def read_hybrid_observables_dump(channels, filename): bytes_shift += double_size_bytes f.seek(bytes_shift, 0) - except: + except struct.error: + # Reached a partial record at end of file: stop reading. break - # Delete last Channel: - RX_time = [row for i, row in enumerate(RX_time) if i != 5] + # 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 != 5] + d_TOW_at_current_symbol)if i != channels] Carrier_Doppler_hz = [row for i, row in enumerate( - Carrier_Doppler_hz) if i != 5] + Carrier_Doppler_hz) if i != channels] Carrier_phase_hz = [row for i, row in enumerate( - Carrier_phase_hz) if i != 5] - Pseudorange_m = [row for i, row in enumerate(Pseudorange_m) if i != 5] - PRN = [row for i, row in enumerate(PRN) if i != 5] - valid = [row for i, row in enumerate(valid) if i != 5] + 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, diff --git a/utils/python/plot_acq_grid.py b/utils/python/plot_acq_grid.py old mode 100644 new mode 100755 index a1ede156f..165cf9513 --- a/utils/python/plot_acq_grid.py +++ b/utils/python/plot_acq_grid.py @@ -1,28 +1,15 @@ +#!/usr/bin/env python3 """ plot_acq_grid.py - Reads GNSS-SDR Acquisition dump .mat file using the provided function and - plots acquisition grid of acquisition statistic of PRN sat + Reads GNSS-SDR Acquisition dump .mat files and plots acquisition grids. - Irene Pérez Riega, 2023. iperrie@inta.es + Irene Perez Riega, 2023. iperrie@inta.es Minhaj Uddin Ahmad, 2026. mahmad12@crimson.ua.edu - Modifiable in the file: - path - Path to folder which contains raw file - fig_path - Path where plots will be save - plot_all_files - Plot all the files in a folder (True/False) - plot_positive_acqs - Only plot the positive acquisition cases - ---- - file_prefix - Fixed part in files names. In our case: acquisition - sat - Satellite. In our case: 1 - channel - Channel. In our case: 1 - execution - In our case: 0 - signal_type - In our case: 1 - ---- - lite_view - True for light grid representation - File format: - {path}/{file}_ch_{system}_{signal}_ch_{channel}_{execution}_sat_{sat}.mat + {path}/{file_prefix}_{system}_{signal}_ch_{channel}_{execution}_sat_{sat}.mat + ----------------------------------------------------------------------------- GNSS-SDR is a Global Navigation Satellite System software-defined receiver. @@ -34,211 +21,444 @@ ----------------------------------------------------------------------------- """ -import os +import argparse from pathlib import Path -import numpy as np -import matplotlib.pyplot as plt -from scipy.interpolate import CubicSpline -import h5py -# ---------- CHANGE HERE: -path = '../../../out/' # path to acquisition .mat dumps -fig_path = '../../../PLOTS/Acquisition/' # path to store figures -plot_all_files = True -plot_positive_acqs = True -file_prefix = "acquisition" +from lib.dump_filename import resolve_dump_prefix +from lib.plot_format import add_output_format_argument, apply_publication_style -signal_types = { - 1: ('G', '1C', 1023), # GPS L1 - 2: ('G', '2S', 10230), # GPS L2M - 3: ('G', 'L5', 10230), # GPS L5 - 4: ('E', '1B', 4092), # Galileo E1B - 5: ('E', '5X', 10230), # Galileo E5 - 6: ('R', '1G', 511), # Glonass 1G - 7: ('R', '2G', 511), # Glonass 2G - 8: ('C', 'B1', 2048), # Beidou B1 - 9: ('C', 'B3', 10230), # Beidou B3 - 10: ('C', '5C', 10230) # Beidou B2a +# Heavy numerical and plotting libraries (h5py, numpy, matplotlib, scipy) are +# imported lazily inside the functions that need them. This keeps --help and +# argument-error invocations fast and free of matplotlib font-cache warnings. + + +# GNSS-SDR signal nomenclature: each two-character signal code maps to its +# system character (as written in the dump filename) and primary code length +# in chips. The signal code matches the value used in GNSS-SDR configuration +# (e.g. Galileo.implementation signals) and in the dump filename. +SIGNAL_TYPES = { + "1C": ("G", 1023), # GPS L1 C/A + "2S": ("G", 10230), # GPS L2C + "L5": ("G", 10230), # GPS L5 + "1B": ("E", 4092), # Galileo E1B + "5X": ("E", 10230), # Galileo E5a + "7X": ("E", 10230), # Galileo E5b + "E6": ("E", 5115), # Galileo E6 + "1G": ("R", 511), # GLONASS L1 + "2G": ("R", 511), # GLONASS L2 + "B1": ("C", 2046), # BeiDou B1I + "B3": ("C", 10230), # BeiDou B3I + "J1": ("J", 1023), # QZSS L1 C/A + "J5": ("J", 10230), # QZSS L5 } -# (system, signal) -> n_chips, derived from the table above -N_CHIPS = {(sys_, sig): nc for sys_, sig, nc in signal_types.values()} +N_CHIPS = {(system, code): nc for code, (system, nc) in SIGNAL_TYPES.items()} +DEFAULT_SAT = 1 +DEFAULT_CHANNEL = 0 +DEFAULT_EXECUTION = 1 +DEFAULT_SIGNAL_TYPE = "1C" -if not os.path.exists(fig_path): - os.makedirs(fig_path) -if not plot_all_files: - # ---------- CHANGE HERE: - sat = 1 - channel = 0 - execution = 1 - signal_type = 1 - lite_view = True +def parse_args(): + parser = argparse.ArgumentParser( + description="Plot GNSS-SDR acquisition dump grids." + ) + parser.add_argument( + "-i", + "--input-path", + type=Path, + default=Path("."), + help="Directory containing acquisition .mat dumps (default: .).", + ) + parser.add_argument( + "-o", + "--fig-path", + type=Path, + default=Path("plots/acquisition"), + help="Directory where plots are saved (default: ./plots/acquisition).", + ) + parser.add_argument( + "--file-prefix", + default="acquisition", + help="GNSS-SDR Acquisition.dump_filename value (default: " + "acquisition). May include a directory and extension; matching " + "___ch___sat_.mat files are read, " + "resolved against --input-path.", + ) + parser.add_argument( + "--output-basename", + "--output-prefix", + dest="output_basename", + default=None, + help="Base name used for saved plot files (default: the --file-prefix " + "basename, without directory or extension).", + ) - # If lite_view -> sets the number of samples per chip in the graphical - # representation - n_samples_per_chip = 3 - d_samples_per_code = 25000 + mode = parser.add_mutually_exclusive_group() + mode.add_argument( + "--all-files", + dest="plot_all_files", + action="store_true", + help="Plot all matching dump files (default).", + ) + mode.add_argument( + "--single-file", + dest="plot_all_files", + action="store_false", + help="Plot one dump selected by --sat, --channel, --execution, and " + "--signal-type.", + ) + parser.set_defaults(plot_all_files=True) - system, signal, n_chips = signal_types.get(signal_type) + acq_filter = parser.add_mutually_exclusive_group() + acq_filter.add_argument( + "--positive-only", + dest="plot_positive_acqs", + action="store_true", + help="Only plot positive acquisitions (default).", + ) + acq_filter.add_argument( + "--include-negative", + dest="plot_positive_acqs", + action="store_false", + help="Plot positive and negative acquisition dumps.", + ) + parser.set_defaults(plot_positive_acqs=True) - # Load data - filename = (f'{path}{file_prefix}_ch_{system}_{signal}_ch_{channel}_{execution}' - f'_sat_{sat}.mat') - img_name_root = (f'{fig_path}{file_prefix}_ch_{system}_{signal}_ch_{channel}_' - f'{execution}_sat_{sat}') + view = parser.add_mutually_exclusive_group() + view.add_argument( + "--lite-view", + dest="lite_view", + action="store_true", + help="Use interpolated light grid representation (default).", + ) + view.add_argument( + "--full-view", + dest="lite_view", + action="store_false", + help="Plot the raw acquisition grid.", + ) + parser.set_defaults(lite_view=True) - with h5py.File(filename, 'r') as data: - acq_grid = data['acq_grid'][:] + parser.add_argument( + "--sat", + type=int, + default=None, + help="Satellite PRN. Filters all-files mode; defaults to 1 in " + "single-file mode.", + ) + parser.add_argument( + "--channel", + type=int, + default=None, + help="Acquisition channel number. Filters all-files mode; defaults " + "to 0 in single-file mode.", + ) + parser.add_argument( + "--execution", + type=int, + default=None, + help="Acquisition dump execution. Filters all-files mode; defaults " + "to 1 in single-file mode.", + ) + parser.add_argument( + "--signal-type", + type=str.upper, + choices=sorted(SIGNAL_TYPES), + default=None, + metavar="CODE", + help="GNSS-SDR signal code (e.g. 1C, 5X, 7X, E6, J1). Filters " + "all-files mode; defaults to 1C (GPS L1 C/A) in single-file mode. " + "Choices: " + ", ".join(sorted(SIGNAL_TYPES)) + ".", + ) + parser.add_argument( + "--samples-per-chip", + type=int, + default=3, + help="Samples per chip used by the light grid interpolation.", + ) + parser.add_argument( + "--samples-per-code", + type=int, + default=25000, + help="Samples per code used for the 2D normalization.", + ) + parser.add_argument( + "--input-power", + type=float, + default=100.0, + help="Input power used for the 2D normalization.", + ) + parser.add_argument( + "--show", + action="store_true", + help="Display figures interactively after saving them.", + ) + add_output_format_argument(parser) + + args = parser.parse_args() + if args.output_basename is None: + # Default the saved-plot basename to the dump basename, dropping any + # directory and extension carried by --file-prefix. + _, args.output_basename = resolve_dump_prefix(args.file_prefix) + return args + + +def positive_acq_value(data): + for key in ("positive_acq", "d_positive_acq"): + value = data.get(key) + if value is not None: + return value[0] + return None + + +def parse_dump_name(file_path): + parts = file_path.stem.split("_") + if len(parts) < 8: + return None + try: + if parts[-5] != "ch" or parts[-2] != "sat": + return None + return { + "file_prefix": "_".join(parts[:-7]), + "system": parts[-7], + "signal": parts[-6], + "channel": int(parts[-4]), + "execution": int(parts[-3]), + "sat": int(parts[-1]), + } + except ValueError: + return None + + +def matching_dump_prefixes(input_path): + prefixes = set() + for file_path in input_path.glob("*.mat"): + metadata = parse_dump_name(file_path) + if metadata is not None: + prefixes.add(metadata["file_prefix"]) + return sorted(prefixes) + + +def read_acquisition_dump(file_path, n_chips, positive_only): + import h5py + import numpy as np + + with h5py.File(file_path, "r") as data: + positive_acq = positive_acq_value(data) + if positive_only: + if positive_acq is None: + print(f"Skipping {file_path.name}: no positive_acq variable") + return None + if positive_acq != 1: + return None + + acq_grid = data["acq_grid"][:] n_fft, n_dop_bins = acq_grid.shape d_max, f_max = np.unravel_index(np.argmax(acq_grid), acq_grid.shape) - doppler_step = data['doppler_step'][0] - doppler_max = data['doppler_max'][0] + doppler_step = data["doppler_step"][0] + doppler_max = data["doppler_max"][0] freq = np.arange(n_dop_bins) * doppler_step - doppler_max delay = np.arange(n_fft) / n_fft * n_chips - # Plot data - # --- Acquisition grid (3D) + return acq_grid, n_fft, d_max, f_max, doppler_step, doppler_max, freq, delay + + +def plot_dump(file_path, fig_path, image_name_root, n_chips, args): + import matplotlib.pyplot as plt + import numpy as np + from scipy.interpolate import CubicSpline + + dump = read_acquisition_dump(file_path, n_chips, args.plot_positive_acqs) + if dump is None: + return False + + ( + acq_grid, + n_fft, + d_max, + f_max, + doppler_step, + doppler_max, + freq, + delay, + ) = dump + fig = plt.figure() - plt.gcf().canvas.manager.set_window_title(filename) - if not lite_view: - ax = fig.add_subplot(111, projection='3d') - X, Y = np.meshgrid(freq, delay) - ax.plot_surface(X, Y, acq_grid, cmap='viridis') + plt.gcf().canvas.manager.set_window_title(str(file_path)) + if not args.lite_view: + ax = fig.add_subplot(111, projection="3d") + x_axis, y_axis = np.meshgrid(freq, delay) + ax.plot_surface(x_axis, y_axis, acq_grid, cmap="viridis") ax.set_ylim([min(delay), max(delay)]) else: - delay_interp = (np.arange(n_samples_per_chip * n_chips) - / n_samples_per_chip) + delay_interp = ( + np.arange(args.samples_per_chip * n_chips) / args.samples_per_chip + ) spline = CubicSpline(delay, acq_grid) grid_interp = spline(delay_interp) - ax = fig.add_subplot(111, projection='3d') - X, Y = np.meshgrid(freq, delay_interp) - ax.plot_surface(X, Y, grid_interp, cmap='inferno') + ax = fig.add_subplot(111, projection="3d") + x_axis, y_axis = np.meshgrid(freq, delay_interp) + ax.plot_surface(x_axis, y_axis, grid_interp, cmap="inferno") ax.set_ylim([min(delay_interp), max(delay_interp)]) - ax.set_xlabel('Doppler shift (Hz)') + ax.set_xlabel("Doppler shift (Hz)") ax.set_xlim([min(freq), max(freq)]) - ax.set_ylabel('Code delay (chips)') - ax.set_zlabel('Test Statistics') - + ax.set_ylabel("Code delay (chips)") + ax.set_zlabel("Test Statistics") plt.tight_layout() - plt.savefig(img_name_root + '_sample_3D.png') - plt.show() - - # --- Acquisition grid (2D) - input_power = 100 # Change Test statistics in Doppler wipe-off plot + plt.savefig(fig_path / f"{image_name_root}_3D.{args.output_format}") + if not args.show: + plt.close(fig) fig2, axes = plt.subplots(2, 1, figsize=(8, 6)) - plt.gcf().canvas.manager.set_window_title(filename) + plt.gcf().canvas.manager.set_window_title(str(file_path)) axes[0].plot(freq, acq_grid[d_max, :]) axes[0].set_xlim([min(freq), max(freq)]) - axes[0].set_xlabel('Doppler shift (Hz)') - axes[0].set_ylabel('Test statistics') - axes[0].set_title(f'Fixed code delay to {(d_max - 1) / n_fft * n_chips} ' - f'chips') + axes[0].set_xlabel("Doppler shift (Hz)") + axes[0].set_ylabel("Test statistics") + axes[0].set_title( + f"Fixed code delay to {(d_max - 1) / n_fft * n_chips} chips" + ) - normalization = (d_samples_per_code**4) * input_power + normalization = (args.samples_per_code ** 4) * args.input_power axes[1].plot(delay, acq_grid[:, f_max] / normalization) axes[1].set_xlim([min(delay), max(delay)]) - axes[1].set_xlabel('Code delay (chips)') - axes[1].set_ylabel('Test statistics') - axes[1].set_title(f'Doppler wipe-off = ' - f'{str((f_max-1) * doppler_step - doppler_max)} Hz') + axes[1].set_xlabel("Code delay (chips)") + axes[1].set_ylabel("Test statistics") + axes[1].set_title( + f"Doppler wipe-off = {(f_max - 1) * doppler_step - doppler_max} Hz" + ) plt.tight_layout() - plt.savefig(img_name_root + '_sample_2D.png') + plt.savefig(fig_path / f"{image_name_root}_2D.{args.output_format}") + if not args.show: + plt.close(fig2) + return True + + +def single_file_path(args): + sat = args.sat if args.sat is not None else DEFAULT_SAT + channel = args.channel if args.channel is not None else DEFAULT_CHANNEL + execution = ( + args.execution if args.execution is not None else DEFAULT_EXECUTION + ) + signal_type = ( + args.signal_type + if args.signal_type is not None + else DEFAULT_SIGNAL_TYPE + ) + signal = signal_type + system, n_chips = SIGNAL_TYPES[signal_type] + directory, base = resolve_dump_prefix(args.file_prefix, args.input_path) + filename = ( + f"{base}_{system}_{signal}_ch_{channel}_" + f"{execution}_sat_{sat}.mat" + ) + image_name_root = image_name_from_metadata( + args.output_basename, + { + "system": system, + "signal": signal, + "channel": channel, + "execution": execution, + "sat": sat, + }, + ) + return directory / filename, image_name_root, n_chips + + +def image_name_from_metadata(output_basename, metadata): + prefix = f"{output_basename}_" if output_basename else "" + return ( + f"{prefix}{metadata['system']}_{metadata['signal']}_ch_" + f"{metadata['channel']}_{metadata['execution']}_sat_{metadata['sat']}" + ) + + +def metadata_matches_filters(metadata, args): + if args.sat is not None and metadata["sat"] != args.sat: + return False + if args.channel is not None and metadata["channel"] != args.channel: + return False + if args.execution is not None and metadata["execution"] != args.execution: + return False + if args.signal_type is not None: + system, _ = SIGNAL_TYPES[args.signal_type] + if metadata["system"] != system or metadata["signal"] != args.signal_type: + return False + return True + + +def show_figures(args): + # Display every saved figure with a single plt.show() call. Showing once + # (instead of once per figure) keeps the GUI event loop to a single entry + # and avoids closing figures afterwards, which prevents the segfault some + # macOS matplotlib backends hit on repeated show()/close() cycles. The + # figures stay open until the user closes all windows, then the process + # exits and the interpreter tears matplotlib down cleanly. + if not args.show: + return + import matplotlib.pyplot as plt + plt.show() -else: - # ---------- CHANGE HERE: - lite_view = True - # If lite_view -> sets the number of samples per chip in the graphical - # representation - n_samples_per_chip = 3 - d_samples_per_code = 25000 - file_paths = Path(path).glob(f'{file_prefix}*.mat') +def main(): + args = parse_args() + args.fig_path.mkdir(parents=True, exist_ok=True) + apply_publication_style() + + if not args.plot_all_files: + file_path, image_name_root, n_chips = single_file_path(args) + if not file_path.exists(): + raise FileNotFoundError(f"dump file not found: {file_path}") + plot_dump(file_path, args.fig_path, image_name_root, n_chips, args) + show_figures(args) + return + + plotted = 0 + directory, base = resolve_dump_prefix(args.file_prefix, args.input_path) + file_paths = sorted(directory.glob(f"{base}*.mat")) for file_path in file_paths: - # {prefix}_{system}_{signal}_ch_{channel}_{execution}_sat_{PRN} - # -7 -6 -5 -4 -3 -2 -1 - parts = file_path.stem.split('_') - - system = parts[-7] - signal = parts[-6] - channel = int(parts[-4]) - execution = int(parts[-3]) - sat = int(parts[-1]) - - filename = file_path.name - - n_chips = N_CHIPS.get((system, signal)) - if n_chips is None: - print(f"Unknown system/signal {system}/{signal}, skipping: " - f"{filename}") + metadata = parse_dump_name(file_path) + if metadata is None: + print(f"Skipping {file_path.name}: unexpected filename format") continue - with h5py.File(file_path, 'r') as data: - if plot_positive_acqs and data.get('positive_acq')[0] != 1: - continue + if not metadata_matches_filters(metadata, args): + continue - acq_grid = data['acq_grid'][:] - n_fft, n_dop_bins = acq_grid.shape - d_max, f_max = np.unravel_index(np.argmax(acq_grid), - acq_grid.shape) - doppler_step = data['doppler_step'][0] - doppler_max = data['doppler_max'][0] - freq = np.arange(n_dop_bins) * doppler_step - doppler_max - delay = np.arange(n_fft) / n_fft * n_chips + n_chips = N_CHIPS.get((metadata["system"], metadata["signal"])) + if n_chips is None: + print( + f"Unknown system/signal {metadata['system']}/" + f"{metadata['signal']}, skipping: {file_path.name}" + ) + continue - # Plot data - # --- Acquisition grid (3D) - fig = plt.figure() - plt.gcf().canvas.manager.set_window_title(filename) - if not lite_view: - ax = fig.add_subplot(111, projection='3d') - X, Y = np.meshgrid(freq, delay) - ax.plot_surface(X, Y, acq_grid, cmap='viridis') - ax.set_ylim([min(delay), max(delay)]) - else: - delay_interp = (np.arange(n_samples_per_chip * n_chips) - / n_samples_per_chip) - spline = CubicSpline(delay, acq_grid) - grid_interp = spline(delay_interp) - ax = fig.add_subplot(111, projection='3d') - X, Y = np.meshgrid(freq, delay_interp) - ax.plot_surface(X, Y, grid_interp, cmap='inferno') - ax.set_ylim([min(delay_interp), max(delay_interp)]) + image_name_root = image_name_from_metadata(args.output_basename, metadata) + if plot_dump(file_path, args.fig_path, image_name_root, n_chips, args): + plotted += 1 - ax.set_xlabel('Doppler shift (Hz)') - ax.set_xlim([min(freq), max(freq)]) - ax.set_ylabel('Code delay (chips)') - ax.set_zlabel('Test Statistics') + if not file_paths: + print(f"No files matching {base}*.mat in {directory}") + prefixes = matching_dump_prefixes(directory) + if prefixes: + print("Available acquisition dump prefixes:") + for prefix in prefixes: + print(f" {prefix}") + print( + "Use --file-prefix to select the input dump prefix. " + "--output-basename only changes saved plot filenames." + ) + elif plotted == 0: + print("No acquisition dumps were plotted.") - plt.savefig(os.path.join(fig_path, filename[:-4]) + '_3D.png') + show_figures(args) - plt.close() - # --- Acquisition grid (2D) - input_power = 100 # Change Test statistics in Doppler wipe-off plot - - fig2, axes = plt.subplots(2, 1, figsize=(8, 6)) - plt.gcf().canvas.manager.set_window_title(filename) - axes[0].plot(freq, acq_grid[d_max, :]) - axes[0].set_xlim([min(freq), max(freq)]) - axes[0].set_xlabel('Doppler shift (Hz)') - axes[0].set_ylabel('Test statistics') - axes[0].set_title(f'Fixed code delay to ' - f'{(d_max - 1) / n_fft * n_chips} chips') - - normalization = (d_samples_per_code ** 4) * input_power - axes[1].plot(delay, acq_grid[:, f_max] / normalization) - axes[1].set_xlim([min(delay), max(delay)]) - axes[1].set_xlabel('Code delay (chips)') - axes[1].set_ylabel('Test statistics') - axes[1].set_title(f'Doppler wipe-off = ' - f'{str((f_max - 1) * doppler_step - doppler_max)} ' - f'Hz') - - plt.tight_layout() - plt.savefig(os.path.join(fig_path, filename[:-4]) + '_2D.png') - # plt.show() - plt.close() +if __name__ == "__main__": + try: + main() + except OSError as exc: + raise SystemExit(f"Error: {exc}") diff --git a/utils/python/plot_tracking_quality_indicators.py b/utils/python/plot_tracking_quality_indicators.py old mode 100644 new mode 100755 index 0078dea6d..994937b38 --- a/utils/python/plot_tracking_quality_indicators.py +++ b/utils/python/plot_tracking_quality_indicators.py @@ -1,15 +1,12 @@ +#!/usr/bin/env python3 """ plot_tracking_quality_indicators.py - Irene Pérez Riega, 2023. iperrie@inta.es - Minhaj Uddin Ahmad, 2026. mahmad12@crimson.ua.edu + Reads GNSS-SDR tracking dump binary files and plots C/N0 plus carrier-lock + quality indicators across channels. - Modifiable in the file: - channels - Number of channels - first_channel - Number of the first channel - path - Path to folder which contains raw files - fig_path - Path where doppler plots will be save - file_prefix - Fixed part of the tracking dump files names + File format: + {input_path}/{file_prefix}{channel}.dat ----------------------------------------------------------------------------- @@ -22,55 +19,126 @@ ----------------------------------------------------------------------------- """ +import argparse +from pathlib import Path import matplotlib.pyplot as plt import numpy as np -import os + from lib.dll_pll_veml_read_tracking_dump import dll_pll_veml_read_tracking_dump +from lib.dump_filename import resolve_dump_prefix +from lib.plot_format import add_output_format_argument, apply_publication_style -GNSS_tracking = [] -plot_names = [] -# ---------- CHANGE HERE: -channels = 5 -first_channel = 0 -path = '../../../out/' -fig_path = '../../../PLOTS/TrackingQualityIndicator' -file_prefix = "track_ch" +def parse_args(): + parser = argparse.ArgumentParser( + description="Plot tracking quality indicators from GNSS-SDR dumps." + ) + parser.add_argument( + "-i", + "--input-path", + type=Path, + default=Path("."), + help="Directory containing tracking .dat dumps (default: .).", + ) + parser.add_argument( + "-o", + "--fig-path", + type=Path, + default=Path("plots/tracking-quality"), + help="Directory where plots are saved.", + ) + parser.add_argument( + "--file-prefix", + default="track_ch", + help="GNSS-SDR Tracking.dump_filename value (default: track_ch). May " + "include a directory and extension; the matching .dat " + "files are read, resolved against --input-path.", + ) + parser.add_argument( + "--channels", + type=int, + default=5, + help="Number of channels to read.", + ) + parser.add_argument( + "--first-channel", + type=int, + default=0, + help="First channel number in the dump filenames.", + ) + parser.add_argument( + "--signal-type", + type=str.upper, + default="1C", + metavar="CODE", + help="GNSS-SDR signal code (e.g. 1C, 5X, 7X, E6, J1) used as the " + "C/N0 legend label. The tracking dump stores only the PRN, so the " + "signal is supplied here (default: 1C).", + ) + parser.add_argument( + "--show", + action="store_true", + help="Display figures interactively after saving them.", + ) + add_output_format_argument(parser) + return parser.parse_args() -for N in range(1, channels + 1): - tracking_log_path = os.path.join(path, - f'{file_prefix}{N-1+first_channel}.dat') - GNSS_tracking.append(dll_pll_veml_read_tracking_dump(tracking_log_path)) -if not os.path.exists(fig_path): - os.makedirs(fig_path) +def read_tracking_dumps(args): + directory, base = resolve_dump_prefix(args.file_prefix, args.input_path) + dumps = [] + for channel in range(args.first_channel, args.first_channel + args.channels): + tracking_log_path = directory / f"{base}{channel}.dat" + dumps.append(dll_pll_veml_read_tracking_dump(tracking_log_path)) + return dumps -# Plot tracking quality indicators -# First plot -plt.figure() -plt.gcf().canvas.manager.set_window_title('Carrier lock test output for all ' - 'the channels') -plt.title('Carrier lock test output for all the channels') -for n in range(len(GNSS_tracking)): - plt.plot(GNSS_tracking[n]['carrier_lock_test']) - plot_names.append(f'SV {str(round(np.mean(GNSS_tracking[n]["PRN"])))}') -plt.legend(plot_names) -plt.savefig(os.path.join(fig_path, - f'carrier_lock_test ' - f'{str(round(np.mean(GNSS_tracking[n]["PRN"])))}')) -# plt.show() -# Second plot -plt.figure() -plt.gcf().canvas.manager.set_window_title('Carrier CN0 output for all the ' - 'channels') -plt.title('Carrier CN0 output for all the channels') -for n in range(len(GNSS_tracking)): - plt.plot(GNSS_tracking[n]['CN0_SNV_dB_Hz']) - plot_names.append(f'CN0_SNV_dB_Hz ' - f'{str(round(np.mean(GNSS_tracking[n]["PRN"])))}') -plt.legend(plot_names) -plt.savefig(os.path.join( - fig_path, f'SV {str(round(np.mean(GNSS_tracking[n]["PRN"])))}')) -# plt.show() +def mean_prn_label(tracking): + return str(round(np.mean(tracking["PRN"]))) + + +def main(): + args = parse_args() + args.fig_path.mkdir(parents=True, exist_ok=True) + apply_publication_style() + + gnss_tracking = read_tracking_dumps(args) + + plt.figure() + plt.gcf().canvas.manager.set_window_title( + "Carrier lock test output for all the channels" + ) + plt.title("Carrier lock test output for all the channels") + carrier_lock_names = [] + for tracking in gnss_tracking: + plt.plot(tracking["carrier_lock_test"]) + carrier_lock_names.append(f"SV {mean_prn_label(tracking)}") + plt.legend(carrier_lock_names) + plt.savefig(args.fig_path / f"carrier_lock_test.{args.output_format}") + if not args.show: + plt.close() + + plt.figure() + plt.gcf().canvas.manager.set_window_title("Carrier CN0 output for all channels") + plt.title("Carrier CN0 output for all channels") + cn0_names = [] + for tracking in gnss_tracking: + plt.plot(tracking["CN0_SNV_dB_Hz"]) + cn0_names.append(f"{args.signal_type} PRN {mean_prn_label(tracking)}") + plt.legend(cn0_names) + plt.savefig(args.fig_path / f"CN0_SNV_dB_Hz.{args.output_format}") + if not args.show: + plt.close() + + # Show all saved figures with a single plt.show() to avoid the repeated + # show()/close() cycle that can crash interactive backends on macOS. + if args.show: + plt.show() + + +if __name__ == "__main__": + try: + main() + except OSError as exc: + raise SystemExit(f"Error: {exc}")