GIS3 exercise2_solution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# K-Means Clustering: A simple introduction\n",
    "Clustering describes a set of problems where the goal is to automatically identify groups from a set of observations (data). These groups are called clusters and can be described by the following general definition: \"Observation from the same cluster are more similar than observations from different clusters.\" \n",
    "\n",
    "A more formal definition of clustering is given by the K-Means problem. Given $n$ datapoints $x_i$, that have to be devided into $k$ clusters, each described by a cluster center $\\mu_j$ which is the mean of all datapoints that belong to the cluster $j$. Find the cluster centers that minimize the following loss function: \n",
    "\n",
    "$$ J = \\sum_{j=1}^{k} \\sum_{x_i \\in S_{j}} {\\| x_i - \\mu_j \\|^2}$$\n",
    "where $S_j$ is the set of all datapoints that are mapped to the cluster $j$ with cluster center $\\mu_j$.\n",
    "\n",
    "Finding an optimal solution for the k-means problem is np-hard, therefore there are several algorithms that find _good_ approximate (locally optimal) solutions. The most well-known algorithm is the _Loyds algorithm_ or _the k-means algorithm_.\n",
    "\n",
    "After the initialization (e.g., randomly choose cluster center among the data points), this algorithm consists of only two steps: \n",
    "- __Assigment:__\n",
    "Assign every point to the cluster center (centroid) that is closest.\n",
    "- __Update:__\n",
    "Calculate new centroids as the average position of all points assigned to a cluster.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate example data\n",
    "We will now try our own implementation of the k-means algorithm.\n",
    "At first we generate some artificial training data. The data will consist of two features (e.g., x-y coordinates). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "\n",
    "from sklearn.datasets import make_blobs\n",
    "\n",
    "n_samples = 1500\n",
    "random_state = 170\n",
    "X, y = make_blobs(n_samples=n_samples, n_features=2, random_state=random_state, cluster_std=1.2)\n",
    "\n",
    "print(X[0:10,:])\n",
    "\n",
    "plt.figure(figsize=(8, 8))\n",
    "plt.scatter(X[:, 0], X[:, 1], s=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: Initialize\n",
    "The first step of the k-means algorithm is to choose an initial set of clustering centers. Here we implement this in a very simple fashion, we simply choose a random datapoint and make sure that we don't choose the same datapoint twice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def init_centroids(k, X):\n",
    "    \"\"\"Choose random datapoints as centers for k-means\n",
    "    \n",
    "    k = number of clusters\n",
    "    X = data [numpy array]\n",
    "    \n",
    "    return: \n",
    "    centroids [tuple]: A tuple of arrays with the coordinates of the \n",
    "                       new cluster centers\"\"\"\n",
    "    \n",
    "    \n",
    "    centroids = []\n",
    "    for i in range(k):\n",
    "        rand = np.random.randint(0,len(X)-1)\n",
    "        randVal = tuple(X[rand,:])\n",
    "        while randVal in centroids:\n",
    "            # the while loop makes sure that you don't choose the same points twice\n",
    "            rand = np.random.randint(0,len(X)-1)\n",
    "            randVal = tuple(X[rand,:])\n",
    "        else:\n",
    "            centroids.append(randVal)\n",
    "    return centroids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(init_centroids(4,X))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: Assign points to cluster\n",
    "The assign step _assings_ every data point to the closest cluster center. In the following we create a function that does this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calcDist(a,b):\n",
    "    \"\"\"Calculate euclidean distance between datapoint and cluster center\"\"\"\n",
    "    return math.sqrt(sum((np.array(a)-np.array(b))**2))\n",
    "\n",
    "def assignToClusters(k, X, centroids):\n",
    "    \"\"\"Assign datapoints to closest cluster\n",
    "    k [] = number of clusters\n",
    "    X = data\n",
    "    centroids [list of tuples] = cluster centers\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    clusters: dictionary\n",
    "        has the coordinates of the centroids as key and all assigned points (via index) as values.\n",
    "    mapping : numpy array\n",
    "        mapping of size nx1, assigns each observation to 1 cluster (enumerated from 0:k).\n",
    "    \"\"\"\n",
    "    clusters = {}\n",
    "    mapping = []\n",
    "    for tup in centroids:\n",
    "        clusters[tup] = []\n",
    "    for i in range(len(X)):\n",
    "        pointDists = {}\n",
    "        for tup in centroids:\n",
    "            dist = calcDist(tuple(X[i]),tup)\n",
    "            pointDists[dist] = tup\n",
    "        ncp = pointDists.get(min(pointDists)) \n",
    "        clusters[ncp].append(i) #or i\n",
    "        mapping.append(centroids.index(ncp))\n",
    "        \n",
    "    return clusters, mapping\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Update centroids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def calcNewCentroids(clusters):\n",
    "    \"\"\"Calculate the new centroids as the mean of all points that are assigned to a cluster.\"\"\"\n",
    "    newcentroids = []\n",
    "    for old_centroid, membership in clusters.items():       \n",
    "        newcentroids.append(tuple(np.mean(X[membership,:], axis=0)))\n",
    "    return newcentroids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import time\n",
    "import pylab as pl\n",
    "from IPython import display\n",
    "\n",
    "def plot_kmeans(title=\"\"):\n",
    "    plt.cla()\n",
    "    plt.scatter(X[:, 0], X[:, 1], s=50, c=mapping)\n",
    "    plt.scatter(*zip(*centroids), s=500, c=range(k), marker='^',\n",
    "            edgecolors='k', linewidth=5)\n",
    "    \n",
    "    plt.title(title)\n",
    "    \n",
    "    display.clear_output(wait=True)\n",
    "    display.display(pl.gcf())\n",
    "    time.sleep(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Testing K-Means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 3 # Number of clusters\n",
    "steps = 2 # Number of steps to run k-means\n",
    "plt.figure(figsize=(8, 8))\n",
    "\n",
    "# initialize\n",
    "centroids = init_centroids(k,X)\n",
    "mapping = np.zeros((len(X))) # the initialization of the mapping is only for colors of the plot\n",
    "\n",
    "plot_kmeans(title=\"Init\")\n",
    "\n",
    "for i in range(steps):\n",
    "    \n",
    "    clusters, mapping = assignToClusters(k, X, centroids)\n",
    "    \n",
    "    plot_kmeans(title=\"Assign\")\n",
    "    centroids = calcNewCentroids(clusters)\n",
    "    plot_kmeans(title=\"Update\")\n",
    "plt.close()\n",
    "    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Clustering Trajectories\n",
    "An important application clustering algorithms in the mobility domain is the clustering of movement trajectories. A movement trajectory is a set of recorded locations (points) that a person _visited_ or passed by and that can be used to describe the movement of a person. A very common type of movement trajectories are GPS-trajectories. Here the movement of a person is continuously recorded by a GPS device usually with a fixed temporal sampling temporal sampling rate of several seconds. \n",
    "\n",
    "When the movement of a person is recorded using a GPS device we call the raw gps trackpoints __positionfixes__, places where a person was stationary for a longer time are called __stay points__ and all positionfixes that are recorded between two stay points are aggregated as a __trip leg__ (a GPS movement trajectory)\n",
    "\n",
    "In the following we look at examplary data of a publicly available GPS tracking set: [The geolife dataset](https://www.microsoft.com/en-us/download/details.aspx?id=52367&from=https%3A%2F%2Fresearch.microsoft.com%2Fen-us%2Fdownloads%2Fb16d359d-d164-469e-9fd4-daa38f2b2e13%2F)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Load modules and some helpful constants/functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import geopandas as gpd\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import os\n",
    "nice_extent = (116.2, 40, 116.3, 40.1)\n",
    "xmin, ymin, xmax, ymax = nice_extent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def write_shapefile(gdf, name):\n",
    "    if not os.path.isdir('./shapefiles'):\n",
    "        os.mkdir('./shapefiles')\n",
    "    \n",
    "\n",
    "    gdf.to_file(driver='ESRI Shapefile', filename= os.path.join('shapefiles', name))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "crs_wgs = {'init' :'epsg:4326'}\n",
    "crs_2D = {'init': 'epsg:32650'}\n",
    "spts_all = gpd.read_file(os.path.join(\"geolife\", \"shp\", \"geolife_staypoints.shp\"))\n",
    "tpls_all = gpd.read_file(os.path.join(\"geolife\", \"shp\", \"geolife_triplegs.shp\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### subsample data for speed\n",
    "Some calculations can be slow, use the subsampled version to try things out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmin = 1242856800 #+10*24*60*60\n",
    "tmax = 1242943199 #+10*24*60*60\n",
    "\n",
    "tpls = tpls_all[(tpls_all['started_at'] > tmin) & (tpls_all['started_at'] < tmax)].copy()\n",
    "spts = spts_all[(spts_all['started_at'] > tmin) & (spts_all['started_at'] < tmax)].copy()\n",
    "\n",
    "print(\"There are {:.0f} triplegs and {:.0f} staypoints\".format(tpls.shape[0], spts.shape[0]))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Clustering trajectories using descriptive features\n",
    "In their usual form, most clustering algorithms are not compatible with trajectories as an input. One way to cluster trajectoires is to use descriptive (scalar) features as an input instead of the trajectory itself.\n",
    "\n",
    "To be able to use the k-means clustering algorithm we create three features that describe each movement trajectory. Duration, distance and speed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate duration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tpls_all[\"duration\"] = tpls_all[\"finished_a\"] - tpls_all[\"started_at\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate distance\n",
    "To calculate a meaningful distance, be aware that you have to use a projected coordinate system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tpls_all2D = tpls_all.to_crs(crs_2D)\n",
    "tpls_all[\"length\"] = tpls_all2D[\"geometry\"].length\n",
    "tpls_all[\"speed\"] = tpls_all[\"length\"]/tpls_all[\"duration\"] *3.6\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspect the data visually\n",
    "If you plot the data using the script below, you will notice that your data is distributed in a strange way with many outliers. Think of a good data transformation that helps to have the data distributed in a more evenly manner! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.scatter(tpls_all.length, tpls_all.duration, marker=\"o\", alpha=0.1)\n",
    "ax.set_xlabel('length in s')\n",
    "ax.set_ylabel('duration in m')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot transformation\n",
    "Plot your data after your transformation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.scatter(np.log10(tpls_all.length), np.log10(tpls_all.duration), marker=\"o\", alpha=0.1)\n",
    "ax.set_xlabel('log length in s')\n",
    "ax.set_ylabel('log duration in m')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## clustering\n",
    "Try clustering the data using the KMeans algorithm. Create two different clustering using two sets of features:\n",
    "- Clustering 1: Use distance and duration as features\n",
    "- Clustering 2: Use only speed as features\n",
    "\n",
    "What could be a meaningful number of clusters if you think of mobility?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.cluster import KMeans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "km = KMeans(n_clusters=8)\n",
    "km.fit(np.log10(tpls_all[['length', 'duration']]))\n",
    "tpls_all['clst_durdist'] = km.labels_\n",
    "\n",
    "km.fit(np.log10(tpls_all[\"speed\"].values.reshape(-1, 1)))\n",
    "tpls_all['clst_speed'] = km.labels_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(1,2, figsize=(16,8))\n",
    "ax[0].scatter(x=np.log10(tpls_all['length']), y=np.log10(tpls_all['duration']), c=tpls_all['clst_durdist'])\n",
    "ax[1].scatter(x=np.log10(tpls_all['length']), y=np.log10(tpls_all['duration']), c=tpls_all['clst_speed'])\n",
    "for a in ax:\n",
    "    a.set_xlabel('log length in s')\n",
    "    a.set_ylabel('log duration in m')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Understanding the results\n",
    "In unsupervised learning it is often hard to estimate if what you did makes any sense. Therefore it is important to carefully look at the results carefully. \n",
    "- Print the cluster center of your clusterings and think about their meaning for your cluster. Does the clustering make sense?\n",
    "- Use the `write_shapefile` function that is defined about to save your geodataframe as a shapefile and examine the clustering result with your favorite GIS software"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "buff = list(km.cluster_centers_**10)\n",
    "buff.sort()\n",
    "print(buff)\n",
    "write_shapefile(tpls_all, 'trajectories')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cluster trajectories\n",
    "Another way to cluster trajectories that takes into account the geometry of the trajectory is based on similarity metrics such as the frechet distance. The common approach is to use the distance metric to calculate the pairwise distance of all trajectories and store it in a so-called distance matrix. Most of the common clustering implementations can cluster the trajectories directly based on this matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import calculate_distance_matrix, plot_nb_dists, norm_coords\n",
    "from sklearn.preprocessing import StandardScaler, MinMaxScaler\n",
    "import shapely"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## calcualte the distance matrix using the frechet distance. \n",
    "- Use the prepared function `utils\\calculate_distance_matrix` to calculate the distance matrix.\n",
    "- Visualize the distance matrix\n",
    "- Normalize your distance matrix\n",
    "- Cluster the trajectories based on the distance matrix\n",
    "- save the clustering into a shapefile and have a look with your favorite GIS software\n",
    "\n",
    "**Important:** Calculating the pairwise distance is very expensive in terms of computation. Best is to only do it with 10 Trajectories at a time. \n",
    "**Important II:** Use the data in a projected coordinate system. \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tpls2D = tpls_all2D.iloc[10:20]\n",
    "gdf = tpls2D.copy()\n",
    "D_frechet = StandardScaler().fit_transform(calculate_distance_matrix(gdf, distance='frechet'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#D_frechet = np.log(D_frechet+1)\n",
    "plt.figure()\n",
    "plt.imshow(D_frechet)\n",
    "plt.colorbar()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This function plots clusters of trajectories \n",
    "\n",
    "def plot_colored_lines(gdf, labels, ax=None):\n",
    "    c = labels\n",
    "    colors = plt.cm.jet(np.linspace(0,1,max(c)+1))\n",
    "    gdf.plot(color=colors[c], ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# kmeans D_frechet\n",
    "km = KMeans(3)\n",
    "km.fit(D_frechet)\n",
    "gdf['D_fre'] =  km.labels_\n",
    "try:\n",
    "    plot_colored_lines(gdf, km.labels_)\n",
    "except:\n",
    "    print(\"Plot not possbile on binder\")\n",
    "    \n",
    "\n",
    "write_shapefile(gdf, 'clustered_trajectories')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bonus: Use dtw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#D_dtw = calculate_distance_matrix(gdf, distance='dtw')\n",
    "\n",
    "#D_dtw_norm = calculate_distance_matrix(gdf_norm, distance='dtw')\n",
    "#D_dtw = StandardScaler().fit_transform(D_dtw)"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}