Cleaned up some files

Author: kanhereashwin

gnss_lib_py/parsers/rinex_nav.py

@@ -135,7 +135,6 @@ def preprocess(self, rinex_paths, satellites):
 
         data = pd.DataFrame()
         self.iono_params = {}
-        #TODO: Convert this to a dictionary
         for rinex_path in rinex_paths:
             new_data, rinex_header = self._get_ephemeris_dataframe(rinex_path,
                                                                    constellations)
@@ -242,7 +241,6 @@ def _get_ephemeris_dataframe(self, rinex_path, constellations=None):
         data['t_oc']  = 1e-9 * data['t_oc'] - consts.WEEKSEC * np.floor(1e-9 * data['t_oc'] / consts.WEEKSEC)
         data['time'] = data['time'].dt.tz_localize('UTC')
         # Rename Keplerian orbital parameters to match a GLP standard
-        # TODO: Change these to a standard nomenclature
         data.rename(columns={'M0': 'M_0', 'Eccentricity': 'e', 'Toe': 't_oe', 'DeltaN': 'deltaN', 'Cuc': 'C_uc', 'Cus': 'C_us',
                              'Cic': 'C_ic', 'Crc': 'C_rc', 'Cis': 'C_is', 'Crs': 'C_rs', 'Io': 'i_0', 'Omega0': 'Omega_0'}, inplace=True)
         data.rename(columns={'X': 'sv_x_m', 'dX': 'sv_dx_mps', 'dX2': 'sv_dx2_mps2',
@@ -362,211 +360,6 @@ def load_leapseconds(self, rinex_header):
                 leap_seconds = np.nan
         return leap_seconds
 
-# def rinex_to_sv_states(gps_millis, rinex_nav):
-#     # keplerian_rinex =
-#     # num_int_rinex =
-#     # keplerian_sv_states =
-#     # num_int_sv_states =
-#     # sv_states = keplerian_sv_states.concat(num_int_sv_states)
-#     # sv_states.sort('gps_millis', inplace=True)
-#     # rotate the states
-#     # return sv_states
-#     pass
-#
-# def orbit_params_to_sv_states(gps_millis, ephem_orbit_params):
-#     """Compute position and velocities for all satellites in ephemeris file
-#     given time of clock.
-#
-#     `ephem` contains broadcast ephemeris parameters (similar in form to GPS
-#     broadcast parameters).
-#
-#     Must contain the following rows (description in [1]_):
-#     * :code:`gnss_id`
-#     * :code:`sv_id`
-#     * :code:`gps_week`
-#     * :code:`t_oe`
-#     * :code:`e`
-#     * :code:`omega`
-#     * :code:`Omega_0`
-#     * :code:`OmegaDot`
-#     * :code:`sqrtA`
-#     * :code:`deltaN`
-#     * :code:`IDOT`
-#     * :code:`i_0`
-#     * :code:`C_is`
-#     * :code:`C_ic`
-#     * :code:`C_rs`
-#     * :code:`C_rc`
-#     * :code:`C_uc`
-#     * :code:`C_us`
-#
-#     Parameters
-#     ----------
-#     gps_millis : int
-#         Time at which measurements are needed, measured in milliseconds
-#         since start of GPS epoch [ms].
-#     ephem : gnss_lib_py.parsers.navdata.NavData
-#         NavData instance containing ephemeris parameters of satellites
-#         for which states are required.
-#
-#     Returns
-#     -------
-#     sv_posvel : gnss_lib_py.parsers.navdata.NavData
-#         NavData containing satellite positions, velocities, corresponding
-#         time with GNSS ID and SV number.
-#
-#     Notes
-#     -----
-#     Based on code written by J. Makela.
-#     AE 456, Global Navigation Sat Systems, University of Illinois
-#     Urbana-Champaign. Fall 2017
-#
-#     More details on the algorithm used to compute satellite positions
-#     from broadcast navigation message can be found in [1]_.
-#
-#     Satellite velocity calculations based on algorithms introduced in [2]_.
-#
-#     References
-#     ----------
-#     ..  [1] Misra, P. and Enge, P,
-#         "Global Positioning System: Signals, Measurements, and Performance."
-#         2nd Edition, Ganga-Jamuna Press, 2006.
-#     ..  [2] B. F. Thompson, S. W. Lewis, S. A. Brown, and T. M. Scott,
-#         “Computing GPS satellite velocity and acceleration from the broadcast
-#         navigation message,” NAVIGATION, vol. 66, no. 4, pp. 769–779, Dec. 2019,
-#         doi: 10.1002/navi.342.
-#
-#     """
-#
-#     # Convert time from GPS millis to TOW
-#     gps_week, gps_tow = gps_millis_to_tow(gps_millis)
-#     # Extract parameters
-#
-#     c_is = ephem_orbit_params['C_is']
-#     c_ic = ephem_orbit_params['C_ic']
-#     c_rs = ephem_orbit_params['C_rs']
-#     c_rc = ephem_orbit_params['C_rc']
-#     c_uc = ephem_orbit_params['C_uc']
-#     c_us = ephem_orbit_params['C_us']
-#     delta_n   = ephem_orbit_params['deltaN']
-#
-#     ecc        = ephem_orbit_params['e']     # eccentricity
-#     omega    = ephem_orbit_params['omega'] # argument of perigee
-#     omega_0  = ephem_orbit_params['Omega_0']
-#     sqrt_sma = ephem_orbit_params['sqrtA'] # sqrt of semi-major axis
-#     sma      = sqrt_sma**2      # semi-major axis
-#
-#     sqrt_mu_a = np.sqrt(consts.MU_EARTH) * sqrt_sma**-3 # mean angular motion
-#     gpsweek_diff = (np.mod(gps_week,1024) - np.mod(ephem_orbit_params['gps_week'],1024))*604800.
-#     sv_posvel = NavData()
-#     sv_posvel['gnss_id'] = ephem_orbit_params['gnss_id']
-#     sv_posvel['sv_id'] = ephem_orbit_params['sv_id']
-#     sv_posvel['gps_millis'] = gps_millis
-#
-#     delta_t = gps_tow - ephem_orbit_params['t_oe'] + gpsweek_diff
-#
-#     # Calculate the mean anomaly with corrections
-#     ecc_anom = _compute_eccentric_anomaly(gps_week, gps_tow,
-#                                             ephem_orbit_params)
-#
-#     cos_e   = np.cos(ecc_anom)
-#     sin_e   = np.sin(ecc_anom)
-#     e_cos_e = (1 - ecc*cos_e)
-#
-#     # Calculate the true anomaly from the eccentric anomaly
-#     sin_nu = np.sqrt(1 - ecc**2) * (sin_e/e_cos_e)
-#     cos_nu = (cos_e-ecc) / e_cos_e
-#     nu_rad     = np.arctan2(sin_nu, cos_nu)
-#
-#     # Calcualte the argument of latitude iteratively
-#     phi_0 = nu_rad + omega
-#     phi   = phi_0
-#     for incl in range(5):
-#         cos_to_phi = np.cos(2.*phi)
-#         sin_to_phi = np.sin(2.*phi)
-#         phi_corr = c_uc * cos_to_phi + c_us * sin_to_phi
-#         phi = phi_0 + phi_corr
-#
-#     # Calculate the longitude of ascending node with correction
-#     omega_corr = ephem_orbit_params['OmegaDot'] * delta_t
-#
-#     # Also correct for the rotation since the beginning of the GPS week for which the Omega0 is
-#     # defined.  Correct for GPS week rollovers.
-#
-#     # Also correct for the rotation since the beginning of the GPS week for
-#     # which the Omega0 is defined.  Correct for GPS week rollovers.
-#     omega = omega_0 - (consts.OMEGA_E_DOT*(gps_tow + gpsweek_diff)) + omega_corr
-#
-#     # Calculate orbital radius with correction
-#     r_corr = c_rc * cos_to_phi + c_rs * sin_to_phi
-#     orb_radius      = sma*e_cos_e + r_corr
-#
-#     ############################################
-#     ######  Lines added for velocity (1)  ######
-#     ############################################
-#     delta_e   = (sqrt_mu_a + delta_n) / e_cos_e
-#     dphi = np.sqrt(1 - ecc**2)*delta_e / e_cos_e
-#     # Changed from the paper
-#     delta_r   = (sma * ecc * delta_e * sin_e) + 2*(c_rs*cos_to_phi - c_rc*sin_to_phi)*dphi
-#
-#     # Calculate the inclination with correction
-#     i_corr = c_ic*cos_to_phi + c_is*sin_to_phi + ephem_orbit_params['IDOT']*delta_t
-#     incl = ephem_orbit_params['i_0'] + i_corr
-#
-#     ############################################
-#     ######  Lines added for velocity (2)  ######
-#     ############################################
-#     delta_i = 2*(c_is*cos_to_phi - c_ic*sin_to_phi)*dphi + ephem_orbit_params['IDOT']
-#
-#     # Find the position in the orbital plane
-#     x_plane = orb_radius*np.cos(phi)
-#     y_plane = orb_radius*np.sin(phi)
-#
-#     ############################################
-#     ######  Lines added for velocity (3)  ######
-#     ############################################
-#     delta_u = (1 + 2*(c_us * cos_to_phi - c_uc*sin_to_phi))*dphi
-#     dxp = delta_r*np.cos(phi) - orb_radius*np.sin(phi)*delta_u
-#     dyp = delta_r*np.sin(phi) + orb_radius*np.cos(phi)*delta_u
-#     # Find satellite position in ECEF coordinates
-#     cos_omega = np.cos(omega)
-#     sin_omega = np.sin(omega)
-#     cos_i = np.cos(incl)
-#     sin_i = np.sin(incl)
-#
-#     sv_posvel['x_sv_m'] = x_plane*cos_omega - y_plane*cos_i*sin_omega
-#     sv_posvel['y_sv_m'] = x_plane*sin_omega + y_plane*cos_i*cos_omega
-#     sv_posvel['z_sv_m'] = y_plane*sin_i
-#
-#     ############################################
-#     ######  Lines added for velocity (4)  ######
-#     ############################################
-#     omega_dot = ephem_orbit_params['OmegaDot'] - consts.OMEGA_E_DOT
-#     sv_posvel['vx_sv_mps'] = (dxp * cos_omega
-#                         - dyp * cos_i*sin_omega
-#                         + y_plane  * sin_omega*sin_i*delta_i
-#                         - (x_plane * sin_omega + y_plane*cos_i*cos_omega)*omega_dot)
-#
-#     sv_posvel['vy_sv_mps'] = (dxp * sin_omega
-#                         + dyp * cos_i * cos_omega
-#                         - y_plane  * sin_i * cos_omega * delta_i
-#                         + (x_plane * cos_omega - (y_plane*cos_i*sin_omega)) * omega_dot)
-#
-#     sv_posvel['vz_sv_mps'] = dyp*sin_i + y_plane*cos_i*delta_i
-#
-#     # Estimate SV clock corrections, including polynomial and relativistic
-#     # clock corrections
-#     clock_corr, _, _ = _estimate_sv_clock_corr(gps_millis,
-#                                                 ephem_orbit_params)
-#
-#     sv_posvel['b_sv_m'] = clock_corr
-#
-#     return sv_posvel
-#
-#
-# def numerical_integration_for_sv_states(self, gps_millis):
-#     pass
-
 
 def _compute_eccentric_anomaly(gps_week, gps_tow, ephem, tol=1e-5, max_iter=10):
     """Compute the eccentric anomaly from ephemeris parameters.

gnss_lib_py/utils/gnss_models.py

@@ -21,7 +21,6 @@
                         find_sv_location, find_sv_states, \
                         find_visible_sv_posvel, _sort_ephem_measures, \
                         _filter_ephemeris_measurements
-from gnss_lib_py.parsers.rinex_nav import _compute_eccentric_anomaly
 from gnss_lib_py.utils.ephemeris_downloader import DEFAULT_EPHEM_PATH
 
 
@@ -339,16 +338,6 @@ def expected_measures(gps_millis, state, ephem=None, sv_posvel=None):
     return measurements, sv_posvel
 
 
-# def _matching_state_for_measure(measurements, state_estimate,
-#                                 delta_t_decimals=2):
-#     measurement_times = np.sort(np.unique(np.around(times,
-#                                          decimals=delta_t_decimals)))
-#     for
-
-
-#     return measure_time_state_row
-
-
 def _extract_state_variables(state):
     """Extract position, velocity and clock bias terms from state.
 
@@ -450,14 +439,6 @@ def calculate_pseudorange_corr(gps_millis, state=None, ephem=None, sv_posvel=Non
                 "SV states must be given when ephemeris isn't"
         satellites = len(sv_posvel)
 
-    # # calculate clock pseudorange correction
-    # if ephem is not None:
-    #     clock_corr, _, _ = _calculate_clock_delay(gps_millis, ephem)
-    # else:
-    #     warnings.warn("Broadcast ephemeris not given, returning 0 "\
-    #                 + "clock corrections", RuntimeWarning)
-    #     clock_corr = np.zeros(satellites)
-
     if rx_ecef is not None:
         # Calculate the tropospheric delays
         tropo_delay = _calculate_tropo_delay(gps_millis, rx_ecef, ephem, sv_posvel)
@@ -475,69 +456,9 @@ def calculate_pseudorange_corr(gps_millis, state=None, ephem=None, sv_posvel=Non
                         RuntimeWarning)
         iono_delay = np.zeros(satellites)
 
-    # return clock_corr, tropo_delay, iono_delay
     return tropo_delay, iono_delay
 
 
-# def _calculate_clock_delay(gps_millis, ephem):
-#     """Calculate the modelled satellite clock delay
-
-#     Parameters
-#     ---------
-#     gps_millis : int
-#         Time at which measurements are needed, measured in milliseconds
-#         since start of GPS epoch [ms].
-#     ephem : gnss_lib_py.parsers.navdata.NavData
-#         Satellite ephemeris parameters for measurement SVs.
-
-#     Returns
-#     -------
-#     clock_corr : np.ndarray
-#         Satellite clock corrections containing all terms [m].
-#     corr_polynomial : np.ndarray
-#         Polynomial clock perturbation terms [m].
-#     clock_relativistic : np.ndarray
-#         Relativistic clock correction terms [m].
-
-#     """
-#     # Extract required GPS constants
-#     ecc        = ephem['e']     # eccentricity
-#     sqrt_sma = ephem['sqrtA'] # sqrt of semi-major axis
-
-#     # if np.abs(delta_t).any() > 302400:
-#     #     delta_t = delta_t - np.sign(delta_t)*604800
-
-#     gps_week, gps_tow = gps_millis_to_tow(gps_millis)
-
-#     # Compute Eccentric Anomaly
-#     ecc_anom = _compute_eccentric_anomaly(gps_week, gps_tow, ephem)
-
-#     # Determine pseudorange corrections due to satellite clock corrections.
-#     # Calculate time offset from satellite reference time
-#     t_offset = gps_tow - ephem['t_oc']
-#     if np.abs(t_offset).any() > 302400:  # pragma: no cover
-#         t_offset = t_offset-np.sign(t_offset)*604800
-
-#     # Calculate clock corrections from the polynomial corrections in
-#     # broadcast message
-#     corr_polynomial = (ephem['SVclockBias']
-#                      + ephem['SVclockDrift']*t_offset
-#                      + ephem['SVclockDriftRate']*t_offset**2)
-
-#     # Calcualte the relativistic clock correction
-#     corr_relativistic = consts.F * ecc * sqrt_sma * np.sin(ecc_anom)
-
-#     # Calculate the total clock correction including the Tgd term
-#     clk_corr = (corr_polynomial - ephem['TGD'] + corr_relativistic)
-
-#     #Convert values to equivalent meters from seconds
-#     clk_corr = np.array(consts.C*clk_corr, ndmin=1)
-#     corr_polynomial = np.array(consts.C*corr_polynomial, ndmin=1)
-#     corr_relativistic = np.array(consts.C*corr_relativistic, ndmin=1)
-
-#     return clk_corr, corr_polynomial, corr_relativistic
-
-
 def _calculate_tropo_delay(gps_millis, rx_ecef, ephem=None, sv_posvel=None):
     """Calculate tropospheric delay