3_1_Filtering.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filtering and Preprocessing\n",
    "We have some GPS data of a vehicle driving around munich for a few weeks. The data is recorded with a smartphone app, that should detect if the vehicle is driving to trigger a record. This dosn't work every time, therefore we need to find and remove points where the vehicle is not driving. When driving in cities the GPS signal isn't that good all the time as well. So we need to find out when a signal is bad and remove this points.\n",
    "\n",
    "Let's start importing the modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:20:53.278861Z",
     "start_time": "2019-11-18T09:20:52.412499Z"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd, geopandas as gpd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from shapely import wkb\n",
    "import mplleaflet\n",
    "import sys\n",
    "import pickle\n",
    "# import additional functions\n",
    "%run ./data/custom_functions.py\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:20:54.860860Z",
     "start_time": "2019-11-18T09:20:54.044205Z"
    }
   },
   "outputs": [],
   "source": [
    "# get the raw data from the smartphone. This is prepared in a geopandas dataframe with all the gps Points converted as shapely points\n",
    "data = pd.read_pickle('data/gps_data.p')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:20:55.950812Z",
     "start_time": "2019-11-18T09:20:55.916882Z"
    }
   },
   "outputs": [],
   "source": [
    "# First display the data we got. \n",
    "display(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparing data for filtering\n",
    "The first task will be to find and remove points indicationg a error in measurement.    \n",
    "Find some aspects that can indicate these errors if you have a time series of GPS points as database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:20:57.664194Z",
     "start_time": "2019-11-18T09:20:57.659519Z"
    }
   },
   "outputs": [],
   "source": [
    "# Aspects incdication measurment error\n",
    "#<<solution>>\n",
    "# - Big change in Speed over short time period = Acceleration\n",
    "# - Big change of location/distance over short time period = Acceleration\n",
    "# - Bad GPS signal indication from GPS module (HDOP, VDOP)\n",
    "#<</solution>>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "No we precalculate some values between two conscutive data points, that we will need several times.    \n",
    "- Time difference\n",
    "- Distance\n",
    "- Speed difference = Acceleration\n",
    "- Speed from Distance difference (for comparison)\n",
    "\n",
    "For distance calculation between two points there aready is a function given in 'custom_functiions.py'. It Return the distance in meter.\n",
    "\n",
    "    lat_lon_2_m(latitude_1, longitude_1, latitude_2, longitude_2)    \n",
    "    \n",
    "Write a function to calculate the distance of all datapoints in the geopandas dataframe. \n",
    "Hint: \n",
    "\n",
    "    You can access the longitude/latitude with .x /.y on the geometry object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:22:50.705039Z",
     "start_time": "2019-11-18T09:22:50.694838Z"
    }
   },
   "outputs": [],
   "source": [
    "def calculate_distance_points(points):\n",
    "    # initialize distance array\n",
    "    distance = np.zeros(len(points))\n",
    "    # loop all points\n",
    "    for i in range(0, len(points) - 1):\n",
    "        # calculate distance between two consecutive points using lat_lon_2_km function\n",
    "        d = lat_lon_2_m(\n",
    "            points[i].x,\n",
    "            points[i].y,\n",
    "            points[i+1].x,\n",
    "            points[i+1].y,)\n",
    "        # append distance to array\n",
    "        distance[i + 1] = d\n",
    "    return distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:23:10.139704Z",
     "start_time": "2019-11-18T09:22:52.647226Z"
    }
   },
   "outputs": [],
   "source": [
    "# calculate time difference between points using pandas.shift() function\n",
    "data['time_diff'] = data.time-data.time.shift(1,fill_value = 0)\n",
    "data['distance_diff'] = calculate_distance_points(data.geom)\n",
    "# calculate distance between a point and the one following in meter\n",
    "data['acceleration'] = (data.speed-data.speed.shift())/(data.time_diff.dt.total_seconds())\n",
    "# for comparison get the speed and acceleration values from the positions and thime between them\n",
    "data['speed_calc'] = (data.distance_diff/(data.time_diff.dt.total_seconds()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:23:12.707637Z",
     "start_time": "2019-11-18T09:23:12.663756Z"
    }
   },
   "outputs": [],
   "source": [
    "# display rows where the time difference is bigger than one minute. The time difference needs to be a timeDelta object.\n",
    "display(data[data.time_diff > np.timedelta64(1, 'm')])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:23:23.573540Z",
     "start_time": "2019-11-18T09:23:15.307336Z"
    }
   },
   "outputs": [],
   "source": [
    "# compare speed from GPS vs speed claculated from distance/time\n",
    "def rsme(predictions,targets):\n",
    "    return np.sqrt(np.mean((predictions-targets)**2))\n",
    "\n",
    "# Inser your code here...\n",
    "#<<solution>>\n",
    "data['gps_vs_calc'] = data.speed - data.speed_calc\n",
    "data['speed_error'] = data.apply( lambda x: rsme(x.speed_calc, x.speed).astype(float), axis =1)\n",
    "#<</solution>>\n",
    "\n",
    "#<<solution>>\n",
    "data.gps_vs_calc.plot(ylim = (0,10), figsize = (30,7))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:23:32.601363Z",
     "start_time": "2019-11-18T09:23:32.120508Z"
    }
   },
   "outputs": [],
   "source": [
    "data.speed_error.plot(ylim = (0,20), figsize = (30,7))\n",
    "print(\"Total RSME: \" + str(data.speed_error.mean(axis = 0)))\n",
    "#<</solution>>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filter values with low speed\n",
    "Remove rows from data where speed is near zero. But let's visualize the current speed distribution by a KDE Plot.    \n",
    "Hint: \n",
    "\n",
    "    You can use the plotting function integrated in pandas for a quick plot. \n",
    "    Parameters are explained in the soure.\n",
    "    https://github.com/pandas-dev/pandas/blob/v0.25.1/pandas/plotting/_core.py#L504-L1533"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:24:09.802785Z",
     "start_time": "2019-11-18T09:24:06.173214Z"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "# Plot the KDE before removing the near zero values\n",
    "data['speed'].plot.kde(figsize = (30,7))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:24:18.676663Z",
     "start_time": "2019-11-18T09:24:18.265850Z"
    }
   },
   "outputs": [],
   "source": [
    "data.speed.plot(figsize = (30,7))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot some points with very low speed on a map. How do they behave?    \n",
    "Hint:\n",
    "\n",
    "    Use the prepared function 'map_folium' for plotting a dataframe.\n",
    "    map_folium(pos_data, colorvalue, geometry = 'geom', c_min = -1, c_max = -1, line = False, cm_type='jet')\n",
    "    take a subset of the data to make it printable (max. 6500 points) by conditional filtering (see link below for examples)\n",
    "    https://chrisalbon.com/python/data_wrangling/pandas_selecting_rows_on_conditions/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:24:37.514711Z",
     "start_time": "2019-11-18T09:24:33.408882Z"
    }
   },
   "outputs": [],
   "source": [
    "condition = (data.time > '2015-01-14 14:05:29.748') & (data.time < '2015-01-14 16:28:26.000') & (data.speed < 1)\n",
    "map_folium(data[condition], 'speed')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:24:43.605954Z",
     "start_time": "2019-11-18T09:24:43.387551Z"
    }
   },
   "outputs": [],
   "source": [
    "# Filter Data\n",
    "# Select a threshold value for clipping points by speed\n",
    "speed_thr = 0.1 # m/s\n",
    "# First calculate a filtered speed with a rolling window of 10 seconds. This prevents short stops to be removed\n",
    "data['speed_median']= data['speed'].rolling(3).median()\n",
    "\n",
    "#apply filter\n",
    "data = data[data.speed_median > speed_thr]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:24:46.365434Z",
     "start_time": "2019-11-18T09:24:45.982056Z"
    }
   },
   "outputs": [],
   "source": [
    "data.speed.plot(figsize = (30,7))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:24:50.466875Z",
     "start_time": "2019-11-18T09:24:49.804722Z"
    }
   },
   "outputs": [],
   "source": [
    "# Plot the KDE after removing the near zero values\n",
    "#<<solution>>\n",
    "data['speed'].plot.kde()\n",
    "#<</solution>>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Filter GPS outliers\n",
    "Filter out points with bad signal quality index (Horizontal Dilution of precision = HDOP)    \n",
    "\n",
    "| DOP Value | Rating    | Description                                                                                                                                                                              |   |   |\n",
    "|-----------|-----------|------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|---|---|\n",
    "| 1         | Ideal     | Highest possible confidence level to be used for applications demanding the highest possible precision at all times.                                                                     |   |   |\n",
    "| 1-2       | Excellent | At this confidence level, positional measurements are considered accurate enough to meet all but the most sensitive applications.                                                        |   |   |\n",
    "| 2-5       | Good      | Represents a level that marks the minimum appropriate for making accurate decisions. Positional measurements could be used to make reliable in-route navigation suggestions to the user. |   |   |\n",
    "| 5-10      | Moderate  | Positional measurements could be used for calculations, but the fix quality could still be improved. A more open view of the sky is recommended.                                         |   |   |\n",
    "| 10-20     | Fair      | Represents a low confidence level. Positional measurements should be discarded or used only to indicate a very rough estimate of the current location.                                   |   |   |\n",
    "| >20       | Poor      | At this level, measurements are inaccurate by as much as 300 meters with a 6-meter accurate device (50 DOP × 6 meters) and should be discarded.                                          |   |   |    \n",
    "https://en.wikipedia.org/wiki/Dilution_of_precision_(navigation)    \n",
    "\n",
    "Visualize the HDOP Values indicating a poor quality. Check out the locations on the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:24:58.043245Z",
     "start_time": "2019-11-18T09:24:57.665272Z"
    }
   },
   "outputs": [],
   "source": [
    "condition = data.hdop > 8\n",
    "map_folium(data[condition], 'hdop')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:25:03.063359Z",
     "start_time": "2019-11-18T09:25:03.028468Z"
    }
   },
   "outputs": [],
   "source": [
    "# display rows where the speed is very high\n",
    "#display(data[ #condition])\n",
    "#<<solution>>\n",
    "display(data[(data.acceleration > 5) & (data.time_diff < np.timedelta64(5, 's'))])\n",
    "#<</solution>>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:25:07.032031Z",
     "start_time": "2019-11-18T09:25:06.993240Z"
    }
   },
   "outputs": [],
   "source": [
    "# display rows where hdop is very high\n",
    "display(data[data.hdop > 20])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:25:23.239156Z",
     "start_time": "2019-11-18T09:25:23.226629Z"
    }
   },
   "outputs": [],
   "source": [
    "# Remove Points with very bad (high) HDOP\n",
    "#hdop_thr = select_a_trehold_for_bad_hdop_value\n",
    "#filter data\n",
    "#<<solution>>\n",
    "hdop_thr = 20\n",
    "data = data[data.hdop < hdop_thr]\n",
    "#<<solution>>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:25:51.368533Z",
     "start_time": "2019-11-18T09:25:48.626906Z"
    }
   },
   "outputs": [],
   "source": [
    "# clean up\n",
    "# first reset the index to comensate the removed datapoints\n",
    "data =data.reset_index(drop=True)\n",
    "# Recalculate all values to their neighbours usinf the functions we already had above. This hels us finding the start end of a trip\n",
    "# recalculate distance between a point and the one following in meter\n",
    "data['distance_diff'] = calculate_distance_points(data.geom)\n",
    "data['time_diff'] = data.time-data.time.shift(1,fill_value = 0)\n",
    "# for comparison get the speed and acceleration values from the positions and thime between them\n",
    "data['speed_calc'] = (data.distance_diff/(data.time_diff.dt.total_seconds()))\n",
    "data['acceleration_calc'] = (data.speed_calc-data.speed_calc.shift(1))/(data.time_diff.dt.total_seconds())\n",
    "#display(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analyze Tracks\n",
    "After some basic filtering we will analyze the raw bulk of data to find single Trips / Tracks. This can be done finding a time difference between two consecutive data points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:25:53.751459Z",
     "start_time": "2019-11-18T09:25:53.723482Z"
    }
   },
   "outputs": [],
   "source": [
    "# First define a new Dataframe for storing the tracks we extract\n",
    "tracks = pd.DataFrame(columns=['time_start','time_stop'])\n",
    "# For the Start of a track simply take the time difference we calculated before. If its biffger than x minutes (you can play around with this value), we will declare a new track.\n",
    "tracks['time_start'] = data[data.time_diff > np.timedelta64(5, 'm')].time\n",
    "\n",
    "# For the track ent we will take every point before the one of a new track.\n",
    "tracks['time_stop'] = data.iloc[tracks.index-1].time.sort_values().values\n",
    "# Lets calculate our track durations.\n",
    "tracks['duration'] = tracks.time_stop- tracks.time_start\n",
    "# We can filter short tracks by duration. This are the leftovers of our prevoiouse filtering\n",
    "tracks = tracks[tracks.duration > np.timedelta64(30, 's')]\n",
    "\n",
    "# reindex the Tracks to Track number instead of first Datapoint in original Dataframe\n",
    "tracks.reset_index(drop = True, inplace=True)\n",
    "# Add a date column for grouping tracks by date for export\n",
    "tracks['date'] = tracks.time_start.apply(lambda x: x.date())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:25:57.309991Z",
     "start_time": "2019-11-18T09:25:57.283952Z"
    }
   },
   "outputs": [],
   "source": [
    "display(tracks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-05T06:16:33.080641Z",
     "start_time": "2019-09-05T06:16:33.078141Z"
    }
   },
   "outputs": [],
   "source": [
    "# reset track view counter\n",
    "track_no = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-09-05T06:19:14.631656Z",
     "start_time": "2019-09-05T06:19:13.557549Z"
    }
   },
   "outputs": [],
   "source": [
    "# Plot a single Track on a map, show next track when executing again\n",
    "if track_no >= tracks.shape[0]: \n",
    "    track_no=0\n",
    "    print(\"Restart\")\n",
    "my_track = data[(data.time > tracks.time_start[track_no]) & (data.time < tracks.time_stop[track_no])]\n",
    "print(\"Track No. \" + str(track_no) + \" Duration: \" + str(tracks.duration[track_no]))\n",
    "track_no += 1\n",
    "map_folium(my_track, 'speed' , line=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:26:04.698758Z",
     "start_time": "2019-11-18T09:26:04.682445Z"
    }
   },
   "outputs": [],
   "source": [
    "# Display Data of this track for manual analysis\n",
    "display(my_track)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:26:08.386521Z",
     "start_time": "2019-11-18T09:26:07.968360Z"
    }
   },
   "outputs": [],
   "source": [
    "# Now we store the data in an array of dataframes, one for each day\n",
    "def generate_daily_driving(data, tracks):\n",
    "    trackdata =[]\n",
    "    for i, row in tracks.iterrows():\n",
    "        trackdata.insert(i,[tracks.iloc[i].date, [data[['time', 'geom', 'speed']][(data.time >= tracks.iloc[i].time_start)&(data.time <= tracks.iloc[i].time_stop)]]])\n",
    "    return trackdata\n",
    "    \n",
    "trackdata = generate_daily_driving(data, tracks)\n",
    "pickle.dump( trackdata, open( \"trackdata.p\", \"wb\" ))\n",
    "#load with: trackdata = pickle.load( open( \"trackdata.p\", \"rb\" ))\n",
    "#display(trackdata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extra Task 1\n",
    "Calculate the distance and average speed for each track by summing up the point-to-point distances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-11-18T09:26:12.415210Z",
     "start_time": "2019-11-18T09:26:12.357409Z"
    }
   },
   "outputs": [],
   "source": [
    "#<<solution>>\n",
    "for i, row in tracks.iterrows():\n",
    "        tracks['distance'].iloc[i] = data[(data.time > tracks.iloc[i].time_start)&(data.time <= tracks.iloc[i].time_stop)].distance_diff.sum()/1000\n",
    "tracks['avg_speed'] = tracks.distance/(tracks.duration.dt.total_seconds()/3600)\n",
    "display(tracks)\n",
    "#<</solution>>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Extra Task 2\n",
    " Try to optimize the filtering to get better results regarding track separation. Optimize the low-speed filtering to avoid deleting traffic jams."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}