SusanEAllen/WaterMask.ipynb

## WaterMask.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy\n",
    "from pathlib import Path\n",
    "import rasterio\n",
    "import xarray as xr\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GeoTiff File"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "geotiffs_dir = Path('/Users/sallen/Documents/MIDOSS/ShipTrackDensityGeoTIFFs/')\n",
    "dataset = rasterio.open(geotiffs_dir/'all_2018_09.tif')\n",
    "density = dataset.read()  # just for shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SalishSeaCast Mesh File"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "mesh = xr.open_dataset('/Users/sallen/Documents/MEOPAR/grid/mesh_mask201702.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find the WaterMask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "20\n",
      "40\n",
      "60\n",
      "80\n",
      "100\n",
      "120\n",
      "140\n",
      "160\n",
      "180\n",
      "200\n",
      "220\n",
      "240\n",
      "260\n",
      "280\n",
      "300\n",
      "320\n",
      "340\n",
      "360\n",
      "380\n",
      "400\n"
     ]
    }
   ],
   "source": [
    "coarsen = 20 # for fast testing, production this needs to be 1\n",
    "countpts = numpy.empty_like(density[0])\n",
    "for xx in range(0, density.shape[1], coarsen):\n",
    "    if int(xx/10)*10 == xx:\n",
    "        print(xx)\n",
    "    for yy in range(0, density.shape[2], coarsen):\n",
    "        (llx, lly) = rasterio.transform.xy(dataset.transform, \n",
    "                                       xx+0.5, yy-0.5)\n",
    "        (urx, ury) = rasterio.transform.xy(dataset.transform, \n",
    "                                       xx-0.5, yy+0.5)\n",
    "        inner_points = (numpy.where(mesh.glamt[0] > llx, 1, 0) * \n",
    "                numpy.where(mesh.glamt[0] < urx, 1, 0) *\n",
    "               numpy.where(mesh.gphit[0] > lly, 1, 0) * \n",
    "                numpy.where(mesh.gphit[0] < ury, 1, 0) *\n",
    "                numpy.where(mesh.tmask[0, 0] == 1, 1, 0))\n",
    "        countpts[xx, yy] = inner_points.sum()\n",
    "watermask = numpy.where(countpts > 0, True, False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAARbElEQVR4nO3df6xkZX3H8fen6wIVIYAIwrIV2mxIkSiazaqhbfAHdNkQUWMtm0bxR7JqJIHEP6Q2UdumiWmrbSwEspYNmCBIRZTEVViJDZL4gwtZYHFBtgTLeslukMqPYIDVb/+4Z+P1MrN37szcvXce3q9kMuc85zlznsPJfu7hmXmek6pCktSuP1jqBkiSFpdBL0mNM+glqXEGvSQ1zqCXpMa9bKkb0MshObQO4/ClboYkTYyn+b/Hq+pVvbYty6A/jMN5U96+1M2QpInxvfr6z/tts+tGkhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNW7eoE+yOsn3k+xMcn+Si7vyY5JsS/JQ9350n/3XJ3kwya4kl477BCRJBzbIHf0+4JNV9afAm4FPJDkNuBS4rarWALd1678nyQrgcuBc4DRgY7evJOkgmTfoq+qxqrq7W34a2AmsAs4HrumqXQO8q8fu64BdVfVwVT0PXN/tJ0k6SBbUR5/kZOANwI+B46vqMZj5YwAc12OXVcCjs9Z3d2W9PntTkqkkUy/w3EKaJUk6gIGDPskrgBuBS6rqqUF361HW85FWVbW5qtZW1dqVHDposyRJ8xgo6JOsZCbkr62qb3TFe5Kc0G0/AdjbY9fdwOpZ6ycB08M3V5K0UIP86ibAVcDOqvrirE03Axd2yxcC3+qx+53AmiSnJDkEuKDbT5J0kAxyR38m8H7gbUm2d68NwOeBs5M8BJzdrZPkxCRbAapqH3ARcAszX+LeUFX3L8J5SJL6mHea4qq6g9597QAvmku4qqaBDbPWtwJbh22gJGk0joyVpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcfPORy9JLbll+p6B6/7lia9fxJYcPPMGfZItwHnA3qo6vSv7GnBqV+Uo4FdVdUaPfR8BngZ+A+yrqrVjarckaUCD3NFfDVwGfGV/QVX99f7lJF8AnjzA/m+tqseHbaAkaTSDPErw9iQn99rWPTj8fcDbxtssSdK4jPpl7J8De6rqoT7bC7g1yV1JNo14LEnSEEb9MnYjcN0Btp9ZVdNJjgO2JXmgqm7vVbH7Q7AJ4DBePmKzJEn7DX1Hn+RlwHuAr/WrU1XT3fte4CZg3QHqbq6qtVW1diWHDtssSdIco3TdvAN4oKp299qY5PAkR+xfBs4BdoxwPEnSEOYN+iTXAT8ETk2yO8lHuk0XMKfbJsmJSbZ2q8cDdyS5B/gJ8O2q+u74mi5JGsQgv7rZ2Kf8gz3KpoEN3fLDQBujDSRpgjkyVtJLykJGuw46ina5j6B1rhtJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMGeZTgliR7k+yYVfa5JL9Isr17beiz7/okDybZleTScTZckjSYQe7orwbW9yj/t6o6o3ttnbsxyQrgcuBc4DRgY5LTRmmsJGnh5g36qrodeGKIz14H7Kqqh6vqeeB64PwhPkeSNIJR+ugvSnJv17VzdI/tq4BHZ63v7sp6SrIpyVSSqRd4boRmSZJmG/bh4FcA/whU9/4F4MNz6qTHftXvA6tqM7AZ4Mgc07eepOVn0IdoLwfL/UHei2GoO/qq2lNVv6mq3wJfZqabZq7dwOpZ6ycB08McT5I0vKGCPskJs1bfDezoUe1OYE2SU5IcAlwA3DzM8SRJw5u36ybJdcBZwLFJdgOfBc5KcgYzXTGPAB/t6p4I/GdVbaiqfUkuAm4BVgBbqur+RTkLSVJf8wZ9VW3sUXxVn7rTwIZZ61uBF/30UpJ08DgyVpIaZ9BLUuMMeklqnEEvSY0z6CWpccOOjJU0gRYygnUhI0gXY7TpYo22Xaz/BsuZd/SS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGucUCNJLyFIP6V+s6Qcm6eHkS2HeO/okW5LsTbJjVtm/JHkgyb1JbkpyVJ99H0lyX5LtSabG2XBJ0mAG6bq5Glg/p2wbcHpVvQ74GfC3B9j/rVV1RlWtHa6JkqRRzBv0VXU78MScslural+3+iPgpEVomyRpDMbxZeyHge/02VbArUnuSrLpQB+SZFOSqSRTL/DcGJolSYIRv4xN8nfAPuDaPlXOrKrpJMcB25I80P0fwotU1WZgM8CROaZGaZck6XeGvqNPciFwHvA3VdUzmKtqunvfC9wErBv2eJKk4QwV9EnWA58C3llVz/apc3iSI/YvA+cAO3rVlSQtnkF+Xnkd8EPg1CS7k3wEuAw4gpnumO1Jruzqnphka7fr8cAdSe4BfgJ8u6q+uyhnIUnqa94++qra2KP4qj51p4EN3fLDQBsPXJSkCebIWEnL0kvxId6LxbluJKlxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcY6MlV5CfLbqS5N39JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxgzxKcEuSvUl2zCo7Jsm2JA9170f32Xd9kgeT7Epy6TgbLkkazCB39FcD6+eUXQrcVlVrgNu69d+TZAVwOXAucBqwMclpI7VWkrRg8wZ9Vd0OPDGn+Hzgmm75GuBdPXZdB+yqqoer6nng+m4/SdJBNGwf/fFV9RhA935cjzqrgEdnre/uynpKsinJVJKpF3huyGZJkuZazCkQ0qOs+lWuqs3AZoAjc0zfepJ+32I9RNuHc7dj2Dv6PUlOAOje9/aosxtYPWv9JGB6yONJkoY0bNDfDFzYLV8IfKtHnTuBNUlOSXIIcEG3nyTpIBrk55XXAT8ETk2yO8lHgM8DZyd5CDi7WyfJiUm2AlTVPuAi4BZgJ3BDVd2/OKchSepn3j76qtrYZ9Pbe9SdBjbMWt8KbB26dZKkkTkyVpIaZ9BLUuMMeklqnEEvSY0z6CWpcT4cXFqGlsNDvBdjFK0jaJeGd/SS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4R8ZKB9FyGPG6GBzxurx5Ry9JjRs66JOcmmT7rNdTSS6ZU+esJE/OqvOZ0ZssSVqIobtuqupB4AyAJCuAXwA39aj6g6o6b9jjSJJGM66um7cD/1NVPx/T50mSxmRcQX8BcF2fbW9Jck+S7yR5bb8PSLIpyVSSqRd4bkzNkiSNHPRJDgHeCfxXj813A6+pqtcD/wF8s9/nVNXmqlpbVWtXcuiozZIkdcZxR38ucHdV7Zm7oaqeqqpnuuWtwMokx47hmJKkAY0j6DfSp9smyauTpFte1x3vl2M4piRpQCMNmErycuBs4KOzyj4GUFVXAu8FPp5kH/Br4IKqqlGOKUlamJGCvqqeBV45p+zKWcuXAZeNcgxJ0micAkE6iJwqQEvBKRAkqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxjoyVJtxCHjjuyNyXJu/oJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuNGCvokjyS5L8n2JFM9tifJl5LsSnJvkjeOcjxJ0sKN43f0b62qx/tsOxdY073eBFzRvUuSDpLF7ro5H/hKzfgRcFSSExb5mJKkWUYN+gJuTXJXkk09tq8CHp21vrsre5Ekm5JMJZl6gedGbJYkab9Ru27OrKrpJMcB25I8UFW3z9qeHvtUrw+qqs3AZoAjc0zPOpJezGkNNJ+R7uirarp73wvcBKybU2U3sHrW+knA9CjHlCQtzNBBn+TwJEfsXwbOAXbMqXYz8IHu1zdvBp6sqseGbq0kacFG6bo5Hrgpyf7P+WpVfTfJxwCq6kpgK7AB2AU8C3xotOZKkhZq6KCvqoeBF3UOdgG/f7mATwx7DEnS6BwZK0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0b5Zmxq5N8P8nOJPcnubhHnbOSPJlke/f6zGjNlSQt1CjPjN0HfLKq7u4eEn5Xkm1V9dM59X5QVeeNcBxJ0giGvqOvqseq6u5u+WlgJ7BqXA2TJI3HWProk5wMvAH4cY/Nb0lyT5LvJHntAT5jU5KpJFMv8Nw4miVJYrSuGwCSvAK4Ebikqp6as/lu4DVV9UySDcA3gTW9PqeqNgObAY7MMTVquyRJM0a6o0+ykpmQv7aqvjF3e1U9VVXPdMtbgZVJjh3lmJKkhRnlVzcBrgJ2VtUX+9R5dVePJOu64/1y2GNKkhZulK6bM4H3A/cl2d6VfRr4I4CquhJ4L/DxJPuAXwMXVJXdMpJ0EA0d9FV1B5B56lwGXDbsMSRJo3NkrCQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDVu1IeDr0/yYJJdSS7tsT1JvtRtvzfJG0c5niRp4UZ5OPgK4HLgXOA0YGOS0+ZUOxdY0702AVcMezxJ0nBGuaNfB+yqqoer6nngeuD8OXXOB75SM34EHJXkhBGOKUlaoKEfDg6sAh6dtb4beNMAdVYBj839sCSbmLnrB3jue/X1HSO0bbk6Fnh8qRuxCDyvyeJ5TZZBz+s1/TaMEvTpUVZD1JkprNoMbAZIMlVVa0do27LkeU0Wz2uyeF79jdJ1sxtYPWv9JGB6iDqSpEU0StDfCaxJckqSQ4ALgJvn1LkZ+ED365s3A09W1Yu6bSRJi2forpuq2pfkIuAWYAWwparuT/KxbvuVwFZgA7ALeBb40IAfv3nYdi1zntdk8bwmi+fVR6p6dplLkhrhyFhJapxBL0mNW1ZBP9+UCpMqySNJ7kuyPcnUUrdnWEm2JNmbZMessmOSbEvyUPd+9FK2cRh9zutzSX7RXbPtSTYsZRuHkWR1ku8n2Znk/iQXd+UTfc0OcF4Tfc2SHJbkJ0nu6c7r77vyka/Xsumj76ZU+BlwNjM/y7wT2FhVP13Sho1BkkeAtVU10YM5kvwF8Awzo51P78r+GXiiqj7f/XE+uqo+tZTtXKg+5/U54Jmq+telbNsoulHoJ1TV3UmOAO4C3gV8kAm+Zgc4r/cxwdcsSYDDq+qZJCuBO4CLgfcw4vVaTnf0g0ypoCVUVbcDT8wpPh+4plu+hpl/cBOlz3lNvKp6rKru7pafBnYyMzJ9oq/ZAc5ronVTxTzTra7sXsUYrtdyCvp+0yW0oIBbk9zVTfXQkuP3j43o3o9b4vaM00XdrKtbJq17Y64kJwNvAH5MQ9dsznnBhF+zJCuSbAf2AtuqaizXazkF/cDTJUygM6vqjczM5vmJrqtAy9sVwJ8AZzAzN9MXlrY5w0vyCuBG4JKqemqp2zMuPc5r4q9ZVf2mqs5gZhaBdUlOH8fnLqegb3a6hKqa7t73Ajcx003Vij37ZyTt3vcucXvGoqr2dP/ofgt8mQm9Zl1f743AtVX1ja544q9Zr/Nq5ZoBVNWvgP8G1jOG67Wcgn6QKRUmTpLDuy+MSHI4cA7Q0sycNwMXdssXAt9awraMzZzptN/NBF6z7su9q4CdVfXFWZsm+pr1O69Jv2ZJXpXkqG75D4F3AA8whuu1bH51A9D9HOrf+d2UCv+0xE0aWZI/ZuYuHmamnPjqpJ5XkuuAs5iZNnUP8Fngm8ANwB8B/wv8VVVN1Bebfc7rLGa6AAp4BPjopM3TlOTPgB8A9wG/7Yo/zUx/9sReswOc10Ym+JoleR0zX7auYOYm/Iaq+ockr2TE67Wsgl6SNH7LqetGkrQIDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUuP8HbrNWbRB1t54AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.pcolormesh(watermask[::coarsen, ::coarsen]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Image when run with coarsen=1"
   ]
  },
  {
   "attachments": {
    "watermask.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD5CAYAAADcDXXiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAY3UlEQVR4nO3dXcwcV33H8e8Px7EhL0rcvNSxrdogU8lBjYMshyoVSgkQExAOF0FGgrpSKnORCNIigQ1SgYtIaVVeetFEMpDilhdjBdJYUURwDBE3aYwNTrDz0hjiJk/s2gWCCL0wifPvxZ6FyXp2Z3Z35tmd2d9HerSzZ2dmz+Hlt8dnzpxRRGBmZu3ymklXwMzMqudwNzNrIYe7mVkLOdzNzFrI4W5m1kIOdzOzFjqr7I6SFgD7gecj4j2SlgDfAlYCR4H3R8QLad9twE3AaeAjEfHAoHOfrUWxmHNGaoCZ2ax6kRd+EREX531WOtyBjwJPAOen91uBvRFxu6St6f0nJK0BNgGXA5cBD0p6Y0Sc7nfixZzDVbp2iKqYmdmDcfd/9/us1LCMpOXAu4EvZ4o3AjvS9g7ghkz5zog4FRHPAEeA9cNW2szMRld2zP2LwMeBVzJll0bEcYD0ekkqXwY8l9lvLpWZmdk8KQx3Se8BTkbEgZLnVE7ZGWscSNoiab+k/S9xquSpzcysjDJj7lcD75V0PbAYOF/S14ATkpZGxHFJS4GTaf85YEXm+OXAsd6TRsR2YDvA+VriBW7MzCpU2HOPiG0RsTwiVtK5UPr9iPggsBvYnHbbDNybtncDmyQtkrQKWA3sq7zmZmbW1zCzZXrdDuySdBPwLHAjQEQclrQLeBx4Gbh50EwZMzOrnqZhyd/ztSQ8FdLMbDgPxt0HImJd3me+Q9XMrIUc7mZmLeRwNzNrIYe7mVkLOdzNzFrI4W5m1kIOdzOzFnK4m5m1kMPdzKyFHO5mZi3kcDczayGHu5lZCznczcxayOFuZtZCDnczsxZyuJuZtZDD3cyshRzuZmYtVBjukhZL2ifpUUmHJX02lX9G0vOSDqa/6zPHbJN0RNJTkq6rswFmZnamMj33U8DbIuIKYC2wQdJb0mdfiIi16e9+AElrgE3A5cAG4A5JC2qou5nNsAeOPZq7bR1nFe0QnSdo/za9XZj+Bj1VeyOwMyJOAc9IOgKsBx4es65mZr8P8usuu6Iw1LP7zppSY+6SFkg6CJwE9kTEI+mjWyQ9JukuSRemsmXAc5nD51KZmdnYskGd3e4Gee9rd7v7NytKhXtEnI6ItcByYL2kNwF3Am+gM1RzHPhc2l15p+gtkLRF0n5J+1/i1EiVN7PZ1O2194Z1XrB39+/dp+2Gmi0TEb8GHgI2RMSJFPqvAF+iM/QCnZ76isxhy4FjOefaHhHrImLdQhaNVHkzm13XXXbFGcMtRb36WVJmtszFki5I268F3g48KWlpZrf3AYfS9m5gk6RFklYBq4F91VbbzKyjG+LZMfhsWd7+sxD2hRdUgaXAjjTj5TXAroi4T9K/S1pLZ8jlKPBhgIg4LGkX8DjwMnBzRJyupfZmZkPKDt20+UJrmdkyjwFX5pR/aMAxtwG3jVc1M7Px9OuhZ8fs83r8beA7VM2s8XpDOS+su2P0/fZtG4e7mTVCmRAet+fdpimT6tyjNFnna0lcpWsnXQ0zm0KDgjYvzPuNpef15nuHZcqcf5o8GHcfiIh1eZ+VuaBqZjYxg8bEs6FcFMRFxxYFfdM43M1s6vUL3by57GV723nnHObHYtp5zN3MGqNoiKaqQG56sIPD3cwaohu4VQ+b5J2vDUMzDncza4yqe9T9zueeu5nZBFTVs85baKwNwQ4OdzNrmHGGZ9oyh70Mz3M3s0Yap7c9HzdEzYdB89zdczezRhp3Sd8qZ9dMI89zN7PGyg7RDHpAR9ag1SDbFPbuuZtZ4/VbbqBfj37Y8iZyuJtZ65Tpgbc94B3uZjZTsmPtbQnyPA53M2uF3tAuCvB+67u3hcPdzFpjlKGWtvbeyzwge7GkfZIelXRY0mdT+RJJeyQ9nV4vzByzTdIRSU9Juq7OBpiZ9dNvCKatgZ5Vpud+CnhbRFwBrAU2SHoLsBXYGxGrgb3pPZLWAJuAy4ENwB3p4dpmZrXqN9Qy7Bh7G+5kLQz36Phtersw/QWwEdiRyncAN6TtjcDOiDgVEc8AR4D1ldbazGxIs3ARNavUmLukBZIOAieBPRHxCHBpRBwHSK+XpN2XAc9lDp9LZWZmU2EWgr7UHaoRcRpYK+kC4B5Jbxqwu/JOccZO0hZgC8BiXlemGmZmI+m9K7XsU52abKjZMhHxa+AhOmPpJyQtBUivJ9Nuc8CKzGHLgWM559oeEesiYt1CFo1QdTOzcoZ9NuqgZQ2aosxsmYtTjx1JrwXeDjwJ7AY2p902A/em7d3AJkmLJK0CVgP7qq64mVlZg+a6Z1/77d/EgC/Tc18K/EDSY8CP6Iy53wfcDrxD0tPAO9J7IuIwsAt4HPgucHMa1jEzq9Uwj8zr9ubL3sTUtIAvHHOPiMeAK3PKfwnkLsIeEbcBt41dOzOzHN2gzRsn7wZ27z5F678XPXy7aeHuO1TNrJHygj1vjnuZu1aLgr33+5rA67mbWaPkPfe07DF574vG45vK4W5mE1EUzqOEa3YMfdBQSt5nTRx6GcThbmYjK7ogOU5YjnpsmePyxt/bFOzgcDezAkWh15RQHHSTUlPaMAxfUDWzXE2+gadXv2Bv6sXSMtxzN7Pfa0uYd3XH0cu2a9i579PM4W42w9oW5r2KnsjU5uEZRZyxpte8O19L4irl3g9lZhVqU3gNo9+c+Oz7QcdMqwfj7gMRsS7vM/fczWbApEN90tMMe6dIZst7NSHUy3C4m7XYJAN1mkOyyb31sjxbxqxluhcQ5zPYy8w6mXRwFi0xMOn6Vc1j7mYtMF9BXlUATnqYqKvpge4xd7MWqzso6wjAQVMO62pP2560VMThbtZgVYfUfPZkq/6uoou2vT8mvfs3vRffy+Fu1lBVBfu0hFqZ9hTdaVrm2LxVJaflP4MqOdzNZtAkwmzQnZ/jrhBZ5q7SWRiKySrzDNUVkn4g6QlJhyV9NJV/RtLzkg6mv+szx2yTdETSU5Kuq7MBZrNm2JkwvaE3yV7qfM3i6dfGUX9cmqhMz/1l4GMR8WNJ5wEHJO1Jn30hIv4pu7OkNcAm4HLgMuBBSW/0c1TN5s80DjNkx7jHWb+l33K9bXtM3rjKPEP1OHA8bb8o6Qlg2YBDNgI7I+IU8IykI8B64OEK6ms200a9+WYaw76r7PrreQHdG/T9ntJU5ulNbVkwrGuom5gkraTzsOxHUtEtkh6TdJekC1PZMuC5zGFzDP4xMLMS+t0q3+RAGudReVm9s2BG0eT/HPOUDndJ5wLfBm6NiN8AdwJvANbS6dl/rrtrzuFn3CklaYuk/ZL2v8SpoStuNkvasAZKXm963DZ0f9y6wzLZO2Xzzj3sWHyTlQp3SQvpBPvXI+I7ABFxIiJOR8QrwJfoDL1Ap6e+InP4cuBY7zkjYntErIuIdQtZNE4bzFot71mfTQyjvAu7o46DFy3jmzcGP2vj7oXLD0gSsAP4VUTcmilfmsbjkfS3wFURsUnS5cA36IT9ZcBeYPWgC6pefsCs/aY1WJv4Q9k17vIDVwMfAn4q6WAq+yTwAUlr6Qy5HAU+DBARhyXtAh6nM9PmZs+UMZtNeb3nfp/1qqOn3e+ibJMDvh8vHGZm86LsTUzZAO7df9yw7x3OaXqoe+EwM5sKvQFfFNb9pkCOq+mhXobD3cxql9dr7xfYdTzXdBbCvJcf1mFmtSrzYGqrnnvuZlaLYde/qSPwZ7HH3uVwN7PKDTvHvK6efFtnwpThYRkzq0y/BbzmeximzDNd287hbmaVGHdN9iqVWSis7RzuZja2+XpYxrBLL7jnbmY2hmlct2WWgx0c7mZWkWFWYayzDk1dWK1qni1jZpWYzwup/W6ImuXZMb3cczezsQ1aIKzO78l+n3vsr+aeu5kNbVDIVvkQ7LywbsuiX3XzqpBmNpR+vfS6L6gWrU0zi2E/aFVIh7uZlTaJm5GKnp3au1zwLHG4m1llJjXlsRvkZdaE7+7fdl7P3cxGNskwL1M26LNZnj3jcDezV6nrARmj1CNr2JCedP0nrTDcJa0A/g34Y+AVYHtE/LOkJcC3gJV0nqH6/oh4IR2zDbgJOA18JCIeqKX2ZlaJOh6QMWmz2mPvKtNzfxn4WET8WNJ5wAFJe4C/BvZGxO2StgJbgU9IWgNsAi4HLgMelPRGPyTbbLo0JcBHCelZHo7pKgz3iDgOHE/bL0p6AlgGbASuSbvtAB4CPpHKd0bEKeAZSUeA9cDDVVfezIpNc4j3e57qOMHsYO8Y6g5VSSuBK4FHgEtT8Hd/AC5Juy0DnsscNpfKzGyedG8kmqZg73fBs8plCxzsf1D6gqqkc4FvA7dGxG8k9d01p+yM+ZaStgBbABbzurLVMLM+pinIu3qDtt8NT1XU3cH+aqV67pIW0gn2r0fEd1LxCUlL0+dLgZOpfA5YkTl8OXCs95wRsT0i1kXEuoUsGrX+ZjalJrEipP1BmdkyAr4CPBERn898tBvYDNyeXu/NlH9D0ufpXFBdDeyrstJmVl1Pvei2/nHPNx/caz9TmWGZq4EPAT+VdDCVfZJOqO+SdBPwLHAjQEQclrQLeJzOTJubPVPGrFrZxbPytov0u5CZ/Xw+HptXVI8iXkSsPy8/YNYgvTNK8sKtip530Xou4ypaVbKq72k7Lz9g1iKDlsEd93xFM2yqCtx+P0bDtMNDMYO5527WEFXPhhl2rH0+gnQWF/8ah3vuZvZ7VV9ArdJ8rQ0/CxzuZlOuyqAbJdgn2Xv20Mvo/AxVsxkx6jCMe9HN5J672ZQbdybMsD3faegpZ39YPHtmNL6gatYQwwT7uDNq8nr5Dtfp4wuqZg02bG992OGUvHnt/ebPz/ozS5vE4W42xcZZHTFPvzAuu6CXh0iaw8MyZlOozhky49TDPffpMmhYxrNlzKbItK3BnseB3gwOd7MpUUeoVxnEVS97YPVyuJtNgWkN9qJAdy9+ejnczSasjjVj6uqx1/EEJauHZ8uYTVBdD9yompfibR733M0mZNybkrrlDlvL43A3m7Ay4dzvaUnzqXdJAA/JTDcPy5hNwCiPxsua1DDJoAd7+F8Q06Ww5y7pLkknJR3KlH1G0vOSDqa/6zOfbZN0RNJTkq6rq+JmTVUUyr3lg2asTDJQe4eE3JOfLmWGZb4KbMgp/0JErE1/9wNIWgNsAi5Px9whaUFVlTVrujKzTbKLdPW+5u07TaE6TXWZdYXDMhHxQ0krS55vI7AzIk4Bz0g6AqwHHh65hmYtkDd8UeVj7aZh1cZJf7+92jhj7rdI+itgP/CxiHgBWAb8Z2afuVRmNvP6Dbf0670POibvuPkOeI+3T7dRZ8vcCbwBWAscBz6XypWzb+7KZJK2SNovaf9LnBqxGmbNUbRS46Cx9ezroPH2uodpemfKeCrm9Bqp5x4RJ7rbkr4E3JfezgErMrsuB471Ocd2YDt0VoUcpR5mTVI2BHtnwuT1yPv1mudjzNth3gwj9dwlLc28fR/QnUmzG9gkaZGkVcBqYN94VTRrtt7e96DPh3kCUplhnao52JujsOcu6ZvANcBFkuaATwPXSFpLZ8jlKPBhgIg4LGkX8DjwMnBzRJyup+pm029QyPZ72tEwsr18j4Fblh/WYVaT3rDtDfOicB/nBiebDX5Yh9kU6b0g2W/YZpTZLw526/LyA2YVKtPTzpve2K/HPijg8x6BZ9blcDcrYdwVHIc5Phvo3V5+3jnL3gxls8nhbpYxSs+799iywyllQ7noXN2eftE8epstDnebecMup9sN796hlN7edvY8vb3sMjNlhv3XQr9zevhmNnm2jM2kUYZZ+t0ZWnTH6KBz9tu3jgD2VMn2GTRbxuFuM613fHvQnaBQPCQz6POyPxK9qg7iaVgu2KrhcDer2TDhPsz5uhzElsfz3M1qVvXMFYe5jcvhbjYPRgnrogujZoM43M0qUkdv24+ys1E53M0aoOhxe2a9HO5mDVC0bLBZL4e7WYW8JIBNC4e7WQ1GXa530Pm6/MNhZXj5AbOaVDU+7jnvNgr33M3MWqgw3CXdJemkpEOZsiWS9kh6Or1emPlsm6Qjkp6SdF1dFTebBf0WJzMrUqbn/lVgQ0/ZVmBvRKwG9qb3SFoDbAIuT8fcIWlBZbU1m2G+kcmGURjuEfFD4Fc9xRuBHWl7B3BDpnxnRJyKiGeAI8D6iupqNnO6vXX32G1Yo465XxoRxwHS6yWpfBnwXGa/uVRmZhVxD97KqPqCqnLKcpedlLRF0n5J+1/iVMXVMGsf995tGKNOhTwhaWlEHJe0FDiZyueAFZn9lgPH8k4QEduB7dBZ8nfEepjNFAe8lTVqz303sDltbwbuzZRvkrRI0ipgNbBvvCqamdmwCnvukr4JXANcJGkO+DRwO7BL0k3As8CNABFxWNIu4HHgZeDmiDhdU93NzKyPwnCPiA/0+Sj30UkRcRtw2ziVMjOz8fgOVTOzFnK4m5m1kMPdzKyFHO5mZi3kcDczayGHu5lZCznczcxayOFuZtZCDnczsxZyuJuZtZDD3cyshRzuZmYt5HA3M2shh7uZWQs53M3MWsjhbmbWQg53M7MWcribmbVQ4WP2BpF0FHgROA28HBHrJC0BvgWsBI4C74+IF8arppmZDaOKnvtfRsTaiFiX3m8F9kbEamBvem9mZvOojmGZjcCOtL0DuKGG7zAzswHGDfcAvifpgKQtqezSiDgOkF4vGfM7zMxsSGONuQNXR8QxSZcAeyQ9WfbA9GOwBWAxrxuzGmZmljVWzz0ijqXXk8A9wHrghKSlAOn1ZJ9jt0fEuohYt5BF41TDzMx6jBzuks6RdF53G3gncAjYDWxOu20G7h23kmZmNpxxhmUuBe6R1D3PNyLiu5J+BOySdBPwLHDj+NU0M7NhjBzuEfFz4Iqc8l8C145TKTMzG4/vUDUzayGHu5lZCznczcxayOFuZtZCDnczsxZyuJuZtZDD3cyshRzuZmYt5HA3M2shh7uZWQs53M3MWsjhbmbWQg53M7MWcribmbWQw93MrIUc7mZmLeRwNzNrIYe7mVkL1RbukjZIekrSEUlb6/oeMzM7Uy3hLmkB8C/Au4A1wAckranju8zM7Ex19dzXA0ci4ucR8TtgJ7Cxpu8yM7MedYX7MuC5zPu5VGZmZvPgrJrOq5yyeNUO0hZgS3p76sG4+1BNdZmki4BfTLoSNXC7msXtapZh2vUn/T6oK9zngBWZ98uBY9kdImI7sB1A0v6IWFdTXSbG7WoWt6tZ3K7B6hqW+RGwWtIqSWcDm4DdNX2XmZn1qKXnHhEvS7oFeABYANwVEYfr+C4zMztTXcMyRMT9wP0ld99eVz0mzO1qFrerWdyuARQRxXuZmVmjePkBM7MWmni4N3mZAkl3STop6VCmbImkPZKeTq8XZj7bltr5lKTrJlPrYpJWSPqBpCckHZb00VTe2LZJWixpn6RHU5s+m8ob26YsSQsk/UTSfel9W9p1VNJPJR2UtD+VNb5tki6QdLekJ9P/z/688nZFxMT+6Fxs/RnweuBs4FFgzSTrNGT93wq8GTiUKftHYGva3gr8Q9pek9q3CFiV2r1g0m3o066lwJvT9nnAf6X6N7ZtdO69ODdtLwQeAd7S5Db1tO/vgG8A97Xlf4epvkeBi3rKGt82YAfwN2n7bOCCqts16Z57o5cpiIgfAr/qKd5I57840usNmfKdEXEqIp4BjtBp/9SJiOMR8eO0/SLwBJ07jBvbtuj4bXq7MP0FDW5Tl6TlwLuBL2eKG9+uARrdNknn0+kYfgUgIn4XEb+m4nZNOtzbuEzBpRFxHDohCVySyhvZVkkrgSvp9HQb3bY0dHEQOAnsiYjGtyn5IvBx4JVMWRvaBZ0f4O9JOpDuaofmt+31wP8C/5qG0r4s6Rwqbtekw71wmYIWaVxbJZ0LfBu4NSJ+M2jXnLKpa1tEnI6ItXTumF4v6U0Ddm9EmyS9BzgZEQfKHpJTNnXtyrg6It5MZ4XZmyW9dcC+TWnbWXSGc++MiCuB/6MzDNPPSO2adLgXLlPQQCckLQVIrydTeaPaKmkhnWD/ekR8JxW3om3pn8APARtofpuuBt4r6SidYc23SfoazW8XABFxLL2eBO6hMxzR9LbNAXPpX44Ad9MJ+0rbNelwb+MyBbuBzWl7M3BvpnyTpEWSVgGrgX0TqF8hSaIzHvhERHw+81Fj2ybpYkkXpO3XAm8HnqTBbQKIiG0RsTwiVtL5/8/3I+KDNLxdAJLOkXRedxt4J3CIhrctIv4HeE7Sn6aia4HHqbpdU3DV+Ho6szF+Bnxq0vUZsu7fBI4DL9H5db0J+CNgL/B0el2S2f9TqZ1PAe+adP0HtOsv6Pyz7zHgYPq7vsltA/4M+Elq0yHg71N5Y9uU08Zr+MNsmca3i87Y9KPp73A3H1rStrXA/vS/x/8ALqy6Xb5D1cyshSY9LGNmZjVwuJuZtZDD3cyshRzuZmYt5HA3M2shh7uZWQs53M3MWsjhbmbWQv8PQlOHMOK5OscAAAAASUVORK5CYII="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![watermask.png](attachment:watermask.png)"
   ]
  },
  {
   "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
}
	{
	"cells": [
	{
	"cell_type": "code",
	"execution_count": 1,
	"metadata": {},
	"outputs": [],
	"source": [
	"import matplotlib.pyplot as plt\n",
	"import numpy\n",
	"from pathlib import Path\n",
	"import rasterio\n",
	"import xarray as xr\n",
	"\n",
	"%matplotlib inline"
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"### GeoTiff File"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 5,
	"metadata": {},
	"outputs": [],
	"source": [
	"geotiffs_dir = Path('/Users/sallen/Documents/MIDOSS/ShipTrackDensityGeoTIFFs/')\n",
	"dataset = rasterio.open(geotiffs_dir/'all_2018_09.tif')\n",
	"density = dataset.read() # just for shape"
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"### SalishSeaCast Mesh File"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 6,
	"metadata": {},
	"outputs": [],
	"source": [
	"mesh = xr.open_dataset('/Users/sallen/Documents/MEOPAR/grid/mesh_mask201702.nc')"
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"### Find the WaterMask"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 7,
	"metadata": {},
	"outputs": [
	{
	"name": "stdout",
	"output_type": "stream",
	"text": [
	"0\n",
	"20\n",
	"40\n",
	"60\n",
	"80\n",
	"100\n",
	"120\n",
	"140\n",
	"160\n",
	"180\n",
	"200\n",
	"220\n",
	"240\n",
	"260\n",
	"280\n",
	"300\n",
	"320\n",
	"340\n",
	"360\n",
	"380\n",
	"400\n"
	]
	}
	],
	"source": [
	"coarsen = 20 # for fast testing, production this needs to be 1\n",
	"countpts = numpy.empty_like(density[0])\n",
	"for xx in range(0, density.shape[1], coarsen):\n",
	" if int(xx/10)*10 == xx:\n",
	" print(xx)\n",
	" for yy in range(0, density.shape[2], coarsen):\n",
	" (llx, lly) = rasterio.transform.xy(dataset.transform, \n",
	" xx+0.5, yy-0.5)\n",
	" (urx, ury) = rasterio.transform.xy(dataset.transform, \n",
	" xx-0.5, yy+0.5)\n",
	" inner_points = (numpy.where(mesh.glamt[0] > llx, 1, 0) * \n",
	" numpy.where(mesh.glamt[0] < urx, 1, 0) *\n",
	" numpy.where(mesh.gphit[0] > lly, 1, 0) * \n",
	" numpy.where(mesh.gphit[0] < ury, 1, 0) *\n",
	" numpy.where(mesh.tmask[0, 0] == 1, 1, 0))\n",
	" countpts[xx, yy] = inner_points.sum()\n",
	"watermask = numpy.where(countpts > 0, True, False)"
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"### Plot"
	]
	},
	{
	"cell_type": "code",
	"execution_count": 8,
	"metadata": {},
	"outputs": [
	{
	"data": {
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAARbElEQVR4nO3df6xkZX3H8fen6wIVIYAIwrIV2mxIkSiazaqhbfAHdNkQUWMtm0bxR7JqJIHEP6Q2UdumiWmrbSwEspYNmCBIRZTEVViJDZL4gwtZYHFBtgTLeslukMqPYIDVb/+4Z+P1MrN37szcvXce3q9kMuc85zlznsPJfu7hmXmek6pCktSuP1jqBkiSFpdBL0mNM+glqXEGvSQ1zqCXpMa9bKkb0MshObQO4/ClboYkTYyn+b/Hq+pVvbYty6A/jMN5U96+1M2QpInxvfr6z/tts+tGkhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNW7eoE+yOsn3k+xMcn+Si7vyY5JsS/JQ9350n/3XJ3kwya4kl477BCRJBzbIHf0+4JNV9afAm4FPJDkNuBS4rarWALd1678nyQrgcuBc4DRgY7evJOkgmTfoq+qxqrq7W34a2AmsAs4HrumqXQO8q8fu64BdVfVwVT0PXN/tJ0k6SBbUR5/kZOANwI+B46vqMZj5YwAc12OXVcCjs9Z3d2W9PntTkqkkUy/w3EKaJUk6gIGDPskrgBuBS6rqqUF361HW85FWVbW5qtZW1dqVHDposyRJ8xgo6JOsZCbkr62qb3TFe5Kc0G0/AdjbY9fdwOpZ6ycB08M3V5K0UIP86ibAVcDOqvrirE03Axd2yxcC3+qx+53AmiSnJDkEuKDbT5J0kAxyR38m8H7gbUm2d68NwOeBs5M8BJzdrZPkxCRbAapqH3ARcAszX+LeUFX3L8J5SJL6mHea4qq6g9597QAvmku4qqaBDbPWtwJbh22gJGk0joyVpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcfPORy9JLbll+p6B6/7lia9fxJYcPPMGfZItwHnA3qo6vSv7GnBqV+Uo4FdVdUaPfR8BngZ+A+yrqrVjarckaUCD3NFfDVwGfGV/QVX99f7lJF8AnjzA/m+tqseHbaAkaTSDPErw9iQn99rWPTj8fcDbxtssSdK4jPpl7J8De6rqoT7bC7g1yV1JNo14LEnSEEb9MnYjcN0Btp9ZVdNJjgO2JXmgqm7vVbH7Q7AJ4DBePmKzJEn7DX1Hn+RlwHuAr/WrU1XT3fte4CZg3QHqbq6qtVW1diWHDtssSdIco3TdvAN4oKp299qY5PAkR+xfBs4BdoxwPEnSEOYN+iTXAT8ETk2yO8lHuk0XMKfbJsmJSbZ2q8cDdyS5B/gJ8O2q+u74mi5JGsQgv7rZ2Kf8gz3KpoEN3fLDQBujDSRpgjkyVtJLykJGuw46ina5j6B1rhtJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMGeZTgliR7k+yYVfa5JL9Isr17beiz7/okDybZleTScTZckjSYQe7orwbW9yj/t6o6o3ttnbsxyQrgcuBc4DRgY5LTRmmsJGnh5g36qrodeGKIz14H7Kqqh6vqeeB64PwhPkeSNIJR+ugvSnJv17VzdI/tq4BHZ63v7sp6SrIpyVSSqRd4boRmSZJmG/bh4FcA/whU9/4F4MNz6qTHftXvA6tqM7AZ4Mgc07eepOVn0IdoLwfL/UHei2GoO/qq2lNVv6mq3wJfZqabZq7dwOpZ6ycB08McT5I0vKGCPskJs1bfDezoUe1OYE2SU5IcAlwA3DzM8SRJw5u36ybJdcBZwLFJdgOfBc5KcgYzXTGPAB/t6p4I/GdVbaiqfUkuAm4BVgBbqur+RTkLSVJf8wZ9VW3sUXxVn7rTwIZZ61uBF/30UpJ08DgyVpIaZ9BLUuMMeklqnEEvSY0z6CWpccOOjJU0gRYygnUhI0gXY7TpYo22Xaz/BsuZd/SS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4g16SGucUCNJLyFIP6V+s6Qcm6eHkS2HeO/okW5LsTbJjVtm/JHkgyb1JbkpyVJ99H0lyX5LtSabG2XBJ0mAG6bq5Glg/p2wbcHpVvQ74GfC3B9j/rVV1RlWtHa6JkqRRzBv0VXU78MScslural+3+iPgpEVomyRpDMbxZeyHge/02VbArUnuSrLpQB+SZFOSqSRTL/DcGJolSYIRv4xN8nfAPuDaPlXOrKrpJMcB25I80P0fwotU1WZgM8CROaZGaZck6XeGvqNPciFwHvA3VdUzmKtqunvfC9wErBv2eJKk4QwV9EnWA58C3llVz/apc3iSI/YvA+cAO3rVlSQtnkF+Xnkd8EPg1CS7k3wEuAw4gpnumO1Jruzqnphka7fr8cAdSe4BfgJ8u6q+uyhnIUnqa94++qra2KP4qj51p4EN3fLDQBsPXJSkCebIWEnL0kvxId6LxbluJKlxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcY6MlV5CfLbqS5N39JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxgzxKcEuSvUl2zCo7Jsm2JA9170f32Xd9kgeT7Epy6TgbLkkazCB39FcD6+eUXQrcVlVrgNu69d+TZAVwOXAucBqwMclpI7VWkrRg8wZ9Vd0OPDGn+Hzgmm75GuBdPXZdB+yqqoer6nng+m4/SdJBNGwf/fFV9RhA935cjzqrgEdnre/uynpKsinJVJKpF3huyGZJkuZazCkQ0qOs+lWuqs3AZoAjc0zfepJ+32I9RNuHc7dj2Dv6PUlOAOje9/aosxtYPWv9JGB6yONJkoY0bNDfDFzYLV8IfKtHnTuBNUlOSXIIcEG3nyTpIBrk55XXAT8ETk2yO8lHgM8DZyd5CDi7WyfJiUm2AlTVPuAi4BZgJ3BDVd2/OKchSepn3j76qtrYZ9Pbe9SdBjbMWt8KbB26dZKkkTkyVpIaZ9BLUuMMeklqnEEvSY0z6CWpcT4cXFqGlsNDvBdjFK0jaJeGd/SS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDXOoJekxhn0ktQ4R8ZKB9FyGPG6GBzxurx5Ry9JjRs66JOcmmT7rNdTSS6ZU+esJE/OqvOZ0ZssSVqIobtuqupB4AyAJCuAXwA39aj6g6o6b9jjSJJGM66um7cD/1NVPx/T50mSxmRcQX8BcF2fbW9Jck+S7yR5bb8PSLIpyVSSqRd4bkzNkiSNHPRJDgHeCfxXj813A6+pqtcD/wF8s9/nVNXmqlpbVWtXcuiozZIkdcZxR38ucHdV7Zm7oaqeqqpnuuWtwMokx47hmJKkAY0j6DfSp9smyauTpFte1x3vl2M4piRpQCMNmErycuBs4KOzyj4GUFVXAu8FPp5kH/Br4IKqqlGOKUlamJGCvqqeBV45p+zKWcuXAZeNcgxJ0micAkE6iJwqQEvBKRAkqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxjoyVJtxCHjjuyNyXJu/oJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuNGCvokjyS5L8n2JFM9tifJl5LsSnJvkjeOcjxJ0sKN43f0b62qx/tsOxdY073eBFzRvUuSDpLF7ro5H/hKzfgRcFSSExb5mJKkWUYN+gJuTXJXkk09tq8CHp21vrsre5Ekm5JMJZl6gedGbJYkab9Ru27OrKrpJMcB25I8UFW3z9qeHvtUrw+qqs3AZoAjc0zPOpJezGkNNJ+R7uirarp73wvcBKybU2U3sHrW+knA9CjHlCQtzNBBn+TwJEfsXwbOAXbMqXYz8IHu1zdvBp6sqseGbq0kacFG6bo5Hrgpyf7P+WpVfTfJxwCq6kpgK7AB2AU8C3xotOZKkhZq6KCvqoeBF3UOdgG/f7mATwx7DEnS6BwZK0mNM+glqXEGvSQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0b5Zmxq5N8P8nOJPcnubhHnbOSPJlke/f6zGjNlSQt1CjPjN0HfLKq7u4eEn5Xkm1V9dM59X5QVeeNcBxJ0giGvqOvqseq6u5u+WlgJ7BqXA2TJI3HWProk5wMvAH4cY/Nb0lyT5LvJHntAT5jU5KpJFMv8Nw4miVJYrSuGwCSvAK4Ebikqp6as/lu4DVV9UySDcA3gTW9PqeqNgObAY7MMTVquyRJM0a6o0+ykpmQv7aqvjF3e1U9VVXPdMtbgZVJjh3lmJKkhRnlVzcBrgJ2VtUX+9R5dVePJOu64/1y2GNKkhZulK6bM4H3A/cl2d6VfRr4I4CquhJ4L/DxJPuAXwMXVJXdMpJ0EA0d9FV1B5B56lwGXDbsMSRJo3NkrCQ1zqCXpMYZ9JLUOINekhpn0EtS4wx6SWqcQS9JjTPoJalxBr0kNc6gl6TGGfSS1DiDXpIaZ9BLUuMMeklqnEEvSY0z6CWpcQa9JDVu1IeDr0/yYJJdSS7tsT1JvtRtvzfJG0c5niRp4UZ5OPgK4HLgXOA0YGOS0+ZUOxdY0702AVcMezxJ0nBGuaNfB+yqqoer6nngeuD8OXXOB75SM34EHJXkhBGOKUlaoKEfDg6sAh6dtb4beNMAdVYBj839sCSbmLnrB3jue/X1HSO0bbk6Fnh8qRuxCDyvyeJ5TZZBz+s1/TaMEvTpUVZD1JkprNoMbAZIMlVVa0do27LkeU0Wz2uyeF79jdJ1sxtYPWv9JGB6iDqSpEU0StDfCaxJckqSQ4ALgJvn1LkZ+ED365s3A09W1Yu6bSRJi2forpuq2pfkIuAWYAWwparuT/KxbvuVwFZgA7ALeBb40IAfv3nYdi1zntdk8bwmi+fVR6p6dplLkhrhyFhJapxBL0mNW1ZBP9+UCpMqySNJ7kuyPcnUUrdnWEm2JNmbZMessmOSbEvyUPd+9FK2cRh9zutzSX7RXbPtSTYsZRuHkWR1ku8n2Znk/iQXd+UTfc0OcF4Tfc2SHJbkJ0nu6c7r77vyka/Xsumj76ZU+BlwNjM/y7wT2FhVP13Sho1BkkeAtVU10YM5kvwF8Awzo51P78r+GXiiqj7f/XE+uqo+tZTtXKg+5/U54Jmq+telbNsoulHoJ1TV3UmOAO4C3gV8kAm+Zgc4r/cxwdcsSYDDq+qZJCuBO4CLgfcw4vVaTnf0g0ypoCVUVbcDT8wpPh+4plu+hpl/cBOlz3lNvKp6rKru7pafBnYyMzJ9oq/ZAc5ronVTxTzTra7sXsUYrtdyCvp+0yW0oIBbk9zVTfXQkuP3j43o3o9b4vaM00XdrKtbJq17Y64kJwNvAH5MQ9dsznnBhF+zJCuSbAf2AtuqaizXazkF/cDTJUygM6vqjczM5vmJrqtAy9sVwJ8AZzAzN9MXlrY5w0vyCuBG4JKqemqp2zMuPc5r4q9ZVf2mqs5gZhaBdUlOH8fnLqegb3a6hKqa7t73Ajcx003Vij37ZyTt3vcucXvGoqr2dP/ofgt8mQm9Zl1f743AtVX1ja544q9Zr/Nq5ZoBVNWvgP8G1jOG67Wcgn6QKRUmTpLDuy+MSHI4cA7Q0sycNwMXdssXAt9awraMzZzptN/NBF6z7su9q4CdVfXFWZsm+pr1O69Jv2ZJXpXkqG75D4F3AA8whuu1bH51A9D9HOrf+d2UCv+0xE0aWZI/ZuYuHmamnPjqpJ5XkuuAs5iZNnUP8Fngm8ANwB8B/wv8VVVN1Bebfc7rLGa6AAp4BPjopM3TlOTPgB8A9wG/7Yo/zUx/9sReswOc10Ym+JoleR0zX7auYOYm/Iaq+ockr2TE67Wsgl6SNH7LqetGkrQIDHpJapxBL0mNM+glqXEGvSQ1zqCXpMYZ9JLUuP8HbrNWbRB1t54AAAAASUVORK5CYII=\n",
	"text/plain": [
	"<Figure size 432x288 with 1 Axes>"
	]
	},
	"metadata": {
	"needs_background": "light"
	},
	"output_type": "display_data"
	}
	],
	"source": [
	"plt.pcolormesh(watermask[::coarsen, ::coarsen]);"
	]
	},
	{
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"Image when run with coarsen=1"
	]
	},
	{
	"attachments": {
	"watermask.png": {
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD5CAYAAADcDXXiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAY3UlEQVR4nO3dXcwcV33H8e8Px7EhL0rcvNSxrdogU8lBjYMshyoVSgkQExAOF0FGgrpSKnORCNIigQ1SgYtIaVVeetFEMpDilhdjBdJYUURwDBE3aYwNTrDz0hjiJk/s2gWCCL0wifPvxZ6FyXp2Z3Z35tmd2d9HerSzZ2dmz+Hlt8dnzpxRRGBmZu3ymklXwMzMqudwNzNrIYe7mVkLOdzNzFrI4W5m1kIOdzOzFjqr7I6SFgD7gecj4j2SlgDfAlYCR4H3R8QLad9twE3AaeAjEfHAoHOfrUWxmHNGaoCZ2ax6kRd+EREX531WOtyBjwJPAOen91uBvRFxu6St6f0nJK0BNgGXA5cBD0p6Y0Sc7nfixZzDVbp2iKqYmdmDcfd/9/us1LCMpOXAu4EvZ4o3AjvS9g7ghkz5zog4FRHPAEeA9cNW2szMRld2zP2LwMeBVzJll0bEcYD0ekkqXwY8l9lvLpWZmdk8KQx3Se8BTkbEgZLnVE7ZGWscSNoiab+k/S9xquSpzcysjDJj7lcD75V0PbAYOF/S14ATkpZGxHFJS4GTaf85YEXm+OXAsd6TRsR2YDvA+VriBW7MzCpU2HOPiG0RsTwiVtK5UPr9iPggsBvYnHbbDNybtncDmyQtkrQKWA3sq7zmZmbW1zCzZXrdDuySdBPwLHAjQEQclrQLeBx4Gbh50EwZMzOrnqZhyd/ztSQ8FdLMbDgPxt0HImJd3me+Q9XMrIUc7mZmLeRwNzNrIYe7mVkLOdzNzFrI4W5m1kIOdzOzFnK4m5m1kMPdzKyFHO5mZi3kcDczayGHu5lZCznczcxayOFuZtZCDnczsxZyuJuZtZDD3cyshRzuZmYtVBjukhZL2ifpUUmHJX02lX9G0vOSDqa/6zPHbJN0RNJTkq6rswFmZnamMj33U8DbIuIKYC2wQdJb0mdfiIi16e9+AElrgE3A5cAG4A5JC2qou5nNsAeOPZq7bR1nFe0QnSdo/za9XZj+Bj1VeyOwMyJOAc9IOgKsBx4es65mZr8P8usuu6Iw1LP7zppSY+6SFkg6CJwE9kTEI+mjWyQ9JukuSRemsmXAc5nD51KZmdnYskGd3e4Gee9rd7v7NytKhXtEnI6ItcByYL2kNwF3Am+gM1RzHPhc2l15p+gtkLRF0n5J+1/i1EiVN7PZ1O2194Z1XrB39+/dp+2Gmi0TEb8GHgI2RMSJFPqvAF+iM/QCnZ76isxhy4FjOefaHhHrImLdQhaNVHkzm13XXXbFGcMtRb36WVJmtszFki5I268F3g48KWlpZrf3AYfS9m5gk6RFklYBq4F91VbbzKyjG+LZMfhsWd7+sxD2hRdUgaXAjjTj5TXAroi4T9K/S1pLZ8jlKPBhgIg4LGkX8DjwMnBzRJyupfZmZkPKDt20+UJrmdkyjwFX5pR/aMAxtwG3jVc1M7Px9OuhZ8fs83r8beA7VM2s8XpDOS+su2P0/fZtG4e7mTVCmRAet+fdpimT6tyjNFnna0lcpWsnXQ0zm0KDgjYvzPuNpef15nuHZcqcf5o8GHcfiIh1eZ+VuaBqZjYxg8bEs6FcFMRFxxYFfdM43M1s6vUL3by57GV723nnHObHYtp5zN3MGqNoiKaqQG56sIPD3cwaohu4VQ+b5J2vDUMzDncza4yqe9T9zueeu5nZBFTVs85baKwNwQ4OdzNrmHGGZ9oyh70Mz3M3s0Yap7c9HzdEzYdB89zdczezRhp3Sd8qZ9dMI89zN7PGyg7RDHpAR9ag1SDbFPbuuZtZ4/VbbqBfj37Y8iZyuJtZ65Tpgbc94B3uZjZTsmPtbQnyPA53M2uF3tAuCvB+67u3hcPdzFpjlKGWtvbeyzwge7GkfZIelXRY0mdT+RJJeyQ9nV4vzByzTdIRSU9Juq7OBpiZ9dNvCKatgZ5Vpud+CnhbRFwBrAU2SHoLsBXYGxGrgb3pPZLWAJuAy4ENwB3p4dpmZrXqN9Qy7Bh7G+5kLQz36Phtersw/QWwEdiRyncAN6TtjcDOiDgVEc8AR4D1ldbazGxIs3ARNavUmLukBZIOAieBPRHxCHBpRBwHSK+XpN2XAc9lDp9LZWZmU2EWgr7UHaoRcRpYK+kC4B5Jbxqwu/JOccZO0hZgC8BiXlemGmZmI+m9K7XsU52abKjZMhHxa+AhOmPpJyQtBUivJ9Nuc8CKzGHLgWM559oeEesiYt1CFo1QdTOzcoZ9NuqgZQ2aosxsmYtTjx1JrwXeDjwJ7AY2p902A/em7d3AJkmLJK0CVgP7qq64mVlZg+a6Z1/77d/EgC/Tc18K/EDSY8CP6Iy53wfcDrxD0tPAO9J7IuIwsAt4HPgucHMa1jEzq9Uwj8zr9ubL3sTUtIAvHHOPiMeAK3PKfwnkLsIeEbcBt41dOzOzHN2gzRsn7wZ27z5F678XPXy7aeHuO1TNrJHygj1vjnuZu1aLgr33+5rA67mbWaPkPfe07DF574vG45vK4W5mE1EUzqOEa3YMfdBQSt5nTRx6GcThbmYjK7ogOU5YjnpsmePyxt/bFOzgcDezAkWh15RQHHSTUlPaMAxfUDWzXE2+gadXv2Bv6sXSMtxzN7Pfa0uYd3XH0cu2a9i579PM4W42w9oW5r2KnsjU5uEZRZyxpte8O19L4irl3g9lZhVqU3gNo9+c+Oz7QcdMqwfj7gMRsS7vM/fczWbApEN90tMMe6dIZst7NSHUy3C4m7XYJAN1mkOyyb31sjxbxqxluhcQ5zPYy8w6mXRwFi0xMOn6Vc1j7mYtMF9BXlUATnqYqKvpge4xd7MWqzso6wjAQVMO62pP2560VMThbtZgVYfUfPZkq/6uoou2vT8mvfs3vRffy+Fu1lBVBfu0hFqZ9hTdaVrm2LxVJaflP4MqOdzNZtAkwmzQnZ/jrhBZ5q7SWRiKySrzDNUVkn4g6QlJhyV9NJV/RtLzkg6mv+szx2yTdETSU5Kuq7MBZrNm2JkwvaE3yV7qfM3i6dfGUX9cmqhMz/1l4GMR8WNJ5wEHJO1Jn30hIv4pu7OkNcAm4HLgMuBBSW/0c1TN5s80DjNkx7jHWb+l33K9bXtM3rjKPEP1OHA8bb8o6Qlg2YBDNgI7I+IU8IykI8B64OEK6ms200a9+WYaw76r7PrreQHdG/T9ntJU5ulNbVkwrGuom5gkraTzsOxHUtEtkh6TdJekC1PZMuC5zGFzDP4xMLMS+t0q3+RAGudReVm9s2BG0eT/HPOUDndJ5wLfBm6NiN8AdwJvANbS6dl/rrtrzuFn3CklaYuk/ZL2v8SpoStuNkvasAZKXm963DZ0f9y6wzLZO2Xzzj3sWHyTlQp3SQvpBPvXI+I7ABFxIiJOR8QrwJfoDL1Ap6e+InP4cuBY7zkjYntErIuIdQtZNE4bzFot71mfTQyjvAu7o46DFy3jmzcGP2vj7oXLD0gSsAP4VUTcmilfmsbjkfS3wFURsUnS5cA36IT9ZcBeYPWgC6pefsCs/aY1WJv4Q9k17vIDVwMfAn4q6WAq+yTwAUlr6Qy5HAU+DBARhyXtAh6nM9PmZs+UMZtNeb3nfp/1qqOn3e+ibJMDvh8vHGZm86LsTUzZAO7df9yw7x3OaXqoe+EwM5sKvQFfFNb9pkCOq+mhXobD3cxql9dr7xfYdTzXdBbCvJcf1mFmtSrzYGqrnnvuZlaLYde/qSPwZ7HH3uVwN7PKDTvHvK6efFtnwpThYRkzq0y/BbzmeximzDNd287hbmaVGHdN9iqVWSis7RzuZja2+XpYxrBLL7jnbmY2hmlct2WWgx0c7mZWkWFWYayzDk1dWK1qni1jZpWYzwup/W6ImuXZMb3cczezsQ1aIKzO78l+n3vsr+aeu5kNbVDIVvkQ7LywbsuiX3XzqpBmNpR+vfS6L6gWrU0zi2E/aFVIh7uZlTaJm5GKnp3au1zwLHG4m1llJjXlsRvkZdaE7+7fdl7P3cxGNskwL1M26LNZnj3jcDezV6nrARmj1CNr2JCedP0nrTDcJa0A/g34Y+AVYHtE/LOkJcC3gJV0nqH6/oh4IR2zDbgJOA18JCIeqKX2ZlaJOh6QMWmz2mPvKtNzfxn4WET8WNJ5wAFJe4C/BvZGxO2StgJbgU9IWgNsAi4HLgMelPRGPyTbbLo0JcBHCelZHo7pKgz3iDgOHE/bL0p6AlgGbASuSbvtAB4CPpHKd0bEKeAZSUeA9cDDVVfezIpNc4j3e57qOMHsYO8Y6g5VSSuBK4FHgEtT8Hd/AC5Juy0DnsscNpfKzGyedG8kmqZg73fBs8plCxzsf1D6gqqkc4FvA7dGxG8k9d01p+yM+ZaStgBbABbzurLVMLM+pinIu3qDtt8NT1XU3cH+aqV67pIW0gn2r0fEd1LxCUlL0+dLgZOpfA5YkTl8OXCs95wRsT0i1kXEuoUsGrX+ZjalJrEipP1BmdkyAr4CPBERn898tBvYDNyeXu/NlH9D0ufpXFBdDeyrstJmVl1Pvei2/nHPNx/caz9TmWGZq4EPAT+VdDCVfZJOqO+SdBPwLHAjQEQclrQLeJzOTJubPVPGrFrZxbPytov0u5CZ/Xw+HptXVI8iXkSsPy8/YNYgvTNK8sKtip530Xou4ypaVbKq72k7Lz9g1iKDlsEd93xFM2yqCtx+P0bDtMNDMYO5527WEFXPhhl2rH0+gnQWF/8ah3vuZvZ7VV9ArdJ8rQ0/CxzuZlOuyqAbJdgn2Xv20Mvo/AxVsxkx6jCMe9HN5J672ZQbdybMsD3faegpZ39YPHtmNL6gatYQwwT7uDNq8nr5Dtfp4wuqZg02bG992OGUvHnt/ebPz/ozS5vE4W42xcZZHTFPvzAuu6CXh0iaw8MyZlOozhky49TDPffpMmhYxrNlzKbItK3BnseB3gwOd7MpUUeoVxnEVS97YPVyuJtNgWkN9qJAdy9+ejnczSasjjVj6uqx1/EEJauHZ8uYTVBdD9yompfibR733M0mZNybkrrlDlvL43A3m7Ay4dzvaUnzqXdJAA/JTDcPy5hNwCiPxsua1DDJoAd7+F8Q06Ww5y7pLkknJR3KlH1G0vOSDqa/6zOfbZN0RNJTkq6rq+JmTVUUyr3lg2asTDJQe4eE3JOfLmWGZb4KbMgp/0JErE1/9wNIWgNsAi5Px9whaUFVlTVrujKzTbKLdPW+5u07TaE6TXWZdYXDMhHxQ0krS55vI7AzIk4Bz0g6AqwHHh65hmYtkDd8UeVj7aZh1cZJf7+92jhj7rdI+itgP/CxiHgBWAb8Z2afuVRmNvP6Dbf0670POibvuPkOeI+3T7dRZ8vcCbwBWAscBz6XypWzb+7KZJK2SNovaf9LnBqxGmbNUbRS46Cx9ezroPH2uodpemfKeCrm9Bqp5x4RJ7rbkr4E3JfezgErMrsuB471Ocd2YDt0VoUcpR5mTVI2BHtnwuT1yPv1mudjzNth3gwj9dwlLc28fR/QnUmzG9gkaZGkVcBqYN94VTRrtt7e96DPh3kCUplhnao52JujsOcu6ZvANcBFkuaATwPXSFpLZ8jlKPBhgIg4LGkX8DjwMnBzRJyup+pm029QyPZ72tEwsr18j4Fblh/WYVaT3rDtDfOicB/nBiebDX5Yh9kU6b0g2W/YZpTZLw526/LyA2YVKtPTzpve2K/HPijg8x6BZ9blcDcrYdwVHIc5Phvo3V5+3jnL3gxls8nhbpYxSs+799iywyllQ7noXN2eftE8epstDnebecMup9sN796hlN7edvY8vb3sMjNlhv3XQr9zevhmNnm2jM2kUYZZ+t0ZWnTH6KBz9tu3jgD2VMn2GTRbxuFuM613fHvQnaBQPCQz6POyPxK9qg7iaVgu2KrhcDer2TDhPsz5uhzElsfz3M1qVvXMFYe5jcvhbjYPRgnrogujZoM43M0qUkdv24+ys1E53M0aoOhxe2a9HO5mDVC0bLBZL4e7WYW8JIBNC4e7WQ1GXa530Pm6/MNhZXj5AbOaVDU+7jnvNgr33M3MWqgw3CXdJemkpEOZsiWS9kh6Or1emPlsm6Qjkp6SdF1dFTebBf0WJzMrUqbn/lVgQ0/ZVmBvRKwG9qb3SFoDbAIuT8fcIWlBZbU1m2G+kcmGURjuEfFD4Fc9xRuBHWl7B3BDpnxnRJyKiGeAI8D6iupqNnO6vXX32G1Yo465XxoRxwHS6yWpfBnwXGa/uVRmZhVxD97KqPqCqnLKcpedlLRF0n5J+1/iVMXVMGsf995tGKNOhTwhaWlEHJe0FDiZyueAFZn9lgPH8k4QEduB7dBZ8nfEepjNFAe8lTVqz303sDltbwbuzZRvkrRI0ipgNbBvvCqamdmwCnvukr4JXANcJGkO+DRwO7BL0k3As8CNABFxWNIu4HHgZeDmiDhdU93NzKyPwnCPiA/0+Sj30UkRcRtw2ziVMjOz8fgOVTOzFnK4m5m1kMPdzKyFHO5mZi3kcDczayGHu5lZCznczcxayOFuZtZCDnczsxZyuJuZtZDD3cyshRzuZmYt5HA3M2shh7uZWQs53M3MWsjhbmbWQg53M7MWcribmbVQ4WP2BpF0FHgROA28HBHrJC0BvgWsBI4C74+IF8arppmZDaOKnvtfRsTaiFiX3m8F9kbEamBvem9mZvOojmGZjcCOtL0DuKGG7zAzswHGDfcAvifpgKQtqezSiDgOkF4vGfM7zMxsSGONuQNXR8QxSZcAeyQ9WfbA9GOwBWAxrxuzGmZmljVWzz0ijqXXk8A9wHrghKSlAOn1ZJ9jt0fEuohYt5BF41TDzMx6jBzuks6RdF53G3gncAjYDWxOu20G7h23kmZmNpxxhmUuBe6R1D3PNyLiu5J+BOySdBPwLHDj+NU0M7NhjBzuEfFz4Iqc8l8C145TKTMzG4/vUDUzayGHu5lZCznczcxayOFuZtZCDnczsxZyuJuZtZDD3cyshRzuZmYt5HA3M2shh7uZWQs53M3MWsjhbmbWQg53M7MWcribmbWQw93MrIUc7mZmLeRwNzNrIYe7mVkL1RbukjZIekrSEUlb6/oeMzM7Uy3hLmkB8C/Au4A1wAckranju8zM7Ex19dzXA0ci4ucR8TtgJ7Cxpu8yM7MedYX7MuC5zPu5VGZmZvPgrJrOq5yyeNUO0hZgS3p76sG4+1BNdZmki4BfTLoSNXC7msXtapZh2vUn/T6oK9zngBWZ98uBY9kdImI7sB1A0v6IWFdTXSbG7WoWt6tZ3K7B6hqW+RGwWtIqSWcDm4DdNX2XmZn1qKXnHhEvS7oFeABYANwVEYfr+C4zMztTXcMyRMT9wP0ld99eVz0mzO1qFrerWdyuARQRxXuZmVmjePkBM7MWmni4N3mZAkl3STop6VCmbImkPZKeTq8XZj7bltr5lKTrJlPrYpJWSPqBpCckHZb00VTe2LZJWixpn6RHU5s+m8ob26YsSQsk/UTSfel9W9p1VNJPJR2UtD+VNb5tki6QdLekJ9P/z/688nZFxMT+6Fxs/RnweuBs4FFgzSTrNGT93wq8GTiUKftHYGva3gr8Q9pek9q3CFiV2r1g0m3o066lwJvT9nnAf6X6N7ZtdO69ODdtLwQeAd7S5Db1tO/vgG8A97Xlf4epvkeBi3rKGt82YAfwN2n7bOCCqts16Z57o5cpiIgfAr/qKd5I57840usNmfKdEXEqIp4BjtBp/9SJiOMR8eO0/SLwBJ07jBvbtuj4bXq7MP0FDW5Tl6TlwLuBL2eKG9+uARrdNknn0+kYfgUgIn4XEb+m4nZNOtzbuEzBpRFxHDohCVySyhvZVkkrgSvp9HQb3bY0dHEQOAnsiYjGtyn5IvBx4JVMWRvaBZ0f4O9JOpDuaofmt+31wP8C/5qG0r4s6Rwqbtekw71wmYIWaVxbJZ0LfBu4NSJ+M2jXnLKpa1tEnI6ItXTumF4v6U0Ddm9EmyS9BzgZEQfKHpJTNnXtyrg6It5MZ4XZmyW9dcC+TWnbWXSGc++MiCuB/6MzDNPPSO2adLgXLlPQQCckLQVIrydTeaPaKmkhnWD/ekR8JxW3om3pn8APARtofpuuBt4r6SidYc23SfoazW8XABFxLL2eBO6hMxzR9LbNAXPpX44Ad9MJ+0rbNelwb+MyBbuBzWl7M3BvpnyTpEWSVgGrgX0TqF8hSaIzHvhERHw+81Fj2ybpYkkXpO3XAm8HnqTBbQKIiG0RsTwiVtL5/8/3I+KDNLxdAJLOkXRedxt4J3CIhrctIv4HeE7Sn6aia4HHqbpdU3DV+Ho6szF+Bnxq0vUZsu7fBI4DL9H5db0J+CNgL/B0el2S2f9TqZ1PAe+adP0HtOsv6Pyz7zHgYPq7vsltA/4M+Elq0yHg71N5Y9uU08Zr+MNsmca3i87Y9KPp73A3H1rStrXA/vS/x/8ALqy6Xb5D1cyshSY9LGNmZjVwuJuZtZDD3cyshRzuZmYt5HA3M2shh7uZWQs53M3MWsjhbmbWQv8PQlOHMOK5OscAAAAASUVORK5CYII="
	}
	},
	"cell_type": "markdown",
	"metadata": {},
	"source": [
	"![watermask.png](attachment:watermask.png)"
	]
	},
	{
	"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
	}