{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"Baseline submission tried from [saitodevel01's notebook](https://www.kaggle.com/code/saitodevel01/gsdc2-baseline-submission) and EDA from [nayuts's notebook](https://www.kaggle.com/code/nayuts/let-s-visualize-dataset-to-understand)\n\n**Goal**\nThe goal of this competition is to compute location (latDeg and lngDeg) down to decimeter or even centimeter resolution, if possible for each test phone and time.\n\n**About dataset**\nGoogle releases 60+ datasets collected from phones in the Android GPS team, together with corrections from SwiftNavigation Inc. and Verizon Inc. These datasets were collected on highways in the US San Francisco Bay Area in the summer of 2020. We can see the video for dataset.\n\nWe are given data from actual runs with android devices installed in cars, see following.\n\n![](https://raw.githubusercontent.com/tasotasoso/kaggle_media/main/Android_smartphones_high_accuracy_GNSS_datasets/fig3_fig4.JPG)\n\nThe figures come from Fu, Guoyu (Michael), Khider, Mohammed, van Diggelen, Frank, \"Android Raw GNSS Measurement Datasets for Precise Positioning,\" Proceedings of the 33rd International Technical Meeting of the Satellite Division of The Institute of Navigation (ION GNSS+ 2020), September 2020, pp. 1925-1937. [https://doi.org/10.33012/2020.17628](https://www.ion.org/publications/abstract.cfm?articleID=17628)\n\nWe can see more detail of data collection process at Android smartphones high accuracy GNSS datasets.\n","metadata":{}},{"cell_type":"code","source":"import glob\nfrom dataclasses import dataclass\nimport numpy as np\nimport pandas as pd\nfrom tqdm.notebook import tqdm\nfrom scipy.interpolate import InterpolatedUnivariateSpline\n\nINPUT_PATH = '/kaggle/input/smartphone-decimeter-2023'\n\nWGS84_SEMI_MAJOR_AXIS = 6378137.0\nWGS84_SEMI_MINOR_AXIS = 6356752.314245\nWGS84_SQUARED_FIRST_ECCENTRICITY  = 6.69437999013e-3\nWGS84_SQUARED_SECOND_ECCENTRICITY = 6.73949674226e-3\n\nHAVERSINE_RADIUS = 6_371_000","metadata":{"execution":{"iopub.status.busy":"2023-09-15T10:56:08.355533Z","iopub.execute_input":"2023-09-15T10:56:08.355976Z","iopub.status.idle":"2023-09-15T10:56:08.363568Z","shell.execute_reply.started":"2023-09-15T10:56:08.355942Z","shell.execute_reply":"2023-09-15T10:56:08.362157Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sample_trail_gnss = pd.read_csv(\"/kaggle/input/smartphone-decimeter-2023/sdc2023/train/2020-06-25-00-34-us-ca-mtv-sb-101/pixel4/device_gnss.csv\")\ndf_sample_trail_gt = pd.read_csv(\"/kaggle/input/smartphone-decimeter-2023/sdc2023/train/2020-06-25-00-34-us-ca-mtv-sb-101/pixel4/ground_truth.csv\")","metadata":{"execution":{"iopub.status.busy":"2023-09-15T11:00:50.641253Z","iopub.execute_input":"2023-09-15T11:00:50.641731Z","iopub.status.idle":"2023-09-15T11:00:51.200047Z","shell.execute_reply.started":"2023-09-15T11:00:50.641694Z","shell.execute_reply":"2023-09-15T11:00:51.198943Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sample_trail_gnss.head()","metadata":{"execution":{"iopub.status.busy":"2023-09-15T11:00:53.171442Z","iopub.execute_input":"2023-09-15T11:00:53.171860Z","iopub.status.idle":"2023-09-15T11:00:53.198879Z","shell.execute_reply.started":"2023-09-15T11:00:53.171829Z","shell.execute_reply":"2023-09-15T11:00:53.197781Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sample_trail_gt.head()","metadata":{"execution":{"iopub.status.busy":"2023-09-15T11:00:53.442158Z","iopub.execute_input":"2023-09-15T11:00:53.442579Z","iopub.status.idle":"2023-09-15T11:00:53.463618Z","shell.execute_reply.started":"2023-09-15T11:00:53.442546Z","shell.execute_reply":"2023-09-15T11:00:53.462414Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"@dataclass\nclass ECEF:\n    x: np.array\n    y: np.array\n    z: np.array\n\n    def to_numpy(self):\n        return np.stack([self.x, self.y, self.z], axis=0)\n\n    @staticmethod\n    def from_numpy(pos):\n        x, y, z = [np.squeeze(w) for w in np.split(pos, 3, axis=-1)]\n        return ECEF(x=x, y=y, z=z)\n\n@dataclass\nclass BLH:\n    lat : np.array\n    lng : np.array\n    hgt : np.array\n\ndef ECEF_to_BLH(ecef):\n    a = WGS84_SEMI_MAJOR_AXIS\n    b = WGS84_SEMI_MINOR_AXIS\n    e2  = WGS84_SQUARED_FIRST_ECCENTRICITY\n    e2_ = WGS84_SQUARED_SECOND_ECCENTRICITY\n    x = ecef.x\n    y = ecef.y\n    z = ecef.z\n    r = np.sqrt(x**2 + y**2)\n    t = np.arctan2(z * (a/b), r)\n    B = np.arctan2(z + (e2_*b)*np.sin(t)**3, r - (e2*a)*np.cos(t)**3)\n    L = np.arctan2(y, x)\n    n = a / np.sqrt(1 - e2*np.sin(B)**2)\n    H = (r / np.cos(B)) - n\n    return BLH(lat=B, lng=L, hgt=H)\n\ndef haversine_distance(blh_1, blh_2):\n    dlat = blh_2.lat - blh_1.lat\n    dlng = blh_2.lng - blh_1.lng\n    a = np.sin(dlat/2)**2 + np.cos(blh_1.lat) * np.cos(blh_2.lat) * np.sin(dlng/2)**2\n    dist = 2 * HAVERSINE_RADIUS * np.arcsin(np.sqrt(a))\n    return dist\n\ndef pandas_haversine_distance(df1, df2):\n    blh1 = BLH(\n        lat=np.deg2rad(df1['LatitudeDegrees'].to_numpy()),\n        lng=np.deg2rad(df1['LongitudeDegrees'].to_numpy()),\n        hgt=0,\n    )\n    blh2 = BLH(\n        lat=np.deg2rad(df2['LatitudeDegrees'].to_numpy()),\n        lng=np.deg2rad(df2['LongitudeDegrees'].to_numpy()),\n        hgt=0,\n    )\n    return haversine_distance(blh1, blh2)\n","metadata":{"execution":{"iopub.status.busy":"2023-09-15T10:13:22.965325Z","iopub.execute_input":"2023-09-15T10:13:22.966288Z","iopub.status.idle":"2023-09-15T10:13:22.984372Z","shell.execute_reply.started":"2023-09-15T10:13:22.966256Z","shell.execute_reply":"2023-09-15T10:13:22.982894Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def ecef_to_lat_lng(tripID, gnss_df, UnixTimeMillis):\n    ecef_columns = ['WlsPositionXEcefMeters', 'WlsPositionYEcefMeters', 'WlsPositionZEcefMeters']\n    columns = ['utcTimeMillis'] + ecef_columns\n    ecef_df = (gnss_df.drop_duplicates(subset='utcTimeMillis')[columns]\n               .dropna().reset_index(drop=True))\n    ecef = ECEF.from_numpy(ecef_df[ecef_columns].to_numpy())\n    blh  = ECEF_to_BLH(ecef)\n\n    TIME = ecef_df['utcTimeMillis'].to_numpy()\n    lat = InterpolatedUnivariateSpline(TIME, blh.lat, ext=3)(UnixTimeMillis)\n    lng = InterpolatedUnivariateSpline(TIME, blh.lng, ext=3)(UnixTimeMillis)\n    return pd.DataFrame({\n        'tripId' : tripID,\n        'UnixTimeMillis'   : UnixTimeMillis,\n        'LatitudeDegrees'  : np.degrees(lat),\n        'LongitudeDegrees' : np.degrees(lng),\n    })\n\ndef calc_score(tripID, pred_df, gt_df):\n    d = pandas_haversine_distance(pred_df, gt_df)\n    score = np.mean([np.quantile(d, 0.50), np.quantile(d, 0.95)])    \n    return score\n","metadata":{"execution":{"iopub.status.busy":"2023-09-15T10:13:22.985833Z","iopub.execute_input":"2023-09-15T10:13:22.986243Z","iopub.status.idle":"2023-09-15T10:13:22.999957Z","shell.execute_reply.started":"2023-09-15T10:13:22.986204Z","shell.execute_reply":"2023-09-15T10:13:22.998834Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%capture --no-stdout\n\npred_dfs  = []\nscore_list = []\nfor dirname in sorted(glob.glob(f'/kaggle/input/smartphone-decimeter-2023/sdc2023/train/*/*')):\n    drive, phone = dirname.split('/')[-2:]\n    tripID  = f'{drive}/{phone}'\n    gnss_df = pd.read_csv(f'{dirname}/device_gnss.csv')\n    gt_df   = pd.read_csv(f'{dirname}/ground_truth.csv')\n    pred_df = ecef_to_lat_lng(tripID, gnss_df, gt_df['UnixTimeMillis'].to_numpy())\n    pred_dfs.append(pred_df)\n    score = calc_score(tripID, pred_df, gt_df)\n    print(f'{tripID:<45}: score = {score:.3f}')\n    score_list.append(score)","metadata":{"execution":{"iopub.status.busy":"2023-09-15T10:17:39.097239Z","iopub.execute_input":"2023-09-15T10:17:39.097684Z","iopub.status.idle":"2023-09-15T10:20:39.594043Z","shell.execute_reply.started":"2023-09-15T10:17:39.097647Z","shell.execute_reply":"2023-09-15T10:20:39.592674Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"baseline_train_df = pd.concat(pred_dfs)\nbaseline_train_df.to_csv('baseline_train.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2023-09-15T10:21:29.659929Z","iopub.execute_input":"2023-09-15T10:21:29.660326Z","iopub.status.idle":"2023-09-15T10:21:31.909877Z","shell.execute_reply.started":"2023-09-15T10:21:29.660295Z","shell.execute_reply":"2023-09-15T10:21:31.908545Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mean_score = np.mean(score_list)\nprint(f'mean_score = {mean_score:.3f}')","metadata":{"execution":{"iopub.status.busy":"2023-09-15T10:21:35.314796Z","iopub.execute_input":"2023-09-15T10:21:35.315239Z","iopub.status.idle":"2023-09-15T10:21:35.322159Z","shell.execute_reply.started":"2023-09-15T10:21:35.315204Z","shell.execute_reply":"2023-09-15T10:21:35.320663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sample_df = pd.read_csv(f'/kaggle/input/smartphone-decimeter-2023/sdc2023/sample_submission.csv')\npred_dfs  = []\nfor dirname in tqdm(sorted(glob.glob(f'/kaggle/input/smartphone-decimeter-2023/sdc2023/test/*/*'))):\n    drive, phone = dirname.split('/')[-2:]\n    tripID  = f'{drive}/{phone}'\n    gnss_df = pd.read_csv(f'{dirname}/device_gnss.csv')\n    UnixTimeMillis = sample_df[sample_df['tripId'] == tripID]['UnixTimeMillis'].to_numpy()\n    pred_dfs.append(ecef_to_lat_lng(tripID, gnss_df, UnixTimeMillis))\nbaseline_test_df = pd.concat(pred_dfs)\nbaseline_test_df.to_csv('baseline_test.csv', index=False)\nbaseline_test_df.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2023-09-15T10:30:31.488467Z","iopub.execute_input":"2023-09-15T10:30:31.489222Z","iopub.status.idle":"2023-09-15T10:31:17.330893Z","shell.execute_reply.started":"2023-09-15T10:30:31.489183Z","shell.execute_reply":"2023-09-15T10:31:17.329932Z"},"_kg_hide-output":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}