{"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"}],"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Preparation","metadata":{}},{"cell_type":"code","source":"!pip install pymap3d\n\n# Import functions first\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","metadata":{"execution":{"iopub.status.busy":"2024-04-17T09:27:44.110452Z","iopub.execute_input":"2024-04-17T09:27:44.111132Z","iopub.status.idle":"2024-04-17T09:28:04.886577Z","shell.execute_reply.started":"2024-04-17T09:27:44.111094Z","shell.execute_reply":"2024-04-17T09:28:04.885154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Set up features","metadata":{}},{"cell_type":"code","source":"# Constants: speed of light (m/s), earth semimajor axis (m), earth angular velocity (rad/s)\nCLIGHT = 299_792_458\nRE_WGS84 = 6_378_137\nOMGE = 7.2921151467E-5","metadata":{"execution":{"iopub.status.busy":"2024-04-17T09:37:05.404093Z","iopub.execute_input":"2024-04-17T09:37:05.404506Z","iopub.status.idle":"2024-04-17T09:37:05.409915Z","shell.execute_reply.started":"2024-04-17T09:37:05.404476Z","shell.execute_reply":"2024-04-17T09:37:05.408705Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Adjust satellite direction 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\n    idx &= df['SvElevationDegrees'] > 10.0\n    idx &= df['Cn0DbHz'] > 15.0\n    idx &= df['MultipathIndicator'] == 0\n    \n    return df[idx]","metadata":{"execution":{"iopub.status.busy":"2024-04-17T09:37:06.780274Z","iopub.execute_input":"2024-04-17T09:37:06.780739Z","iopub.status.idle":"2024-04-17T09:37:06.786345Z","shell.execute_reply.started":"2024-04-17T09:37:06.780705Z","shell.execute_reply":"2024-04-17T09:37:06.785429Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Compute line-of-sight vector (user to satellite)\ndef los_vector(xusr, xsat):\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 primary Jacobian matrix\ndef jac_pr_residuals(x, xsat, pr, W):\n    u, _ = los_vector(x[:3], xsat)\n    J = np.hstack([-u, np.ones([len(pr), 1])])\n    \n    return W @ J\n\n# Compute pseudo-range residuals\ndef pr_residuals(x, xsat, pr, W):\n    u, rng = los_vector(x[:3], xsat)\n    # Approximate correction of the earth rotation (Sagnac effect): GNSS positioning\n    rng += OMGE * (xsat[:, 0] * x[1] - xsat[:, 1] * x[0]) / CLIGHT\n    # Enable GPS L1 clock offset\n    residuals = rng - (pr - x[3])\n    \n    return residuals @ W\n\n# Compute secondary Jacobian matrix\ndef jac_prr_residuals(v, vsat, prr, x, xsat, W):\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 pseudo-range rate residuals\ndef prr_residuals(v, vsat, prr, x, xsat, W):\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    residuals = rate - (prr - v[3])\n    \n    return residuals @ W","metadata":{"execution":{"iopub.status.busy":"2024-04-17T09:37:09.442110Z","iopub.execute_input":"2024-04-17T09:37:09.442801Z","iopub.status.idle":"2024-04-17T09:37:09.456491Z","shell.execute_reply.started":"2024-04-17T09:37:09.442766Z","shell.execute_reply":"2024-04-17T09:37:09.455140Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Advanced setting, 1st step","metadata":{}},{"cell_type":"code","source":"# Carrier smoothing of pseudarange\ndef carrier_smoothing(gnss_df):\n    carr_th = 1.0\n    pr_th =  15.0\n    prsmooth = np.full_like(gnss_df['RawPseudorangeMeters'], np.nan)\n    \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}})\n        \n        # Compare time difference between pseudo-range/carrier using Doppler's method\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\n        slip2 = (df['AccumulatedDeltaRangeState'].to_numpy() & 2**2) != 0\n        slip3 = np.fabs(drng1.to_numpy()) > carr_th\n        slip4 = np.fabs(drng2.to_numpy()) > pr_th\n        idx_slip = slip1 | slip2 | slip3 | slip4\n        idx_slip[0] = True\n        \n        # Groups w/continuous carrier phase tracking\n        df['group_slip'] = np.cumsum(idx_slip)\n        \n        # 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 w/bias\n        prsmooth[idx] = df['AccumulatedDeltaRangeMeters'] + df['dpc_Mean']\n        \n    # If carrier smoothing is not possible, use original pseudo-range\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","metadata":{"execution":{"iopub.status.busy":"2024-04-17T09:37:12.153166Z","iopub.execute_input":"2024-04-17T09:37:12.153622Z","iopub.status.idle":"2024-04-17T09:37:12.167460Z","shell.execute_reply.started":"2024-04-17T09:37:12.153587Z","shell.execute_reply":"2024-04-17T09:37:12.165952Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Compute distance using Vincenty's formula\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","metadata":{"execution":{"iopub.status.busy":"2024-04-17T09:37:15.773540Z","iopub.execute_input":"2024-04-17T09:37:15.773957Z","iopub.status.idle":"2024-04-17T09:37:15.781072Z","shell.execute_reply.started":"2024-04-17T09:37:15.773926Z","shell.execute_reply":"2024-04-17T09:37:15.779880Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# GNSS single-point positioning using pseudo-range\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)\n    v0 = np.zeros(4)\n    x_wls = np.full([nepoch, 3], np.nan)\n    v_wls = np.full([nepoch, 3], np.nan)\n    \n    # Loop for epochs\n    for i, (t_utc, df) in enumerate(tqdm(gnss_df.groupby('utcTimeMillis'), total=nepoch)):\n        \n        # Valid satellite selection\n        df_pr = satellite_selection(df, 'pr_smooth')\n        df_prr = satellite_selection(df, 'PseudorangeRateMetersPerSecond')\n        \n        # Corrected pseudo-range/pseudo-range 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 pseudo-range/pseudo-range 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            \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            \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            \n            # Velocity estimation\n            if len(df_prr) >= 4:\n                if np.all(v0 == 0):\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            \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","metadata":{"execution":{"iopub.status.busy":"2024-04-17T09:37:19.108272Z","iopub.execute_input":"2024-04-17T09:37:19.108689Z","iopub.status.idle":"2024-04-17T09:37:19.131062Z","shell.execute_reply.started":"2024-04-17T09:37:19.108658Z","shell.execute_reply":"2024-04-17T09:37:19.129623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# 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    # Velocity jump detection\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()","metadata":{"execution":{"iopub.status.busy":"2024-04-17T09:37:23.175884Z","iopub.execute_input":"2024-04-17T09:37:23.176522Z","iopub.status.idle":"2024-04-17T09:37:23.186188Z","shell.execute_reply.started":"2024-04-17T09:37:23.176489Z","shell.execute_reply":"2024-04-17T09:37:23.184741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Define Kalman filter function\ndef Kalman_filter(zs, us, phone):\n    # Set up parameters\n    sigma_v = 0.6 if phone == 'OppoF3' else 0.1\n    sigma_x = 5.0\n    sigma_mahalanobis = 30.0\n    \n    n, dim_x = zs.shape\n    F = np.eye(3)\n    Q = sigma_v**2 * np.eye(3)\n\n    H = np.eye(3)\n    R = sigma_x**2 * np.eye(3)\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    # Initiate filtering process\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 and 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)","metadata":{"execution":{"iopub.status.busy":"2024-04-17T09:37:27.081480Z","iopub.execute_input":"2024-04-17T09:37:27.082333Z","iopub.status.idle":"2024-04-17T09:37:27.104554Z","shell.execute_reply.started":"2024-04-17T09:37:27.082288Z","shell.execute_reply":"2024-04-17T09:37:27.103284Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Advanced setting, 2nd step","metadata":{}},{"cell_type":"code","source":"# Set the file destination\npath = '/kaggle/input/smartphone-decimeter-2023/sdc2023/train/2023-09-07-22-48-us-ca-routebc2/pixel4xl'\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# Assign scores\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# Visualize the distance error rate\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# Visualize the velocity error rate\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()","metadata":{"execution":{"iopub.status.busy":"2024-04-17T09:37:32.660043Z","iopub.execute_input":"2024-04-17T09:37:32.661313Z","iopub.status.idle":"2024-04-17T09:38:08.783657Z","shell.execute_reply.started":"2024-04-17T09:37:32.661271Z","shell.execute_reply":"2024-04-17T09:38:08.782305Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create a submission","metadata":{}},{"cell_type":"code","source":"# Set the file destination\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# Save as CSV submission file\ntest_df = pd.concat(test_dfs)\ntest_df.to_csv('submission.csv', index=False)\nprint(\"Successfully saved as CSV file\")","metadata":{"execution":{"iopub.status.busy":"2024-04-17T09:38:32.736822Z","iopub.execute_input":"2024-04-17T09:38:32.737300Z","iopub.status.idle":"2024-04-17T10:01:21.203521Z","shell.execute_reply.started":"2024-04-17T09:38:32.737264Z","shell.execute_reply":"2024-04-17T10:01:21.201934Z"},"trusted":true},"execution_count":null,"outputs":[]}]}