{"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":"In this notebook, I introduce how to detect road from open street map and how to create grid points from road data. \nThe generated grid points may be used for \"snap to grid\". Please refer [original notebook](https://www.kaggle.com/robikscube/indoor-navigation-snap-to-grid-post-processing ) to know the detail of \"snap to grid\".\n\nActually, I haven't applied these grids to \"Snap to Grid\" well yet by some problem, and I'm still trying to figure it out.\nPlese comment if there are my mistakes or any idea.\n\nReference sites:  \nhttps://medium.com/@brendan_ward/how-to-leverage-geopandas-for-faster-snapping-of-points-to-lines-6113c94e59aa","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport glob\nimport os\nimport matplotlib.pyplot as plt\nfrom tqdm.notebook import tqdm\nfrom pathlib import Path\nimport plotly.express as px","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-06-09T13:12:06.36565Z","iopub.execute_input":"2021-06-09T13:12:06.366061Z","iopub.status.idle":"2021-06-09T13:12:07.743994Z","shell.execute_reply.started":"2021-06-09T13:12:06.365973Z","shell.execute_reply":"2021-06-09T13:12:07.742902Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_dir = Path(\"../input/google-smartphone-decimeter-challenge\")\ntrain_df = pd.read_csv(data_dir / \"baseline_locations_train.csv\")","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:12:08.532491Z","iopub.execute_input":"2021-06-09T13:12:08.532872Z","iopub.status.idle":"2021-06-09T13:12:08.877856Z","shell.execute_reply.started":"2021-06-09T13:12:08.532837Z","shell.execute_reply":"2021-06-09T13:12:08.876736Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get all ground truth dataframe\ngt_df = pd.DataFrame()\nfor (collection_name, phone_name), df in tqdm(train_df.groupby([\"collectionName\", \"phoneName\"])):\n    path = data_dir / f\"train/{collection_name}/{phone_name}/ground_truth.csv\"\n    df = pd.read_csv(path)  \n    gt_df = pd.concat([gt_df, df]).reset_index(drop=True)   \ngt_df.head()","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:12:33.905574Z","iopub.execute_input":"2021-06-09T13:12:33.905939Z","iopub.status.idle":"2021-06-09T13:12:35.02635Z","shell.execute_reply.started":"2021-06-09T13:12:33.905901Z","shell.execute_reply":"2021-06-09T13:12:35.02571Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = px.scatter_mapbox(gt_df,\n\n                    # Here, plotly gets, (x,y) coordinates\n                    lat=\"latDeg\",\n                    lon=\"lngDeg\",\n                    text='phoneName',\n\n                    #Here, plotly detects color of series\n                    color=\"collectionName\",\n                    labels=\"collectionName\",\n\n                    zoom=9,\n                    center={\"lat\":37.423576, \"lon\":-122.094132},\n                    height=600,\n                    width=800)\nfig.update_layout(mapbox_style='stamen-terrain')\nfig.update_layout(margin={\"r\": 0, \"t\": 0, \"l\": 0, \"b\": 0})\nfig.update_layout(title_text=\"GPS trafic\")\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:12:37.109868Z","iopub.execute_input":"2021-06-09T13:12:37.110363Z","iopub.status.idle":"2021-06-09T13:12:39.53995Z","shell.execute_reply.started":"2021-06-09T13:12:37.11033Z","shell.execute_reply":"2021-06-09T13:12:39.538319Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Target place\nLet's take '2021-04-29-US-SJC-2' as an example.","metadata":{}},{"cell_type":"code","source":"target_collection = '2021-04-29-US-SJC-2'\ntarget_gt_df = gt_df[gt_df[\"collectionName\"]==target_collection].reset_index(drop=True)\ntarget_collection","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:13:42.128139Z","iopub.execute_input":"2021-06-09T13:13:42.12865Z","iopub.status.idle":"2021-06-09T13:13:42.146214Z","shell.execute_reply.started":"2021-06-09T13:13:42.128615Z","shell.execute_reply":"2021-06-09T13:13:42.145118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = px.scatter_mapbox(target_gt_df,\n\n                    # Here, plotly gets, (x,y) coordinates\n                    lat=\"latDeg\",\n                    lon=\"lngDeg\",\n                    text='phoneName',\n\n                    #Here, plotly detects color of series\n                    color=\"collectionName\",\n                    labels=\"collectionName\",\n\n                    zoom=15,\n                    center={\"lat\":37.33351, \"lon\":-121.8906},\n                    height=600,\n                    width=800)\nfig.update_layout(mapbox_style='stamen-terrain')\nfig.update_layout(margin={\"r\": 0, \"t\": 0, \"l\": 0, \"b\": 0})\nfig.update_layout(title_text=\"GPS trafic\")\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:14:13.15847Z","iopub.execute_input":"2021-06-09T13:14:13.158837Z","iopub.status.idle":"2021-06-09T13:14:13.252183Z","shell.execute_reply.started":"2021-06-09T13:14:13.158806Z","shell.execute_reply":"2021-06-09T13:14:13.251441Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Import geographical library","metadata":{}},{"cell_type":"code","source":"!pip install osmnx momepy geopandas","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:15:14.575449Z","iopub.execute_input":"2021-06-09T13:15:14.576021Z","iopub.status.idle":"2021-06-09T13:15:27.877773Z","shell.execute_reply.started":"2021-06-09T13:15:14.575971Z","shell.execute_reply":"2021-06-09T13:15:27.87671Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from shapely.geometry import Point\nimport osmnx as ox\nimport momepy\nimport geopandas as gpd","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:15:27.879477Z","iopub.execute_input":"2021-06-09T13:15:27.879806Z","iopub.status.idle":"2021-06-09T13:15:29.462567Z","shell.execute_reply.started":"2021-06-09T13:15:27.879772Z","shell.execute_reply":"2021-06-09T13:15:29.461355Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# change pd.DataFrame -> gpd.GeoDataFrame\ntarget_gt_df[\"geometry\"] = [Point(p) for p in target_gt_df[[\"lngDeg\", \"latDeg\"]].to_numpy()]\ntarget_gt_gdf = gpd.GeoDataFrame(target_gt_df, geometry=target_gt_df[\"geometry\"])\ntarget_gt_gdf.head(5)","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:15:31.224499Z","iopub.execute_input":"2021-06-09T13:15:31.225007Z","iopub.status.idle":"2021-06-09T13:15:31.288833Z","shell.execute_reply.started":"2021-06-09T13:15:31.224973Z","shell.execute_reply":"2021-06-09T13:15:31.288131Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"target_gt_gdf.plot()","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:15:32.325412Z","iopub.execute_input":"2021-06-09T13:15:32.325763Z","iopub.status.idle":"2021-06-09T13:15:32.746123Z","shell.execute_reply.started":"2021-06-09T13:15:32.325733Z","shell.execute_reply":"2021-06-09T13:15:32.745033Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can get road data from open street map by creating bounding box. ","metadata":{}},{"cell_type":"code","source":"# get road data from open street map by osmnx\noffset = 0.1**5\nbbox = target_gt_gdf.bounds + [-offset, -offset, offset, offset]\neast = bbox[\"minx\"].min()\nwest = bbox[\"maxx\"].max()\nsouth = bbox[\"miny\"].min()\nnorth = bbox[\"maxy\"].max()\nG = ox.graph.graph_from_bbox(north, south, east, west, network_type='drive')","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:15:42.119619Z","iopub.execute_input":"2021-06-09T13:15:42.119953Z","iopub.status.idle":"2021-06-09T13:15:56.336407Z","shell.execute_reply.started":"2021-06-09T13:15:42.119924Z","shell.execute_reply":"2021-06-09T13:15:56.335362Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ox.folium.plot_graph_folium(G)","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:15:56.338758Z","iopub.execute_input":"2021-06-09T13:15:56.339191Z","iopub.status.idle":"2021-06-09T13:15:56.595793Z","shell.execute_reply.started":"2021-06-09T13:15:56.339146Z","shell.execute_reply":"2021-06-09T13:15:56.594748Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"nodes, edges = momepy.nx_to_gdf(G)","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:17:55.514978Z","iopub.execute_input":"2021-06-09T13:17:55.51533Z","iopub.status.idle":"2021-06-09T13:17:55.549062Z","shell.execute_reply.started":"2021-06-09T13:17:55.515287Z","shell.execute_reply":"2021-06-09T13:17:55.548228Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"nodes.head()","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:17:57.450649Z","iopub.execute_input":"2021-06-09T13:17:57.450982Z","iopub.status.idle":"2021-06-09T13:17:57.466115Z","shell.execute_reply.started":"2021-06-09T13:17:57.450953Z","shell.execute_reply":"2021-06-09T13:17:57.465014Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"edges.head()","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:17:58.429407Z","iopub.execute_input":"2021-06-09T13:17:58.429761Z","iopub.status.idle":"2021-06-09T13:17:58.449213Z","shell.execute_reply.started":"2021-06-09T13:17:58.429729Z","shell.execute_reply":"2021-06-09T13:17:58.448187Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In this notebook, I use only edges data.","metadata":{}},{"cell_type":"code","source":"edges.plot()","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:17:59.886838Z","iopub.execute_input":"2021-06-09T13:17:59.887192Z","iopub.status.idle":"2021-06-09T13:18:00.09174Z","shell.execute_reply.started":"2021-06-09T13:17:59.88716Z","shell.execute_reply":"2021-06-09T13:18:00.090562Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Since it still contains extra roads, we will leave only the relevant roads.","metadata":{}},{"cell_type":"code","source":"edges = edges.dropna(subset=[\"geometry\"]).reset_index(drop=True)\nhits = bbox.apply(lambda row: list(edges.sindex.intersection(row)), axis=1)\ntmp = pd.DataFrame({\n    # index of points table\n    \"pt_idx\": np.repeat(hits.index, hits.apply(len)),\n    # ordinal position of line - access via iloc later\n    \"line_i\": np.concatenate(hits.values)\n})\n# Join back to the lines on line_i; we use reset_index() to \n# give us the ordinal position of each line\ntmp = tmp.join(edges.reset_index(drop=True), on=\"line_i\")\n# Join back to the original points to get their geometry\n# rename the point geometry as \"point\"\ntmp = tmp.join(target_gt_gdf.geometry.rename(\"point\"), on=\"pt_idx\")\n# Convert back to a GeoDataFrame, so we can do spatial ops\ntmp = gpd.GeoDataFrame(tmp, geometry=\"geometry\", crs=target_gt_gdf.crs)","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:18:43.410879Z","iopub.execute_input":"2021-06-09T13:18:43.411208Z","iopub.status.idle":"2021-06-09T13:18:43.785628Z","shell.execute_reply.started":"2021-06-09T13:18:43.411177Z","shell.execute_reply":"2021-06-09T13:18:43.784648Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tmp.head()","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:18:49.21246Z","iopub.execute_input":"2021-06-09T13:18:49.212857Z","iopub.status.idle":"2021-06-09T13:18:49.238134Z","shell.execute_reply.started":"2021-06-09T13:18:49.212826Z","shell.execute_reply":"2021-06-09T13:18:49.237197Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Find closest road","metadata":{}},{"cell_type":"code","source":"tmp[\"snap_dist\"] = tmp.geometry.distance(gpd.GeoSeries(tmp.point))\n\n# Discard any lines that are greater than tolerance from points\ntolerance = 0.0005  \ntmp = tmp.loc[tmp.snap_dist <= tolerance]\n# Sort on ascending snap distance, so that closest goes to top\ntmp = tmp.sort_values(by=[\"snap_dist\"])\n\n# group by the index of the points and take the first, which is the\n# closest line \nclosest = tmp.groupby(\"pt_idx\").first()\n# construct a GeoDataFrame of the closest lines\nclosest = gpd.GeoDataFrame(closest, geometry=\"geometry\")\nclosest = closest.drop_duplicates(\"line_i\").reset_index(drop=True)","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:19:00.535729Z","iopub.execute_input":"2021-06-09T13:19:00.536064Z","iopub.status.idle":"2021-06-09T13:19:01.22324Z","shell.execute_reply.started":"2021-06-09T13:19:00.536034Z","shell.execute_reply":"2021-06-09T13:19:01.221768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"closest.plot()","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:19:04.781096Z","iopub.execute_input":"2021-06-09T13:19:04.781445Z","iopub.status.idle":"2021-06-09T13:19:04.943133Z","shell.execute_reply.started":"2021-06-09T13:19:04.781414Z","shell.execute_reply":"2021-06-09T13:19:04.941919Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"closest.head()","metadata":{"execution":{"iopub.status.busy":"2021-06-09T13:19:13.331925Z","iopub.execute_input":"2021-06-09T13:19:13.332265Z","iopub.status.idle":"2021-06-09T13:19:13.359461Z","shell.execute_reply.started":"2021-06-09T13:19:13.332234Z","shell.execute_reply":"2021-06-09T13:19:13.358353Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Then, we can obtain the road data corresponding to the target data. These features may be used for modeling.\n  \nNext, I generate grid points from road data.","metadata":{}},{"cell_type":"markdown","source":"## Generate road points","metadata":{}},{"cell_type":"code","source":"line_points_list = []\nsplit = 50  # param: number of split in each LineString\nfor dist in range(0, split, 1):\n    dist = dist/split\n    line_points = closest[\"geometry\"].interpolate(dist, normalized=True)\n    line_points_list.append(line_points)\nline_points = pd.concat(line_points_list).reset_index(drop=True)\nline_points = line_points.reset_index().rename(columns={0:\"geometry\"})\nline_points[\"lngDeg\"] = line_points[\"geometry\"].x\nline_points[\"latDeg\"] = line_points[\"geometry\"].y","metadata":{"execution":{"iopub.status.busy":"2021-06-09T11:22:12.734587Z","iopub.execute_input":"2021-06-09T11:22:12.735072Z","iopub.status.idle":"2021-06-09T11:22:12.776598Z","shell.execute_reply.started":"2021-06-09T11:22:12.735038Z","shell.execute_reply":"2021-06-09T11:22:12.775913Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = px.scatter_mapbox(line_points,\n\n                    # Here, plotly gets, (x,y) coordinates\n                    lat=\"latDeg\",\n                    lon=\"lngDeg\",\n\n                    zoom=15,\n                    center={\"lat\":37.33351, \"lon\":-121.8906},\n                    height=600,\n                    width=800)\nfig.update_layout(mapbox_style='stamen-terrain')\nfig.update_layout(margin={\"r\": 0, \"t\": 0, \"l\": 0, \"b\": 0})\nfig.update_layout(title_text=\"GPS trafic\")\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2021-06-09T11:24:17.68717Z","iopub.execute_input":"2021-06-09T11:24:17.687468Z","iopub.status.idle":"2021-06-09T11:24:17.755112Z","shell.execute_reply.started":"2021-06-09T11:24:17.68744Z","shell.execute_reply":"2021-06-09T11:24:17.75437Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I shared the road detection and creating grid points in this notebook.\n\nI hope it helps. Thanks!","metadata":{}}]}