RouteClustering.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8626fd9d-dd67-42fd-8a7d-d4d67047ebaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "from bs4 import BeautifulSoup\n",
    "import os\n",
    "import requests\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn import svm\n",
    "from sklearn.metrics import r2_score, median_absolute_error, mean_absolute_error\n",
    "from pandas import read_csv\n",
    "from keras.models import Sequential\n",
    "from keras.layers import Dense\n",
    "from keras.optimizers import SGD, Adam\n",
    "from keras.wrappers.scikit_learn import KerasRegressor\n",
    "from sklearn.model_selection import cross_val_score\n",
    "from sklearn.model_selection import KFold\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.pipeline import Pipeline\n",
    "from keras.layers import Dropout\n",
    "from tensorflow.keras.callbacks import ModelCheckpoint, EarlyStopping\n",
    "from sklearn.model_selection import RepeatedKFold\n",
    "from sklearn.metrics import classification_report, accuracy_score\n",
    "import numpy as np\n",
    "from sklearn.ensemble import GradientBoostingClassifier\n",
    "import matplotlib.pyplot as plt\n",
    "import plotly.express as px\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import seaborn as sns\n",
    "from datetime import datetime, timedelta\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from pathlib import Path\n",
    "import pandas as pd\n",
    "import glob\n",
    "from io import StringIO\n",
    "import csv\n",
    "from sys import argv"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07547e1e-dc4b-45b9-a55b-f0caa9c2bc18",
   "metadata": {},
   "source": [
    "Step 1: Clustering should be applied for each ship"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "321e3e8f-19a3-43a9-8ae4-51a42ed834eb",
   "metadata": {},
   "source": [
    "define a function that can cluster time THEN space for a ship filter by its id. Using extended DBSCAN, one can extract multiple routes for the same ship, and extract more specifiecd corelations between weather and speed of this ship."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5dbf634-9cf5-43f9-a9d3-e0a87d641956",
   "metadata": {},
   "outputs": [],
   "source": [
    "from geopy.distance import great_circle\n",
    "def ST_DBSCAN(df, spatial_threshold, temporal_threshold, min_neighbors):\n",
    "    cluster_label = 0\n",
    "    NOISE = -1\n",
    "    UNMARKED = 777777\n",
    "    stack = []\n",
    "\n",
    "    # initialize each point with unmarked\n",
    "    df['cluster'] = UNMARKED\n",
    "    \n",
    "    # for each point in database\n",
    "    for index, point in df.iterrows():\n",
    "        if df.loc[index]['cluster'] == UNMARKED:\n",
    "            neighborhood = retrieve_neighbors(index, df, spatial_threshold, temporal_threshold)\n",
    "            \n",
    "            if len(neighborhood) < min_neighbors:\n",
    "                df.at[index, 'cluster'] = NOISE\n",
    "\n",
    "            else: # found a core point\n",
    "                cluster_label = cluster_label + 1\n",
    "                df.at[index, 'cluster'] = cluster_label# assign a label to core point\n",
    "\n",
    "                for neig_index in neighborhood: # assign core's label to its neighborhood\n",
    "                    df.at[neig_index, 'cluster'] = cluster_label\n",
    "                    stack.append(neig_index) # append neighborhood to stack\n",
    "                \n",
    "                while len(stack) > 0: # find new neighbors from core point neighborhood\n",
    "                    current_point_index = stack.pop()\n",
    "                    new_neighborhood = retrieve_neighbors(current_point_index, df, spatial_threshold, temporal_threshold)\n",
    "                    \n",
    "                    if len(new_neighborhood) >= min_neighbors: # current_point is a new core\n",
    "                        for neig_index in new_neighborhood:\n",
    "                            neig_cluster = df.loc[neig_index]['cluster']\n",
    "                            if (neig_cluster != NOISE) & (neig_cluster == UNMARKED): \n",
    "                                # TODO: verify cluster average before add new point\n",
    "                                df.at[neig_index, 'cluster'] = cluster_label\n",
    "                                stack.append(neig_index)\n",
    "    return df\n",
    "\n",
    "\n",
    "def retrieve_neighbors(index_center, df, spatial_threshold, temporal_threshold):\n",
    "    neigborhood = []\n",
    "\n",
    "    center_point = df.loc[index_center]\n",
    "\n",
    "    # filter by time \n",
    "    min_time = center_point['BaseDateTime'] - timedelta(minutes = temporal_threshold)\n",
    "    max_time = center_point['BaseDateTime'] + timedelta(minutes = temporal_threshold)\n",
    "    df = df[(df['BaseDateTime'] >= min_time) & (df['BaseDateTime'] <= max_time)]\n",
    "\n",
    "    # filter by distance\n",
    "    for index, point in df.iterrows():\n",
    "        if index != index_center:\n",
    "            distance = great_circle((center_point['LAT'], center_point['LON']), (point['LAT'], point['LON'])).meters\n",
    "            if distance <= spatial_threshold:\n",
    "                neigborhood.append(index)\n",
    "\n",
    "    return neigborhood\n",
    "\n",
    "def route_cluster(df, dtw_threshold):\n",
    "    route_label = 0\n",
    "    NOISE = -1\n",
    "    UNMARKED = 777777\n",
    "    stack = []\n",
    "\n",
    "    # initialize each point with unmarked\n",
    "    df['route'] = UNMARKED\n",
    "    df = df.groupby('cluster').agg(lambda x: list(x))\n",
    "\n",
    "        \n",
    "    # for each trip cluster\n",
    "    for index, point in df.iterrows():\n",
    "        xxx.iloc[index][['COG','LAT','LON']].values\n",
    "        if df.loc[index]['route'] == UNMARKED:\n",
    "            route_label = route_label + 1\n",
    "            df.at[index, 'route'] = route_label # assign a label to core point\n",
    "\n",
    "            for neighbor_index, neighbor_point in df.iterrows():\n",
    "                dtw_distance = dtw_ndim.distance(series1, series2)\n",
    "    return df\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9eacd4c-5e75-4ba6-9909-7f5e011160dc",
   "metadata": {},
   "source": [
    "Read a raw table of AIS data merged with weather information. \"Raw\" in sense of no subsampling is done so that all ship AIS-points are perserved."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d1e2021-1e0c-4bdf-a75b-0ec79243d2fa",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('dna_raw19.csv').dropna()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "356dfcdb-1409-4f30-a0bb-7009d6591e6d",
   "metadata": {},
   "source": [
    "Extract unique IDs "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53c40272-6945-49f4-a46a-780f257cac3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.MMSI.value_counts()[20:30]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c25beb7-70fb-49cb-8e3b-c854120156a8",
   "metadata": {},
   "source": [
    "Select a ship with relatively enough points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf1c4e19-c8a4-44a9-8f9f-289d6977a5eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "ship = df[df.MMSI == 367133160]\n",
    "ship.BaseDateTime = pd.to_datetime(ship.BaseDateTime)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c99b9a2-71a4-4ca1-a66c-2002a33814ae",
   "metadata": {},
   "source": [
    "Apply clustering Algorithm --> Spatial Temporal DBSCAN by using threshold of 12 km and 4 hours and #min_neighbors is 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5ff529b-0e98-42dd-b75b-ef1cdb003ad5",
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_threshold = 12000 # meters\n",
    "temporal_threshold = 240 # minutes\n",
    "min_neighbors = 3\n",
    "df_clustering = ST_DBSCAN(ship, spatial_threshold, temporal_threshold, min_neighbors)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "609e43b1-ef87-4bde-b0e1-8b919f7006b6",
   "metadata": {},
   "source": [
    "Trying to Group by cluster then show only the ones with the longest distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c54ef8d1-46e4-4beb-bd83-c2eab5b5b4cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.DataFrame(df_clustering.groupby('cluster').apply(lambda x : (great_circle((min(x['LAT']), min(x['LON'])), (max(x['LAT']), max(x['LON']))).meters, min(x['LAT']), max(x['LAT']),  min(x['LON']), max(x['LON']) )))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fbea715-efad-494b-b2ad-77bee113c5fe",
   "metadata": {},
   "source": [
    "Clean the clusters using the mean of ckusters as a threshold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ce6a592-ea11-4601-a8c1-d3a9bfa818ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "counts = df_clustering.cluster.value_counts() > df_clustering.cluster.value_counts().mean()\n",
    "colusters = counts[counts].index\n",
    "clean_clusters = df_clustering[df_clustering.cluster.isin(colusters)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c03ed5ee-a807-4113-aeb6-2823bd94cb46",
   "metadata": {},
   "source": [
    "Visualize either all clusters as lines or one cluster as scattered points with different colors which illustrates speed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b5a3ff8-3463-49a6-bb1a-47acdc3dc028",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "c = clean_clusters #[clean_clusters.cluster == 13]\n",
    "\n",
    "fig = px.line_geo(c, lat=c.LAT, lon=c.LON, hover_data=[c.BaseDateTime, c.COG, c.SOG], color=c.cluster)\n",
    "# fig = px.scatter_geo(c, lat=c.LAT, lon=c.LON, hover_data=[c.BaseDateTime, c.COG], color=c.SOG)\n",
    "\n",
    "fig.update_geos(projection_type=\"natural earth\")\n",
    "fig.update_layout( autosize=False, template='plotly_dark',\n",
    "        width=1200,\n",
    "        height=900,\n",
    "        hovermode=\"closest\"   \n",
    "        )\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cae7ad3-4281-4d2c-a481-583e69eaadd3",
   "metadata": {},
   "source": [
    "Step 2 Extracting similar routes using Dynamic-Time-Warping (DTW): "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80cb6aca-8ac5-4571-9c53-12fac7b7c478",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# all_c = c.groupby('cluster').agg(lambda x: list(x)).reset_index(drop=True)\n",
    "# dtw_i = dtw.distance_matrix(all_c['LAT'].values) + dtw.distance_matrix(all_c['LON'].values) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2dd707ec-e76e-4fb3-a8e0-adf24bea67ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import numpy as np\n",
    "# import seaborn as sns\n",
    "# import matplotlib.pylab as plt\n",
    "# plt.Figure(figsize=(20,2)) \n",
    "# ax = sns.heatmap(dtw_i, linewidth=0.5)\n",
    "# plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73429f4e-74d3-4317-8a46-c123759b49a2",
   "metadata": {},
   "source": [
    "Step 3: apply time serise analysis and predict the speed of the ship using the last couple AIS data points as input\n",
    "This can be achived by e.g.  LSTM  Model"
   ]
  }
 ],
 "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}