FleetGrouping-Blog-Support.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fleet Clustering\n",
    "\n",
    "### Tim Hochberg, 2019-01-16\n",
    "\n",
    "## Blog Support Edition\n",
    "\n",
    "We cluster vessel using HDBSCAN and a custom metric to derive fleets\n",
    "that are related in the sense that they spend a lot of time in the same\n",
    "location while at sea.\n",
    "\n",
    "## See Also\n",
    "\n",
    "* Other notebooks in https://github.com/GlobalFishingWatch/fleet-clustering for \n",
    "examples of clustering Squid Jiggers, etc.\n",
    "* This workspace that Nate put together: https://globalfishingwatch.org/map/workspace/udw-v2-85ff8c4f-fbfe-4126-b067-4d94cdd2b737"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "from __future__ import division\n",
    "from collections import Counter, OrderedDict\n",
    "import datetime as dt\n",
    "import hdbscan\n",
    "import logging\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as mpl_animation\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from skimage import color\n",
    "from IPython.display import HTML\n",
    "from fleet_clustering import bq\n",
    "from fleet_clustering import filters\n",
    "from fleet_clustering import distances\n",
    "from fleet_clustering import animation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load AIS Clustering Data\n",
    "\n",
    "Load the AIS data that we use for clustering. Note that it onlyu includes vessels away\n",
    "from shores so as to exclude clustering on ports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2016-01-01\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "//anaconda/envs/ml2/lib/python2.7/site-packages/google/auth/_default.py:66: UserWarning: Your application has authenticated using end user credentials from Google Cloud SDK. We recommend that most server applications use service accounts instead. If your application continues to use end user credentials from Cloud SDK, you might receive a \"quota exceeded\" or \"API not enabled\" error. For more information about service accounts, see https://cloud.google.com/docs/authentication/.\n",
      "  warnings.warn(_CLOUD_SDK_CREDENTIALS_WARNING)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2016-07-03\n",
      "2017-01-03\n",
      "2017-07-06\n",
      "2018-01-06\n",
      "2018-07-09\n"
     ]
    }
   ],
   "source": [
    "all_by_date = bq.load_ais_by_date('drifting_longlines', dt.date(2016, 1, 1), dt.date(2018, 12, 31),\n",
    "                                 fishing_only=False, min_km_from_shore=0)    \n",
    "pruned_by_date = {k : filters.remove_near_shore(10,\n",
    "                            filters.remove_chinese_coast(v)) for (k, v) in all_by_date.items()}\n",
    "valid_ssvid = sorted(filters.find_valid_ssvid(pruned_by_date))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple, Daily Clustering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create Distance Metrics\n",
    "\n",
    "Create an array of distance metrics. The details are still evolving, but in general\n",
    "we want to deal with two things.  Days on which a boat is missing and days where the\n",
    "boat is away from the fleet.\n",
    "\n",
    "* Distances to/from a boat on days when it is missing are represented by $\\infty$ in \n",
    "  the distance matrix. HDBSCAN ignores these values.\n",
    "* Only the closest N days are kept for each boat pair, allowing boats to leave the fleet\n",
    "  for up to half the year without penalty.\n",
    "  \n",
    "In addition, distances have a floor of 1 km to prevent overclustering when boats tie up\n",
    "up together, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dists_by_date = {}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "computing distance for 20170101 20170101\n",
      "computing distance for 20170107 20170107\n"
     ]
    }
   ],
   "source": [
    "for i in [1, 7]:\n",
    "    start_date = end_date = '201701{:02d}'.format(i)\n",
    "    key = start_date + '_daily'\n",
    "    if key in dists_by_date:\n",
    "        continue\n",
    "    print(\"computing distance for\", start_date, end_date)\n",
    "    subset_by_date = {k : v for (k, v) in pruned_by_date.items() if start_date <= k <= end_date}\n",
    "    C = distances.create_composite_lonlat_array(subset_by_date, valid_ssvid)\n",
    "    dists = distances.compute_distances_RMS(C)\n",
    "    dists_by_date[key] = dists"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit the Clusterer\n",
    "\n",
    "This is pretty straightforward -- all the complicated stuff is\n",
    "embedded in the matrix computations. Fleet size can be tweaked\n",
    "using `min_cluster_size` and `min_sample_size`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "//anaconda/envs/ml2/lib/python2.7/site-packages/hdbscan/hdbscan_.py:101: UserWarning: The minimum spanning tree contains edge weights with value infinity. Potentially, you are missing too many distances in the initial distance matrix for the given neighborhood size.\n",
      "  'size.', UserWarning)\n"
     ]
    }
   ],
   "source": [
    "raw_clusterers = {}\n",
    "for key in ['20170101_daily', '20170107_daily']:\n",
    "    clusterer = hdbscan.HDBSCAN(metric='precomputed', \n",
    "                                min_cluster_size=9,\n",
    "                               )\n",
    "    dists = dists_by_date[key]\n",
    "    clusterer.fit(dists)\n",
    "    raw_clusterers[key] = clusterer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up Fleets\n",
    "\n",
    "Set up the fleets for viewing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def to_rgb(string):\n",
    "    string = string.strip('#')\n",
    "    r = string[:2]\n",
    "    g = string[2:4]\n",
    "    b = string[4:]\n",
    "    return [int(x, 16) / 225.0 for x in (r, g, b)]\n",
    "\n",
    "\n",
    "def find_labels(dists):\n",
    "    clusterer = hdbscan.HDBSCAN(metric='precomputed', \n",
    "                                min_cluster_size=9).fit(dists)\n",
    "    \n",
    "    all_fleet_ssvid_set = set([s for (s, f) in zip(valid_ssvid, clusterer.labels_) if f >= 0])\n",
    "    valid_ssvid_set = set(valid_ssvid)\n",
    "\n",
    "    valid_ssvid_set = set(valid_ssvid)\n",
    "    joint_ssvid = valid_ssvid\n",
    "    labels = list(clusterer.labels_) \n",
    "\n",
    "    # Remove vessels that have no connection to other vessels\n",
    "    for i, ssvid in enumerate(valid_ssvid):\n",
    "        connections = (~np.isinf(dists[i])).sum()\n",
    "        if connections == 0:\n",
    "            labels[i] = -1\n",
    "            \n",
    "    return joint_ssvid, labels\n",
    "\n",
    "\n",
    "def create_fleet_mapping(labels, include_carriers=False):\n",
    "    counts = []\n",
    "    skip = []\n",
    "    for i in range(max(labels) + 1):\n",
    "        if i in skip:\n",
    "            counts.append(0)\n",
    "        else:\n",
    "            counts.append((np.array(labels) == i).sum())\n",
    "\n",
    "    fleet_ids = [x for x in np.argsort(counts)[::-1] if counts[x] > 0]\n",
    "    fleet_ids_without_carriers = [x for x in fleet_ids if x != max(labels)]\n",
    "\n",
    "    fleets = OrderedDict()\n",
    "    n_hues = int(np.ceil(len(fleet_ids) / 4.0))\n",
    "    used = set()\n",
    "    for i, fid in enumerate(fleet_ids_without_carriers):\n",
    "        b = (i // (2 * n_hues)) % 2\n",
    "        c = (i // 2)% n_hues\n",
    "        d = i  % 2\n",
    "        symbol = 'H^'[d]\n",
    "        assert (b, c, d) not in used, (i, b, c, d)\n",
    "        used.add((b, c, d))\n",
    "        sat = 1\n",
    "        val = 1\n",
    "        raw_hue = c / float(n_hues)\n",
    "        # We remap the raw hue in order to avoid the 60 degree segment around blue\n",
    "        hue = 5. / 6. * raw_hue\n",
    "        if hue > 7. / 12.:\n",
    "            hue += 1. / 6.\n",
    "        assert 0 <= hue < 1, hue\n",
    "        [[clr]] = color.hsv2rgb([[(hue, sat, val)]])\n",
    "        fg = [[0.1511111111111111, 0.2, 0.3333333333333333], clr][b]\n",
    "        bg = [clr, [0.1511111111111111, 0.2, 0.3333333333333333]][b]\n",
    "        w = [1, 2][b]\n",
    "        sz = [9, 7][b]\n",
    "        fleets[fid] = (symbol, tuple(fg), tuple(bg), sz, w,  str(i + 1))\n",
    "    if include_carriers:\n",
    "        fleets[max(labels)] = ('1', 'k', 'k', 8, 2, 'Carrier Vessel')\n",
    "    print(len(set([x for x in fleets if x != -1])), \"fleets\")\n",
    "    return fleets\n",
    "    \n",
    "    \n",
    "def iou(a, b):\n",
    "    a = set(a)\n",
    "    b = set(b)\n",
    "    return len(a & b) / len(a | b)\n",
    "\n",
    "def best_match(a, bs):\n",
    "    ious = [iou(a, b) for b in bs]\n",
    "    i = np.argmax(ious)\n",
    "    if ious[i] == 0:\n",
    "        return None\n",
    "    return i\n",
    "\n",
    "    \n",
    "def adapt_fleet_mapping(base_fleets, base_ssvid, base_labels, new_ssivd, new_labels):\n",
    "    new_labels = np.array(new_labels)\n",
    "    ssvid_base = []\n",
    "    base_fleet_ids = sorted(base_fleets)\n",
    "    for fid in base_fleet_ids:\n",
    "        mask = (base_labels == fid)\n",
    "        ssvid_base.append(np.array(base_ssvid)[mask])\n",
    "\n",
    "    ssvid_new = []\n",
    "    new_fleet_ids = sorted(set([x for x in new_labels if x != -1]))\n",
    "    for fid in new_fleet_ids:\n",
    "        mask = (new_labels == fid)\n",
    "        ssvid_new.append(np.array(new_ssivd)[mask])\n",
    "\n",
    "    rev_mapping = {}\n",
    "    for fid, ssvid_list in zip(new_fleet_ids, ssvid_new):\n",
    "        i = best_match(ssvid_list, ssvid_base)\n",
    "        if i is None:\n",
    "            rev_mapping[fid] = None\n",
    "        else:\n",
    "            rev_mapping[fid] = base_fleet_ids[i]\n",
    "            \n",
    "    mapping = {}\n",
    "    for k, v in rev_mapping.items():\n",
    "        if v in mapping:\n",
    "            mask = (new_labels == k)\n",
    "            new_labels[mask] = mapping[v]\n",
    "        else:\n",
    "            mapping[v] = k\n",
    "                         \n",
    "    new_fleets = OrderedDict()\n",
    "    for i, fid in enumerate(base_fleets):\n",
    "        if fid in mapping and mapping[fid] is not None:\n",
    "            k = mapping[fid]\n",
    "            if k in new_fleets:\n",
    "                print(\"Skipping\", k, fid, \"because of double match\")\n",
    "                new_fleets[i + max(base_fleets)] = base_fleets[fid]\n",
    "            else:\n",
    "                new_fleets[mapping[fid]] = base_fleets[fid]\n",
    "        else:\n",
    "            new_fleets[i + max(base_fleets)] = base_fleets[fid]\n",
    "            \n",
    "    return new_fleets, new_labels\n",
    "                         "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "26 fleets\n"
     ]
    }
   ],
   "source": [
    "ssvid_20170101, labels_20170101 = find_labels(dists_by_date['20170101_daily'])\n",
    "fleets_20170101 = create_fleet_mapping(labels_20170101)\n",
    "all_by_date_20170101 = {k : v for (k, v) in all_by_date.items() if '20170101' <= k <= '20170101'}\n",
    "\n",
    "animation.make_anim(ssvid_20170101, \n",
    "                           labels_20170101, \n",
    "                           all_by_date_20170101, \n",
    "                           interval=100,\n",
    "                           fleets=fleets_20170101, \n",
    "                           show_ungrouped=True,\n",
    "                           alpha=1,\n",
    "                           legend_cols=12,\n",
    "                           plot_frame=0,\n",
    "                           text_color='white',\n",
    "                           ungrouped_legend=\"Ungrouped\")\n",
    "\n",
    "plt.savefig('fig1.png', dpi=450, facecolor='#222D4B')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "24 fleets\n"
     ]
    }
   ],
   "source": [
    "ssvid_20170107, labels_20170107 = find_labels(dists_by_date['20170107_daily'])\n",
    "fleets_20170107 = create_fleet_mapping(labels_20170107)\n",
    "all_by_date_20170107 = {k : v for (k, v) in all_by_date.items() if '20170107' <= k <= '20170107'}\n",
    "\n",
    "animation.make_anim(ssvid_20170107, \n",
    "                           labels_20170107, \n",
    "                           all_by_date_20170107, \n",
    "                           interval=100,\n",
    "                           fleets=fleets_20170107, \n",
    "                           show_ungrouped=True,\n",
    "                           alpha=1,\n",
    "                           legend_cols=12,\n",
    "                           plot_frame=0,\n",
    "                           text_color='white',\n",
    "                           ungrouped_legend=\"Ungrouped\")\n",
    "\n",
    "plt.savefig('fig2.png', dpi=450, facecolor='#222D4B')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Clustering Across a Year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "computing distance for 20170101 20171231\n"
     ]
    }
   ],
   "source": [
    "start_date = '20170101'\n",
    "end_date = '20171231'\n",
    "key = start_date\n",
    "if key not in dists_by_date:\n",
    "    print(\"computing distance for\", start_date, end_date)\n",
    "    subset_by_date = {k : v for (k, v) in pruned_by_date.items() if start_date <= k <= end_date}\n",
    "    C = distances.create_composite_lonlat_array(subset_by_date, valid_ssvid)\n",
    "    dists = distances.compute_distances_RMS(C)\n",
    "    dists_by_date[key] = dists"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "raw_clusterers = {}\n",
    "for key in ['20170101']:\n",
    "    clusterer = hdbscan.HDBSCAN(metric='precomputed', \n",
    "                                min_cluster_size=9,\n",
    "                               )\n",
    "    dists = dists_by_date[key]\n",
    "    clusterer.fit(dists)\n",
    "    raw_clusterers[key] = clusterer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "28 fleets\n"
     ]
    }
   ],
   "source": [
    "ssvid_2017, labels_2017 = find_labels(dists_by_date['20170101'])\n",
    "fleets_2017 = create_fleet_mapping(labels_2017)\n",
    "\n",
    "animation.make_anim(ssvid_2017, \n",
    "                           labels_2017, \n",
    "                           all_by_date_20170101, \n",
    "                           interval=100,\n",
    "                           fleets=fleets_2017, \n",
    "                           show_ungrouped=True,\n",
    "                           alpha=1,\n",
    "                           legend_cols=12,\n",
    "                           plot_frame=0,\n",
    "                           text_color='white',\n",
    "                           ungrouped_legend=\"Ungrouped\")\n",
    "plt.savefig('fig3.png', dpi=450, facecolor='#222D4B')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "animation.make_anim(ssvid_2017, \n",
    "                           labels_2017, \n",
    "                           all_by_date_20170107, \n",
    "                           interval=100,\n",
    "                           fleets=fleets_2017, \n",
    "                           show_ungrouped=True,\n",
    "                           alpha=1,\n",
    "                           legend_cols=12,\n",
    "                           plot_frame=0,\n",
    "                           text_color='white',\n",
    "                           ungrouped_legend=\"Ungrouped\")\n",
    "plt.savefig('fig4.png', dpi=450, facecolor='#222D4B')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tweaked Distance Matrix\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "computing distance for 20170101 20171231\n"
     ]
    }
   ],
   "source": [
    "start_date = '20170101'\n",
    "end_date = '20171231'\n",
    "key = start_date + '_tweaked'\n",
    "if key not in dists_by_date:\n",
    "    print(\"computing distance for\", start_date, end_date)\n",
    "    subset_by_date = {k : v for (k, v) in pruned_by_date.items() if start_date <= k <= end_date}\n",
    "    C = distances.create_composite_lonlat_array(subset_by_date, valid_ssvid)\n",
    "    dists = distances.compute_distances_tweaked(C)\n",
    "    dists_by_date[key] = dists"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "raw_clusterers = {}\n",
    "for key in ['20170101_tweaked']:\n",
    "    clusterer = hdbscan.HDBSCAN(metric='precomputed', \n",
    "                                min_cluster_size=9,\n",
    "                               )\n",
    "    dists = dists_by_date[key]\n",
    "    clusterer.fit(dists)\n",
    "    raw_clusterers[key] = clusterer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "39 fleets\n"
     ]
    }
   ],
   "source": [
    "ssvid_2017_tweaked, labels_2017_tweaked = find_labels(dists_by_date['20170101_tweaked'])\n",
    "fleets_2017_tweaked = create_fleet_mapping(labels_2017_tweaked)\n",
    "\n",
    "animation.make_anim(ssvid_2017_tweaked, \n",
    "                           labels_2017_tweaked, \n",
    "                           all_by_date_20170101, \n",
    "                           interval=100,\n",
    "                           fleets=fleets_2017_tweaked, \n",
    "                           show_ungrouped=True,\n",
    "                           alpha=1,\n",
    "                           legend_cols=12,\n",
    "                           plot_frame=0,\n",
    "                           text_color='white',\n",
    "                           ungrouped_legend=\"Ungrouped\")\n",
    "plt.savefig('fig5.png', dpi=450, facecolor='#222D4B')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "39 fleets\n"
     ]
    }
   ],
   "source": [
    "ssvid_2017_tweaked, labels_2017_tweaked = find_labels(dists_by_date['20170101_tweaked'])\n",
    "fleets_2017_tweaked = create_fleet_mapping(labels_2017_tweaked)\n",
    "\n",
    "animation.make_anim(ssvid_2017_tweaked, \n",
    "                           labels_2017_tweaked, \n",
    "                           all_by_date_20170107, \n",
    "                           interval=100,\n",
    "                           fleets=fleets_2017_tweaked, \n",
    "                           show_ungrouped=True,\n",
    "                           alpha=1,\n",
    "                           legend_cols=12,\n",
    "                           plot_frame=0,\n",
    "                           text_color='white',\n",
    "                           ungrouped_legend=\"Ungrouped\")\n",
    "plt.savefig('fig6.png', dpi=450, facecolor='#222D4B')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}