{"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\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns","metadata":{"execution":{"iopub.status.busy":"2022-10-24T11:47:02.829978Z","iopub.execute_input":"2022-10-24T11:47:02.830397Z","iopub.status.idle":"2022-10-24T11:47:02.835765Z","shell.execute_reply.started":"2022-10-24T11:47:02.830359Z","shell.execute_reply":"2022-10-24T11:47:02.834681Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"This notebook explains calculating speed, accleration, \n\nand menger curvature using latitude and longitude.\n\nThese feature will be an important feature when considering the race strategy.","metadata":{}},{"cell_type":"code","source":"df = pd.read_csv(\"../input/big-data-derby-2022/nyra_tracking_table.csv\")\ndf = df.sort_values(by=['track_id','race_date','race_number',\n                        'program_number','trakus_index']).reset_index(drop=True)\ndf = df.groupby(['track_id','race_date','race_number',\n                 'program_number'])[['latitude','longitude']].agg(lambda x: x.tolist())\ndf","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-10-24T11:47:03.197379Z","iopub.execute_input":"2022-10-24T11:47:03.198090Z","iopub.status.idle":"2022-10-24T11:47:15.213937Z","shell.execute_reply.started":"2022-10-24T11:47:03.198031Z","shell.execute_reply":"2022-10-24T11:47:15.212893Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# sample_record\nN = len(df)\nsample_idx = np.random.choice(N)\nsample_record = df.iloc[sample_idx]\nlnglat = np.vstack([sample_record['longitude'], sample_record['latitude']]).T\n\ndisplay(df.iloc[sample_idx])\nplt.plot(lnglat[:, 0], lnglat[:, 1])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-24T11:47:15.217110Z","iopub.execute_input":"2022-10-24T11:47:15.217655Z","iopub.status.idle":"2022-10-24T11:47:15.426279Z","shell.execute_reply.started":"2022-10-24T11:47:15.217619Z","shell.execute_reply":"2022-10-24T11:47:15.424934Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Speed","metadata":{}},{"cell_type":"code","source":"def calculate_speed(lnglat, dt=0.25):\n    lnglat_diff = lnglat[2:] - lnglat[:-2]\n    mag = np.linalg.norm(lnglat_diff, axis=1)\n    speed = mag/(2*dt)\n    speed = np.concatenate([np.array([speed[0]/2]),\n                            speed,\n                            np.array([2*speed[-1] - speed[-2]])])\n    speed = 3600*speed\n    return speed","metadata":{"execution":{"iopub.status.busy":"2022-10-24T11:47:15.427759Z","iopub.execute_input":"2022-10-24T11:47:15.428124Z","iopub.status.idle":"2022-10-24T11:47:15.435506Z","shell.execute_reply.started":"2022-10-24T11:47:15.428090Z","shell.execute_reply":"2022-10-24T11:47:15.434143Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"speed = calculate_speed(lnglat)\nplt.plot(speed)\nplt.title('speed  (lnglat_unit/hour)')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-24T11:47:15.438261Z","iopub.execute_input":"2022-10-24T11:47:15.438628Z","iopub.status.idle":"2022-10-24T11:47:15.617657Z","shell.execute_reply.started":"2022-10-24T11:47:15.438586Z","shell.execute_reply":"2022-10-24T11:47:15.616251Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Accelration","metadata":{}},{"cell_type":"code","source":"def calculate_acc(lnglat, dt=0.25):\n    lnglat_diff = (-lnglat[4:] + 16*lnglat[3:-1] - 30*lnglat[2:-2] \\\n                    + 16*lnglat[1:-3] - lnglat[:-4])\n    mag = np.linalg.norm(lnglat_diff, axis=1)\n    acc = mag/(12*dt*dt)\n    \n    v1 = np.linalg.norm(lnglat[2] - 2*lnglat[1] + lnglat[0])/(dt*dt)\n    v0 = 2*v1 - acc[0]\n    v2 =  np.linalg.norm(lnglat[-1] - 2*lnglat[-2] + lnglat[-3])/(dt*dt)\n    v3 = 2*v2 - acc[-1]\n    acc = np.concatenate([np.array([v0, v1]),\n                          acc,\n                          np.array([v2, v3])])\n    acc = 3600*acc\n    return acc","metadata":{"execution":{"iopub.status.busy":"2022-10-24T11:47:15.619076Z","iopub.execute_input":"2022-10-24T11:47:15.619453Z","iopub.status.idle":"2022-10-24T11:47:15.628804Z","shell.execute_reply.started":"2022-10-24T11:47:15.619413Z","shell.execute_reply":"2022-10-24T11:47:15.627756Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"acc = calculate_acc(lnglat)\nplt.plot(acc)\nplt.title('acceleration  (lnglat_unit/hour/sec)')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-24T11:47:15.630234Z","iopub.execute_input":"2022-10-24T11:47:15.630725Z","iopub.status.idle":"2022-10-24T11:47:15.829334Z","shell.execute_reply.started":"2022-10-24T11:47:15.630688Z","shell.execute_reply":"2022-10-24T11:47:15.826741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Menger curvature\n- https://en.wikipedia.org/wiki/Menger_curvature","metadata":{}},{"cell_type":"code","source":"def calculate_menger_curvature(lnglat, step=20):\n    p0 = lnglat[:-2*step]\n    p1 = lnglat[step:-step]\n    p2 = lnglat[2*step:]\n    \n    vec1 = p0 - p1\n    vec2 = p2 - p1\n    inp_value = np.abs(np.einsum('ij, ij->i', vec1, vec2))\n    d1 = np.linalg.norm(vec1, axis=1)\n    d2 = np.linalg.norm(vec2, axis=1)\n    d3 = np.linalg.norm(p2-p0, axis=1)\n    d = d1*d2*d3\n\n    curv = np.divide(inp_value, d, out=np.zeros_like(inp_value), where=d!=0)\n    curv = np.pad(curv, (step,step), mode='edge')\n    return curv","metadata":{"execution":{"iopub.status.busy":"2022-10-24T11:47:15.830780Z","iopub.execute_input":"2022-10-24T11:47:15.831313Z","iopub.status.idle":"2022-10-24T11:47:15.841487Z","shell.execute_reply.started":"2022-10-24T11:47:15.831280Z","shell.execute_reply":"2022-10-24T11:47:15.839908Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"curv = calculate_menger_curvature(lnglat)\nplt.plot(curv)\nplt.title('curvature (trakus_index_step = 20)')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-24T11:47:15.843910Z","iopub.execute_input":"2022-10-24T11:47:15.844262Z","iopub.status.idle":"2022-10-24T11:47:16.051731Z","shell.execute_reply.started":"2022-10-24T11:47:15.844230Z","shell.execute_reply":"2022-10-24T11:47:16.050512Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualization","metadata":{}},{"cell_type":"code","source":"def get_marker_size(feat):\n    percentile = np.zeros_like(feat)\n    percentile[np.argsort(feat)] = np.arange(1, len(feat)+1)/len(feat)\n    \n    marker_size = 50* np.ones_like(feat)\n    marker_size[percentile >= 0.75] = 200\n    marker_size[(percentile > 0.5)&(percentile <= 0.75)] = 100\n    marker_size[(percentile > 0.25)&(percentile <= 0.5)] = 50\n    marker_size[percentile <= 0.25] = 15\n    \n    return marker_size","metadata":{"execution":{"iopub.status.busy":"2022-10-24T11:47:16.053663Z","iopub.execute_input":"2022-10-24T11:47:16.054024Z","iopub.status.idle":"2022-10-24T11:47:16.060334Z","shell.execute_reply.started":"2022-10-24T11:47:16.053993Z","shell.execute_reply":"2022-10-24T11:47:16.059522Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"speed = calculate_speed(lnglat)\nacc = calculate_acc(lnglat)\ncurv = calculate_menger_curvature(lnglat)\n\nspeed_marker_size = get_marker_size(speed)\nacc_marker_size = get_marker_size(acc)\n\n#plot\nplt.figure(figsize=(16,6))\nplt.subplot(1,2,1)\nplt.title('latitude and longitude (lnglat scale). Speed (size), Curvature (color)')\nplt.scatter(lnglat[::5, 0], lnglat[::5, 1], s=speed_marker_size[::5], c=curv[::5])\nplt.subplot(1,2,2)\nplt.title('latitude and longitude (lnglat scale). Accleration (size), Curvature (color)')\nplt.scatter(lnglat[::5, 0], lnglat[::5, 1], s=acc_marker_size[::5], c=curv[::5])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-24T11:47:16.062923Z","iopub.execute_input":"2022-10-24T11:47:16.063864Z","iopub.status.idle":"2022-10-24T11:47:16.406166Z","shell.execute_reply.started":"2022-10-24T11:47:16.063830Z","shell.execute_reply":"2022-10-24T11:47:16.405058Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sns.histplot(x=speed, y=curv)\nplt.xlabel('speed')\nplt.ylabel('menger curvature')\nplt.title('2d histogram of speed and menger curvature')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-10-24T11:47:16.407492Z","iopub.execute_input":"2022-10-24T11:47:16.407941Z","iopub.status.idle":"2022-10-24T11:47:16.607897Z","shell.execute_reply.started":"2022-10-24T11:47:16.407907Z","shell.execute_reply":"2022-10-24T11:47:16.606770Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The curvature value mainly decreases in the corner section.\n\nI'm using Curvature to check if the horse is slowing down in a corner.\n","metadata":{}},{"cell_type":"markdown","source":"This notebook covered how to calculate features (speed, acceleration, and curvature).\n\nAfter calculating a few more features and indexes, I plan to conduct a case study. \n\nThanks.","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}