Recipe: ROTI

Recipe: ROTI#

In this recipe we will replicate Figure 2b of [Juan et al., 2017] (for the SEY1 station and GPS PRN26) using the methodology described in Algorithm .

This recipe includes the following steps:

Loading a Rinex file into a pandas DataFrame for later processing
Computing the LI for a satellite
Computing the \(\Delta STEC\) for a satellite
Computing the ROTI

First we will import some necessary modules:

Load a Rinex#

The roktools module contains a utility method that loads a RINEX file into a convenient data structure (pandas DataFrame) that eases data analysis.

rinex_file = '../assets/SEYG00SYC_R_20140581500_05H_01S_MO.rnx'

df = rinex.to_dataframe(rinex_file)

The contents of the DataFrame are layout as a data table, with various columns (pseudorange, phase, Doppler, …). To peek the first contents of the RINEX file and check how the data is organized and which are the differents columns of the table, use the head function.

df.head()

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

Within each column, data is stored as Python classes. For instance, the column sat contains the satellite, but it is not stored as a string, but as a Satellite object. To check the class type, use the following command

df['sat'].dtype

dtype('O')

Unfortunately, the O identifier tells us that the sat is an object, but we do not have information in the specific class. To know the specific class type. Do this instead:

df['sat'].apply(type).head()

  <class 'roktools.gnss.types.Satellite'>
  <class 'roktools.gnss.types.Satellite'>
  <class 'roktools.gnss.types.Satellite'>
  <class 'roktools.gnss.types.Satellite'>
  <class 'roktools.gnss.types.Satellite'>
Name: sat, dtype: object

Now that the whole RINEX is loaded into the DataFrame, you can perform some basics checks such as getting the e.g. list of the satellites contained in the file:

df['sat'].unique()

array([E12, E19, E20, G01, G04, G05, G07, G08, G09, G10, G11, G13, G17,
       G20, G23, G28, R03, R04, R05, R17, R18, R19, S26, S27, S28, E11,
       G26, R06, G02, R15, R16, R20, G15, R21, S20, G24, R09, R10],
      dtype=object)

You can also select the data for just one satellite

# Create a Satellite object, that will be used for the DataFrame indexing
sat = Satellite.from_string('G26')

# Create a "sub-DataFrame" with the contents for this satellite only
df_sat = df[df['sat'] == sat]
df_sat.head()

	epoch	constellation	sat	channel	signal	range	phase	doppler	snr
71704	2014-02-27 15:30:00	ConstellationId.GPS	G26	1C	G261C	2.431778e+07	1.277908e+08	3446.301	30.5
71705	2014-02-27 15:30:00	ConstellationId.GPS	G26	2W	G262W	NaN	NaN	NaN	NaN
71706	2014-02-27 15:30:00	ConstellationId.GPS	G26	2X	G262X	NaN	NaN	NaN	NaN
71707	2014-02-27 15:30:00	ConstellationId.GPS	G26	5X	G265X	NaN	NaN	NaN	NaN
71787	2014-02-27 15:30:01	ConstellationId.GPS	G26	1C	G261C	2.431713e+07	1.277874e+08	3446.324	31.0

Compute the ionospheric combination#

In this section we will compute the ionospheric (or geometry-free) combination (see equation (1)) for a specific satellite and data combination. To do this, we will use pandas merge function. The merge function allow us to join two DataFrames based on the contents of a column in an efficient, vectorized manner, without cumbersome for loops, which typically slow down Python code (and in general should be avoided)

First let’s check which are the different tracking channels available for a the previous satellite:

df_sat['channel'].unique()

array(['1C', '2W', '2X', '5X'], dtype=object)

To compute the ionospheric combination we need to pick two channels, for this example we will pick the observables generated with the C/A tracking loop at the L1 frequency (RINEX code 1C) and the encrypted code at the L2 (RINEX code 2W, usually obtained with proprietary techniques such as semi-codeless tracking)

ch_a = TrackingChannel.from_string('1C')
ch_b = TrackingChannel.from_string('2W')

# Select the channels that will be used for this recipe
df_sat_ch_a = df_sat[df_sat['channel'] == str(ch_a)]
df_sat_ch_b = df_sat[df_sat['channel'] == str(ch_b)]

# Check the results
df_sat_ch_a.head()

	epoch	constellation	sat	channel	signal	range	phase	doppler	snr
71704	2014-02-27 15:30:00	ConstellationId.GPS	G26	1C	G261C	2.431778e+07	1.277908e+08	3446.301	30.5
71787	2014-02-27 15:30:01	ConstellationId.GPS	G26	1C	G261C	2.431713e+07	1.277874e+08	3446.324	31.0
71870	2014-02-27 15:30:02	ConstellationId.GPS	G26	1C	G261C	2.431647e+07	1.277839e+08	3446.621	30.8
71953	2014-02-27 15:30:03	ConstellationId.GPS	G26	1C	G261C	2.431582e+07	1.277805e+08	3446.117	30.1
72036	2014-02-27 15:30:04	ConstellationId.GPS	G26	1C	G261C	2.431516e+07	1.277770e+08	3446.406	30.8

Now use the merge method to compute the LI combination

df_sat_li = pd.merge(df_sat_ch_a, df_sat_ch_b, on='epoch', how='inner', suffixes=('_a', '_b'))
df_sat_li.tail()

	epoch	constellation_a	sat_a	channel_a	signal_a	range_a	phase_a	doppler_a	snr_a	...	constellation_b	sat_b	channel_b	signal_b	range_b	phase_b	doppler_b	snr_b
7190	2014-02-27 19:59:55	ConstellationId.GPS	G26	1C	G261C	2.301577e+07	1.209488e+08	-1662.426	46.7	...	ConstellationId.GPS	G26	2W	G262W	2.301578e+07	9.424584e+07	-1295.403	31.5
7191	2014-02-27 19:59:56	ConstellationId.GPS	G26	1C	G261C	2.301609e+07	1.209505e+08	-1662.098	46.7	...	ConstellationId.GPS	G26	2W	G262W	2.301610e+07	9.424713e+07	-1295.284	31.1
7192	2014-02-27 19:59:57	ConstellationId.GPS	G26	1C	G261C	2.301640e+07	1.209521e+08	-1661.891	46.8	...	ConstellationId.GPS	G26	2W	G262W	2.301641e+07	9.424843e+07	-1295.184	31.6
7193	2014-02-27 19:59:58	ConstellationId.GPS	G26	1C	G261C	2.301672e+07	1.209538e+08	-1662.270	46.6	...	ConstellationId.GPS	G26	2W	G262W	2.301673e+07	9.424972e+07	-1295.126	31.8
7194	2014-02-27 19:59:59	ConstellationId.GPS	G26	1C	G261C	2.301703e+07	1.209555e+08	-1661.457	46.6	...	ConstellationId.GPS	G26	2W	G262W	2.301705e+07	9.425102e+07	-1294.960	31.8

5 rows × 21 columns

The RINEX format specifies the carrier phase as cycles, and therefore to compute the ionospheric combination we will need to convert the data to meters. To do this we first need to compute the wavelength of the data for each tracking channel

# Get the wavelength for each tracking channel
wl_a = ch_a.get_wavelength(sat.constellation)
wl_b = ch_b.get_wavelength(sat.constellation)

Now computing the ionospheric combination is straightforward. We will create the LI combination as a new column into the DataFrame (to align the data with the epochs for later analysis).

df_sat_li['li'] = df_sat_li['phase_a'] * wl_a - df_sat_li['phase_b'] * wl_b

We can now plot the LI combination using matplotlib plotting functions:

# Close all previous figures
plt.close()

# Plot the LI against the time
plt.plot(df_sat_li['epoch'], df_sat_li['li'], '.')
plt.title(f'LI combination for {sat}')
plt.xlabel('Time [hour of 2014 doy 058]')
plt.ylabel('LI [m]')

Text(0, 0.5, 'LI [m]')

Compute \(\Delta STEC\)#

We now will create another cell that contain the \(\Delta STEC\), computed as the time derivative of the geometric free combination (\(LI\)). As shown before (see equation (1), besides the \(STEC\), the geometry free combination of phases contain the phase ambiguities and the uncalibrated phase biases, which vanish with the time differentiation thanks to its constant nature. Therefore, the only term that will accompany the \(\Delta STEC\) is the phase measurement noise (\(\sqrt{2}\cdot \varepsilon_L\)), which can be considered in the millimeter range.

To compute the time difference of the geometry free combination, we can simply use the numpy.diff function, which, again, will save us from using for loops thanks to its vectorized nature.

d_li = np.diff(df_sat_li['li'])

We could now add this new time series into the dataframe. However, we must be careful due to a mismatch in the number of values: because of the difference operator, there is one less value in the d_li array, therefore, we need to add a sample manually.

# Number of samples in the difference time series
len(d_li)

# Number of samples in the original geometry combination time series
len(df_sat_li['li'])

We can insert (in this case prepend) a new sample in a numpy array using the insert method

d_li = np.insert(d_li, 0, np.nan, axis=0)
len(d_li)

Now that we have homogenized sizes, we can create the column, but first we will perform some data editing, discarding (setting to NaN) all those LI values larger than a threshold (due to time gaps and cycle slips in the data). We can do this first by simple array indexing, to find those values larger than a threshold and then assign those rows to NaN.

In order to have an idea on which are the elements to remove, we can plot the histogram to have an idea on the distribution of samples and have an idea on where to place the thresholds

plt.close()
_ = plt.hist(d_li)

From the histogram, it seems that samples below 2.5 could be safely discarded

li_thresold = 1.0

d_li[d_li > +li_thresold] = np.nan
d_li[d_li < -li_thresold] = np.nan

The \(\Delta STEC\) is simply the \(\Delta LI\) scaled by the \(\alpha\) factor, that depends on the frequencies used to compute the ionospheric combination

f_a = ch_a.get_frequency(sat.constellation)
f_b = ch_b.get_frequency(sat.constellation)

alpha = 40.3 / (f_b * f_b) - 40.3 / (f_a * f_a)

df_sat_li['d_stec_tecu_per_s'] = d_li / alpha / 1.0e16

Now we can plot the STEC

plt.close()
plt.plot(df_sat_li['epoch'], df_sat_li['d_stec_tecu_per_s'] * 60, marker='.')
plt.ylabel('Rate of STEC (TECU/min)')
plt.xlabel('Epoch')
plt.title('Rate of STEC')

Text(0.5, 1.0, 'Rate of STEC')

Compute ROTI#

In order to compute the ROTI, we need to group the data in batches of a certain time period (e.g. 5 minutes) and compute the standard deviation for each of these batches. To do this we will use the resample method, but to use this we first need to set the epoch column as index of the DataFrame

# Set the epoch as DataFrame time index. The resample method will use
# this index as basis for computation
df_sat_li.set_index('epoch', inplace=True)

# Resample the dSTEC data every 5 minutes
df_sat_dstec_sampled = df_sat_li['d_stec_tecu_per_s'].resample('5T').std()
df_sat_dstec_sampled.head()

epoch
2014-02-27 15:30:00         NaN
2014-02-27 15:35:00    0.034045
2014-02-27 15:40:00    0.021636
2014-02-27 15:45:00    0.022799
2014-02-27 15:50:00    0.018959
Freq: 5T, Name: d_stec_tecu_per_s, dtype: float64

Now we can finally plot the data

# Convert the pandas index (DateTimeIndex) to seconds for plotting
datetime_array = df_sat_dstec_sampled.index.to_numpy()

# Calculate the seconds of the day as a vectorized operation
seconds_of_day = (datetime_array - datetime_array.astype('datetime64[D]')) 

# Convert the timedelta values to seconds
seconds_of_day = seconds_of_day / np.timedelta64(1, 's')

plt.close()
plt.ylim(0,65)
plt.xlim(0,86400)
plt.xticks(np.linspace(0,86400, 7))
plt.yticks(np.arange(0,65, step=5))
plt.grid()
plt.plot(seconds_of_day, df_sat_dstec_sampled * 60)
plt.ylabel('Rate of TEC Index (TECU/min)')
plt.xlabel('Seconds of day 2014, doy 58')
plt.title('Rate of TEC Index for SEY1 station, GPS PRN26')

Text(0.5, 1.0, 'Rate of TEC Index for SEY1 station, GPS PRN26')

ROTI with affordable receivers#

The process can be repeated using an affordable receiver. In this case a MEDEA computer (based on the u-blox ZED-F9P and a Talysmann antenna)

rinex_file = '../assets/MTIC00ESP_R_20191221131_05H_01S_MO.rnx'

df = rinex.to_dataframe(rinex_file)

After loading the RINEX, let’s perform some basic data editing to flag phase breaks due to time gaps

df = mark_time_gap(df)

Selection of the channels for which the ionospheric combination needs to be performed

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

Once the channels have been selected, the ionospheric (or geometry-free) combination can be computed

df_li_gps = compute_geometry_free(df, constellation, channel_a, channel_b)

# Preview the LI values
df_li_gps.head()

	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	2019-05-02 11:31:51.400	ConstellationId.GPS	G02	1C	G021C	2.073340e+07	1.089547e+08	False	ConstellationId.GPS	2L	G022L	NaN	NaN	False	NaN	NaN
1	2019-05-02 11:31:51.400	ConstellationId.GPS	G06	1C	G061C	2.165911e+07	1.138194e+08	False	ConstellationId.GPS	2L	G062L	NaN	NaN	False	NaN	NaN
2	2019-05-02 11:31:51.400	ConstellationId.GPS	G12	1C	G121C	1.844391e+07	9.692342e+07	False	ConstellationId.GPS	2L	G122L	NaN	NaN	False	NaN	NaN
3	2019-05-02 11:31:51.400	ConstellationId.GPS	G14	1C	G141C	2.002953e+07	1.052559e+08	False	ConstellationId.GPS	2L	G142L	NaN	NaN	False	NaN	NaN
4	2019-05-02 11:31:51.400	ConstellationId.GPS	G24	1C	G241C	2.019349e+07	1.061175e+08	False	ConstellationId.GPS	2L	G242L	NaN	NaN	False	NaN	NaN

To compute the Slant Total Electron Content, we will need the \(\alpha_{LI}\) coefficient (that transforms LI to Slant Total Electron Content), which can be computed using the frequency associated to the channel bands:

f_a = channel_a.get_frequency(constellation)
f_b = channel_b.get_frequency(constellation)

alpha = 40.3 / (f_b * f_b) - 40.3 / (f_a * f_a)

sat = Satellite.from_string('G29')
df_sat = df_li_gps[df_li_gps['sat'] == sat]

t = compute_decimal_hours(df_sat['epoch'])
plt.close()
plt.title(f"LI combination for {sat}")
plt.plot(t, df_sat['li_m'], '.')
plt.xlabel(f"Time [ hour of {df_sat.iloc[0]['epoch'].date()} ]")
plt.ylabel("LI [m]")

Text(0, 0.5, 'LI [m]')

Once the LI combination is obtained, we are in the position of computing the \(\Delta STEC\) as follows:

df_li_gps['d_stec_tecu_per_s'] = df_li_gps.groupby('sat')['li_m'].diff() / alpha / 1.0e16

df_sat = df_li_gps[df_li_gps['sat'] == sat]

t = compute_decimal_hours(df_sat['epoch'])
plt.close()
plt.title(f"Time difference of STEC for {sat}")
plt.plot(t, df_sat['d_stec_tecu_per_s'], '.')
plt.xlabel(f"Time [ hour of {df_sat.iloc[0]['epoch'].date()} ]")
plt.ylabel("DSTEC [TECU/s]")
plt.ylim(-1, 1)

(-1.0, 1.0)

Once the time difference of the STEC has been computed we can now proceed to compute the standard deviation \(\sigma\) for intervals of 5 minutes on a per-satellite basis (i.e. definition of ROTI)

# Set 'epoch' as the DataFrame index
df_li_gps.set_index('epoch', inplace=True)

# Group the samples in 5 minute intervals and compute the sigma (i.e. ROTI)
df_roti = df_li_gps.groupby('sat').resample('5T').std(numeric_only=True).reset_index()

Plot the ROTI for a satellite

df_sat = df_roti[df_roti['sat'] == sat]

plt.close()

t = compute_decimal_hours(df_sat['epoch'])

plt.close()
plt.ylim(0,65)
plt.yticks(np.arange(0,65, step=5))
plt.grid()
plt.ylabel('Rate of TEC Index (TECU/min)')
plt.xlabel(f"Time [ hour of {df_sat.iloc[0]['epoch'].date()} ]")
plt.title(f'ROTI for affordable receiver, satellite {sat}')
plt.plot(t, df_sat['d_stec_tecu_per_s'] *60)

[<matplotlib.lines.Line2D at 0x7f012d67dc90>]

In this particular case, no scintillation event was detected (also likely due to the fact that the take was performed in mid latitude, in a quiet ionospheric period). Note the peaks observed in the data. These are artifacts due to cycle slips.

Higher ROTI values have been also detected in other cases at low elevations. Care must be exercised in affordable receivers (i.e. atennas) since low elevations might include multipath. Therefore, a conservative elevation mask (e.g. \(>15^\circ\)) is recommended when processing affordable receivers