{"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":"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":"2023-09-21T09:24:10.247723Z","iopub.execute_input":"2023-09-21T09:24:10.248270Z","iopub.status.idle":"2023-09-21T09:24:13.646034Z","shell.execute_reply.started":"2023-09-21T09:24:10.248223Z","shell.execute_reply":"2023-09-21T09:24:13.644600Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import glob\nfrom dataclasses import dataclass\nimport seaborn as sns\nimport matplotlib.pyplot as plt\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","metadata":{"execution":{"iopub.status.busy":"2023-09-21T09:24:13.648684Z","iopub.execute_input":"2023-09-21T09:24:13.649445Z","iopub.status.idle":"2023-09-21T09:24:13.654845Z","shell.execute_reply.started":"2023-09-21T09:24:13.649408Z","shell.execute_reply":"2023-09-21T09:24:13.653973Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"INPUT_PATH","metadata":{"execution":{"iopub.status.busy":"2023-09-21T09:24:13.656520Z","iopub.execute_input":"2023-09-21T09:24:13.657176Z","iopub.status.idle":"2023-09-21T09:24:13.676364Z","shell.execute_reply.started":"2023-09-21T09:24:13.657141Z","shell.execute_reply":"2023-09-21T09:24:13.674617Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"WGS84_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-21T09:24:13.679311Z","iopub.execute_input":"2023-09-21T09:24:13.680250Z","iopub.status.idle":"2023-09-21T09:24:13.688936Z","shell.execute_reply.started":"2023-09-21T09:24:13.680212Z","shell.execute_reply":"2023-09-21T09:24:13.687678Z"},"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\")\ndf_sample_trail_gnss.head()","metadata":{"execution":{"iopub.status.busy":"2023-09-21T09:24:13.690431Z","iopub.execute_input":"2023-09-21T09:24:13.691636Z","iopub.status.idle":"2023-09-21T09:24:14.828469Z","shell.execute_reply.started":"2023-09-21T09:24:13.691591Z","shell.execute_reply":"2023-09-21T09:24:14.827165Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gnss_null=df_sample_trail_gnss.head().isnull()","metadata":{"execution":{"iopub.status.busy":"2023-09-21T09:24:14.830001Z","iopub.execute_input":"2023-09-21T09:24:14.830556Z","iopub.status.idle":"2023-09-21T09:24:14.837560Z","shell.execute_reply.started":"2023-09-21T09:24:14.830504Z","shell.execute_reply":"2023-09-21T09:24:14.836217Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(10, 6))\nsns.heatmap(gnss_null, cmap='viridis', cbar=False)\nplt.title('Missing Values Heatmap')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-09-21T09:24:14.839307Z","iopub.execute_input":"2023-09-21T09:24:14.839688Z","iopub.status.idle":"2023-09-21T09:24:15.439005Z","shell.execute_reply.started":"2023-09-21T09:24:14.839658Z","shell.execute_reply":"2023-09-21T09:24:15.437562Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gnss_null_counts=df_sample_trail_gnss.head().isnull().sum()","metadata":{"execution":{"iopub.status.busy":"2023-09-21T09:24:15.440400Z","iopub.execute_input":"2023-09-21T09:24:15.440726Z","iopub.status.idle":"2023-09-21T09:24:15.448759Z","shell.execute_reply.started":"2023-09-21T09:24:15.440697Z","shell.execute_reply":"2023-09-21T09:24:15.447341Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gnss_null_counts","metadata":{"execution":{"iopub.status.busy":"2023-09-21T09:24:15.450607Z","iopub.execute_input":"2023-09-21T09:24:15.450957Z","iopub.status.idle":"2023-09-21T09:24:15.465096Z","shell.execute_reply.started":"2023-09-21T09:24:15.450927Z","shell.execute_reply":"2023-09-21T09:24:15.463501Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gt_null_values_count=df_sample_trail_gt.head().isnull()","metadata":{"execution":{"iopub.status.busy":"2023-09-21T09:24:15.470296Z","iopub.execute_input":"2023-09-21T09:24:15.471391Z","iopub.status.idle":"2023-09-21T09:24:15.479189Z","shell.execute_reply.started":"2023-09-21T09:24:15.471309Z","shell.execute_reply":"2023-09-21T09:24:15.477291Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gt_null_values_count","metadata":{"execution":{"iopub.status.busy":"2023-09-21T09:24:15.480935Z","iopub.execute_input":"2023-09-21T09:24:15.481501Z","iopub.status.idle":"2023-09-21T09:24:15.507615Z","shell.execute_reply.started":"2023-09-21T09:24:15.481447Z","shell.execute_reply":"2023-09-21T09:24:15.506343Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(10, 6))\nsns.heatmap(gt_null_values_count, cmap='viridis', cbar=False)\nplt.title('Missing Values Heatmap')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-09-21T09:24:15.509239Z","iopub.execute_input":"2023-09-21T09:24:15.509618Z","iopub.status.idle":"2023-09-21T09:24:15.863748Z","shell.execute_reply.started":"2023-09-21T09:24:15.509588Z","shell.execute_reply":"2023-09-21T09:24:15.862408Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gt_null_values_counts=df_sample_trail_gt.head().isnull().sum()","metadata":{"execution":{"iopub.status.busy":"2023-09-21T09:24:15.865253Z","iopub.execute_input":"2023-09-21T09:24:15.865713Z","iopub.status.idle":"2023-09-21T09:24:15.871628Z","shell.execute_reply.started":"2023-09-21T09:24:15.865675Z","shell.execute_reply":"2023-09-21T09:24:15.870673Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gt_null_values_counts","metadata":{"execution":{"iopub.status.busy":"2023-09-21T09:24:15.872759Z","iopub.execute_input":"2023-09-21T09:24:15.873379Z","iopub.status.idle":"2023-09-21T09:24:15.890192Z","shell.execute_reply.started":"2023-09-21T09:24:15.873344Z","shell.execute_reply":"2023-09-21T09:24:15.888389Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sample_trail_gt.head()","metadata":{"execution":{"iopub.status.busy":"2023-09-21T09:24:15.891697Z","iopub.execute_input":"2023-09-21T09:24:15.892042Z","iopub.status.idle":"2023-09-21T09:24:15.916249Z","shell.execute_reply.started":"2023-09-21T09:24:15.892005Z","shell.execute_reply":"2023-09-21T09:24:15.915209Z"},"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)","metadata":{"execution":{"iopub.status.busy":"2023-09-21T09:24:15.918268Z","iopub.execute_input":"2023-09-21T09:24:15.918660Z","iopub.status.idle":"2023-09-21T09:24:15.938145Z","shell.execute_reply.started":"2023-09-21T09:24:15.918629Z","shell.execute_reply":"2023-09-21T09:24:15.937034Z"},"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","metadata":{"execution":{"iopub.status.busy":"2023-09-21T09:24:15.940042Z","iopub.execute_input":"2023-09-21T09:24:15.940520Z","iopub.status.idle":"2023-09-21T09:24:15.957101Z","shell.execute_reply.started":"2023-09-21T09:24:15.940478Z","shell.execute_reply":"2023-09-21T09:24:15.955903Z"},"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-21T09:24:15.959311Z","iopub.execute_input":"2023-09-21T09:24:15.959807Z","iopub.status.idle":"2023-09-21T09:27:48.823968Z","shell.execute_reply.started":"2023-09-21T09:24:15.959766Z","shell.execute_reply":"2023-09-21T09:27:48.822690Z"},"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-21T09:27:48.826097Z","iopub.execute_input":"2023-09-21T09:27:48.826890Z","iopub.status.idle":"2023-09-21T09:27:51.232688Z","shell.execute_reply.started":"2023-09-21T09:27:48.826851Z","shell.execute_reply":"2023-09-21T09:27:51.231329Z"},"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-21T09:27:51.234335Z","iopub.execute_input":"2023-09-21T09:27:51.235140Z","iopub.status.idle":"2023-09-21T09:27:51.244204Z","shell.execute_reply.started":"2023-09-21T09:27:51.235093Z","shell.execute_reply":"2023-09-21T09:27:51.242083Z"},"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-21T09:27:51.246882Z","iopub.execute_input":"2023-09-21T09:27:51.248157Z","iopub.status.idle":"2023-09-21T09:28:52.851614Z","shell.execute_reply.started":"2023-09-21T09:27:51.248112Z","shell.execute_reply":"2023-09-21T09:28:52.850261Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}