{"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":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os\nfrom pathlib import Path\nimport pyproj\nfrom pyproj import Proj, transform\n\ndef calc_haversine(lat1, lon1, lat2, lon2):\n    RADIUS = 6_367_000\n    lat1, lon1, lat2, lon2 = map(np.radians, [lat1, lon1, lat2, lon2])\n    dlat = lat2 - lat1\n    dlon = lon2 - lon1\n    a = np.sin(dlat/2)**2 + \\\n        np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2\n    dist = 2 * RADIUS * np.arcsin(a**0.5)\n    return dist\n \ndef compute_dist(fname, fname2 = 'gt.csv'):\n    btrain = pd.read_csv(fname)\n    gt = pd.read_csv(fname2)\n    df = btrain.merge(gt, on = ['phone','millisSinceGpsEpoch'])\n    dst_btrain = calc_haversine(df.latDeg_x,df.lngDeg_x, df.latDeg_y, df.lngDeg_y)\n    scores = pd.DataFrame({'phone': df.phone,'dst': dst_btrain})\n    scores_grp = scores.groupby('phone')\n    d50 = scores_grp.quantile(.50).reset_index()\n    d50.columns = ['phone','q50']\n    d95 = scores_grp.quantile(.95).reset_index()\n    d95.columns = ['phone','q95']\n    return (scores_grp.quantile(.50).mean() + scores_grp.quantile(.95).mean())/2, d50.merge(d95)\n\ndef WGS84_to_ECEF(lat, lon, alt):\n    # convert to radians\n    rad_lat = lat * (np.pi / 180.0)\n    rad_lon = lon * (np.pi / 180.0)\n    a    = 6378137.0\n    # f is the flattening factor\n    finv = 298.257223563\n    f = 1 / finv   \n    # e is the eccentricity\n    e2 = 1 - (1 - f) * (1 - f)    \n    # N is the radius of curvature in the prime vertical\n    N = a / np.sqrt(1 - e2 * np.sin(rad_lat) * np.sin(rad_lat))\n    x = (N + alt) * np.cos(rad_lat) * np.cos(rad_lon)\n    y = (N + alt) * np.cos(rad_lat) * np.sin(rad_lon)\n    z = (N * (1 - e2) + alt)        * np.sin(rad_lat)\n    return x, y, z\n\ntransformer = pyproj.Transformer.from_crs(\n    {\"proj\":'geocent', \"ellps\":'WGS84', \"datum\":'WGS84'},\n    {\"proj\":'latlong', \"ellps\":'WGS84', \"datum\":'WGS84'},)\ndef ECEF_to_WGS84(x,y,z):\n    lon, lat, alt = transformer.transform(x,y,z,radians=False)\n    return lon, lat, alt\n\ndatadir = Path('/kaggle/input/google-smartphone-decimeter-challenge/')\ntestdir = datadir / 'test'\ntraindir = datadir / 'train'\n\nsample_sub = pd.read_csv(datadir/'sample_submission.csv')\nsub_columns = sample_sub.columns\n\nbaseline_train = pd.read_csv(datadir / 'baseline_locations_train.csv')\nbaseline_train[sub_columns].to_csv('btrain.csv',index = False)\nbaseline_test = pd.read_csv(datadir / 'baseline_locations_test.csv')\nbaseline_test[sub_columns].to_csv('btest.csv',index = False)\n\nmsge = 'millisSinceGpsEpoch'\n\ngt = pd.DataFrame()\nfor d in os.listdir(traindir):\n    for p in os.listdir(traindir/d):\n        gt = gt.append(pd.read_csv(traindir/d/p/'ground_truth.csv'))\n\ngt['phone'] = gt['collectionName'] + '_' + gt['phoneName']\ngt[sub_columns].to_csv('gt.csv', index = False)\ngt['heightAboveWgs84EllipsoidM'].describe()\ndf = pd.read_csv('btrain.csv').merge(gt, on = ['phone','millisSinceGpsEpoch'])\nprint(df)\nscore, scores = compute_dist('btrain.csv','gt.csv')\nprint(score)\nscores","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-12-24T10:39:57.081062Z","iopub.execute_input":"2021-12-24T10:39:57.081474Z","iopub.status.idle":"2021-12-24T10:40:00.511004Z","shell.execute_reply.started":"2021-12-24T10:39:57.081441Z","shell.execute_reply":"2021-12-24T10:40:00.510336Z"},"trusted":true},"execution_count":null,"outputs":[]}]}