In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import georinex as gr
import xarray as xr
from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()
In many ephemeris entries, the SISA and health fields are swapped (probably a bug in some receiver). The function below is used to detect when this happens and undo the swap.
In [2]:
def fix_SISA_health_swap(eph):
for sv in eph.coords['sv']:
swapped = ((np.fmod(eph.sel(sv = sv)['health'], 1) != 0) & ~np.isnan(eph.sel(sv = sv)['health'])) \
| (eph.sel(sv = sv)['health'] == -1)
sisa = np.copy(eph.sel(sv = sv)['SISA'][swapped])
eph.sel(sv = sv)['SISA'][swapped] = eph.sel(sv = sv)['health'][swapped]
eph.sel(sv = sv)['health'][swapped] = sisa
In [3]:
def plot_ephtimes(ephs):
plt.figure(figsize = (10,12), facecolor = 'w')
for sv in ephs.coords['sv']:
if sv.item()[0] == 'E':
eph_times = ephs.sel(sv = sv)['IODnav'].dropna('time').coords['time'].values
plt.plot(eph_times, int(sv.item()[1:3])*np.ones(eph_times.size), '.', color = 'C0')
plt.yticks(np.arange(1,37), [f'E{svn:02d}' for svn in range(1,37)])
plt.grid()
plt.ylabel('SVN')
plt.xlabel('GST time')
plt.title('Ephemeris availability')
def plot_SISA(ephs, freq_pair = 'E5b,E1', timesel = None):
sisa_table = [-1, 0, 3.12, 3.28, 3.44, 3.6]
freq_pair_bit = {'E5a,E1' : 1 << 8, 'E5b,E1' : 1 << 9}
plt.figure(figsize = (10,12), facecolor = 'w')
for sv in ephs.coords['sv']:
if sv.item()[0] == 'E':
for j,sisa_val in enumerate(sisa_table):
sisa = ephs.sel(sv = sv)['SISA'].dropna('time')
datasrc = np.int32(ephs.sel(sv = sv)['DataSrc'].dropna('time').values) & freq_pair_bit[freq_pair] != 0
sisa = sisa[datasrc]
sisa = sisa.coords['time'][sisa == sisa_val]
plt.plot(sisa, int(sv.item()[1:3])*np.ones(sisa.size), '.', color = f'C{j}')
plt.yticks(np.arange(1,37), [f'E{svn:02d}' for svn in range(1,37)])
plt.grid()
plt.ylabel('SVN')
plt.xlabel('GST time')
plt.title(f'SISA {freq_pair} (blue = undefined/unknown, orange = 0, green = 3.12m)')
def plot_DVS(ephs, band = 'E1B'):
dvs_bit = {'E1B' : 1 << 0, 'E5a' : 1 << 3, 'E5b' : 1 << 6}
dvs_src = {'E1B' : 1 << 0 | 1 << 2, 'E5a' : 1 << 1, 'E5b' : 1 << 0 | 1 << 2}
plt.figure(figsize = (10,12), facecolor = 'w')
for sv in ephs.coords['sv']:
if sv.item()[0] == 'E':
dvs = ephs.sel(sv = sv)['health'].dropna('time')
datasrc = np.int32(ephs.sel(sv = sv)['DataSrc'].dropna('time').values) & dvs_src[band] != 0
dvs = dvs[datasrc]
dvs_set = np.int32(dvs) & dvs_bit[band] != 0
dvs_1 = dvs.coords['time'][dvs_set].values
dvs_0 = dvs.coords['time'][~dvs_set].values
plt.plot(dvs_0, int(sv.item()[1:3])*np.ones(dvs_0.size), '.', color = 'green')
plt.plot(dvs_1, int(sv.item()[1:3])*np.ones(dvs_1.size), '.', color = 'red')
plt.yticks(np.arange(1,37), [f'E{svn:02d}' for svn in range(1,37)])
plt.grid()
plt.ylabel('SVN')
plt.xlabel('GST time')
plt.title(f'{band} DVS (green = nav data valid, red = working w/o guarantee)')
def plot_HS(ephs, band = 'E1B'):
hs_shift = {'E1B' : 1, 'E5a' : 4, 'E5b' : 7}
hs_src = {'E1B' : 1 << 0 | 1 << 2, 'E5a' : 1 << 1, 'E5b' : 1 << 0 | 1 << 2}
plt.figure(figsize = (10,12), facecolor = 'w')
for sv in ephs.coords['sv']:
if sv.item()[0] == 'E':
hs = ephs.sel(sv = sv)['health'].dropna('time')
datasrc = np.int32(ephs.sel(sv = sv)['DataSrc'].dropna('time').values) & hs_src[band] != 0
hs = hs[datasrc]
hs_val = (np.int32(hs) >> hs_shift[band]) & 0x3
color = ['green', 'red', 'orange', 'blue']
for h in range(4):
hh = hs.coords['time'][hs_val == h].values
plt.plot(hh, int(sv.item()[1:3])*np.ones(hh.size), '.', color = color[h])
plt.yticks(np.arange(1,37), [f'E{svn:02d}' for svn in range(1,37)])
plt.grid()
plt.ylabel('SVN')
plt.xlabel('GST time')
plt.title(f'{band} HS (green = OK, red = out of service, orange = will be out of service, blue = in test)')
def all_plots(eph):
plot_ephtimes(eph)
plot_SISA(eph, freq_pair = 'E5a,E1')
plot_SISA(eph, freq_pair = 'E5b,E1')
plot_DVS(eph, band = 'E5a')
plot_DVS(eph, band = 'E5b')
plot_DVS(eph, band = 'E1B')
plot_HS(eph, band = 'E5a')
plot_HS(eph, band = 'E5b')
plot_HS(eph, band = 'E1B')
In [4]:
control = gr.load('BRDC00IGS_R_20191800000_01D_MN.rnx.gz')
fix_SISA_health_swap(control)
recovery = xr.merge((gr.load('BRDC00IGS_R_20191970000_01D_MN.rnx.gz'), gr.load('BRDC00IGS_R_20191980000_01D_MN.rnx.gz')))
fix_SISA_health_swap(recovery)
List all the different SISA values that appear in the ephemeris. Usually -1 (NAPA), 0 and 3.12m appear.
In [5]:
x = np.ravel(xr.merge((control, recovery))['SISA'].values)
x = x[~np.isnan(x)]
np.unique(x)
Out[5]:
In [6]:
all_plots(control)
In [7]:
all_plots(recovery)
In [8]:
plot_ephtimes(recovery.sel(time = slice('2019-07-16 18:00', '2019-07-16 21:00')))
In [9]:
plot_SISA(recovery.sel(time = slice('2019-07-17 17:00', '2019-07-17 19:00')))
Below there is some information about the gaps in the ephemeris. The IODnav is increased every 10 minutes, and usually every sequential IODnav appears in the RINEX, but for some reason some IODnav's are skipped (either in the RINEX or in the signal in space).
In [10]:
plt.figure(figsize = (12, 7))
control.sel(sv = 'E36')['IODnav'].dropna('time').plot(marker = '.')
plt.title('IODnav evolution (E36)')
plt.xlabel('GST time')
plt.ylabel('IODnav')
Out[10]:
In [11]:
control.sel(sv = 'E36')['IODnav'].dropna('time')[:20]
Out[11]:
In [12]:
control.sel(sv = 'E36')['IODnav'].dropna('time').coords['time'][:20]
Out[12]:
A short comparison with the facts stated in this video.
In [13]:
failure = xr.merge((gr.load('BRDC00IGS_R_20191920000_01D_MN.rnx.gz'),
gr.load('BRDC00IGS_R_20191910000_01D_MN.rnx.gz')))
fix_SISA_health_swap(failure)
In [14]:
x = np.ravel(failure['SISA'].values)
x = x[~np.isnan(x)]
np.unique(x)
Out[14]:
In [15]:
def get_iodnavs(ephs, sv, freq_pair):
freq_pair_bit = {'E5a,E1' : 1 << 8, 'E5b,E1' : 1 << 9}
iodnavs = ephs.sel(sv = sv)['IODnav'].dropna('time')
datasrc = np.int32(ephs.sel(sv = sv)['DataSrc'].dropna('time').values) & freq_pair_bit[freq_pair] != 0
return iodnavs[datasrc]
def plot_iodnav(ephs, freq_pair = 'E5b,E1'):
plt.figure(figsize = (12,8), facecolor = 'w')
svs = set([sv.item()[:3] for sv in ephs.coords['sv'] if sv.item()[0] == 'E'])
for sv in sorted(svs):
svlabels = [s.item() for s in ephs.coords['sv'] if s.item()[:3] == sv]
iodnavs = xr.merge((get_iodnavs(ephs, s, freq_pair).drop('sv') for s in svlabels))['IODnav']
iodnavs.plot(marker = '.', linestyle = '', label = sv)
plt.legend()
plt.ylabel('IODnav')
plt.xlabel('GST time')
plt.title(f'IODnav {freq_pair}')
plt.grid()
In [16]:
plot_iodnav(control)
In [17]:
plot_iodnav(failure)
In [18]:
plot_iodnav(recovery)
In [19]:
bogus_iodnav = np.where(recovery['IODnav'] > 127)
recovery['IODnav'].isel(time = bogus_iodnav[0], sv = bogus_iodnav[1])
Out[19]:
In [20]:
plt.figure(figsize = (12, 7))
recovery.sel(sv = 'E07')['IODnav'].dropna('time').plot(marker = '.')
plt.title('IODnav evolution (E07)')
plt.xlabel('GST time')
plt.ylabel('IODnav')
Out[20]:
In [21]:
all_plots(failure)