Skip to content

Commit 804d192

Browse files
boriodaborioda
authored
Plotting routines
Simple scripts to plot the different corrections. The ADEV computation is also included.
1 parent 1793f0e commit 804d192

File tree

5 files changed

+462
-0
lines changed

5 files changed

+462
-0
lines changed

‎plot/clk_adev.py

+107
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,107 @@
1+
# -*- coding: utf-8 -*-
2+
"""
3+
Created on Mon Jan 30 17:05:49 2023
4+
5+
@author: Daniele
6+
"""
7+
8+
import pandas as pd
9+
import matplotlib.pyplot as plt
10+
import numpy as np
11+
import allantools as at
12+
13+
14+
import matplotlib as mp
15+
fsize = 16
16+
mp.rc('xtick', labelsize=fsize)
17+
mp.rc('ytick', labelsize=fsize)
18+
19+
20+
plt.style.use('ggplot')
21+
gnss_ids = ["G", "", "E"]
22+
gnss_spell = ["GPS", "", "Galileo"]
23+
24+
# load the csv file with the clock corrections
25+
# filename = "SEPT1010.22_has_clk.csv"
26+
filename = "SEPT293_GALRawCNAV.zip_has_clk.csv"
27+
fname = filename[:filename.rfind(".")]
28+
df = pd.read_csv(filename)
29+
30+
# determine which GNSS is present in the correction file
31+
GNSS = np.unique(df["gnssID"].values)
32+
33+
###################### SOME INFO #####################
34+
# light speed [m/sec]
35+
v_light = 299792458
36+
37+
# # Frequency E1
38+
# f1 = 1575.42e6
39+
40+
# # Frequency E5a
41+
# f2 = 1176.45e6
42+
43+
# # Wavelength E1
44+
# lambda1 = v_light/ f1
45+
46+
# # Wavelength E5a
47+
# lambda2 = v_light / f2
48+
49+
# # Coefficient linear combiantion
50+
# a = f1**2 / (f1**2 + f2**2)
51+
# b = f2**2 / (f1**2 - f2**2)
52+
53+
# # Equivalent wavelength
54+
# lambda_eq = lambda1 * lambda2 / (a * lambda1 + b * lambda2)
55+
56+
# Tau for ADEV
57+
tau = np.logspace(1, 5, 40)
58+
59+
#######################################################
60+
61+
# make different plots for each gnss
62+
for gnss in GNSS :
63+
# filter by gnss
64+
gnss_df = df[df["gnssID"] == gnss]
65+
66+
# Create a new plot
67+
fig, ax = plt.subplots()
68+
fig.set_size_inches((12, 8))
69+
# Determine the list of prn
70+
sats = np.unique(gnss_df["PRN"].values)
71+
72+
# Make a plot for each satellite
73+
for sat in sats :
74+
sat_df = gnss_df[gnss_df["PRN"] == sat]
75+
76+
# build the time index
77+
time = np.floor( sat_df["ToW"].values / 3600 ) * 3600 + sat_df["ToH"].values
78+
79+
# build the clock correction
80+
clk_corr = sat_df["delta_clock_c0"].values * sat_df["multiplier"].values
81+
82+
# remove duplicated values
83+
time, ind = np.unique(time, return_index = True)
84+
clk_corr = clk_corr[ind] / v_light
85+
86+
# measurement rate
87+
Ts = np.median(np.diff(time))
88+
rate = 1 / Ts
89+
clk_freq = np.diff(clk_corr) / Ts
90+
91+
t, ad, ade, adn = at.oadev( clk_corr, rate = rate, \
92+
data_type = "phase", taus = tau[:-3])
93+
94+
# finally plot the result
95+
ax.loglog(t, ad, "-.", label = f"{gnss_ids[gnss]}{sat}")
96+
97+
# Some cosmetics for the plots
98+
# ax.tick_params(axis='x', labelrotation = 45)
99+
ax.set_xlabel("Averaging Interval [s]", fontsize = fsize)
100+
ax.set_ylabel("ADEV", fontsize = fsize)
101+
ax.legend(loc=(1.05, 0), fontsize = 14, ncol=2)
102+
ax.autoscale(tight=True)
103+
ax.set_title(f"{gnss_spell[gnss]}", fontsize = fsize)
104+
plt.tight_layout()
105+
106+
# Save as pdf
107+
plt.savefig(f"{fname}_{gnss_spell[gnss]}_adev.png", bbox_inches="tight")

‎plot/plot_cb.py

+101
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,101 @@
1+
#!/usr/bin/env python3
2+
# -*- coding: utf-8 -*-
3+
"""
4+
Created on Mon Jan 23 13:10:31 2023
5+
6+
@author: daniele
7+
"""
8+
9+
import pandas as pd
10+
import matplotlib.pyplot as plt
11+
import numpy as np
12+
13+
import matplotlib as mp
14+
fsize = 16
15+
mp.rc('xtick', labelsize=fsize)
16+
mp.rc('ytick', labelsize=fsize)
17+
18+
19+
plt.style.use('ggplot')
20+
gnss_ids = ["G", "", "E"]
21+
gnss_spell = ["GPS", "", "Galileo"]
22+
gnss_ind = [1, -1, 0]
23+
24+
# Signal table - indentifies the different signals
25+
signal_table = [["E1-B I/NAV OS", "L1 C/A"],
26+
["E1-C", ""],
27+
["E1-B + E1-C", ""],
28+
["E5a-I F/NAV OS", "L1C(D)"],
29+
["E5a-Q", "L1C(P)"],
30+
["E5a-I+E5a-Q", "L1C(D+P)"],
31+
["E5b-I I/NAV OS", "L2 CM"],
32+
["E5b-Q", "L2 CL"],
33+
["E5b-I+E5b-Q", "L2 CM+CL"],
34+
["E5-I", "L2 P"],
35+
["E5-Q", "Reserved"],
36+
["E5-I + E5-Q", "L5 I"],
37+
["E6-B C/NAV HAS", "L5 Q"],
38+
["E6-C", "L5 I + L5 Q"],
39+
["E6-B + E6-C", ""],
40+
["", ""]]
41+
42+
# load the csv file with the code bias corrections
43+
filename = "SEPT293_GALRawCNAV.zip_has_cb.csv"
44+
fname = filename[:filename.rfind(".")]
45+
df = pd.read_csv(filename)
46+
47+
# determine which GNSS is present in the correction file
48+
GNSS = np.unique(df["gnssID"].values)
49+
50+
# make different plots for each gnss
51+
for gnss in GNSS :
52+
# filter by gnss
53+
gnss_df = df[df["gnssID"] == gnss]
54+
55+
# determine the number of unique signals in the file
56+
signals = np.unique(gnss_df["signal"].values)
57+
58+
# Create a new plot with as many subplots as signals
59+
fig, ax = plt.subplots(nrows=len(signals), ncols = 1)
60+
fig.set_size_inches((8, 8))
61+
62+
# Loop on the different signals
63+
for ii, signal in enumerate(signals) :
64+
signal_df = gnss_df[gnss_df["signal"] == signal]
65+
66+
# Determine the list of prn
67+
sats = np.unique(signal_df["PRN"].values)
68+
69+
# Make a plot for each satellite
70+
for sat in sats :
71+
sat_df = signal_df[signal_df["PRN"] == sat]
72+
73+
# build the time index
74+
time = np.floor( sat_df["ToW"].values / 3600 ) * 3600 + sat_df["ToH"].values
75+
76+
# build the code bias
77+
code_bias = sat_df["code_bias"].values
78+
79+
# remove duplicated values
80+
time, ind = np.unique(time, return_index = True)
81+
code_bias = code_bias[ind]
82+
83+
# finally plot the result
84+
ax[ii].plot(time, code_bias, ".--", label = f"{gnss_ids[gnss]}{sat}")
85+
86+
ax[ii].set_ylabel("Code bias [m]", fontsize = fsize)
87+
ax[ii].autoscale(tight=True)
88+
ax[ii].set_title(f"{gnss_spell[gnss]} - {signal_table[signal - 1][gnss_ind[gnss]]}", fontsize = fsize)
89+
if ax[ii] != ax[-1] :
90+
ax[ii].set_xticks([])
91+
92+
93+
# Some cosmetics for the plots
94+
ax[-1].tick_params(axis='x', labelrotation = 45)
95+
ax[-1].set_xlabel("Time of Week [s]", fontsize = fsize)
96+
ax[-1].legend(loc=(1.05, 0), fontsize = 14, ncol=2)
97+
98+
# plt.tight_layout()
99+
100+
# Save as pdf
101+
plt.savefig(f"{fname}_{gnss_spell[gnss]}.png", bbox_inches="tight")

‎plot/plot_clk.py

+70
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
#!/usr/bin/env python3
2+
# -*- coding: utf-8 -*-
3+
"""
4+
Created on Mon Jan 23 13:10:31 2023
5+
6+
@author: daniele
7+
"""
8+
9+
import pandas as pd
10+
import matplotlib.pyplot as plt
11+
import numpy as np
12+
13+
import matplotlib as mp
14+
fsize = 16
15+
mp.rc('xtick', labelsize=fsize)
16+
mp.rc('ytick', labelsize=fsize)
17+
18+
19+
plt.style.use('ggplot')
20+
gnss_ids = ["G", "", "E"]
21+
gnss_spell = ["GPS", "", "Galileo"]
22+
23+
# load the csv file with the clock corrections
24+
# filename = "SEPT1010.22_has_clk.csv"
25+
filename = "SEPT293_GALRawCNAV.zip_has_clk.csv"
26+
fname = filename[:filename.rfind(".")]
27+
df = pd.read_csv(filename)
28+
29+
# determine which GNSS is present in the correction file
30+
GNSS = np.unique(df["gnssID"].values)
31+
32+
# make different plots for each gnss
33+
for gnss in GNSS :
34+
# filter by gnss
35+
gnss_df = df[df["gnssID"] == gnss]
36+
37+
# Create a new plot
38+
fig, ax = plt.subplots()
39+
fig.set_size_inches((12, 8))
40+
# Determine the list of prn
41+
sats = np.unique(gnss_df["PRN"].values)
42+
43+
# Make a plot for each satellite
44+
for sat in sats :
45+
sat_df = gnss_df[gnss_df["PRN"] == sat]
46+
47+
# build the time index
48+
time = np.floor( sat_df["ToW"].values / 3600 ) * 3600 + sat_df["ToH"].values
49+
50+
# build the clock correction
51+
clk_corr = sat_df["delta_clock_c0"].values * sat_df["multiplier"].values
52+
53+
# remove duplicated values
54+
time, ind = np.unique(time, return_index = True)
55+
clk_corr = clk_corr[ind]
56+
57+
# finally plot the result
58+
ax.plot(time, clk_corr, ".", label = f"{gnss_ids[gnss]}{sat}")
59+
60+
# Some cosmetics for the plots
61+
ax.tick_params(axis='x', labelrotation = 45)
62+
ax.set_xlabel("Time of Week [s]", fontsize = fsize)
63+
ax.set_ylabel("Clock Corrections [m]", fontsize = fsize)
64+
ax.legend(loc=(1.05, 0), fontsize = 14, ncol=2)
65+
ax.autoscale(tight=True)
66+
ax.set_title(f"{gnss_spell[gnss]}", fontsize = fsize)
67+
plt.tight_layout()
68+
69+
# Save as pdf
70+
plt.savefig(f"{fname}_{gnss_spell[gnss]}.png", bbox_inches="tight")

‎plot/plot_cp.py

+101
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,101 @@
1+
#!/usr/bin/env python3
2+
# -*- coding: utf-8 -*-
3+
"""
4+
Created on Mon Jan 23 13:10:31 2023
5+
6+
@author: daniele
7+
"""
8+
9+
import pandas as pd
10+
import matplotlib.pyplot as plt
11+
import numpy as np
12+
13+
import matplotlib as mp
14+
fsize = 16
15+
mp.rc('xtick', labelsize=fsize)
16+
mp.rc('ytick', labelsize=fsize)
17+
18+
19+
plt.style.use('ggplot')
20+
gnss_ids = ["G", "", "E"]
21+
gnss_spell = ["GPS", "", "Galileo"]
22+
gnss_ind = [1, -1, 0]
23+
24+
# Signal table - indentifies the different signals
25+
signal_table = [["E1-B I/NAV OS", "L1 C/A"],
26+
["E1-C", ""],
27+
["E1-B + E1-C", ""],
28+
["E5a-I F/NAV OS", "L1C(D)"],
29+
["E5a-Q", "L1C(P)"],
30+
["E5a-I+E5a-Q", "L1C(D+P)"],
31+
["E5b-I I/NAV OS", "L2 CM"],
32+
["E5b-Q", "L2 CL"],
33+
["E5b-I+E5b-Q", "L2 CM+CL"],
34+
["E5-I", "L2 P"],
35+
["E5-Q", "Reserved"],
36+
["E5-I + E5-Q", "L5 I"],
37+
["E6-B C/NAV HAS", "L5 Q"],
38+
["E6-C", "L5 I + L5 Q"],
39+
["E6-B + E6-C", ""],
40+
["", ""]]
41+
42+
# load the csv file with the code bias corrections
43+
filename = "SEPT293_GALRawCNAV.zip_has_cp.csv"
44+
fname = filename[:filename.rfind(".")]
45+
df = pd.read_csv(filename)
46+
47+
# determine which GNSS is present in the correction file
48+
GNSS = np.unique(df["gnssID"].values)
49+
50+
# make different plots for each gnss
51+
for gnss in GNSS :
52+
# filter by gnss
53+
gnss_df = df[df["gnssID"] == gnss]
54+
55+
# determine the number of unique signals in the file
56+
signals = np.unique(gnss_df["signal"].values)
57+
58+
# Create a new plot with as many subplots as signals
59+
fig, ax = plt.subplots(nrows=len(signals), ncols = 1)
60+
fig.set_size_inches((8, 8))
61+
62+
# Loop on the different signals
63+
for ii, signal in enumerate(signals) :
64+
signal_df = gnss_df[gnss_df["signal"] == signal]
65+
66+
# Determine the list of prn
67+
sats = np.unique(signal_df["PRN"].values)
68+
69+
# Make a plot for each satellite
70+
for sat in sats :
71+
sat_df = signal_df[signal_df["PRN"] == sat]
72+
73+
# build the time index
74+
time = np.floor( sat_df["ToW"].values / 3600 ) * 3600 + sat_df["ToH"].values
75+
76+
# build the code bias
77+
phase_bias = sat_df["phase_bias"].values
78+
79+
# remove duplicated values
80+
time, ind = np.unique(time, return_index = True)
81+
phase_bias = phase_bias[ind]
82+
83+
# finally plot the result
84+
ax[ii].plot(time, phase_bias, ".--", label = f"{gnss_ids[gnss]}{sat}")
85+
86+
ax[ii].set_ylabel("Phase bias [cycles]", fontsize = fsize)
87+
ax[ii].autoscale(tight=True)
88+
ax[ii].set_title(f"{gnss_spell[gnss]} - {signal_table[signal - 1][gnss_ind[gnss]]}", fontsize = fsize)
89+
if ax[ii] != ax[-1] :
90+
ax[ii].set_xticks([])
91+
92+
93+
# Some cosmetics for the plots
94+
ax[-1].tick_params(axis='x', labelrotation = 45)
95+
ax[-1].set_xlabel("Time of Week [s]", fontsize = fsize)
96+
ax[-1].legend(loc=(1.05, 0), fontsize = 14, ncol=2)
97+
98+
# plt.tight_layout()
99+
100+
# Save as pdf
101+
plt.savefig(f"{fname}_{gnss_spell[gnss]}.png", bbox_inches="tight")

0 commit comments

Comments
 (0)