{"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":"# Google Smartphone Decimeter Challenge 2022\n![myimg](https://thehackernews.com/images/-MWrikkgPuOk/Wk4SJfLxqAI/AAAAAAAAvYw/sPCXIrCBxvQfuUBFg-v_yJvD1wllXdTcgCLcBGAs/s728-e100/gps-location-tracking-device.png)\n\nThis notebook was created during a live stream on twitch:\n- Check out my twitch channel here: [link](https://www.twitch.tv/medallionstallion_)\n- Shameless plug for my youtube channel with videos about data science and machine learning. [Check it out here.](https://www.youtube.com/channel/UCxladMszXan-jfgzyeIMyvw)","metadata":{}},{"cell_type":"markdown","source":"# Data\n\nIn this competition we are tasked with identifying the exact location of a phone using it's GPS data. As the competition description says they would like us to find the... \"location down to decimeter or even centimeter resolution which could enable services that require lane-level accuracy such as HOV lane ETA estimation.\"\n","metadata":{}},{"cell_type":"code","source":"!pip install nb_black > /dev/null","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:46:17.014424Z","iopub.execute_input":"2022-05-03T02:46:17.015082Z","iopub.status.idle":"2022-05-03T02:46:31.617492Z","shell.execute_reply.started":"2022-05-03T02:46:17.014983Z","shell.execute_reply":"2022-05-03T02:46:31.616515Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%load_ext lab_black","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:46:31.619517Z","iopub.execute_input":"2022-05-03T02:46:31.619806Z","iopub.status.idle":"2022-05-03T02:46:31.932557Z","shell.execute_reply.started":"2022-05-03T02:46:31.61977Z","shell.execute_reply":"2022-05-03T02:46:31.931391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Imports\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pylab as plt\nimport plotly.express as px\n\npd.set_option(\"max_columns\", 500)","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:46:31.934259Z","iopub.execute_input":"2022-05-03T02:46:31.935032Z","iopub.status.idle":"2022-05-03T02:46:33.327478Z","shell.execute_reply.started":"2022-05-03T02:46:31.934979Z","shell.execute_reply":"2022-05-03T02:46:33.326046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Training Data\n\nThe data is setup in the `smartphone-decimeter-2022` folder into train and test. Each route is it's own folder like `2020-05-15-US-MTV-1` each phone is a folder under the route like `GooglePixel4XL`.\n\nThe main data sources are:\n- The target: For each phone there is a `ground_truth.csv` with it's location at timestamps.\n- Training data `device_gnss.csv` Each row contains raw GNSS measurements.\n- Training data `device_imu.csv` Readings the phone's accelerometer, gyroscope, and magnetometer. \n","metadata":{}},{"cell_type":"code","source":"trip_id = \"2020-05-15-US-MTV-1/GooglePixel4XL\"","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:46:33.329416Z","iopub.execute_input":"2022-05-03T02:46:33.329689Z","iopub.status.idle":"2022-05-03T02:46:33.337734Z","shell.execute_reply.started":"2022-05-03T02:46:33.329659Z","shell.execute_reply":"2022-05-03T02:46:33.336506Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Pull Data for an example phone\ngt = pd.read_csv(\n    \"../input/smartphone-decimeter-2022/train/2020-05-15-US-MTV-1/GooglePixel4XL/ground_truth.csv\"\n)\ngnss = pd.read_csv(\n    \"../input/smartphone-decimeter-2022/train/2020-05-15-US-MTV-1/GooglePixel4XL/device_gnss.csv\"\n)\nimu = pd.read_csv(\n    \"../input/smartphone-decimeter-2022/train/2020-05-15-US-MTV-1/GooglePixel4XL/device_imu.csv\"\n)","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:46:33.339612Z","iopub.execute_input":"2022-05-03T02:46:33.340243Z","iopub.status.idle":"2022-05-03T02:46:36.198913Z","shell.execute_reply.started":"2022-05-03T02:46:33.340191Z","shell.execute_reply":"2022-05-03T02:46:36.197557Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create Baseline Submission using GNSS Data\nShout out to this notebook for the helper code that will get us started converting GNSS to lat/long\n\nhttps://www.kaggle.com/code/saitodevel01/gsdc2-baseline-submission","metadata":{"execution":{"iopub.status.busy":"2022-05-03T00:49:49.024469Z","iopub.execute_input":"2022-05-03T00:49:49.024794Z","iopub.status.idle":"2022-05-03T00:49:49.090549Z","shell.execute_reply.started":"2022-05-03T00:49:49.02476Z","shell.execute_reply":"2022-05-03T00:49:49.089737Z"}}},{"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 = \"../input/smartphone-decimeter-2022\"\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\n\n\n@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\n@dataclass\nclass BLH:\n    lat: np.array\n    lng: np.array\n    hgt: np.array\n\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\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 = (\n        np.sin(dlat / 2) ** 2\n        + np.cos(blh_1.lat) * np.cos(blh_2.lat) * np.sin(dlng / 2) ** 2\n    )\n    dist = 2 * HAVERSINE_RADIUS * np.arcsin(np.sqrt(a))\n    return dist\n\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\n\ndef ecef_to_lat_lng(tripID, gnss_df, UnixTimeMillis):\n    ecef_columns = [\n        \"WlsPositionXEcefMeters\",\n        \"WlsPositionYEcefMeters\",\n        \"WlsPositionZEcefMeters\",\n    ]\n    columns = [\"utcTimeMillis\"] + ecef_columns\n    ecef_df = (\n        gnss_df.drop_duplicates(subset=\"utcTimeMillis\")[columns]\n        .dropna()\n        .reset_index(drop=True)\n    )\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        {\n            \"tripId\": tripID,\n            \"UnixTimeMillis\": UnixTimeMillis,\n            \"LatitudeDegrees\": np.degrees(lat),\n            \"LongitudeDegrees\": np.degrees(lng),\n        }\n    )\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":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-05-03T02:46:36.200299Z","iopub.execute_input":"2022-05-03T02:46:36.200569Z","iopub.status.idle":"2022-05-03T02:46:36.519233Z","shell.execute_reply.started":"2022-05-03T02:46:36.200531Z","shell.execute_reply":"2022-05-03T02:46:36.518141Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ss = pd.read_csv(\"../input/smartphone-decimeter-2022/sample_submission.csv\")","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:46:36.520839Z","iopub.execute_input":"2022-05-03T02:46:36.521291Z","iopub.status.idle":"2022-05-03T02:46:36.622518Z","shell.execute_reply.started":"2022-05-03T02:46:36.521255Z","shell.execute_reply":"2022-05-03T02:46:36.621369Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trip_id = \"2020-05-15-US-MTV-1/GooglePixel4XL\"\nbaseline = ecef_to_lat_lng(trip_id, gnss, gt[\"UnixTimeMillis\"].values)","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:46:36.62371Z","iopub.execute_input":"2022-05-03T02:46:36.62573Z","iopub.status.idle":"2022-05-03T02:46:36.671299Z","shell.execute_reply.started":"2022-05-03T02:46:36.625685Z","shell.execute_reply":"2022-05-03T02:46:36.670216Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Plotting Function for the Paths","metadata":{"execution":{"iopub.status.busy":"2022-05-03T00:55:46.46941Z","iopub.execute_input":"2022-05-03T00:55:46.469769Z","iopub.status.idle":"2022-05-03T00:55:46.481442Z","shell.execute_reply.started":"2022-05-03T00:55:46.469732Z","shell.execute_reply":"2022-05-03T00:55:46.480822Z"}}},{"cell_type":"code","source":"def visualize_traffic(\n    df,\n    lat_col=\"LatitudeDegrees\",\n    lon_col=\"LongitudeDegrees\",\n    center=None,\n    color_col=\"phone\",\n    label_col=\"tripId\",\n    zoom=9,\n    opacity=1,\n):\n    if center is None:\n        center = {\n            \"lat\": df[lat_col].mean(),\n            \"lon\": df[lon_col].mean(),\n        }\n    fig = px.scatter_mapbox(\n        df,\n        # Here, plotly gets, (x,y) coordinates\n        lat=lat_col,\n        lon=lon_col,\n        # Here, plotly detects color of series\n        color=color_col,\n        labels=label_col,\n        zoom=zoom,\n        center=center,\n        height=600,\n        width=800,\n        opacity=0.5,\n    )\n    fig.update_layout(mapbox_style=\"stamen-terrain\")\n    fig.update_layout(margin={\"r\": 0, \"t\": 0, \"l\": 0, \"b\": 0})\n    fig.update_layout(title_text=\"GPS trafic\")\n    fig.show()\n\n\ndef plot_gt_vs_baseline(tripId):\n    \"\"\"\n    Create a plot of the baseline predictions vs. the ground truth\n    for a given tripId\n    \"\"\"\n    # Pull Data for an example phone\n    gt = pd.read_csv(\n        f\"../input/smartphone-decimeter-2022/train/{tripId}/ground_truth.csv\"\n    )\n    gnss = pd.read_csv(\n        f\"../input/smartphone-decimeter-2022/train/{tripId}/device_gnss.csv\"\n    )\n    imu = pd.read_csv(\n        f\"../input/smartphone-decimeter-2022/train/{tripId}/device_imu.csv\"\n    )\n    baseline = ecef_to_lat_lng(trip_id, gnss, gt[\"UnixTimeMillis\"].values)\n    # Combine ground truth with baseline predictions\n    baseline[\"isGT\"] = False\n    gt[\"isGT\"] = True\n    gt[\"tripId\"] = tripId\n\n    combined = (\n        pd.concat([baseline, gt[baseline.columns]], axis=0)\n        .reset_index(drop=True)\n        .copy()\n    )\n\n    # Plotting the route\n    visualize_traffic(\n        combined,\n        lat_col=\"LatitudeDegrees\",\n        lon_col=\"LongitudeDegrees\",\n        color_col=\"isGT\",\n        zoom=10,\n    )","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:46:36.672764Z","iopub.execute_input":"2022-05-03T02:46:36.673025Z","iopub.status.idle":"2022-05-03T02:46:36.858926Z","shell.execute_reply.started":"2022-05-03T02:46:36.672994Z","shell.execute_reply":"2022-05-03T02:46:36.857853Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Concat the Baseline Predictions with the Ground Truth\n- If you zoom in you can see how the blue and red dots are different.\n- The baseline predictions are very noisy when the car is moving at slow speeds or stopped.","metadata":{}},{"cell_type":"code","source":"plot_gt_vs_baseline(trip_id)","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:46:36.86304Z","iopub.execute_input":"2022-05-03T02:46:36.863305Z","iopub.status.idle":"2022-05-03T02:46:40.013577Z","shell.execute_reply.started":"2022-05-03T02:46:36.863276Z","shell.execute_reply":"2022-05-03T02:46:40.012554Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Plot an Example and Zoom at noisy area","metadata":{}},{"cell_type":"code","source":"from glob import glob\n\ntrain_gts = glob(\"../input/smartphone-decimeter-2022/train/*/*/ground_truth.csv\")\ntrip_ids = [\"/\".join(p.split(\"/\")[-3:-1]) for p in train_gts]","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:46:40.014822Z","iopub.execute_input":"2022-05-03T02:46:40.015844Z","iopub.status.idle":"2022-05-03T02:46:40.482667Z","shell.execute_reply.started":"2022-05-03T02:46:40.015785Z","shell.execute_reply":"2022-05-03T02:46:40.481283Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tripId = trip_ids[10]\n# Pull Data for an example phone\ngt = pd.read_csv(f\"../input/smartphone-decimeter-2022/train/{tripId}/ground_truth.csv\")\ngnss = pd.read_csv(f\"../input/smartphone-decimeter-2022/train/{tripId}/device_gnss.csv\")\nimu = pd.read_csv(f\"../input/smartphone-decimeter-2022/train/{tripId}/device_imu.csv\")\nbaseline = ecef_to_lat_lng(trip_id, gnss, gt[\"UnixTimeMillis\"].values)\n# Combine ground truth with baseline predictions\nbaseline[\"isGT\"] = False\ngt[\"isGT\"] = True\ngt[\"tripId\"] = tripId\n\ncombined = (\n    pd.concat([baseline, gt[baseline.columns]], axis=0).reset_index(drop=True).copy()\n)","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:46:40.484365Z","iopub.execute_input":"2022-05-03T02:46:40.484701Z","iopub.status.idle":"2022-05-03T02:46:41.787227Z","shell.execute_reply.started":"2022-05-03T02:46:40.484656Z","shell.execute_reply":"2022-05-03T02:46:41.786206Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Example of baseline noise\n- Note that the blue (baseline predictions) are very noisy at the start of this trip","metadata":{}},{"cell_type":"code","source":"offset = 150_000\nstart = 1607640760432 + offset\ncombined.query(\"UnixTimeMillis < @start\")\n\nvisualize_traffic(combined.query(\"UnixTimeMillis < @start\"), color_col=\"isGT\", zoom=16)","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:46:41.788625Z","iopub.execute_input":"2022-05-03T02:46:41.789737Z","iopub.status.idle":"2022-05-03T02:46:41.865859Z","shell.execute_reply.started":"2022-05-03T02:46:41.789647Z","shell.execute_reply":"2022-05-03T02:46:41.864891Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create Baselines for Training Set and Test Set\n\nUsing the awesome code from: https://www.kaggle.com/code/saitodevel01/gsdc2-baseline-submission","metadata":{}},{"cell_type":"code","source":"import glob\n\nINPUT_PATH = \"../input/smartphone-decimeter-2022\"\n\nsample_df = pd.read_csv(f\"{INPUT_PATH}/sample_submission.csv\")\npred_dfs = []\nfor dirname in tqdm(sorted(glob.glob(f\"{INPUT_PATH}/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][\n        \"UnixTimeMillis\"\n    ].to_numpy()\n    pred_dfs.append(ecef_to_lat_lng(tripID, gnss_df, UnixTimeMillis))\nsub_df = pd.concat(pred_dfs)\n\nbaselines = []\ngts = []\nfor dirname in tqdm(sorted(glob.glob(f\"{INPUT_PATH}/train/*/*\"))):\n    drive, phone = dirname.split(\"/\")[-2:]\n    tripID = f\"{drive}/{phone}\"\n    gnss_df = pd.read_csv(f\"{dirname}/device_gnss.csv\", low_memory=False)\n    gt_df = pd.read_csv(f\"{dirname}/ground_truth.csv\", low_memory=False)\n    baseline_df = ecef_to_lat_lng(tripID, gnss_df, gt_df[\"UnixTimeMillis\"].to_numpy())\n    baselines.append(baseline_df)\n    gts.append(gt_df)\nbaselines = pd.concat(baselines)\ngts = pd.concat(gts)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2022-05-03T02:46:41.867567Z","iopub.execute_input":"2022-05-03T02:46:41.868091Z","iopub.status.idle":"2022-05-03T02:50:34.479969Z","shell.execute_reply.started":"2022-05-03T02:46:41.868042Z","shell.execute_reply":"2022-05-03T02:50:34.478821Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"baselines[\"group\"] = \"train_baseline\"\nsub_df[\"group\"] = \"submission_baseline\"\ngts[\"group\"] = \"train_ground_truth\"\ncombined = pd.concat([baselines, sub_df, gts]).reset_index(drop=True).copy()","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:50:34.481885Z","iopub.execute_input":"2022-05-03T02:50:34.482199Z","iopub.status.idle":"2022-05-03T02:50:34.65608Z","shell.execute_reply.started":"2022-05-03T02:50:34.482154Z","shell.execute_reply":"2022-05-03T02:50:34.655188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# This Dataset has to main areas\n- Paths from the bay area and LA. Lets split the data into the two areas","metadata":{}},{"cell_type":"code","source":"sf_paths = combined.query(\"LatitudeDegrees > 36\").copy()\nla_paths = combined.query(\"LatitudeDegrees < 36\").copy()","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:51:13.459865Z","iopub.execute_input":"2022-05-03T02:51:13.460261Z","iopub.status.idle":"2022-05-03T02:51:13.573639Z","shell.execute_reply.started":"2022-05-03T02:51:13.460232Z","shell.execute_reply":"2022-05-03T02:51:13.572766Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Plot SF Paths\n- Ground Truth\n- Baseline for training paths\n- Baseline for test paths","metadata":{}},{"cell_type":"code","source":"visualize_traffic(\n    sf_paths.sample(frac=0.2),\n    lat_col=\"LatitudeDegrees\",\n    lon_col=\"LongitudeDegrees\",\n    color_col=\"group\",\n    zoom=9,\n)","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:51:14.758842Z","iopub.execute_input":"2022-05-03T02:51:14.759122Z","iopub.status.idle":"2022-05-03T02:51:15.427146Z","shell.execute_reply.started":"2022-05-03T02:51:14.759093Z","shell.execute_reply":"2022-05-03T02:51:15.426046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"visualize_traffic(\n    la_paths.sample(frac=0.2),\n    lat_col=\"LatitudeDegrees\",\n    lon_col=\"LongitudeDegrees\",\n    color_col=\"group\",\n    zoom=9,\n)","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:51:21.248669Z","iopub.execute_input":"2022-05-03T02:51:21.24896Z","iopub.status.idle":"2022-05-03T02:51:21.361397Z","shell.execute_reply.started":"2022-05-03T02:51:21.248929Z","shell.execute_reply":"2022-05-03T02:51:21.360616Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Post Processing Idea 1\n- Find Stopped Locations and Isolate","metadata":{}},{"cell_type":"code","source":"def calc_haversine(lat1, lon1, lat2, lon2):\n    \"\"\"Calculates the great circle distance between two points\n    on the earth. Inputs are array-like and specified in decimal degrees.\n    \"\"\"\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 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2\n    dist = 2 * RADIUS * np.arcsin(a**0.5)\n    return dist\n\n\ndef add_prev_post_shift(\n    df,\n    lat_col=\"LatitudeDegrees\",\n    lng_col=\"LongitudeDegrees\",\n    dist_suffix=\"\",\n    sortby=[\"tripId\", \"UnixTimeMillis\"],\n):\n    df = df.sort_values(sortby).reset_index(drop=True)\n    df[f\"{lat_col}_shift1\"] = df.groupby([\"tripId\"])[lat_col].shift(1)\n    df[f\"{lng_col}_shift1\"] = df.groupby([\"tripId\"])[lng_col].shift(1)\n    df[f\"{lat_col}_shift-1\"] = df.groupby([\"tripId\"])[lat_col].shift(-1)\n    df[f\"{lng_col}_shift-1\"] = df.groupby([\"tripId\"])[lng_col].shift(-1)\n\n    df[f\"UnixTimeMillis_shift1\"] = df.groupby([\"tripId\"])[\"UnixTimeMillis\"].shift(1)\n    df[f\"UnixTimeMillis_shift-1\"] = df.groupby([\"tripId\"])[\"UnixTimeMillis\"].shift(-1)\n\n    df[f\"dist_prev{dist_suffix}\"] = calc_haversine(\n        df[lat_col], df[lng_col], df[f\"{lat_col}_shift1\"], df[f\"{lng_col}_shift1\"]\n    )\n    df[f\"dist_post{dist_suffix}\"] = calc_haversine(\n        df[lat_col], df[lng_col], df[f\"{lat_col}_shift-1\"], df[f\"{lng_col}_shift-1\"]\n    )\n    return df","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:51:25.906269Z","iopub.execute_input":"2022-05-03T02:51:25.906549Z","iopub.status.idle":"2022-05-03T02:51:25.942762Z","shell.execute_reply.started":"2022-05-03T02:51:25.906522Z","shell.execute_reply":"2022-05-03T02:51:25.94205Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"baselines2 = add_prev_post_shift(baselines)\nbaselines2[\"UnixTimeMillis_prev_diff\"] = (\n    baselines2[\"UnixTimeMillis\"] - baselines2[\"UnixTimeMillis_shift1\"]\n)\nbaselines2[\"speed_calc\"] = (\n    baselines2[\"dist_prev\"] / baselines2[\"UnixTimeMillis_prev_diff\"]\n)\n\n\nsub_df = add_prev_post_shift(sub_df)\nsub_df[\"UnixTimeMillis_prev_diff\"] = (\n    sub_df[\"UnixTimeMillis\"] - sub_df[\"UnixTimeMillis_shift1\"]\n)\nsub_df[\"speed_calc\"] = sub_df[\"dist_prev\"] / sub_df[\"UnixTimeMillis_prev_diff\"]","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:51:28.859844Z","iopub.execute_input":"2022-05-03T02:51:28.86076Z","iopub.status.idle":"2022-05-03T02:51:29.419486Z","shell.execute_reply.started":"2022-05-03T02:51:28.860716Z","shell.execute_reply":"2022-05-03T02:51:29.418449Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Find when the cars are not moving:","metadata":{}},{"cell_type":"code","source":"baselines2.query(\"dist_prev < 5\")[\"dist_prev\"].plot(\n    kind=\"hist\", bins=50, title=\"Distribution of Speeds < 5\"\n)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:51:31.065526Z","iopub.execute_input":"2022-05-03T02:51:31.065802Z","iopub.status.idle":"2022-05-03T02:51:31.392975Z","shell.execute_reply.started":"2022-05-03T02:51:31.065774Z","shell.execute_reply":"2022-05-03T02:51:31.392112Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def do_postprocess(sub_df, thres=1):\n    sub = sub_df.copy()\n    for c, sub_stopped in sub.groupby(\"tripId\"):\n        sub_stopped = sub_stopped.loc[sub_stopped[\"dist_prev\"] < thres].copy()\n        sub_stopped[\"UnixTimeMillis_diff\"] = sub_stopped[\"UnixTimeMillis\"].diff()\n        sub_stopped[\"big_timeshift\"] = sub_stopped[\"UnixTimeMillis_diff\"] > 2_000\n        sub_stopped[\"time_group\"] = sub_stopped[\"big_timeshift\"].astype(\"int\").cumsum()\n\n        for stop_group, d in sub_stopped.groupby(\"time_group\"):\n            tstart, tstop = d[\"UnixTimeMillis\"].min(), d[\"UnixTimeMillis\"].max()\n            stopped_len = len(\n                sub.loc[\n                    (sub[\"UnixTimeMillis\"] >= tstart) & (sub[\"UnixTimeMillis\"] <= tstop)\n                ]\n            )\n            if stopped_len >= 20:\n                buffer = 800\n                latDegmean = sub.loc[\n                    (sub[\"UnixTimeMillis\"] >= (tstart - buffer))\n                    & (sub[\"UnixTimeMillis\"] <= (tstop + buffer))\n                ][\"LatitudeDegrees\"].mean()\n                lngDegmean = sub.loc[\n                    (sub[\"UnixTimeMillis\"] >= (tstart - buffer))\n                    & (sub[\"UnixTimeMillis\"] <= (tstop + buffer))\n                ][\"LongitudeDegrees\"].mean()\n                sub.loc[\n                    (sub[\"UnixTimeMillis\"] >= tstart)\n                    & (sub[\"UnixTimeMillis\"] <= tstop),\n                    \"LatitudeDegrees\",\n                ] = latDegmean\n                sub.loc[\n                    (sub[\"UnixTimeMillis\"] >= tstart)\n                    & (sub[\"UnixTimeMillis\"] <= tstop),\n                    \"LongitudeDegrees\",\n                ] = lngDegmean\n                sub.loc[\n                    (sub[\"UnixTimeMillis\"] >= tstart)\n                    & (sub[\"UnixTimeMillis\"] <= tstop),\n                    \"stopped\",\n                ] = True\n    sub[\"stopped\"] = sub[\"stopped\"].fillna(False)\n    return sub","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub = do_postprocess(sub_df, thres=1.5)","metadata":{"execution":{"iopub.status.busy":"2022-05-03T03:30:44.110903Z","iopub.execute_input":"2022-05-03T03:30:44.111764Z","iopub.status.idle":"2022-05-03T03:30:46.865870Z","shell.execute_reply.started":"2022-05-03T03:30:44.111716Z","shell.execute_reply":"2022-05-03T03:30:46.865059Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"visualize_traffic(\n    sub.query(\"LatitudeDegrees > 36\").sample(frac=0.2),\n    lat_col=\"LatitudeDegrees\",\n    lon_col=\"LongitudeDegrees\",\n    color_col=\"stopped\",\n    zoom=10,\n)","metadata":{"execution":{"iopub.status.busy":"2022-05-03T03:30:46.867504Z","iopub.execute_input":"2022-05-03T03:30:46.867943Z","iopub.status.idle":"2022-05-03T03:30:46.961258Z","shell.execute_reply.started":"2022-05-03T03:30:46.867898Z","shell.execute_reply":"2022-05-03T03:30:46.960418Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Check Post Processing on Training Dataset","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:59:17.420615Z","iopub.execute_input":"2022-05-03T02:59:17.421235Z","iopub.status.idle":"2022-05-03T02:59:17.467528Z","shell.execute_reply.started":"2022-05-03T02:59:17.421187Z","shell.execute_reply":"2022-05-03T02:59:17.466610Z"}}},{"cell_type":"code","source":"baselines2 = add_prev_post_shift(baselines)\nbaselines2[\"UnixTimeMillis_prev_diff\"] = (\n    baselines2[\"UnixTimeMillis\"] - baselines2[\"UnixTimeMillis_shift1\"]\n)\nbaselines2[\"speed_calc\"] = (\n    baselines2[\"dist_prev\"] / baselines2[\"UnixTimeMillis_prev_diff\"]\n)\n\nbaselines2_pp = do_postprocess(baselines2, thres=1.5)\nscores = []\nfor tripID in baselines2_pp[\"tripId\"].unique():\n    score = calc_score(tripID, baselines2_pp, gts)\n    scores.append(score)\n\nmean_score = np.mean(scores)\nprint(f\"mean_score = {mean_score:.3f}\")","metadata":{"execution":{"iopub.status.busy":"2022-05-03T03:30:47.526615Z","iopub.execute_input":"2022-05-03T03:30:47.526876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Save off results","metadata":{}},{"cell_type":"code","source":"sub.reset_index(drop=True)[ss.columns].to_csv(\"submission.csv\", index=False)","metadata":{"execution":{"iopub.status.busy":"2022-05-03T02:53:13.764641Z","iopub.execute_input":"2022-05-03T02:53:13.765329Z","iopub.status.idle":"2022-05-03T02:53:14.088916Z","shell.execute_reply.started":"2022-05-03T02:53:13.765293Z","shell.execute_reply":"2022-05-03T02:53:14.087987Z"},"trusted":true},"execution_count":null,"outputs":[]}]}