2_feature_extraction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c8f57f1f",
   "metadata": {},
   "source": [
    "# Feature Extraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "63456ce6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from tqdm import tqdm\n",
    "from joblib import Parallel, delayed\n",
    "\n",
    "import os\n",
    "os.environ['USE_PYGEOS'] = '0'\n",
    "\n",
    "import geopandas as gpd\n",
    "import trackintel as ti\n",
    "from shapely.geometry import Point\n",
    "from scipy.spatial import KDTree\n",
    "import statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "43157503",
   "metadata": {},
   "outputs": [],
   "source": [
    "# save\n",
    "if not os.path.exists(\"fig\"):\n",
    "    os.makedirs(\"fig\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a09576fd",
   "metadata": {},
   "source": [
    "## Load Preprocessed Data\n",
    "\n",
    "\n",
    "Load stages from 1_preprocess script, example data format:\n",
    "\n",
    "|index | id | user_id | started_at    | finished_at | mode |geom | length | dur_s | user_t_quality | recording_interval | speed_av |\n",
    "|----|----|----|----|----|----|----|----|----|----|----|----|\n",
    "| 0 | 29 | 16 | 2016-04-14 18:17:23.000+00:00 | 2016-04-14 18:20:57.000+00:00 | Mode::Walk | \"LINESTRING (8.0125 47.5064, 8.0129 47.50666)\" | 235.018 | 213.999 | 0.9182 | 7.1333|  3.953 |\n",
    "| 1 | 32 | 17 | 2017-04-30 11:05:40.000+00:00 | 2017-04-30 11:19:45.000+00:00 | Mode::Walk | \"LINESTRING (6.6961 46.8441, 6.6962 46.8446)\"  | 1035.4314 | 845.063 | 0.859 | 16.569 | 4.410 |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7c0e7547",
   "metadata": {},
   "outputs": [],
   "source": [
    "tl = ti.io.file.read_triplegs_csv(\n",
    "    os.path.join(\"data\", \"triplegs_preprocessed.csv\"), index_col=\"index\", geom_col=\"geom\", crs=4326\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "61c4d1ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "# transform coordinates from WGS84 to LV95\n",
    "tl[\"geom_LV95\"] = tl[\"geom\"]\n",
    "tl.set_geometry(\"geom_LV95\", inplace=True)\n",
    "tl.to_crs(2056, inplace=True)\n",
    "\n",
    "tl.rename(columns={\"dur_s\": \"dur\"}, inplace=True)\n",
    "\n",
    "# add list of intervals to each tripleg\n",
    "tl[\"intervals\"] = tl[\"geom_LV95\"].apply(np.array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "e603260a",
   "metadata": {},
   "outputs": [],
   "source": [
    "tl['mode'] = tl['mode'].replace({\n",
    "    'Mode::Car': 'Car',\n",
    "    'Mode::Ecar': 'Car',\n",
    "    'Mode::Bicycle': 'Bicycle',\n",
    "    'Mode::Ebicycle': 'Bicycle',\n",
    "    'Mode::Walk': 'Walk',\n",
    "    'Mode::Train': 'Train',\n",
    "    'Mode::Tram': 'Tram',\n",
    "    'Mode::Bus': 'Bus',\n",
    "    'Mode::Boat': 'Boat'\n",
    "})\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6424a8c0",
   "metadata": {},
   "source": [
    "## Motion Feature Extraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a6c993b1",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Bar: 100%|██████████| 365307/365307 [00:17<00:00, 21356.44it/s]\n",
      "Bar: 100%|██████████| 365307/365307 [03:55<00:00, 1548.61it/s]\n",
      "Bar: 100%|██████████| 365307/365307 [00:07<00:00, 52132.71it/s]\n",
      "Bar: 100%|██████████| 365307/365307 [00:14<00:00, 24952.95it/s]\n",
      "Bar: 100%|██████████| 365307/365307 [00:07<00:00, 52021.63it/s]\n"
     ]
    }
   ],
   "source": [
    "def _calculate_length(points):\n",
    "    \"\"\"Calculate length of intervals\"\"\"\n",
    "    lengths = np.linalg.norm(np.array(points.coords)[1:] - np.array(points.coords)[:-1], axis=1)\n",
    "    return lengths\n",
    "\n",
    "\n",
    "def _calculate_bearing(points):\n",
    "    \"\"\"Calculate bearing (heading direction) of intervals\"\"\"\n",
    "    interval = np.array(points.coords)[1:] - np.array(points.coords)[:-1]\n",
    "    return np.arctan2(interval[:, 0], interval[:, 1])\n",
    "\n",
    "\n",
    "def _calculate_bearing_rate(bearings):\n",
    "    \"\"\"Calculate bearing rate (change in heading direction) between intervals\"\"\"\n",
    "    diff = bearings[1:] - bearings[:-1]\n",
    "    return np.minimum(diff % (2 * np.pi), ((2 * np.pi) - diff) % (2 * np.pi))\n",
    "\n",
    "\n",
    "def _calculate_speed(row):\n",
    "    \"\"\"Calculate speed of intervals\"\"\"\n",
    "    return np.array([length / row[\"recording_interval\"] for length in row[\"lengths_of_intervals\"]])\n",
    "\n",
    "\n",
    "def _calculate_acceleration(row):\n",
    "    \"\"\"Calculate acceleration between intervals\"\"\"\n",
    "    return abs((row[\"speeds_of_intervals\"][1:] - row[\"speeds_of_intervals\"][:-1]) / row[\"recording_interval\"])\n",
    "\n",
    "\n",
    "\n",
    "def motion_features(tl):\n",
    "    # calculate properties of every interval\n",
    "    tqdm.pandas(desc=\"Bar\")\n",
    "    \n",
    "    # calculate lengths [m]\n",
    "    tl[\"lengths_of_intervals\"] = tl[\"intervals\"].progress_apply(_calculate_length)\n",
    "    # calculate speeds [m/s]\n",
    "    tl[\"speeds_of_intervals\"] = tl.progress_apply(_calculate_speed, axis=1)\n",
    "    # calculate accelerations [m/(s^2)]\n",
    "    tl[\"acc_of_intervals\"] = tl.progress_apply(_calculate_acceleration, axis=1)\n",
    "    # calculate bearings [rad] within (-pi,pi]\n",
    "    tl[\"bearings_of_intervals\"] = tl[\"intervals\"].progress_apply(_calculate_bearing)\n",
    "    # calculate absolute value of minimal angle between bearings [rad] within [0,pi]\n",
    "    tl[\"br_of_intervals\"] = tl[\"bearings_of_intervals\"].progress_apply(_calculate_bearing_rate)\n",
    "\n",
    "    # calculate 85th percentile of speed based on speeds of elementar intervals\n",
    "    tl[\"speed_85th\"] = [np.percentile(speeds, 85) * 3.6 for speeds in tl[\"speeds_of_intervals\"]]\n",
    "\n",
    "    tl[\"speed_90th\"] = [np.percentile(speeds, 90) * 3.6 for speeds in tl[\"speeds_of_intervals\"]]\n",
    "    tl[\"speed_95th\"] = [np.percentile(speeds, 95) * 3.6 for speeds in tl[\"speeds_of_intervals\"]]\n",
    "    \n",
    "    # calculate average acceleraton based on acceleration of elementar intervals\n",
    "    tl[\"acc_av\"] = [np.average(acc) for acc in tl[\"acc_of_intervals\"]]\n",
    "    tl[\"acc_85th\"] = [np.percentile(acc, 85) for acc in tl[\"acc_of_intervals\"]]\n",
    "\n",
    "    # calculate average bearing rate based on bearing rates of elementar intervals\n",
    "    tl[\"br_av\"] = [np.average(br) for br in tl[\"br_of_intervals\"]]\n",
    "    tl[\"br_85th\"] = [np.percentile(br, 85) for br in tl[\"br_of_intervals\"]]\n",
    "\n",
    "    return tl\n",
    "\n",
    "tl = motion_features(tl)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a175158c",
   "metadata": {},
   "source": [
    "## GeoSpatial Context Feature Extraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "eea6ef47",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load OSM point layers\n",
    "OSM_transport = gpd.read_file(os.path.join(\"data\", \"osm_switzerland_17-01-01\", \"gis_osm_transport_free_1.shp\"))\n",
    "OSM_traffic = gpd.read_file(os.path.join(\"data\", \"osm_switzerland_17-01-01\", \"gis_osm_traffic_free_1.shp\"))\n",
    "OSM_poi = gpd.read_file(os.path.join(\"data\", \"osm_switzerland_17-01-01\", \"gis_osm_pois_free_1.shp\"))\n",
    "OSM_pofw = gpd.read_file(os.path.join(\"data\", \"osm_switzerland_17-01-01\", \"gis_osm_pofw_free_1.shp\"))\n",
    "\n",
    "# load OSM line layers\n",
    "OSM_roads = gpd.read_file(os.path.join(\"data\", \"osm_switzerland_17-01-01\", \"gis_osm_roads_free_1.shp\"))\n",
    "OSM_railways = gpd.read_file(os.path.join(\"data\", \"osm_switzerland_17-01-01\", \"gis_osm_railways_free_1.shp\"))\n",
    "\n",
    "# load OSM polygon layers\n",
    "OSM_water = gpd.read_file(os.path.join(\"data\", \"osm_switzerland_latest\", \"gis_osm_water_a_free_1.shp\"))\n",
    "OSM_landuse = gpd.read_file(os.path.join(\"data\", \"osm_switzerland_latest\", \"gis_osm_landuse_a_free_1.shp\"))\n",
    "\n",
    "# load swissTLMregio rivers\n",
    "v25_rivers = gpd.read_file(os.path.join(\"data\", \"oberflachengewasser_vector25\", \"Typisierung_LV95\", \"FGT.shp\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "00aecf60",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _apply_parallel(df, func, n=-1, **kwargs):\n",
    "    \"\"\"parallel apply for spending up.\"\"\"\n",
    "    if n is None:\n",
    "        n = -1\n",
    "\n",
    "    cpunum = 20\n",
    "    dflength = len(df)\n",
    "    if dflength < cpunum:\n",
    "        spnum = dflength\n",
    "    if n < 0:\n",
    "        spnum = cpunum + n + 1\n",
    "    else:\n",
    "        spnum = n or 1\n",
    "\n",
    "    sp = list(range(dflength)[:: int(dflength / spnum + 0.5)])\n",
    "    sp.append(dflength)\n",
    "    slice_gen = (slice(*idx) for idx in zip(sp[:-1], sp[1:]))\n",
    "\n",
    "    results = Parallel(n_jobs=n, verbose=0)(delayed(func)(df.iloc[slc], **kwargs) for slc in slice_gen)\n",
    "    return pd.concat(results)\n",
    "\n",
    "def _create_index(df, geom_type):\n",
    "    \"\"\"Create spatial index. For point features return sindex and kdtree, otherwise only sindex.\"\"\"\n",
    "    # transform from WGS84 to LV95\n",
    "    df_LV95 = df.to_crs(2056)\n",
    "    # create sindex\n",
    "    df_LV95_sindex = df_LV95.sindex\n",
    "    if geom_type != \"point\":\n",
    "        return df_LV95, df_LV95_sindex\n",
    "    else:\n",
    "        # create kdtree\n",
    "        df_LV95_kdtree = KDTree(np.squeeze(np.array([point.coords for point in df_LV95.geometry])))\n",
    "        return df_LV95, df_LV95_sindex, df_LV95_kdtree"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7999ac79",
   "metadata": {},
   "source": [
    "### Extract Spatial Context Information Types and Create Spatial Index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "879564c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract railway stations\n",
    "rw_stations = OSM_transport[\n",
    "    ((OSM_transport[\"fclass\"] == \"railway_station\") | (OSM_transport[\"fclass\"] == \"railway_halt\"))\n",
    "].copy()\n",
    "_, _, rw_stations_LV95_kdtree = _create_index(rw_stations, \"point\")\n",
    "\n",
    "# extract tram stops\n",
    "tram_stops = OSM_transport[(OSM_transport[\"fclass\"] == \"tram_stop\")].copy()\n",
    "_, _, tram_stops_LV95_kdtree = _create_index(tram_stops, \"point\")\n",
    "\n",
    "# extract bus stops\n",
    "bus_stops = OSM_transport[\n",
    "    ((OSM_transport[\"fclass\"] == \"bus_stop\") | (OSM_transport[\"fclass\"] == \"bus_station\"))\n",
    "].copy()\n",
    "_, _, bus_stops_LV95_kdtree = _create_index(bus_stops, \"point\")\n",
    "\n",
    "# extract car parkings\n",
    "car_parking_types = [\"parking\", \"parking_site\", \"parking_multistorey\", \"parking_underground\"]\n",
    "car_parkings = OSM_traffic.loc[OSM_traffic[\"fclass\"].isin(car_parking_types)].copy()\n",
    "_, _, car_parkings_LV95_kdtree = _create_index(car_parkings, \"point\")\n",
    "\n",
    "# extract bycicle parkings\n",
    "bicycle_parkings = OSM_traffic[OSM_traffic[\"fclass\"] == \"parking_bicycle\"].copy()\n",
    "_, _, bicycle_parkings_LV95_kdtree = _create_index(bicycle_parkings, \"point\")\n",
    "\n",
    "# extract landing stages\n",
    "ferry_terminals = OSM_transport[OSM_transport[\"fclass\"] == \"ferry_terminal\"]\n",
    "harbours = OSM_traffic[OSM_traffic[\"fclass\"] == \"marina\"]\n",
    "landing_stages = pd.concat([ferry_terminals, harbours])\n",
    "landing_stages.reset_index(inplace=True)\n",
    "_, _, landing_stages_LV95_kdtree = _create_index(landing_stages, \"point\")\n",
    "\n",
    "# extract points of interest\n",
    "non_poi_types = [\n",
    "    \"fire_station\",\n",
    "    \"prison\",\n",
    "    \"hunting_stand\",\n",
    "    \"camera_surveillance\",\n",
    "    \"emergency_phone\",\n",
    "    \"emergency_access\",\n",
    "    \"tower_comms\",\n",
    "    \"water_tower\",\n",
    "    \"tower_observation\",\n",
    "    \"waste_basket\",\n",
    "    \"wastewater_plant\",\n",
    "    \"water_works\",\n",
    "]\n",
    "poi = OSM_poi.loc[~OSM_poi[\"fclass\"].isin(non_poi_types)].copy()\n",
    "poi = pd.concat([poi, OSM_pofw])\n",
    "poi.reset_index(inplace=True)\n",
    "_, _, poi_LV95_kdtree = _create_index(poi, \"point\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "6692cf1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# extract railway network\n",
    "rw_types = [\"rail\", \"light_rail\", \"narrow_gauge\", \"rack\"]\n",
    "rw_network = OSM_railways.loc[OSM_railways[\"fclass\"].isin(rw_types)].copy()\n",
    "_, rw_network_LV95_sindex = _create_index(rw_network, \"line\")\n",
    "\n",
    "# extract tram network\n",
    "tram_types = [\"tram\", \"subway\"]\n",
    "tram_network = OSM_railways.loc[OSM_railways[\"fclass\"].isin(tram_types)].copy()\n",
    "_, tram_network_LV95_sindex = _create_index(tram_network, \"line\")\n",
    "\n",
    "# extract road network\n",
    "road_types = [\n",
    "    \"motorway\",\n",
    "    \"trunk\",\n",
    "    \"primary\",\n",
    "    \"secondary\",\n",
    "    \"tertiary\",\n",
    "    \"unclassified\",\n",
    "    \"residential\",\n",
    "    \"living_street\",\n",
    "    \"motorway_link\",\n",
    "    \"trunk_link\",\n",
    "    \"primary_link\",\n",
    "    \"secondary_link\",\n",
    "    \"tertiary_link\",\n",
    "]\n",
    "road_network = OSM_roads.loc[OSM_roads[\"fclass\"].isin(road_types)].copy()\n",
    "_, road_network_LV95_sindex = _create_index(road_network, \"line\")\n",
    "\n",
    "# extract pedestrian and bicycle network\n",
    "pb_types = [\n",
    "    \"primary\",\n",
    "    \"secondary\",\n",
    "    \"tertiary\",\n",
    "    \"unclassified\",\n",
    "    \"residential\",\n",
    "    \"living_street\",\n",
    "    \"pedestrian\",\n",
    "    \"service\",\n",
    "    \"track\",\n",
    "    \"track_grade1\",\n",
    "    \"track_grade2\",\n",
    "    \"track_grade3\",\n",
    "    \"track_grade4\",\n",
    "    \"track_grade5\",\n",
    "    \"cycleway\",\n",
    "    \"footway\",\n",
    "    \"path\",\n",
    "    \"steps\",\n",
    "]\n",
    "pb_network = OSM_roads.loc[OSM_roads[\"fclass\"].isin(pb_types)].copy()\n",
    "_, pb_network_LV95_sindex = _create_index(pb_network, \"line\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5d83b014",
   "metadata": {},
   "outputs": [],
   "source": [
    "road_network.to_file(os.path.join(\"data\", \"road_network_visualize.shp\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "332b4ad4",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# extract large water bodies\n",
    "# select lakes and large rivers\n",
    "lakes = OSM_water[((OSM_water[\"fclass\"] == \"water\") | (OSM_water[\"fclass\"] == \"reservoir\"))].copy()\n",
    "rivers = v25_rivers[v25_rivers[\"GROSSERFLU\"] != \"NA\"].copy()\n",
    "# filter lakes based on area and simplify borders with Douglas Peuker\n",
    "lakes.to_crs(2056, inplace=True)\n",
    "lakes = lakes[lakes[\"geometry\"].area >= 100000]\n",
    "lakes[\"geometry\"] = lakes[\"geometry\"].simplify(0.5)\n",
    "lakes.to_crs(4326, inplace=True)\n",
    "# buffer rivers (from line to polygon) and simplify borders with Douglas Peuker\n",
    "rivers.rename(columns={\"OBJECTID_G\": \"id\", \"GROSSERFLU\": \"name\"}, inplace=True)\n",
    "rivers = rivers[[\"id\", \"name\", \"geometry\"]].copy()\n",
    "rivers.reset_index(inplace=True)\n",
    "rivers[\"geometry\"] = rivers[\"geometry\"].buffer(40)\n",
    "rivers[\"geometry\"] = rivers[\"geometry\"].simplify(0.5)\n",
    "rivers.to_crs(4326, inplace=True)\n",
    "\n",
    "# merge lakes and rivers\n",
    "water_bodies = pd.concat([lakes, rivers])\n",
    "water_bodies.reset_index(inplace=True)\n",
    "_, water_bodies_LV95_sindex = _create_index(water_bodies, \"polygon\")\n",
    "\n",
    "# extract public green spaces\n",
    "green_space_types = [\"park\", \"cemetery\", \"recreation_ground\"]\n",
    "green_spaces = OSM_landuse.loc[OSM_landuse[\"fclass\"].isin(green_space_types)].copy()\n",
    "_, green_spaces_LV95_sindex = _create_index(green_spaces, \"polygon\")\n",
    "\n",
    "# extract residental areas\n",
    "residental_areas = OSM_landuse[(OSM_landuse[\"fclass\"] == \"residential\")].copy()\n",
    "_, residental_areas_LV95_sindex = _create_index(residental_areas, \"polygon\")\n",
    "\n",
    "# extract forest areas\n",
    "forest_areas = OSM_landuse[(OSM_landuse[\"fclass\"] == \"forest\")].copy()\n",
    "_, forest_areas_LV95_sindex = _create_index(forest_areas, \"polygon\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "2e8774cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "green_spaces.to_file(os.path.join(\"data\", \"green_spaces_visualize.shp\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee415204",
   "metadata": {},
   "source": [
    "### Extract Start/End Point Features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "fb9a0767",
   "metadata": {},
   "outputs": [],
   "source": [
    "# add start and end coordinates to triplegs\n",
    "tl['start_geom'] = [Point(segment.coords[0]) for segment in tl['geom_LV95']]\n",
    "tl['end_geom'] = [Point(segment.coords[-1]) for segment in tl['geom_LV95']]\n",
    "tl['start_point'] = [np.array(segment.coords)[0] for segment in tl['geom_LV95']]\n",
    "tl['end_point'] = [np.array(segment.coords)[-1] for segment in tl['geom_LV95']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1044ba46",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _compute_startEnd_feature(tl, feature_kdtree):\n",
    "    \"\"\"Calculate min/max Euclidean distance from start/end point to closest geospatial object\"\"\"\n",
    "    tl[\"feature_start\"] = feature_kdtree.query(np.array(tl[\"start_point\"].tolist()), k=[1])[0]\n",
    "    tl[\"feature_end\"] = feature_kdtree.query(np.array(tl[\"end_point\"].tolist()), k=[1])[0]\n",
    "    tl[\"feature_min\"] = tl[[\"feature_start\", \"feature_end\"]].values.min(1)\n",
    "    tl[\"feature_max\"] = tl[[\"feature_start\", \"feature_end\"]].values.max(1)\n",
    "    return tl[[\"feature_min\", \"feature_max\"]]\n",
    "\n",
    "tl[[\"rwStation_se_dist_min\", \"rwStation_se_dist_max\"]] = _compute_startEnd_feature(tl, rw_stations_LV95_kdtree)\n",
    "tl[[\"tramStop_se_dist_min\", \"tramStop_se_dist_max\"]] = _compute_startEnd_feature(tl, tram_stops_LV95_kdtree)\n",
    "tl[[\"busStop_se_dist_min\", \"busStop_se_dist_max\"]] = _compute_startEnd_feature(tl, bus_stops_LV95_kdtree)\n",
    "tl[[\"carParking_se_dist_min\", \"carParking_se_dist_max\"]] = _compute_startEnd_feature(tl, car_parkings_LV95_kdtree)\n",
    "tl[[\"bicycleParking_se_dist_min\", \"bicycleParking_se_dist_max\"]] = _compute_startEnd_feature(\n",
    "    tl, bicycle_parkings_LV95_kdtree\n",
    ")\n",
    "tl[[\"landingStage_se_dist_min\", \"landingStage_se_dist_max\"]] = _compute_startEnd_feature(\n",
    "    tl, landing_stages_LV95_kdtree\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33ea2c52",
   "metadata": {},
   "source": [
    "### Extract Point Object Features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "8e6a6876",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _compute_trackPoints_point_feature(tl_geom, feature_kdtree):\n",
    "    \"\"\"Calculate average Euclidean distance of track points to closest geospatial point object\"\"\"\n",
    "\n",
    "    # extract track point of trajectory\n",
    "    points_list = [np.array([x, y]) for x, y in tl_geom.coords]\n",
    "\n",
    "    # calculate distance to closest object for every track point\n",
    "    dists = feature_kdtree.query(np.array(points_list), k=[1])[0]\n",
    "\n",
    "    # calculate average distance\n",
    "    return dists.mean()\n",
    "\n",
    "tl[\"rwStations_dist_av\"] = tl.geom_LV95.apply(\n",
    "    lambda x: _compute_trackPoints_point_feature(x, rw_stations_LV95_kdtree)\n",
    ")\n",
    "tl[\"tramStops_dist_av\"] = tl.geom_LV95.apply(\n",
    "    lambda x: _compute_trackPoints_point_feature(x, tram_stops_LV95_kdtree)\n",
    ")\n",
    "tl[\"busStops_dist_av\"] = tl.geom_LV95.apply(\n",
    "    lambda x: _compute_trackPoints_point_feature(x, bus_stops_LV95_kdtree)\n",
    ")\n",
    "tl[\"poi_dist_av\"] = tl.geom_LV95.apply(lambda x: _compute_trackPoints_point_feature(x, poi_LV95_kdtree))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33b60e8f",
   "metadata": {},
   "source": [
    "### Extract Network Features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "98b37952",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _compute_trackPoints_line_feature(tl_geom, feature_sindex):\n",
    "    \"\"\"Calculate average Euclidean distance of track points to closest point on network\"\"\"\n",
    "\n",
    "    # extract track point of trajectory\n",
    "    points_list = [Point(x, y) for x, y in tl_geom.coords]\n",
    "    points_gs = gpd.GeoSeries(points_list)\n",
    "\n",
    "    # find closest object for every track point and return the distances \n",
    "    _, dists = feature_sindex.nearest(points_gs, return_all=False, return_distance=True)\n",
    "\n",
    "    # calculate average distance\n",
    "    return statistics.mean(dists)\n",
    "\n",
    "def _apply_extract_line(tl, **kwargs):\n",
    "    tqdm.pandas(desc=\"Line\")\n",
    "    return tl.geom_LV95.progress_apply(lambda x: _compute_trackPoints_line_feature(x, **kwargs))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "1ed0efe1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finish rail network feature\n"
     ]
    }
   ],
   "source": [
    "tl[\"rwNetwork_dist_av\"] = _apply_parallel(\n",
    "    tl, _apply_extract_line, n=-2, feature_sindex=rw_network_LV95_sindex\n",
    ")\n",
    "print(\"Finish rail network feature\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "64d22390",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finish tram network feature\n"
     ]
    }
   ],
   "source": [
    "tl[\"tramNetwork_dist_av\"] = _apply_parallel(\n",
    "    tl, _apply_extract_line, n=-2, feature_sindex=tram_network_LV95_sindex\n",
    ")\n",
    "print(\"Finish tram network feature\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "81b76091",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finish road network feature\n"
     ]
    }
   ],
   "source": [
    "tl[\"roadNetwork_dist_av\"] = _apply_parallel(\n",
    "    tl, _apply_extract_line, n=-2, feature_sindex=road_network_LV95_sindex\n",
    ")\n",
    "print(\"Finish road network feature\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "905b86ce",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\ProgramData\\Anaconda3\\envs\\mode\\lib\\site-packages\\joblib\\externals\\loky\\process_executor.py:700: UserWarning: A worker stopped while some jobs were given to the executor. This can be caused by a too short worker timeout or by a memory leak.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finish pb network feature\n"
     ]
    }
   ],
   "source": [
    "tl[\"pbNetwork_dist_av\"] = _apply_parallel(\n",
    "    tl, _apply_extract_line, n=-2, feature_sindex=pb_network_LV95_sindex\n",
    ")\n",
    "print(\"Finish pb network feature\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c7cf26e",
   "metadata": {},
   "source": [
    "### Extract Land Cover Features"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d5810b7-f917-4783-8185-9f06b3381ba9",
   "metadata": {},
   "source": [
    "#### Overlap"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "08539882",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _compute_trackPoints_polygon_feature(tl_geom, feature_sindex):\n",
    "    \"\"\"Calculate proportion of overlap of track points to land cover objects\"\"\"\n",
    "\n",
    "    # extract track point of trajectory\n",
    "    points_list = [Point(x, y) for x, y in tl_geom.coords]\n",
    "    points_gs = gpd.GeoSeries(points_list)\n",
    "\n",
    "    # check if track point lies within object\n",
    "    indices_within = feature_sindex.query_bulk(points_gs, predicate=\"within\")\n",
    "\n",
    "    # calculate proportion of track points lieing within object\n",
    "    n_total = len(points_gs)\n",
    "    n_within = len(np.unique(indices_within[0]))\n",
    "    return n_within / n_total\n",
    "\n",
    "def _apply_extract_polygon_overlap(tl, **kwargs):\n",
    "    tqdm.pandas(desc=\"Polygon\")\n",
    "    return tl.geom_LV95.progress_apply(lambda x: _compute_trackPoints_polygon_feature(x, **kwargs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "bad2e077",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finish water body feature\n"
     ]
    }
   ],
   "source": [
    "tl[\"water_ovlp_prop\"] = _apply_parallel(tl, _apply_extract_polygon_overlap, n=-1, feature_sindex=water_bodies_LV95_sindex)\n",
    "print(\"Finish water body feature\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "a8bf4b8e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finish green space feature\n"
     ]
    }
   ],
   "source": [
    "tl[\"greenSpaces_ovlp_prop\"] = _apply_parallel(tl, _apply_extract_polygon_overlap, n=-2, feature_sindex=green_spaces_LV95_sindex)\n",
    "print(\"Finish green space feature\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "bb12d048",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finish residential space feature\n"
     ]
    }
   ],
   "source": [
    "tl[\"residental_ovlp_prop\"] = _apply_parallel(\n",
    "    tl, _apply_extract_polygon_overlap, n=-2, feature_sindex=residental_areas_LV95_sindex\n",
    ")\n",
    "print(\"Finish residential space feature\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "6ce468a1",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\ProgramData\\Anaconda3\\envs\\mode\\lib\\site-packages\\joblib\\externals\\loky\\process_executor.py:700: UserWarning: A worker stopped while some jobs were given to the executor. This can be caused by a too short worker timeout or by a memory leak.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finish forest overlap feature\n"
     ]
    }
   ],
   "source": [
    "tl[\"forest_ovlp_prop\"] = _apply_parallel(tl, _apply_extract_polygon_overlap, n=-2, feature_sindex=forest_areas_LV95_sindex)\n",
    "print(\"Finish forest overlap feature\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c6e1846-8c5a-4a28-baea-c3bc57b2fcd9",
   "metadata": {},
   "source": [
    "#### Distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "7a667fe1-7c29-4e19-86b5-c5ba30309b65",
   "metadata": {},
   "outputs": [],
   "source": [
    "def _compute_trackPoints_polygon_distance_feature(tl_geom, feature_sindex):\n",
    "    \"\"\"Calculate average Euclidean distance of track points to closest point on the polygon land-cover\"\"\"\n",
    "\n",
    "    # extract every n-th track point of trajectory (simplification with regular time intervals)\n",
    "    points_list = [Point(x, y) for x, y in tl_geom.coords]\n",
    "    points_gs = gpd.GeoSeries(points_list)\n",
    "\n",
    "    # find closest object for every track point and return the distances \n",
    "    _, distances = feature_sindex.nearest(points_gs, return_all=False, return_distance=True)\n",
    "\n",
    "    # calculate average distance\n",
    "    return statistics.mean(distances)\n",
    "\n",
    "def _apply_extract_polygon_dist(tl, **kwargs):\n",
    "    tqdm.pandas(desc=\"Polygon distance\")\n",
    "    return tl.geom_LV95.progress_apply(lambda x: _compute_trackPoints_polygon_distance_feature(x, **kwargs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "9dd4ec23-3bef-450a-86b3-66f27378960c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finish water body feature\n"
     ]
    }
   ],
   "source": [
    "tl[\"water_poly_dist\"] = _apply_parallel(tl, _apply_extract_polygon_dist, n=-1, feature_sindex=water_bodies_LV95_sindex)\n",
    "print(\"Finish water body feature\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "9d12737b-50ad-464d-91ef-fc030218af8d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finish green space feature\n"
     ]
    }
   ],
   "source": [
    "tl[\"greenSpaces_poly_dist\"] = _apply_parallel(tl, _apply_extract_polygon_dist, n=-2, feature_sindex=green_spaces_LV95_sindex)\n",
    "print(\"Finish green space feature\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "59fef908-30f8-4091-ad05-4953aa43c72c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finish residential space feature\n"
     ]
    }
   ],
   "source": [
    "tl[\"residental_poly_dist\"] = _apply_parallel(\n",
    "    tl, _apply_extract_polygon_dist, n=-2, feature_sindex=residental_areas_LV95_sindex\n",
    ")\n",
    "print(\"Finish residential space feature\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "aed668c5-9c10-4021-8d49-a715cbcbc4b2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Finish forest overlap feature\n"
     ]
    }
   ],
   "source": [
    "tl[\"forest_poly_dist\"] = _apply_parallel(tl, _apply_extract_polygon_dist, n=-2, feature_sindex=forest_areas_LV95_sindex)\n",
    "print(\"Finish forest overlap feature\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c896246",
   "metadata": {},
   "source": [
    "## Save Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "2c26250a",
   "metadata": {},
   "outputs": [],
   "source": [
    "tl = tl[[\n",
    "    'user_id',  'mode', 'started_at', 'finished_at', 'user_t_quality', 'recording_interval',\n",
    "    'length', 'dur', 'speed_av', 'speed_85th', 'speed_90th', 'speed_95th', 'acc_av', 'acc_85th', 'br_av', 'br_85th',\n",
    "    'rwStation_se_dist_min', 'rwStation_se_dist_max', 'tramStop_se_dist_min', 'tramStop_se_dist_max',\n",
    "    'busStop_se_dist_min', 'busStop_se_dist_max', 'carParking_se_dist_min', 'carParking_se_dist_max',\n",
    "    'bicycleParking_se_dist_min', 'bicycleParking_se_dist_max', 'landingStage_se_dist_min', 'landingStage_se_dist_max',\n",
    "    'rwStations_dist_av', 'tramStops_dist_av', 'busStops_dist_av', 'poi_dist_av',\n",
    "    'rwNetwork_dist_av', 'tramNetwork_dist_av', 'roadNetwork_dist_av', 'pbNetwork_dist_av',\n",
    "    'water_ovlp_prop', 'greenSpaces_ovlp_prop', 'residental_ovlp_prop', 'forest_ovlp_prop',\n",
    "    'water_poly_dist', 'greenSpaces_poly_dist', 'residental_poly_dist', 'forest_poly_dist'\n",
    "]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "a0921018",
   "metadata": {},
   "outputs": [],
   "source": [
    "tl.to_csv(os.path.join(\"data\", \"triplegs_features.csv\"), index_label=\"index\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8b78065",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.8"
  },
  "toc-autonumbering": true,
  "toc-showcode": false,
  "vscode": {
   "interpreter": {
    "hash": "e3307148b578e6aa85616c61b46cb8e0917f8e20a9b0d5c7552aa4399d77cb8c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}