Effect of SalishSeaCast NEMO water mask for GeoTIFFs on calculation of monthly VTE probability

VTEProbWaterMask.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Water Mask & VTE Probability\n",
    "\n",
    "Effect of SalishSeaCast NEMO water mask for GeoTIFFs on calculation of monthly VTE probability."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy\n",
    "import rasterio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_vte_probability_no_mask(geotiffs_dir):\n",
    "    total_vte_by_month = numpy.empty(12)\n",
    "\n",
    "    for month in range(1, 13):\n",
    "        # The filenames are formatted as \"all_2018_MM.tif\"\n",
    "        f_name = geotiffs_dir / f\"all_2018_{month:02d}.tif\"\n",
    "\n",
    "        with rasterio.open(f_name) as dataset:\n",
    "            total_vte_by_month[month - 1] = dataset.read(1, boundless=True, fill_value=0).sum()\n",
    "\n",
    "    # calculate VTE probability by month based on total traffic for each month\n",
    "    vte_probability = total_vte_by_month / total_vte_by_month.sum()\n",
    "\n",
    "    return vte_probability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_vte_probability_masked(geotiffs_dir, geotiff_watermask):\n",
    "    total_vte_by_month = numpy.empty(12)\n",
    "\n",
    "    for month in range(1, 13):\n",
    "        # The filenames are formatted as \"all_2018_MM.tif\"\n",
    "        f_name = geotiffs_dir / f\"all_2018_{month:02d}.tif\"\n",
    "\n",
    "        with rasterio.open(f_name) as dataset:\n",
    "            total_vte_by_month[month - 1] = dataset.read(1, boundless=True, fill_value=0).sum(where=geotiff_watermask)\n",
    "\n",
    "    # calculate VTE probability by month based on total traffic for each month\n",
    "    vte_probability = total_vte_by_month / total_vte_by_month.sum()\n",
    "\n",
    "    return vte_probability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "geotiffs_dir = Path(\"/media/doug/warehouse/MIDOSS/ShipTrackDensityGeoTIFFs/\")\n",
    "geotiff_watermask = numpy.load(Path(\"/media/doug/warehouse/MIDOSS/ShipTrackDensityGeoTIFFs/geotiff-watermask.npy\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "vte_probability_no_mask = calc_vte_probability_no_mask(geotiffs_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "vte_probability_masked = calc_vte_probability_masked(geotiffs_dir, geotiff_watermask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fb1350c4d60>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXRV5dn38e9lggaRScFHBQS0gDInRiDBWim1gqXghEIHH7SWUus8NW191SJ1waulSGFJFRV9yiMtFoUqrSO+rQYEFMosMihErQIVBAFluN4/zuZ4ctg5OSQ5Q5LfZ62snL3vPVwnhPPLvvfe9zZ3R0REJN5RmS5ARESykwJCRERCKSBERCSUAkJEREIpIEREJFRupguoSS1atPB27dplugwRkVrjrbfe2uruLcPa6lRAtGvXjsWLF2e6DBGRWsPM3q+oTV1MIiISSgEhIiKhUhoQZjbAzN4xs3VmVhLSfoaZzTezL8zstiNZV0REUitl5yDMLAeYDJwPlAGLzGyOu6+KWew/wA3ARVVYNyn79u2jrKyMvXv3VvGdSF2Vl5dH69atadCgQaZLEclKqTxJ3QtY5+4bAMxsBjAEiH7Iu/snwCdm9p0jXTdZZWVlNG7cmHbt2mFmVXsnUue4O9u2baOsrIz27dtnuhyRrJTKLqZWwOaY6bJgXo2ua2YjzWyxmS3esmXLYe179+7lhBNOUDhIOWbGCSecoCNLkQRSGRBhn8jJDh2b9Lru/rC7F7p7YcuWoZfyKhwklH4vRBJLZUCUAW1iplsDH6ZhXRERqQGpPAexCOhgZu2BD4BhwPfSsG5C7Uqer4nNRL03Nv70SXk333wzbdu25aabbgLgggsuoE2bNkydOhWAW2+9lVatWvG1r32NVatWUVJSwrPPPkvHjh3p3LkzAOeddx4PPPAAhYWF1a73vvvu45e//OVh8x988EE2btzIhAkTAPjJT37C+vXrefnllwH4/e9/z7vvvsvEiRNDtxtfc1W89tpr9OvXj6lTp/KjH/0IgCVLllBQUMD999/PbbfdVskWyps2bRqLFy9m0qRJVa5JpD5LWUC4+34zuw54AcgBHnP3lWY2KmifYmYnAYuBJsBBM7sJ6Ozun4Wtm6paU6m4uJiZM2dy0003cfDgQbZu3cpnn30WbS8tLWXChAn07t2bwYMHA5EP20GDBlXrw7YiFQVEcXEx06dPj04vXbqUgwcPcuDAAXJycigtLeWiiy46bL1DqlLz/v37yc0t/yvYrVs3/vSnP0UDYsaMGfTo0SPpbUp5yfxBVNkfOVJ/pfQ+CHef6+4d3f10d/9NMG+Ku08JXv/b3Vu7exN3bxa8/qyidWujvn37UlpaCsDKlSvp2rUrjRs35tNPP+WLL75g9erV5OfnM23aNK677jpKS0uZM2cOt99+Oz179mT9+vUAzJw5k169etGxY0f++c9/ApET8FdddRXdunUjPz+fefPmAUS3dcigQYN47bXXKCkpYc+ePfTs2ZPvf//75erMz89n7dq17Nmzhx07dnDsscfSs2dPli9fDkSCrLi4mEceeYSzzz6bHj16cOmll7J79+7QmtevX8+AAQM466yz+PrXv86aNWsAGDFiBLfccgv9+vXj5z//+WE/r1NPPZW9e/fy8ccf4+78/e9/Z+DAgdH2sP0f+vl07dqVHj16cO655x623eeff56ioiK2bt1atX9IkXqoTo3FlI1OOeUUcnNz2bRpE6WlpRQVFfHBBx8wf/58mjZtSvfu3Tn66KOjyxcXFzN48GAGDRrEZZddFp2/f/9+Fi5cyNy5c/n1r3/Nyy+/zOTJkwFYvnw5a9as4dvf/jZr166tsJaxY8cyadIkli5delhbbm4uPXv2ZNGiRezZs4fevXvToUMHSktLOfHEE3F32rRpwyWXXMKPf/xjAO68804effRRrr/++sNq7t+/P1OmTKFDhw68+eabXHvttbz66qsArF27lpdffpmcnJzQOi+77DJmzpxJfn4+BQUFHHPMMdG2ivY/evRoXnjhBVq1asX27dvLbe+ZZ55h/PjxzJ07l+bNm1f8jyUi5Sgg0uDQUURpaSm33HILH3zwAaWlpTRt2pTi4uKktnHJJZcAcNZZZ/Hee+8B8Prrr3P99dcDcMYZZ9C2bduEAZFsnXv27KGoqIgOHTpw33330bJly2idK1as4M4772T79u3s2rWLCy644LDt7Nq1i9LSUoYOHRqd98UXX0RfDx06tMJwALj88su54oorWLNmDcOHD48egSXaf9++fRkxYgSXX3559GcFMG/ePBYvXsyLL75IkyZNqvyzqdPuaZrEMjtSX4dkHY3FlAbFxcWUlpayfPlyunbtSp8+fZg/fz6lpaX07ds3qW0c+is6JyeH/fv3A5GbvcLk5uZy8ODB6HSy1/ofqnP+/PkUFRVx5plnsmrVqnJ1jhgxgkmTJrF8+XLuvvvu0G0fPHiQZs2asXTp0ujX6tWro+2NGjVKWMdJJ51EgwYNeOmll+jfv3+5tor2P2XKFMaMGcPmzZvp2bMn27ZtA+C0005j586d1QpOkfpKAZEGffv25bnnnuP4448nJyeH448/nu3bt0c/iOM1btyYnTt3Vrrdc889N3piee3atWzatIlOnTrRrl276EnmzZs3s3Dhwug6DRo0YN++faHbKy4uZsGCBWzZsoUTTzwRM6Nly5bMnj07egSxc+dOTj75ZPbt21fupHZszU2aNKF9+/bMnDkTiATZv/71ryR/WhGjR49m3Lhxhx1pVLT/9evX07t3b0aPHk2LFi3YvDlyn2Xbtm2ZNWsWV155JStX1srrHEQypt51MWXiio1u3bqxdetWvve975Wbt2vXLlq0aHHY8sOGDePHP/4xEydO5Omnn65wu9deey2jRo2iW7du5ObmMm3aNI455hj69u1L+/bt6datG127dqWgoCC6zsiRI+nevTsFBQXlPmABmjdvTsuWLenSpUt0XlFREW+88Ub0SqJ7772X3r1707ZtW7p16xYNhfiap0+fzk9/+lPGjBnDvn37GDZs2BFdjVRR11tF+7/99tt59913cXf69+9Pjx49oudaOnXqxPTp0xk6dCh//etfOf3005OuQ6Q+s4q6KWqjwsJCj39g0OrVqznzzDMzVJFku0z9fqTr8tOk9pOXxC1GOgdRZ5nZW+4eepOVuphERCSUAkJEREIpIEREJFS9O0ktkkhWDU2h+xMkwxQQUissK9te6TLdWzdLQyUSJquCVWqMuphERCRU/QuIe5rW7Fclbr755ugQ2hAZ7vuaa66JTt96662MHz+eOXPmMHbsWCAyMuqqVV89XfW8884j/vLdeBdffDHPPvtsdLpTp06MGTMmOn3ppZcya9asCtefMGFCdOC7qrrnnnswM9atWxed97vf/Q4zq7T+MCNGjEh4H4iIpFb9C4g0OzR8BRAd7jv2jt5Dw1gMHjyYkpIS4PCAONL9bNu2jeOOO4758+dH2+fPn59w3KeqBMSBAwcOm9etWzdmzJgRnX766adTMmy5iKSeAiLFUjncd0X7KS0tZdCgQWzZsgV3Z+PGjTRs2JCTTjqJn/70pxQWFtKlSxfuvvtuACZOnMiHH35Iv3796NevHwAvvvgiRUVFFBQUMHToUHbt2gVAu3btGD16NOecc050KI1YF110EbNnzwZgw4YNNG3alNhHwYbtH6CkpITOnTvTvXv30AcDTbr/N/yfm68tN8bUYT5cUvlXTaiBI0uR2kAnqVMslcN9xzrrrLNYsWIFX375JaWlpXzjG99gw4YNrF69miVLlkQH2/vNb37D8ccfz4EDB+jfvz/Lli3jhhtuYPz48cybN48WLVqwdetWxowZw8svv0yjRo0YN24c48eP56677gIgLy+P119/PfT9NmnShDZt2rBixQpmz57NFVdcweOPPx5tD9t/69ateeaZZ1izZg1mdthw3XfccQf/2baF0eMn6znStZmuyqp1dASRBrHDfRcVFVFUVBSdrs5w37GOOeYYunTpwttvv82CBQvo3bt36H7+/Oc/U1BQQH5+PitXrgztylqwYAGrVq2ib9++9OzZkyeeeIL3338/2n7FFVckrHXYsGHMmDGDZ599losvvrhcW9j+mzRpQl5eHtdccw2zZs3i2GOPjS5/7733sn37du4aO0HhIJJmOoJIg/jhvtu0acNvf/tbmjRpwtVXX53UNsKG+w7bzz/+8Q927txJ8+bN6dOnD5MmTWLJkiWMGjWKjRs38sADD7Bo0SKaN2/OiBEjQofrdnfOP/98nnrqqdD9VDZc93e/+11uv/12CgsLyz2DoaL95+bmsnDhQl555RVmzJjBpEmTog8XOvvss3nrrbfY8emnNNXDfuo9XU6bXjqCSINUDfcdtp8//OEP0VFTu3fvzoIFC9i0aRNdunThs88+o1GjRjRt2pSPP/6Yv/3tb6H77NOnD2+88Ub0aqTdu3cf0fMUGjZsyLhx4/jVr35Vbn5F+9+1axc7duzgwgsvZMKECeWeeDdgwABKSkq4bsQVfL7ryH8mIlJ19e8IIgN9nKka7jtecXExGzZs4Be/+AUQeXDQiSeeSJs2bTjqqKPo0aMH+fn5dOnShdNOO63cw4pGjhzJwIEDOfnkk5k3bx7Tpk1j+PDh0SfBjRkzho4dOyZdy7Bhww6bV9H+d+7cyZAhQ9i7dy/uzu9+97ty6w0dOpRV73/MDVd/j8lP/Jm8hg2TrkPqIZ3rqDEa7ltqhaTupD5qY+UbOiW/3GT870e6hsfWfrJ7P/WJhvsWEZEjpoAQEZFQ9SIg6lI3mtQc/V6IJFbnAyIvL49t27bpw0DKcXe2bdtGXl5epksRyVp1/iqm1q1bU1ZWxpYtWzJdilTDx5/uqXSZ1ZbEv/GO1dGXeXl5tG7dujplidRpdT4gGjRoQPv27TNdhlTTQF29IpJ2db6LSUREqkYBISIioRQQIiISSgEhIiKh6vxJahGRI6UhPSJ0BCEiIqFSGhBmNsDM3jGzdWZWEtJuZjYxaF9mZgUxbTeb2UozW2FmT5mZ7mgSEUmjlHUxmVkOMBk4HygDFpnZHHePfYTZQKBD8NUbeAjobWatgBuAzu6+x8z+DAwDpqWq3npJwyKLSAKpPILoBaxz9w3u/iUwAxgSt8wQ4EmPWAA0M7OTg7ZcoKGZ5QLHAh+msFYREYmTyoBoBWyOmS4L5lW6jLt/ADwAbAI+Ana4+4thOzGzkWa22MwWazgNEZGak8qACHvCfPyIeaHLmFlzIkcX7YFTgEZm9oOwnbj7w+5e6O6FLVu2rFbBIiLylVQGRBnQJma6NYd3E1W0zLeAje6+xd33AbOA4hTWKiIicVIZEIuADmbW3syOJnKSeU7cMnOAK4OrmfoQ6Ur6iEjXUh8zO9bMDOgPrEZERNImZVcxuft+M7sOeAHIAR5z95VmNiponwLMBS4E1gG7gauCtjfN7GngbWA/sAR4OFW1iojI4VJ6J7W7zyUSArHzpsS8duBnFax7N3B3Kuury5K7EzQNhYhIraWhNkREslkG71fSUBsiIhJKRxAiIhmS7V3BCghJPQ3pIVIrqYtJRERC6QhCqiXbD5FFpOp0BCEiIqEUECIiEkoBISIioRQQIiISSgEhIiKhFBAiIhJKASEiIqEUECIiEkoBISIioRQQIiISSgEhIiKhFBAiIhKq0oAwswfMrEs6ihERkeyRzBHEGuBhM3vTzEaZWRKD+4uISG1X6XDf7j4VmGpmnYCrgGVm9gbwiLvPS3WBWUUPvhGReiSpcxBmlgOcEXxtBf4F3GJmM1JYm4iIZFClRxBmNh74LvAqcJ+7LwyaxpnZO6ksTkREMieZJ8qtAO50990hbb1quB4BdWWJSFZIpovp+/HhYGavALi7PqVEROqoCo8gzCwPOBZoYWbNAQuamgCnpKE2ERHJoERdTD8BbiISBm/HzP8MmJzKojKhXcnzlS7zXl4aChERyRIVBoS7Pwg8aGbXu/vv01hTnaYgEpHaIlEX0zfd/VXgAzO7JL7d3WeltDIREcmoRF1M3yByaet3Q9ocUECIiNRhibqY7g6+X5W+ckREJFsk6mK6JdGK7j6+5ssREZFskaiLqXHaqhARkayTqIvp19XduJkNAB4EcoCp7j42rt2C9guB3cAId387aGsGTAW6EjnncbW7z69uTSIikpxEXUx3uPv/NbPfE/mALsfdb0i04WCAv8nA+UAZsMjM5rj7qpjFBgIdgq/ewEPBd4gEx9/d/TIzO5rITXsiIpImibqYVgffF1dx272Ade6+ASAY+XUIEBsQQ4An3d2BBWbWzMxOBj4HzgVGALj7l8CXVaxDRESqIFEX01+D708AmFmTyKTvTHLbrYDNMdNlfHV0kGiZVsB+YAvwuJn1AN4CbnT3z+N3YmYjgZEAp556apKliYhIZZJ55GihmS0HlgErzOxfZnZWEtu2kHnxXVUVLZMLFAAPuXs+kSOKkrCduPvD7l7o7oUtW7ZMoiwREUlGMqO5PgZc6+7t3L0t8DPg8STWKwPaxEy3Bj5McpkyoMzd3wzmP00kMEREJE2SCYid7v7PQxPu/jqQTDfTIqCDmbUPTjIPA+bELTMHuNIi+gA73P0jd/83sDl4zClAf8qfuxARkRRLdBXTob/YF5rZH4CniHT/XAG8VtmG3X2/mV0HvEDkMtfH3H2lmY0K2qcAc4lc4rqOyGWusXdtXw9MD8JlQ1ybiIikWKKrmH4bN313zOvDLnsN4+5ziYRA7LwpMa+dSJdV2LpLgcJk9iMiIjUv0VVM/dJZiIiIZJdknkmNmX0H6AJEn1Tg7qNTVZSIiGReMpe5TiFy3uF6IpelDgXaprguERHJsGSuYip29yuBT4PxmYoof2mqiIjUQckExJ7g+24zOwXYB7RPXUkiIpINkjkH8Vwwsur9wNtErmB6JKVViYhIxlUaEO5+b/DyL2b2HJDn7jtSW5aIiGRapQFhZnnAtcA5RI4eXjezh9x9b6qLExGRzEmmi+lJIkNr/D6YHg78D5GrmUREpI5KJiA6uXuPmOl5ZvavVBUkIiLZIZmrmJYEA+kBYGa9gTdSV5KIiGSDRIP1LSdyzqEBkRFXNwVNp6KRVUVE6rxEXUyD0laFiIhknUSD9b1/6HXw2M+vB5P/dHedgxARqeOSGYvpRmA6cGLw9Uczuz7VhYmISGYlcxXTj4De7v45gJmNA+bz1WWvIiJSByVzFZMBB2KmDwTzRESkDkvmCOIx4E0zeyaYvgh4NHUliYhINkgYEGZ2FPAm8P+IDLVhwFXuviQNtYmISAYlDAh3P2hmv3X3IiIjuYqISD2RzDmIF83sUjPTeQcRkXokmXMQtwCNgP1mtpdIN5O7e5OUViYiIhmVzPMgGqejEBERyS4VdjGZWQczm21mK8zsf82sVToLExGRzEp0DuIx4DngUmAJujFORKReSdTF1NjdDz17+n4z01VMIiL1SKKAyDOzfL66a7ph7LS7KzBEROqwRAHxETA+ZvrfMdMOfDNVRYmISOYlGu67XzoLERGR7JLMjXIiIlIPKSBERCSUAkJEREIlulHuBzGv+8a1XZfKokREJPMSHUHcEvM6/ia5q5PZuJkNMLN3zGydmZWEtJuZTQzal5lZQVx7jpktMbPnktmfiIjUnEQBYRW8Dps+fGWzHGAyMBDoDAw3s85xiw0EOgRfI4GH4tpvBFZXti8REal5iQLCK3gdNh2mF7DO3Te4+5fADGBI3DJDgCc9YgHQzMxOBjCz1sB3gKlJ7EtERGpYohvlzjCzZUSOFk4PXhNMn5bEtlsBm2Omy4DeSSzTishNehOAO4CEo8ma2UgiRx+ceuqpSZQlIiLJSBQQZ1Zz22HdUPFHHqHLmNkg4BN3f8vMzku0E3d/GHgYoLCwMJkjGxERSUKigHjE3b9djW2XAW1iplsDHya5zGXAYDO7EMgDmpjZH939B4iISFokOgfRsprbXgR0MLP2ZnY0MAyYE7fMHODK4GqmPsAOd//I3X/h7q3dvV2w3qsKBxGR9Ep0BNHUzC6pqNHdZyXasLvvD+6XeAHIAR5z95VmNiponwLMBS4E1gG7gauOsH4REUmRhAEBDKLicwkJAwLA3ecSCYHYeVNiXjvws0q28RrwWmX7EhGRmpUoIDa5e1I3xImISN2T7H0QIiJSzyQKiO1mVpy2SkREJKskCojZwG/N7D0zG2dmPdNVlIiIZF6FAeHuD7p7EfAN4D/A42a22szuMrOOaatQREQyotLnQbj7++4+zt3zge8BF6MB9ERE6rxKA8LMGpjZd81sOvA3YC1wacorExGRjKrwMlczOx8YTmRE1YVERmMd6e6fp6k2ERHJoET3QfwS+F/gNnf/T5rqERGRLFFhQLh7v3QWIiIi2aXScxAiIlI/KSBERCSUAkJEREIpIEREJJQCQkREQikgREQklAJCRERCKSBERCSUAkJEREIpIEREJJQCQkREQikgREQklAJCRERCKSBERCSUAkJEREIpIEREJJQCQkREQikgREQklAJCRERCKSBERCSUAkJEREIpIEREJJQCQkREQqU0IMxsgJm9Y2brzKwkpN3MbGLQvszMCoL5bcxsnpmtNrOVZnZjKusUEZHDpSwgzCwHmAwMBDoDw82sc9xiA4EOwddI4KFg/n7gVnc/E+gD/CxkXRERSaFUHkH0Ata5+wZ3/xKYAQyJW2YI8KRHLACamdnJ7v6Ru78N4O47gdVAqxTWKiIicVIZEK2AzTHTZRz+IV/pMmbWDsgH3qzxCkVEpEKpDAgLmedHsoyZHQf8BbjJ3T8L3YnZSDNbbGaLt2zZUuViRUSkvFQGRBnQJma6NfBhssuYWQMi4TDd3WdVtBN3f9jdC929sGXLljVSuIiIpDYgFgEdzKy9mR0NDAPmxC0zB7gyuJqpD7DD3T8yMwMeBVa7+/gU1igiIhXITdWG3X2/mV0HvADkAI+5+0ozGxW0TwHmAhcC64DdwFXB6n2BHwLLzWxpMO+X7j43VfWKiEh5KQsIgOADfW7cvCkxrx34Wch6rxN+fkJERNJEd1KLiEgoBYSIiIRSQIiISCgFhIiIhFJAiIhIKAWEiIiEUkCIiEgoBYSIiIRSQIiISCgFhIiIhFJAiIhIKAWEiIiEUkCIiEgoBYSIiIRSQIiISCgFhIiIhFJAiIhIKAWEiIiEUkCIiEgoBYSIiIRSQIiISCgFhIiIhFJAiIhIKAWEiIiEUkCIiEgoBYSIiIRSQIiISCgFhIiIhFJAiIhIKAWEiIiEUkCIiEgoBYSIiIRSQIiISCgFhIiIhEppQJjZADN7x8zWmVlJSLuZ2cSgfZmZFSS7roiIpFbKAsLMcoDJwECgMzDczDrHLTYQ6BB8jQQeOoJ1RUQkhVJ5BNELWOfuG9z9S2AGMCRumSHAkx6xAGhmZicnua6IiKSQuXtqNmx2GTDA3a8Jpn8I9Hb362KWeQ4Y6+6vB9OvAD8H2lW2bsw2RhI5+gDoBLyTkjeUXi2ArZkuogbp/WQ3vZ/slur309bdW4Y15KZwpxYyLz6NKlommXUjM90fBh4+stKym5ktdvfCTNdRU/R+spveT3bL5PtJZUCUAW1iplsDHya5zNFJrCsiIimUynMQi4AOZtbezI4GhgFz4paZA1wZXM3UB9jh7h8lua6IiKRQyo4g3H2/mV0HvADkAI+5+0ozGxW0TwHmAhcC64DdwFWJ1k1VrVmoTnWZofeT7fR+slvG3k/KTlKLiEjtpjupRUQklAJCRERCKSCyhJm1MbN5ZrbazFaa2Y2ZrqkmmFmOmS0J7nmp9cysmZk9bWZrgn+rokzXVFVmdnPwu7bCzJ4ys7xM13SkzOwxM/vEzFbEzDvezF4ys3eD780zWeORqOD93B/8vi0zs2fMrFm66lFAZI/9wK3ufibQB/hZHRle5EZgdaaLqEEPAn939zOAHtTS92ZmrYAbgEJ370rkYpBhma2qSqYBA+LmlQCvuHsH4JVguraYxuHv5yWgq7t3B9YCv0hXMQqILOHuH7n728HrnUQ+eFpltqrqMbPWwHeAqZmupSaYWRPgXOBRAHf/0t23Z7aqaskFGppZLnAstfBeI3f/B/CfuNlDgCeC108AF6W1qGoIez/u/qK77w8mFxC5LywtFBBZyMzaAfnAm5mtpNomAHcABzNdSA05DdgCPB50m001s0aZLqoq3P0D4AFgE/ARkXuQXsxsVTXmv4L7qQi+n5jhemrS1cDf0rUzBUSWMbPjgL8AN7n7Z5mup6rMbBDwibu/lelaalAuUAA85O75wOfUru6LqKBffgjQHjgFaGRmP8hsVZKImf2KSFf09HTtUwGRRcysAZFwmO7uszJdTzX1BQab2XtERuP9ppn9MbMlVVsZUObuh47sniYSGLXRt4CN7r7F3fcBs4DiDNdUUz4ORoUm+P5JhuupNjP7b2AQ8H1P481rCogsYWZGpG97tbuPz3Q91eXuv3D31u7ejsjJz1fdvVb/heru/wY2m1mnYFZ/YFUGS6qOTUAfMzs2+N3rTy094R5iDvDfwev/BmZnsJZqM7MBREa5Huzuu9O5bwVE9ugL/JDIX9pLg68LM12UHOZ6YLqZLQN6AvdluJ4qCY6CngbeBpYT+SyodUNUmNlTwHygk5mVmdmPgLHA+Wb2LnB+MF0rVPB+JgGNgZeCz4UpaatHQ22IiEgYHUGIiEgoBYSIiIRSQIiISCgFhIiIhFJAiIhIKAWE1Elm5mb2PzHTuWa2paqjygajuF4bM31eMtsys2lm9oGZHRNMtwhuHqy2ZGsQqSoFhNRVnwNdzaxhMH0+8EE1ttcMuLbSpcIdIDKGTlYxs5xM1yDZTQEhddnfiIwmCzAceOpQQ/DMgGeDMfYXmFn3YP49wZj8r5nZBjO7IVhlLHB6cKPS/cG842KeDTE9uCM5zATg5mDU1Kj4IwAzm2RmI4LX75nZfWY238wWm1mBmb1gZusPPdc90CR4RsAqM5tiZkcF6387WPdtM5sZjPF1aLt3mdnrwNAj/olKvaKAkLpsBjAseBBOd8qPjvtrYEkwxv4vgSdj2s4ALq8THGEAAAHlSURBVAB6AXcHY2SVAOvdvae73x4slw/cBHQmMtJr3wrq2AS8TuRO+SOx2d2LgH8SeU7AZUSeFTI6ZplewK1AN+B04BIzawHcCXzL3QuAxcAtMevsdfdz3H3GEdYj9Uxu5YuI1E7uviwYOn04MDeu+Rzg0mC5V83sBDNrGrQ97+5fAF+Y2SfAf1Wwi4XuXgZgZkuBdkSCIMx9RMYIev4I3sKc4Pty4LjgOSE7zWxvzFPFFrr7hqCGp4L3tZdIaL0RHNQcTWT4hkP+dAQ1SD2mgJC6bg6R5x6cB5wQMz+sO+jQuDNfxMw7QMX/T5JdDndfF4TI5TGz91P+KD7+kZ+Htn8wbl8HY/YVP1aOE3lvL7n78ArK+byiOkViqYtJ6rrHgNHuvjxu/j+A70PkXACwtZLnb+wkMmBadfwGuC1m+n2gs5kdExy99K/CNnuZWfvg3MMVRI5gFgB9zexrAMGIrR2rWbvUQwoIqdPcvczdHwxpugcoDEZlHctXw0NXtJ1tRLpsVsScpD7SWlYSGT310PRm4M/AMiIPgVlShc3OJ1L/CmAj8Iy7bwFGAE8F728BkfMqIkdEo7mKiEgoHUGIiEgoBYSIiIRSQIiISCgFhIiIhFJAiIhIKAWEiIiEUkCIiEio/w9EfeO0R8nveAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "month_numbers = numpy.array(list(range(vte_probability_no_mask.size + 1))[1:])\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "bar_width = 0.35\n",
    "ax.bar(month_numbers - bar_width/2, vte_probability_no_mask, bar_width, label=\"Without Water Mask\")\n",
    "ax.bar(month_numbers + bar_width/2, vte_probability_masked, bar_width, label=\"With Water Mask\")\n",
    "ax.set_ylabel(\"VTE Probability\")\n",
    "ax.set_xlabel(\"Month Number\")\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}