{"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":"### tl;dr they're all really close\n\nThere has been some discussion about which distance formula is better for our data, haversine or Vincenty. My intuition was that it shouldn't matter for an area as relatively small as a race track. I decided to find out and also throw UTM Euclidean distances into the mix.\n\nMy conclusion is that it doesn't matter if the earth is elliptical, spherical, or flat - at least for my purposes. The difference in distances at each point on the Trakus index are rarely over 2cm from what I've seen. The cumulative difference amounts to 0.2% (0.002) of the total race distance. If I had to choose, I would pick Vincenty or UTM because the ellipsoid-based model is more accurate, and the calculation overhead is neglible.\n\nI used 6369km as the earth's radius for haversine distance. The distance is specific to New York according to the calculations at https://rechneronline.de/earth-radius/. That's a little bit different than the average radius of the earth and seems to yield the best match to Vincenty and UTM.\n\nHere's an example of what you might find for a single path from one horse during one race.","metadata":{}},{"cell_type":"code","source":"pip install hvplot utm vincenty  # thanks to shapital for the tip on UTM!","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2022-10-17T05:14:51.345086Z","iopub.execute_input":"2022-10-17T05:14:51.345422Z","iopub.status.idle":"2022-10-17T05:15:02.656990Z","shell.execute_reply.started":"2022-10-17T05:14:51.345353Z","shell.execute_reply":"2022-10-17T05:15:02.655628Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import hvplot.pandas  # noqa: F401\nimport numpy as np\nimport pandas as pd\nimport utm\nimport vincenty","metadata":{"execution":{"iopub.status.busy":"2022-10-17T05:15:02.660388Z","iopub.execute_input":"2022-10-17T05:15:02.660903Z","iopub.status.idle":"2022-10-17T05:15:05.236307Z","shell.execute_reply.started":"2022-10-17T05:15:02.660843Z","shell.execute_reply":"2022-10-17T05:15:05.235001Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def project_utm(df: pd.DataFrame) -> pd.DataFrame:\n    utms = utm.from_latlon(df.latitude.values, df.longitude.values)\n    df = df.assign(utm_easting=utms[0], utm_northing=utms[1])\n\n    utm_diff = np.diff(df[[\"utm_easting\", \"utm_northing\"]].to_numpy(), axis=0)\n    utm_distances = np.linalg.norm(utm_diff, axis=1)\n    df[\"distance_utm\"] = np.insert(utm_distances, 0, 0)\n    return df\n\n\n# Vectorized haversine: https://stackoverflow.com/a/44682708/5859466\ndef get_haversine(df: pd.DataFrame, earth_radius: float) -> pd.DataFrame:\n    coords_arr = np.deg2rad(df[[\"latitude\", \"longitude\"]].to_numpy())\n    lat_diff = np.diff(coords_arr[:, 0])\n    lon_diff = np.diff(coords_arr[:, 1])\n\n    a = (\n        np.sin(lat_diff / 2.0) ** 2\n        + np.cos(coords_arr[:-1, 0])\n        * np.cos(coords_arr[1:, 0])\n        * np.sin(lon_diff / 2) ** 2\n    )\n    c = 2 * np.arcsin(np.sqrt(a))\n    distances_m = earth_radius * c * 1000\n\n    df[\"distance_haversine\"] = np.insert(distances_m, 0, 0)\n    return df\n\n\ndef get_vincenty(df: pd.DataFrame) -> pd.DataFrame:  # taking the easy way here\n    vincenty_list = []\n    df = df.reset_index()\n    for row_tuple in df.itertuples():\n        if row_tuple.trakus_index == 1:\n            vincenty_list.append(0)\n        else:\n            distance = (\n                vincenty.vincenty(\n                    (\n                        df.loc[row_tuple.trakus_index - 2, \"latitude\"],\n                        df.loc[row_tuple.trakus_index - 2, \"longitude\"],\n                    ),\n                    (\n                        df.loc[row_tuple.trakus_index - 1, \"latitude\"],\n                        df.loc[row_tuple.trakus_index - 1, \"longitude\"],\n                    ),\n                )\n                * 1000\n            )\n            vincenty_list.append(distance)\n    df[\"distance_vincenty\"] = vincenty_list\n    return df\n\n\ndef get_diffs(df: pd.DataFrame) -> pd.DataFrame:\n    df[\"mean_distance\"] = df[\n        [\"distance_haversine\", \"distance_vincenty\", \"distance_utm\"]\n    ].mean(axis=1)\n    for col in [\"distance_haversine\", \"distance_vincenty\", \"distance_utm\"]:\n        df[col + \"_diff\"] = df.apply(lambda df: df[col] - df.mean_distance, axis=1) * 100\n    return df\n\n","metadata":{"execution":{"iopub.status.busy":"2022-10-17T05:15:05.238060Z","iopub.execute_input":"2022-10-17T05:15:05.238410Z","iopub.status.idle":"2022-10-17T05:15:05.254475Z","shell.execute_reply.started":"2022-10-17T05:15:05.238374Z","shell.execute_reply":"2022-10-17T05:15:05.253637Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Main\nEARTH_RADIUS_NY = 6369  # for haversine, https://rechneronline.de/earth-radius/\n\none_race_cols = [\"track_id\", \"race_date\", \"race_number\", \"program_number\"]\ntracking = pd.read_csv(\"../input/big-data-derby-2022/nyra_tracking_table.csv\", index_col=one_race_cols)\n\nhorse_choice = np.random.choice(tracking.index.nunique())\ntrack_one_path = (\n    tracking.sort_index()\n    .loc[tracking.index[horse_choice]]\n    .sort_values(\"trakus_index\")\n    .pipe(get_haversine, EARTH_RADIUS_NY)\n    .pipe(get_vincenty)\n    .pipe(project_utm)\n    .pipe(get_diffs)\n)\n\ntrack_one_path","metadata":{"execution":{"iopub.status.busy":"2022-10-17T05:15:05.257005Z","iopub.execute_input":"2022-10-17T05:15:05.257366Z","iopub.status.idle":"2022-10-17T05:15:21.104781Z","shell.execute_reply.started":"2022-10-17T05:15:05.257332Z","shell.execute_reply":"2022-10-17T05:15:21.103718Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"one_path_id = track_one_path.loc[0, one_race_cols].to_string()\n\nchart_opts = dict(\n    xlabel=\"Trakus Index\",\n    width=800,\n    height=600,\n    fontsize={\"ylabel\": \"20px\", \"xlabel\": \"20px\"},\n)\ndists = track_one_path.hvplot.line(\n    x=\"trakus_index\",\n    y=[\"distance_haversine\", \"distance_vincenty\", \"distance_utm\"],\n    alpha=0.4,\n    title=\"Distance between Trakus points \\n \" + one_path_id,\n    ylabel=\"Distance (m)\",\n    **chart_opts,\n).opts(legend_position=\"top\", show_grid=True,)\n\ndiffs = track_one_path.hvplot.line(\n    x=\"trakus_index\",\n    y=[\"distance_haversine_diff\", \"distance_vincenty_diff\", \"distance_utm_diff\"],\n    title=f\"Differences from mean \\n {one_path_id}\",\n    ylabel=\"Difference (cm)\",\n    **chart_opts,\n).opts(legend_position=\"top\", show_grid=True,)\n\n(dists + diffs).opts(shared_axes=False).cols(1)\n","metadata":{"execution":{"iopub.status.busy":"2022-10-17T05:15:21.106163Z","iopub.execute_input":"2022-10-17T05:15:21.106459Z","iopub.status.idle":"2022-10-17T05:15:21.565745Z","shell.execute_reply.started":"2022-10-17T05:15:21.106431Z","shell.execute_reply":"2022-10-17T05:15:21.564620Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Differences are relatively small between all 3. Haversine differences from the mean are greater because it assumes the earth is spherical. \n","metadata":{}},{"cell_type":"markdown","source":"### Sanity Check\n\nLet's do a rough sanity check to see if the difference between Vincenty and haversine is significant in terms of the overall race, including pre and post-race positions.","metadata":{"execution":{"iopub.status.busy":"2022-10-15T22:15:29.577602Z","iopub.execute_input":"2022-10-15T22:15:29.578118Z","iopub.status.idle":"2022-10-15T22:15:29.583718Z","shell.execute_reply.started":"2022-10-15T22:15:29.578080Z","shell.execute_reply":"2022-10-15T22:15:29.582512Z"}}},{"cell_type":"code","source":"track_one_path = track_one_path.assign(\n    clock_time=track_one_path.trakus_index * 0.25,\n    diff_hv=track_one_path.distance_vincenty-track_one_path.distance_haversine,\n    diff_cume_m=lambda df: df.diff_hv.cumsum(),\n    dist_cume_m=track_one_path.distance_haversine.cumsum(),\n    diff_cume_pct=lambda df: df.diff_cume_m / df.dist_cume_m * 100,\n)","metadata":{"execution":{"iopub.status.busy":"2022-10-17T05:15:21.567485Z","iopub.execute_input":"2022-10-17T05:15:21.567848Z","iopub.status.idle":"2022-10-17T05:15:21.578445Z","shell.execute_reply.started":"2022-10-17T05:15:21.567814Z","shell.execute_reply":"2022-10-17T05:15:21.577169Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"track_one_path.hvplot.line(\n    x=\"clock_time\",\n    y=\"diff_cume_pct\",\n    ylim = (-0.3, 0.3),\n#     yticks = list(np.arange(0,0.33,0.05)),\n    label=\"haversine\",\n    title=\"Cumulative difference as a percent of travel \\n \" + one_path_id,\n    xlabel=\"Time (s)\",\n    ylabel=\"Difference percentage (%)\",\n    yformatter=\"%.3f\",\n    c=\"#1170AA\",\n    grid=True,\n    width=800,\n    height=600,\n    fontsize={\"ylabel\": \"20px\", \"xlabel\": \"20px\"},\n)","metadata":{"execution":{"iopub.status.busy":"2022-10-17T05:15:53.788530Z","iopub.execute_input":"2022-10-17T05:15:53.788978Z","iopub.status.idle":"2022-10-17T05:15:53.925580Z","shell.execute_reply.started":"2022-10-17T05:15:53.788941Z","shell.execute_reply":"2022-10-17T05:15:53.924776Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"So even though the cumulative difference between Vincenty/UTM x_y and haversine might be 1-2 meters, it's a very small percentage of the cumulative race distance at any given time. Your application may find the difference to be meaningful in which case I hope you now have a better sense of the differences.\n\nThanks and best wishes for the competition!","metadata":{}}]}