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]:
array([-1.  ,  0.  ,  3.12])

In [6]:
all_plots(control)



In [7]:
all_plots(recovery)