Measurement quality from affordable GNSS receivers

Measurement quality from affordable GNSS receivers#

In this section we will assess some basic differences in terms of GNSS observables between a geodetic-grade and an affordable receiver.

Hide code cell source 
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from roktools import rinex
from roktools.time import compute_elapsed_seconds
from roktools.stats import rms

from roktools.gnss.types import ConstellationId, TrackingChannel, Satellite
from roktools.gnss.observables import compute_code_minus_carrier, compute_geometry_free
from roktools.gnss.edit import mark_time_gap, detrend, remove_mean

%matplotlib widget

Using the roktools library, the RINEX data from the receivers will be stored into a pandas DataFrame (df for short) for manipulation.

# DataFrame with data from a geodetic receiver
#df_geodetic = rinex.to_dataframe('../assets/SEYG00SYC_R_20140581500_05H_01S_MO.rnx')
#df_geodetic.to_parquet('../assets/SEYG00SYC_R_20140581500_05H_01S_MO.rnx.parquet')
df_geodetic = pd.read_parquet('../assets/SEYG00SYC_R_20140581500_05H_01S_MO.rnx.parquet')

# DataFrame with an affordable receiver
#df_afford = rinex.to_dataframe('../assets/MTIC00ESP_R_20191221131_05H_01S_MO.rnx')
#df_afford.to_parquet('../assets/MTIC00ESP_R_20191221131_05H_01S_MO.rnx.parquet')
df_afford = pd.read_parquet('../assets/MTIC00ESP_R_20191221131_05H_01S_MO.rnx.parquet')

The DataFrame can be consider a CSV file of sorts, where each row has a time tag, satellite, observables and other fields that will be explained and used later in the notebook. A preview of the contents can be obtained with the head method

df_geodetic.head()
epoch constellation sat channel signal range phase doppler snr flag slip
0 2014-02-27 15:00:00 E E12 1X E121X 2.492298e+07 1.309712e+08 2196.344 47.0 00000000 0
1 2014-02-27 15:00:00 E E12 5X E125X 2.492299e+07 9.780326e+07 NaN 47.0 00000000 0
2 2014-02-27 15:00:00 E E12 7X E127X 2.492299e+07 1.003546e+08 NaN 46.3 00000000 0
3 2014-02-27 15:00:00 E E12 8X E128X 2.492299e+07 9.907894e+07 NaN 52.2 00000000 0
4 2014-02-27 15:00:00 E E19 1X E191X 2.500360e+07 1.313950e+08 -2145.535 48.7 00000000 0

Observable types#

Once we have loaded the RINEX files into DataFrames, we can perform some basic checks on the differences between geodetic and affordable GNSS data.

The following example gives you the channels tracked by the receiver for each GNSS constellation.

The channel corresponds to the last two characters (i.e. band and attribute) of the RINEX observation code (see Section 5.2.17). For instance 1C for GPS means the observables obtained with the C/A tracking at the L1 frequency.

To know the observables for constellation, we will use the groupby method of pandas.

# We will group the data using various criteria
columns = ['constellation', 'channel']

# Use groupby() to group by the two columns and apply unique()
unique_combinations = df_geodetic.groupby(columns).size()

# Print the unique combinations (along with the number of samples for each
# tracking channel)
print(unique_combinations)

constellation  channel
E              1X         18910
               5X         18910
               7X         18910
               8X         18910
G              1C         87774
               2W         87774
               2X         87774
               5X         87774
R              1C         54304
               1P         54304
               2C         54304
               2P         54304
S              1C         25134
dtype: int64

We can now perform the same for the affordable receivers

unique_combinations = df_afford.groupby(columns).size()

print(unique_combinations)

constellation  channel
C              2I          95280
E              1C         106261
               7Q         106261
G              1C         138731
               2L         138731
R              1C         111380
               2C         111380
dtype: int64

In these examples, the following basic differences are observed:

  • the geodetic receiver tracks various frequencies (GPS L1/L2/L5, Galileo E1/E5a/E5b/E5, …) whereas the affordable receiver tracks typically two frequencies (GPS L1/L2, Galileo E1/E5b, …)

  • the affordable receiver does not attempt to track encrypted codes (i.e. GPS P code) by means of e.g. z-tracking loops. Only civilian codes (e.g. GPS L2C) are used.

Some other strenghts of affordable receivers:

  • Availability of SNR and Doppler measurements (not always available in 30s or high rate CORS data)

  • High rate up to 0.1s (or even higher) available for affordable measurements

Code noise: Detrended code-minus-carrier#

The observable code noise of a GNSS receiver can be estimated using the code-minus-carrier combination. This section illustrates how to estimate to check some basic differences between receiver types

The steps to be followed to compute the unbiased detrended Code minus carrier are detailed in the following steps (taking as example the Geodetic receiver and then applying it to the Affordable receiver)

  1. Edit the data and find phase breaks such as cycle slips. In this example, since receivers already provide with Loss-of-Lock Indicator (LLI), we will only mark phase breaks due to data time gaps

df_geodetic_cmc = mark_time_gap(df_geodetic)

  1. Proceed to compute the code minus carrier (CMC) observable. Each row of the DataFrame contains the range and phase, so there will be as many CMC observables as rows. A note of caution: the phase is usually expressed in cycles, therefore we will need to get the wavelength from the tracking channel in order to be consistent with the units

df_geodetic_cmc = compute_code_minus_carrier(df_geodetic_cmc)

df_geodetic_cmc.head()
epoch constellation sat channel signal range phase doppler snr flag slip arc_id wl cmc
0 2014-02-27 15:00:00 E E12 1X E121X 2.492298e+07 1.309712e+08 2196.344 47.0 00000000 False 0 0.190294 -13.752709
1 2014-02-27 15:00:00 E E12 5X E125X 2.492299e+07 9.780326e+07 NaN 47.0 00000000 False 0 0.254828 -23.022230
2 2014-02-27 15:00:00 E E12 7X E127X 2.492299e+07 1.003546e+08 NaN 46.3 00000000 False 0 0.248349 -20.960017
3 2014-02-27 15:00:00 E E12 8X E128X 2.492299e+07 9.907894e+07 NaN 52.2 00000000 False 0 0.251547 -22.530810
4 2014-02-27 15:00:00 E E19 1X E191X 2.500360e+07 1.313950e+08 -2145.535 48.7 00000000 False 0 0.190294 -36.599890 
# Plot a CMC for a channel and satellite
signal = 'E121X'
df_sample = df_geodetic_cmc[df_geodetic_cmc['signal'] == signal]
plt.close()
plt.plot(compute_elapsed_seconds(df_sample['epoch'])/3600, df_sample['cmc'], ',k')
plt.xlabel('Elapsed hours')
plt.ylabel('CMC [m]')
_ = plt.title(f'Code-minus-carrier for {signal}')

  1. Because the CMC contains twice the ionosphere, which is a nuisance parameter for us (because we are interested in the noise). We will proceed to remove its contribution with a simple detrending, based on a rolling window

n_samples = 5
df_geodetic_cmc = detrend(df_geodetic_cmc, 'cmc', n_samples)

# Plot the detrened CMC for a channel and satellite
signal = 'E121X'
df_sample = df_geodetic_cmc[df_geodetic_cmc['signal'] == signal]
plt.close()
plt.plot(compute_elapsed_seconds(df_sample['epoch'])/3600, df_sample['cmc_detrended'], ',k')
plt.xlabel('Elapsed hours')
plt.ylabel('Detrended CMC [m]')
plt.ylim(-1, 1)
_ = plt.title(f'Detrended Code-minus-carrier for {signal}')

# Analysis for Ublox receiver
df_afford_cmc = mark_time_gap(df_afford)
df_afford_cmc = compute_code_minus_carrier(df_afford_cmc)
df_afford_cmc = detrend(df_afford_cmc, 'cmc', n_samples)

Let’s compute now an estimate of the code noise for the geodetic grade and the affordable receiver, for a specific band and constellation.

# GLONASS is excluded as the slot number to know the frequency (and hence the wavelength)
# is missing
condition1 = df_geodetic_cmc['constellation'] == ConstellationId.GPS.value
condition2 = df_geodetic_cmc['channel'] == '2W'
#condition2 = df_geodetic_cmc['channel'] == '1C'
df_tmp_g = df_geodetic_cmc[condition1 & condition2]

condition1 = df_afford_cmc['constellation'] == ConstellationId.GPS.value
condition2 = df_afford_cmc['channel'] == '2L'
#condition2 = df_afford_cmc['channel'] == '1C'
df_tmp_a = df_afford_cmc[condition1 & condition2]

noise_samples_geodetic = df_tmp_g['cmc_detrended']
noise_samples_afford = df_tmp_a['cmc_detrended']

rms_geodetic = rms(noise_samples_geodetic)
rms_afford = rms(noise_samples_afford)

print(f'Geodetic   receiver code noise: {rms_geodetic:.2} m')
print(f'Affordable receiver code noise: {rms_afford:.2} m')

Geodetic   receiver code noise: 0.34 m
Affordable receiver code noise: 0.32 m

Observable noise: Detrended LI#

In addition, to the estimation of the code, for multifrequency receivers, the geometry-free (ionospheric) combination can be used to estimate the noise of the phase observables

constellation = ConstellationId.GPS
channel_a = TrackingChannel.from_string('1C')
channel_b = TrackingChannel.from_string('2W')

# Let's reuse the previous dataframe that has been already edited (time gaps marked)
df_geodetic_li = compute_geometry_free(df_geodetic_cmc, constellation, channel_a, channel_b)

The ionospheric combination (LI) contains the ionospheric delay, constant terms due to the phase ambiguities of the carrier phases used to build the combination as well as the phase hardware biases and noise of the phase observable

An example of LI for a particular satellite is shown in the Figure below

df_geodetic_li
epoch constellation_a sat channel_a signal_a range_a phase_a slip_a constellation_b channel_b signal_b range_b phase_b slip_b li_m pi_m
0 2014-02-27 15:00:00 G G01 1C G011C 2.434418e+07 1.279295e+08 False G 2W G012W 2.434421e+07 9.968548e+07 False -36.684348 25.582
1 2014-02-27 15:00:00 G G04 1C G041C 2.382347e+07 1.251933e+08 False G 2W G042W 2.382349e+07 9.755329e+07 False -17.711599 17.414
2 2014-02-27 15:00:00 G G05 1C G051C 2.506722e+07 1.317292e+08 False G 2W G052W 2.506725e+07 1.026462e+08 False -30.585518 26.824
3 2014-02-27 15:00:00 G G07 1C G071C 2.346388e+07 1.233036e+08 False G 2W G072W 2.346390e+07 9.608087e+07 False -40.325221 16.801
4 2014-02-27 15:00:00 G G08 1C G081C 2.418957e+07 1.271171e+08 False G 2W G082W 2.418959e+07 9.905244e+07 False -38.773982 20.829
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
87769 2014-02-27 19:59:59 G G10 1C G101C 2.304162e+07 1.210848e+08 False G 2W G102W 2.304163e+07 9.435189e+07 False -23.682447 11.836
87770 2014-02-27 19:59:59 G G15 1C G151C 2.285345e+07 1.200958e+08 False G 2W G152W 2.285346e+07 9.358119e+07 False -19.520759 12.805
87771 2014-02-27 19:59:59 G G24 1C G241C 2.167144e+07 1.138841e+08 False G 2W G242W 2.167145e+07 8.874101e+07 False -38.376173 15.293
87772 2014-02-27 19:59:59 G G26 1C G261C 2.301703e+07 1.209555e+08 False G 2W G262W 2.301705e+07 9.425102e+07 False -0.104817 11.543
87773 2014-02-27 19:59:59 G G28 1C G281C 2.421823e+07 1.272677e+08 False G 2W G282W 2.421824e+07 9.916970e+07 False -12.247415 11.145

87774 rows × 16 columns

sat = 'G04'

condition1 = df_geodetic_li['sat'] == sat
df_tmp = df_geodetic_li[condition1]

t = compute_elapsed_seconds(df_tmp['epoch'])
li = df_tmp['li_m']

plt.close()
plt.plot(t, li, '.')
plt.xlabel('Time [seconds]')
plt.ylabel('LI [m]')
plt.title(f'LI combination for {sat}, Geoetic receiver (SEY1)')

Text(0.5, 1.0, 'LI combination for G04, Geoetic receiver (SEY1)')

import pandas as pd

def detrend_li_combination(df:pd.DataFrame):
    """
    Perform a detrending of the LI combination (if present in the input dataframe)
    """

    df['slip_li'] =  df['slip_a'] | df['slip_b']
    df['arc_id'] = df.groupby(['signal_a', 'signal_b'])['slip_li'].transform('cumsum')
    trend = df.groupby(['signal_a', 'signal_b', 'arc_id'])['li_m'].transform(lambda x: x.rolling(n_samples).mean())
    df['li_trend'] = trend
    df['li_detrended'] = df['li_m'] - df['li_trend']

detrend_li_combination(df_geodetic_li)

condition1 = df_geodetic_li['signal_a'] == 'G041C'
condition2 = df_geodetic_li['signal_b'] == 'G042W'

df_sample = df_geodetic_li[condition1 & condition2]
plt.close()
plt.plot(df_sample['epoch'], df_sample['li_detrended'], '.')
plt.ylim(-0.05, 0.05)

(-0.05, 0.05)

We can proceed with the same analysis for the affordable receiver, taking into account that the observables are slightly different

channel_a = TrackingChannel.from_string('1C')
channel_b = TrackingChannel.from_string('2L')

df_afford_li = compute_geometry_free(df_afford_cmc, constellation, channel_a, channel_b)

detrend_li_combination(df_afford_li)

condition1 = df_afford_li['signal_a'] == 'G241C'
condition2 = df_afford_li['signal_b'] == 'G242L'

df_sample = df_afford_li[condition1 & condition2]

plt.close()
plt.plot(df_sample['epoch'], df_sample['li_detrended'], '.')
plt.ylim(-0.05, 0.05)

(-0.05, 0.05)

The detrended time series can be used to compute the noise figures of the carrier phase for both types of receivers

# Remove outliers for the statistics
noise_li_afford = df_afford_li['li_detrended'][(df_afford_li['li_detrended'] > -1) & (df_afford_li['li_detrended'] < 1)]
noise_li_geodetic = df_geodetic_li['li_detrended'][(df_geodetic_li['li_detrended'] > -1) & (df_geodetic_li['li_detrended'] < 1)]

print(f"Phase Noise Affordable receiver: {np.std(noise_li_afford)*100:.0f} cm")
print(f"Phase Noise Geodetic   receiver: {np.std(noise_li_geodetic)*100:.0f} cm") 

Phase Noise Affordable receiver: 1 cm
Phase Noise Geodetic   receiver: 2 cm

As it can be seen, the differences between receivers is not substantial. This can be further confirmed with the actual distribution of the noise samples, which follow a very similar Gaussian pattern for both receivers

plt.close()
bins = np.linspace(-0.05, 0.05, 200)
_ = plt.hist(df_geodetic_li['li_detrended'], bins=bins, histtype='step', label='Geodetic')
_ = plt.hist(df_afford_li['li_detrended'], bins=bins, histtype='step', label='Affordable')
plt.legend()
plt.title('Histogram of carrier phase noise (estimated with LI)')
plt.xlabel('Code phase error [m]')
plt.ylabel('Count')

Text(0, 0.5, 'Count')
Contents