{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":60095,"databundleVersionId":6542333,"sourceType":"competition"}],"dockerImageVersionId":30673,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-03-29T09:24:52.576457Z","iopub.execute_input":"2024-03-29T09:24:52.577012Z","iopub.status.idle":"2024-03-29T09:24:53.450616Z","shell.execute_reply.started":"2024-03-29T09:24:52.576965Z","shell.execute_reply":"2024-03-29T09:24:53.449551Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install numpy pandas pymap3d matplotlib tqdm scipy","metadata":{"execution":{"iopub.status.busy":"2024-03-29T09:24:53.452032Z","iopub.execute_input":"2024-03-29T09:24:53.452579Z","iopub.status.idle":"2024-03-29T09:27:23.28976Z","shell.execute_reply.started":"2024-03-29T09:24:53.452549Z","shell.execute_reply":"2024-03-29T09:27:23.288234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install pymap3d\n\nimport numpy as np\nimport pandas as pd\nimport pymap3d as pm\nimport pymap3d.vincenty as pmv\nimport matplotlib.pyplot as plt\nimport glob as gl\nimport scipy.optimize\nfrom tqdm.auto import tqdm\nfrom scipy.interpolate import InterpolatedUnivariateSpline\nfrom scipy.spatial import distance\n# Constants\nCLIGHT = 299_792_458   # speed of light (m/s)\nRE_WGS84 = 6_378_137   # earth semimajor axis (WGS84) (m)\nOMGE = 7.2921151467E-5  # earth angular velocity (IS-GPS) (rad/s)\n# Satellite selection using carrier frequency error, elevation angle, and C/N0\ndef satellite_selection(df, column):\n    idx = df[column].notnull()\n    idx &= df['CarrierErrorHz'] < 2.0e6  # carrier frequency error (Hz)\n    idx &= df['SvElevationDegrees'] > 10.0  # elevation angle (deg)\n    idx &= df['Cn0DbHz'] > 15.0  # C/N0 (dB-Hz)\n    idx &= df['MultipathIndicator'] == 0 # Multipath flag\n\n    return df[idx]\n# Compute line-of-sight vector from user to satellite\ndef los_vector(xusr, xsat):\n    \n    u = xsat - xusr\n    rng = np.linalg.norm(u, axis=1).reshape(-1, 1)\n    u /= rng\n    \n    return u, rng.reshape(-1)\n\n# Compute Jacobian matrix\ndef jac_pr_residuals(x, xsat, pr, W):\n    \n    u, _ = los_vector(x[:3], xsat)\n    J = np.hstack([-u, np.ones([len(pr), 1])])  # J = [-ux -uy -uz 1]\n\n    return W @ J\n\n# Compute pseudorange residuals\ndef pr_residuals(x, xsat, pr, W):\n    \n    u, rng = los_vector(x[:3], xsat)\n\n    # Approximate correction of the earth rotation (Sagnac effect) often used in GNSS positioning\n    rng += OMGE * (xsat[:, 0] * x[1] - xsat[:, 1] * x[0]) / CLIGHT\n\n    # Add GPS L1 clock offset\n    residuals = rng - (pr - x[3])\n\n    return residuals @ W\n\n# Compute Jacobian matrix\ndef jac_prr_residuals(v, vsat, prr, x, xsat, W):\n  \n    u, _ = los_vector(x[:3], xsat)\n    J = np.hstack([-u, np.ones([len(prr), 1])])\n\n    return np.dot(W, J)\n\n# Compute pseudorange rate residuals\ndef prr_residuals(v, vsat, prr, x, xsat, W):\n  \n    u, rng = los_vector(x[:3], xsat)\n    rate = np.sum((vsat-v[:3])*u, axis=1) \\\n          + OMGE / CLIGHT * (vsat[:, 1] * x[0] + xsat[:, 1] * v[0]\n                           - vsat[:, 0] * x[1] - xsat[:, 0] * v[1])\n\n    residuals = rate - (prr - v[3])\n\n    return residuals @ W\n# Carrier smoothing of pseudarange\ndef carrier_smoothing(gnss_df):\n  \n    carr_th = 1.0# carrier phase jump threshold [m] 2->1.5 (best)->1.0\n    pr_th =  15.0 # pseudorange jump threshold [m] 20->15\n\n    prsmooth = np.full_like(gnss_df['RawPseudorangeMeters'], np.nan)\n    # Loop for each signal\n    for (i, (svid_sigtype, df)) in enumerate((gnss_df.groupby(['Svid', 'SignalType']))):\n        df = df.replace(\n            {'AccumulatedDeltaRangeMeters': {0: np.nan}})  # 0 to NaN\n\n        # Compare time difference between pseudorange/carrier with Doppler\n        drng1 = df['AccumulatedDeltaRangeMeters'].diff() - df['PseudorangeRateMetersPerSecond']\n        drng2 = df['RawPseudorangeMeters'].diff() - df['PseudorangeRateMetersPerSecond']\n\n        # Check cycle-slip\n        slip1 = (df['AccumulatedDeltaRangeState'].to_numpy() & 2**1) != 0  # reset flag\n        slip2 = (df['AccumulatedDeltaRangeState'].to_numpy() & 2**2) != 0  # cycle-slip flag\n        slip3 = np.fabs(drng1.to_numpy()) > carr_th # Carrier phase jump\n        slip4 = np.fabs(drng2.to_numpy()) > pr_th # Pseudorange jump\n\n        idx_slip = slip1 | slip2 | slip3 | slip4\n        idx_slip[0] = True\n\n        # groups with continuous carrier phase tracking\n        df['group_slip'] = np.cumsum(idx_slip)\n\n        # Psudorange - carrier phase\n        df['dpc'] = df['RawPseudorangeMeters'] - df['AccumulatedDeltaRangeMeters']\n\n        # Absolute distance bias of carrier phase\n        meandpc = df.groupby('group_slip')['dpc'].mean()\n        df = df.merge(meandpc, on='group_slip', suffixes=('', '_Mean'))\n\n        # Index of original gnss_df\n        idx = (gnss_df['Svid'] == svid_sigtype[0]) & (\n            gnss_df['SignalType'] == svid_sigtype[1])\n\n        # Carrier phase + bias\n        prsmooth[idx] = df['AccumulatedDeltaRangeMeters'] + df['dpc_Mean']\n\n    # If carrier smoothing is not possible, use original pseudorange\n    idx_nan = np.isnan(prsmooth)\n    prsmooth[idx_nan] = gnss_df['RawPseudorangeMeters'][idx_nan]\n    gnss_df['pr_smooth'] = prsmooth\n\n    return gnss_df\n# Compute distance by Vincenty's formulae\ndef vincenty_distance(llh1, llh2):\n  \n    d, az = np.array(pmv.vdist(llh1[:, 0], llh1[:, 1], llh2[:, 0], llh2[:, 1]))\n\n    return d\n\n# Compute score\ndef calc_score(llh, llh_gt):\n  \n    d = vincenty_distance(llh, llh_gt)\n    score = np.mean([np.quantile(d, 0.50), np.quantile(d, 0.95)])\n\n    return score\n# GNSS single point positioning using pseudorange\n\ndef point_positioning(gnss_df):\n    # Add nominal frequency to each signal\n    # Note: GLONASS is an FDMA signal, so each satellite has a different frequency\n    CarrierFrequencyHzRef = gnss_df.groupby(['Svid', 'SignalType'])[\n        'CarrierFrequencyHz'].median()\n    gnss_df = gnss_df.merge(CarrierFrequencyHzRef, how='left', on=[\n                            'Svid', 'SignalType'], suffixes=('', 'Ref'))\n    gnss_df['CarrierErrorHz'] = np.abs(\n        (gnss_df['CarrierFrequencyHz'] - gnss_df['CarrierFrequencyHzRef']))\n\n    # Carrier smoothing\n    gnss_df = carrier_smoothing(gnss_df)\n\n    # GNSS single point positioning\n    utcTimeMillis = gnss_df['utcTimeMillis'].unique()\n    nepoch = len(utcTimeMillis)\n    x0 = np.zeros(4)  # [x,y,z,tGPSL1]\n    v0 = np.zeros(4)  # [vx,vy,vz,dtGPSL1]\n    x_wls = np.full([nepoch, 3], np.nan)  # For saving position\n    v_wls = np.full([nepoch, 3], np.nan)  # For saving velocity\n\n    # Loop for epochs\n    for i, (t_utc, df) in enumerate(tqdm(gnss_df.groupby('utcTimeMillis'), total=nepoch)):\n        # Valid satellite selection\n        df_pr = satellite_selection(df, 'pr_smooth')\n        df_prr = satellite_selection(df, 'PseudorangeRateMetersPerSecond')\n\n        # Corrected pseudorange/pseudorange rate\n        pr = (df_pr['pr_smooth'] + df_pr['SvClockBiasMeters'] - df_pr['IsrbMeters'] -\n              df_pr['IonosphericDelayMeters'] - df_pr['TroposphericDelayMeters']).to_numpy()\n        prr = (df_prr['PseudorangeRateMetersPerSecond'] +\n               df_prr['SvClockDriftMetersPerSecond']).to_numpy()\n\n        # Satellite position/velocity\n        xsat_pr = df_pr[['SvPositionXEcefMeters', 'SvPositionYEcefMeters',\n                         'SvPositionZEcefMeters']].to_numpy()\n        xsat_prr = df_prr[['SvPositionXEcefMeters', 'SvPositionYEcefMeters',\n                           'SvPositionZEcefMeters']].to_numpy()\n        vsat = df_prr[['SvVelocityXEcefMetersPerSecond', 'SvVelocityYEcefMetersPerSecond',\n                       'SvVelocityZEcefMetersPerSecond']].to_numpy()\n\n        # Weight matrix for peseudorange/pseudorange rate\n        Wx = np.diag(1 / df_pr['RawPseudorangeUncertaintyMeters'].to_numpy())\n        Wv = np.diag(1 / df_prr['PseudorangeRateUncertaintyMetersPerSecond'].to_numpy())\n\n        # Robust WLS requires accurate initial values for convergence,\n        # so perform normal WLS for the first time\n        if len(df_pr) >= 4:\n            # Normal WLS\n            if np.all(x0 == 0):\n                opt = scipy.optimize.least_squares(\n                    pr_residuals, x0, jac_pr_residuals, args=(xsat_pr, pr, Wx))\n                x0 = opt.x \n            # Robust WLS for position estimation\n            opt = scipy.optimize.least_squares(\n                 pr_residuals, x0, jac_pr_residuals, args=(xsat_pr, pr, Wx), loss='soft_l1')\n            if opt.status < 1 or opt.status == 2:\n                 print(f'i = {i} position lsq status = {opt.status}')\n            else:\n                 x_wls[i, :] = opt.x[:3]\n                 x0 = opt.x\n                 \n        # Velocity estimation\n        if len(df_prr) >= 4:\n            if np.all(v0 == 0): # Normal WLS\n                opt = scipy.optimize.least_squares(\n                    prr_residuals, v0, jac_prr_residuals, args=(vsat, prr, x0, xsat_prr, Wv))\n                v0 = opt.x\n            # Robust WLS for velocity estimation\n            opt = scipy.optimize.least_squares(\n                prr_residuals, v0, jac_prr_residuals, args=(vsat, prr, x0, xsat_prr, Wv), loss='soft_l1')\n            if opt.status < 1:\n                print(f'i = {i} velocity lsq status = {opt.status}')\n            else:\n                v_wls[i, :] = opt.x[:3]\n                v0 = opt.x\n\n    return utcTimeMillis, x_wls, v_wls\n# Simple outlier detection and interpolation\ndef exclude_interpolate_outlier(x_wls, v_wls):\n    # Up velocity threshold\n    v_up_th = 2.0 # m/s\n\n    # Coordinate conversion\n    x_llh = np.array(pm.ecef2geodetic(x_wls[:, 0], x_wls[:, 1], x_wls[:, 2])).T\n    v_enu = np.array(pm.ecef2enuv(\n        v_wls[:, 0], v_wls[:, 1], v_wls[:, 2], x_llh[0, 0], x_llh[0, 1])).T\n\n    # Up velocity jump detection\n    # Cars don't jump suddenly!\n    idx_v_out = np.abs(v_enu[:, 2]) > v_up_th\n    v_wls[idx_v_out, :] = np.nan\n    \n    # Interpolate NaNs at beginning and end of array\n    x_df = pd.DataFrame({'x': x_wls[:, 0], 'y': x_wls[:, 1], 'z': x_wls[:, 2]})\n    x_df = x_df.interpolate(limit_area='outside', limit_direction='both')\n    \n    # Interpolate all NaN data\n    v_df = pd.DataFrame({'x': v_wls[:, 0], 'y': v_wls[:, 1], 'z': v_wls[:, 2]})\n    v_df = v_df.interpolate(limit_area='outside', limit_direction='both')\n    v_df = v_df.interpolate('spline', order=3)\n\n    return x_df.to_numpy(), v_df.to_numpy()\n# Kalman filter\ndef Kalman_filter(zs, us, phone):\n    # Parameters\n    # I don't know why only XiaomiMi8 seems to be inaccurate ... \n    sigma_v = 0.6 if phone == 'XiaomiMi8' else 0.1 # velocity SD m/s\n    sigma_x = 5.0  # position SD m\n    sigma_mahalanobis = 30.0 # Mahalanobis distance for rejecting innovation\n    \n    n, dim_x = zs.shape\n    F = np.eye(3)  # Transition matrix\n    Q = sigma_v**2 * np.eye(3)  # Process noise\n\n    H = np.eye(3)  # Measurement function\n    R = sigma_x**2 * np.eye(3)  # Measurement noise\n\n    # Initial state and covariance\n    x = zs[0, :3].T  # State\n    P = sigma_x**2 * np.eye(3)  # State covariance\n    I = np.eye(dim_x)\n\n    x_kf = np.zeros([n, dim_x])\n    P_kf = np.zeros([n, dim_x, dim_x])\n\n    # Kalman filtering\n    for i, (u, z) in enumerate(zip(us, zs)):\n        # First step\n        if i == 0:\n            x_kf[i] = x.T\n            P_kf[i] = P\n            continue\n\n        # Prediction step\n        x = F @ x + u.T\n        P = (F @ P) @ F.T + Q\n\n        # Check outliers for observation\n        d = distance.mahalanobis(z, H @ x, np.linalg.pinv(P))\n\n        # Update step\n        if d < sigma_mahalanobis:\n            y = z.T - H @ x\n            S = (H @ P) @ H.T + R\n            K = (P @ H.T) @ np.linalg.inv(S)\n            x = x + K @ y\n            P = (I - (K @ H)) @ P\n        else:\n            # If no observation update is available, increase covariance\n            P += 10**2*Q\n\n        x_kf[i] = x.T\n        P_kf[i] = P\n\n    return x_kf, P_kf\n\n# Forward + backward Kalman filter and smoothing\ndef Kalman_smoothing(x_wls, v_wls, phone):\n    n, dim_x = x_wls.shape\n\n    # Forward\n    v = np.vstack([np.zeros([1, 3]), (v_wls[:-1, :] + v_wls[1:, :])/2])\n    x_f, P_f = Kalman_filter(x_wls, v, phone)\n\n    # Backward\n    v = -np.flipud(v_wls)\n    v = np.vstack([np.zeros([1, 3]), (v[:-1, :] + v[1:, :])/2])\n    x_b, P_b = Kalman_filter(np.flipud(x_wls), v, phone)\n\n    # Smoothing\n    x_fb = np.zeros_like(x_f)\n    P_fb = np.zeros_like(P_f)\n    for (f, b) in zip(range(n), range(n-1, -1, -1)):\n        P_fi = np.linalg.inv(P_f[f])\n        P_bi = np.linalg.inv(P_b[b])\n\n        P_fb[f] = np.linalg.inv(P_fi + P_bi)\n        x_fb[f] = P_fb[f] @ (P_fi @ x_f[f] + P_bi @ x_b[b])\n\n    return x_fb, x_f, np.flipud(x_b)\n# Target course/phone\npath = '/kaggle/input/smartphone-decimeter-2023/sdc2023/train/2023-09-07-22-48-us-ca-routebc2/pixel4xl'\n\ndrive, phone = path.split('/')[-2:]\n\n# Read data\ngnss_df = pd.read_csv(f'{path}/device_gnss.csv')  # GNSS data\ngt_df = pd.read_csv(f'{path}/ground_truth.csv')  # ground truth\n\n# Point positioning\nutc, x_wls, v_wls = point_positioning(gnss_df)\n\n# Exclude velocity outliers\nx_wls, v_wls = exclude_interpolate_outlier(x_wls, v_wls)\n\n# Kalman smoothing\nx_kf, _, _ = Kalman_smoothing(x_wls, v_wls, phone)\n\n# Convert to latitude and longitude\nllh_wls = np.array(pm.ecef2geodetic(x_wls[:, 0], x_wls[:, 1], x_wls[:, 2])).T\nllh_kf = np.array(pm.ecef2geodetic(x_kf[:, 0], x_kf[:, 1], x_kf[:, 2])).T\n\n# Baseline\nx_bl = gnss_df.groupby('TimeNanos')[\n    ['WlsPositionXEcefMeters', 'WlsPositionYEcefMeters', 'WlsPositionZEcefMeters']].mean().to_numpy()\nllh_bl = np.array(pm.ecef2geodetic(x_bl[:, 0], x_bl[:, 1], x_bl[:, 2])).T\n\n# Ground truth\nllh_gt = gt_df[['LatitudeDegrees', 'LongitudeDegrees']].to_numpy()\n\n# Distance from ground truth\nvd_bl = vincenty_distance(llh_bl, llh_gt)\nvd_wls = vincenty_distance(llh_wls, llh_gt)\nvd_kf = vincenty_distance(llh_kf, llh_gt)\n\n# Score\nscore_bl = calc_score(llh_bl, llh_gt)\nscore_wls = calc_score(llh_wls, llh_gt)\nscore_kf = calc_score(llh_kf[:-1, :], llh_gt[:-1, :])\n\nprint(f'Score Baseline   {score_bl:.4f} [m]')\nprint(f'Score Robust WLS {score_wls:.4f} [m]')\nprint(f'Score KF         {score_kf:.4f} [m]')\n\n# Plot distance error\nplt.figure()\nplt.title('Distance error')\nplt.ylabel('Distance error [m]')\nplt.plot(vd_bl, label=f'Baseline, Score: {score_bl:.4f} m')\nplt.plot(vd_wls, label=f'Robust WLS, Score: {score_wls:.4f} m')\nplt.plot(vd_kf, label=f'Robust WLS + KF, Score: {score_kf:.4f} m')\nplt.legend()\nplt.grid()\nplt.ylim([0, 30])\n\n# Compute velocity error\nspeed_wls = np.linalg.norm(v_wls[:, :3], axis=1)\nspeed_gt = gt_df['SpeedMps'].to_numpy()\nspeed_rmse = np.sqrt(np.sum((speed_wls-speed_gt)**2)/len(speed_gt))\n\n# Plot velocity error\nplt.figure()\nplt.title('Speed error')\nplt.ylabel('Speed Error [m/s]')\nplt.plot(speed_wls - speed_gt, label=f'Speed RMSE: {speed_rmse:.4f} m')\nplt.legend()\nplt.grid()\npath = '/kaggle/input/smartphone-decimeter-2023/sdc2023'\nsample_df = pd.read_csv(f'{path}/sample_submission.csv')\ntest_dfs = []\n\n# Loop for each trip\nfor i, dirname in enumerate(tqdm(sorted(gl.glob(f'{path}/test/*/*/')))):\n    drive, phone = dirname.split('/')[-3:-1]\n    tripID = f'{drive}/{phone}'\n    print(tripID)\n\n    # Read data\n    gnss_df = pd.read_csv(f'{dirname}/device_gnss.csv')\n\n    # Point positioning\n    utc, x_wls, v_wls = point_positioning(gnss_df)\n\n    # Exclude velocity outliers\n    x_wls, v_wls = exclude_interpolate_outlier(x_wls, v_wls)\n\n    # Kalman smoothing\n    x_kf, _, _ = Kalman_smoothing(x_wls, v_wls, phone)\n\n    # Convert to latitude and longitude\n    llh_kf = np.array(pm.ecef2geodetic(x_kf[:, 0], x_kf[:, 1], x_kf[:, 2])).T\n\n    # Interpolation for submission\n    UnixTimeMillis = sample_df[sample_df['tripId'] == tripID]['UnixTimeMillis'].to_numpy()\n    lat = InterpolatedUnivariateSpline(utc, llh_kf[:,0], ext=3)(UnixTimeMillis)\n    lng = InterpolatedUnivariateSpline(utc, llh_kf[:,1], ext=3)(UnixTimeMillis)\n    trip_df = pd.DataFrame({\n        'tripId' : tripID,\n        'UnixTimeMillis': UnixTimeMillis,\n        'LatitudeDegrees': lat,\n        'LongitudeDegrees': lng\n        })\n\n    test_dfs.append(trip_df)\n\n# Write submission.csv\ntest_df = pd.concat(test_dfs)\ntest_df.to_csv('submission.csv', index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Write submission.csv\ntest_df = pd.concat(test_dfs)\ntest_df.to_csv('submission.csv', index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To work on this GNSS (Global Navigation Satellite System) data processing and position estimation project step by step, follow these guidelines:\n\n1. **Understanding the Project Scope**:\n   - Familiarize yourself with the project goal, which is to process GNSS data and estimate the position of a device accurately.\n\n2. **Setting Up Environment**:\n   - Ensure you have Python installed on your system.\n   - Install the required libraries mentioned in the code (`pymap3d`, `numpy`, `pandas`, etc.). You can use `pip install [library_name]`.\n\n3. **Data Acquisition**:\n   - Gather the GNSS data and ground truth data for your specific scenario. Ensure the data is in the expected format (e.g., CSV files).\n\n4. **Data Loading**:\n   - Modify the code to load your GNSS data and ground truth data. Ensure the file paths are correct.\n\n5. **Data Preprocessing**:\n   - Check the structure of your data and modify the code accordingly if there are any differences in column names or data format.\n   - Ensure that the GNSS data contains necessary information like pseudorange, pseudorange rate, satellite positions, etc.\n\n6. **Point Positioning**:\n   - Understand how the `point_positioning` function works. This function performs GNSS single point positioning using pseudorange data.\n   - Execute the `point_positioning` function with your GNSS data to estimate positions.\n\n7. **Outlier Detection and Interpolation**:\n   - Understand the `exclude_interpolate_outlier` function. This function detects and interpolates outliers in position and velocity data.\n   - Apply this function to your estimated positions and velocities if necessary.\n\n8. **Kalman Smoothing**:\n   - Understand the `Kalman_smoothing` function. This function implements forward-backward Kalman smoothing to refine position estimates.\n   - Execute the `Kalman_smoothing` function with your estimated positions and velocities.\n\n9. **Evaluation**:\n   - Compute distance errors and scores for your baseline, robust WLS, and Kalman smoothed estimates.\n   - Visualize the results using the provided plotting functions.\n\n10. **Submission Preparation**:\n    - Prepare submission data by interpolating estimated positions for your specific scenario.\n\n11. **Writing Submission File**:\n    - Write the submission file in CSV format with the interpolated positions for evaluation or competition purposes.\n\n12. **Iterate and Optimize**:\n    - Test the project with different datasets and parameters.\n    - Fine-tune the algorithms and parameters to improve accuracy if necessary.\n\n13. **Documentation and Reporting**:\n    - Document your findings, observations, and any modifications made to the code.\n    - Generate reports or summaries of your results for presentation or further analysis.\n\n14. **Collaboration and Feedback**:\n    - Share your results with colleagues or peers for feedback and collaboration.\n    - Incorporate any suggestions or improvements to enhance the project further.\n","metadata":{}},{"cell_type":"markdown","source":"**Code Report:**\n\n---\n\n### Libraries Used:\n1. `pymap3d`: Used for conversions between Earth-centered, Earth-fixed (ECEF) coordinates and latitude, longitude, and altitude.\n2. `numpy`: Utilized for numerical computations, particularly with arrays.\n3. `pandas`: Employed for data manipulation and analysis, especially for handling tabular data.\n4. `matplotlib.pyplot`: Utilized for creating visualizations such as plots.\n5. `glob`: Used for pathname pattern expansion (i.e., matching filenames with a pattern).\n6. `scipy.optimize`: Utilized for optimization tasks, such as least-squares fitting.\n7. `tqdm.auto`: Utilized for displaying progress bars during iterations.\n8. `scipy.interpolate.InterpolatedUnivariateSpline`: Utilized for interpolating values.\n9. `scipy.spatial.distance`: Utilized for computing distances between points.\n\n### Constants:\n1. `CLIGHT`: Represents the speed of light in meters per second.\n2. `RE_WGS84`: Represents the Earth's semimajor axis in meters (WGS84).\n3. `OMGE`: Represents the Earth's angular velocity in radians per second (IS-GPS).\n\n### Functions Overview:\n1. **Satellite Selection Functions**: These functions filter satellite data based on criteria like carrier frequency error, elevation angle, C/N0, and multipath flag.\n2. **Vector and Matrix Computation Functions**: These functions compute line-of-sight vectors, Jacobian matrices, and residuals for pseudorange and pseudorange rate calculations.\n3. **Carrier Smoothing Function**: Applies carrier smoothing to GNSS pseudorange data.\n4. **Distance Calculation Functions**: Functions to compute distances between points using Vincenty's formulae.\n5. **Score Calculation Function**: Computes a score based on the distance between estimated and ground truth positions.\n6. **Point Positioning Function**: Performs GNSS single point positioning using pseudorange data.\n7. **Outlier Detection and Interpolation Function**: Detects and interpolates outliers in position and velocity data.\n8. **Kalman Filter Functions**: Implements Kalman filtering for state estimation and smoothing.\n\n### Main Execution Steps:\n1. **Data Reading**: Reads GNSS data and ground truth data.\n2. **Point Positioning**: Estimates positions using pseudorange data.\n3. **Outlier Detection and Interpolation**: Excludes velocity outliers and interpolates missing values in position and velocity data.\n4. **Kalman Smoothing**: Applies Kalman smoothing to improve position estimates.\n5. **Evaluation**: Computes distance errors and scores for baseline, robust WLS, and Kalman smoothed estimates.\n6. **Visualization**: Plots distance and speed errors for evaluation.\n7. **Submission File Preparation**: Prepares submission data by interpolating estimated positions.\n8. **Submission File Writing**: Writes the submission file in CSV format.\n\n### Additional Notes:\n- The code is well-structured and modular, with clear separation of concerns between different functions.\n- It demonstrates the use of various numerical and optimization techniques for processing GNSS data and improving position estimation accuracy.\n\n---\n","metadata":{}}]}