{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define a subset domain area"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To launch this notebook interactively in a Jupyter notebook-like browser interface, please click the \"Launch Binder\" button below. Note that Binder may take several minutes to launch.\n",
    "\n",
    "[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/hydroframe/subsettools-binder/HEAD?labpath=subsettools%2Fdefinte_subset.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of the subset functions require you to define a bounding box that defines the extent of your domain.  We use `i`,`j` indices to define this box where `i`  and `j` are index values relative to the lower left hand corner of whatever the reference grid is that will be subset. \n",
    "\n",
    "The [`define_latlon_domain`](https://hydroframesubsettools.readthedocs.io/en/latest/autoapi/subsettools/subsettools/index.html#subsettools.subsettools.define_latlon_domain) function translates a bounding box in lat-lon  coordinates bounds to `i`,`j` indices in whatever grid system we select. It returns a tuple `(imin, jmin, imax, jmax)` of grid indices that define a bounding box containing our region (or point) of interest (Note: `(imin, jmin, imax, jmax)` are the west, south, east and north boundaries of the box respectively) and a mask for that domain.\n",
    "\n",
    "Here we will show how to define a subset extent for:\n",
    "-  A single point of interest \n",
    "-  A user specified bounding box\n",
    "-  A bounding box that surrounds a user specified HUC watershed\n",
    "-  A bounding box that encompasses the upstream area of a collection of points\n",
    "\n",
    "**IMPORTANT NOTE**: *The i,j indices found in this step are based on whatever grid you select (e.g. `conus1` or `conus2`). It's very important that the grid you use in this step is the same as the grid that the data files (static input and forcing) you are subsetting are in or you will end up subsetting a different location than you expect.*\n",
    "\n",
    "The grids are shown below and described in [Yang et al 2023](https://www.sciencedirect.com/science/article/pii/S0022169423012362)* \n",
    "\n",
    "![CONUS domains](CONUS1_2_domain.jpg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.  Setup \n",
    "\n",
    "In all examples you will need to import the following packages and register your PIN in order to have access to the HydroData datasets\n",
    "\n",
    "Refer to the [getting started](https://hydroframesubsettools.readthedocs.io/en/latest/getting_started.html) instructions for creating your pin if you have not done this already."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import subsettools as st\n",
    "import hf_hydrodata as hf\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "hf.register_api_pin(\"your_email\", \"your_pin\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Defining bounds to extract data for a single point\n",
    "To extract data for a single point we use the same bounding box function as we would to extract a larger domain but just repeat the point values as the upper and lower bounds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bounding box: (4057, 1915, 4058, 1916)\n",
      "mask: [[1]]\n"
     ]
    }
   ],
   "source": [
    "lat = 39.8379\n",
    "lon = -74.3791\n",
    "# Since we want to subset only a single location, both lat-lon bounds are defined by this point:\n",
    "latlon_bounds = [[lat, lon],[lat, lon]]\n",
    "ij_column_bounds, mask = st.define_latlon_domain(latlon_bounds=latlon_bounds, grid=\"conus2\")\n",
    "print(f\"bounding box: {ij_column_bounds}\")\n",
    "# The mask contains a single point:\n",
    "print(f\"mask: {mask}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Defining bounds for a box defined by lat-lon bounds\n",
    "To extract a bounding box, provide the upper and lower latitude and longitude bounds respectively for the area of interest as well as the grid system that you would like to use. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bounding box: (2285, 436, 2359, 496)\n"
     ]
    }
   ],
   "source": [
    "ij_box_bounds, mask = st.define_latlon_domain(latlon_bounds=[[37.91, -91.43], [37.34, -90.63]], grid=\"conus1\")\n",
    "print(f\"bounding box: {ij_box_bounds}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Defining bounds for a HUC watershed\n",
    "The subsettools [`define_huc_domain`](https://hydroframesubsettools.readthedocs.io/en/latest/autoapi/subsettools/subsettools/index.html#subsettools.subsettools.define_huc_domain) function returns a tuple `(imin, jmin, imax, jmax)` of grid indices that define a bounding box containing any HUC, and a mask for that domain. You can provide 2, 4, 6, 8 or 10-digit HUCs.  For help finding your HUC you can refer to the [USGS HUC picker](https://water.usgs.gov/wsc/a_api/wbd/index_wbd.html). **Make sure to use the new (up to 12-digit) USGS HUC maps rather than the legacy 8-digit ones.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bounding box: (1225, 1738, 1347, 1811)\n"
     ]
    }
   ],
   "source": [
    "ij_huc_bounds, mask = st.define_huc_domain(hucs=[\"14050002\"], grid=\"conus2\")\n",
    "print(f\"bounding box: {ij_huc_bounds}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now visualize the mask for our HUC within its bounding box. Note that, because the origin of the CONUS grids it the lower-left corner, we need to use `origin=lower` to orient the mask properly. (The origin for numpy arrays is, of course, the upper-left corner.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x152a41138c90>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiYAAAFUCAYAAAAQzPcRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAjBUlEQVR4nO3df3BU1f3/8deGJAsKuyERdklJNLbYoAJiwLBCK8XYlK9jocRfDFYUpow0UCHTqmmL9ocapFOxKGB1HGynpbRMBUQ/6mCUONjwK0gropFWxlDDLmqb3YjmB+R8/2jdugFM9keyJ5vnY+bOsPfevXnnhCSvOfedcx3GGCMAAAALpCW7AAAAgE8RTAAAgDUIJgAAwBoEEwAAYA2CCQAAsAbBBAAAWINgAgAArJGe7AI66+joUGNjo4YMGSKHw5HscgAAQDcYY9Tc3Kzc3FylpcU+72FdMGlsbFReXl6yywAAADE4cuSIRo4cGfP7rQsmQ4YMkSRN0f9TujKSXA0AAKe36e3XI15/64IxSarEDifUrh36v/Dv8VhZF0w+vX2TrgylOwgmAAA7uYZE3q7o97+z/vuAm3jbMGh+BQAA1rBuxgQAgBca9/f4xyjNvSShNXQ+v6vr4/SYMQEAANYgmAAAAGtEFUzOO+88ORyOU7by8nJJUktLi8rLy5WTk6PBgwerrKxMgUCgRwoHAACpx2GMMd09+f3339fJkyfDrw8cOKCrrrpKL7/8sqZOnaqFCxfq2Wef1ZNPPim3261FixYpLS1Nr776arcLCoVCcrvdmqoZdDgDQIrqjR4S26R6z8kJ067t2qJgMCiXyxXzdaJqfh02bFjE6+XLl+uLX/yirrjiCgWDQT3xxBNav369pk2bJklat26dRo8erZ07d2rSpEkxFwkAAPqHmHtM2tra9Lvf/U7z5s2Tw+FQXV2d2tvbVVJSEj6nsLBQ+fn5qq2tPeN1WltbFQqFIjYAANA/xRxMNm/erKamJt1yyy2SJL/fr8zMTGVlZUWc5/F45Pf7z3idqqoqud3u8MZy9AAA9F8xr2PyxBNPaPr06crNzY2rgMrKSlVUVIRfh0IhwgkAWKY/9oQkGuucdE9MweTdd9/Viy++qKeeeiq8z+v1qq2tTU1NTRGzJoFAQF6v94zXcjqdcjqdsZQBAABSTEy3ctatW6fhw4fr6quvDu8rKipSRkaGqqurw/vq6+vV0NAgn88Xf6UAACDlRT1j0tHRoXXr1mnu3LlKT//f291ut+bPn6+KigplZ2fL5XJp8eLF8vl8/EUOAADolqiDyYsvvqiGhgbNmzfvlGMrV65UWlqaysrK1NraqtLSUq1ZsyYhhQIAeg49JL2PnpPTi2qBtd7AAmsA0PsIJsnX14NJohZY41k5AADAGgQTAABgjZjXMUFica8R6D+4bYLT4ffAfzBjAgAArEEwAQAA1iCYAAAAa9BjkiRd3WOO9h50f70XCfQF9JQgFl39v0nVn/vMmAAAAGsQTAAAgDUIJgAAwBr0mCRJ53uD8d6DTsQ97FS9XwkkGj0jQM9hxgQAAFiDYAIAAKxBMAEAANYgmAAAAGvQ/IowHiAFnB7NrrBRqv7MZsYEAABYg2ACAACsQTABAADWoMcEZ5Sq9y+BrtBTAiQPMyYAAMAaBBMAAGANggkAALAGPSZJ0hfvYUdbMz0psFVf/P4DuhLL/2sbf04zYwIAAKxBMAEAANYgmAAAAGvQY5Ikne/rpeI9b9ZB6XmM8alS8XsJ6E+YMQEAANaIOpi89957uummm5STk6NBgwZpzJgx2rt3b/i4MUZ33323RowYoUGDBqmkpESHDh1KaNEAACA1RRVM/v3vf2vy5MnKyMjQc889p4MHD+qXv/ylhg4dGj5nxYoVWrVqlR599FHt2rVLZ599tkpLS9XS0pLw4gEAQGqJqsfkgQceUF5entatWxfeV1BQEP63MUYPPfSQfvzjH2vGjBmSpN/+9rfyeDzavHmzbrzxxgSVjb6IfggAQFeimjF5+umnNWHCBF133XUaPny4xo8fr8cffzx8/PDhw/L7/SopKQnvc7vdKi4uVm1t7Wmv2draqlAoFLEBAID+Kapg8s4772jt2rUaNWqUXnjhBS1cuFDf+9739Jvf/EaS5Pf7JUkejyfifR6PJ3yss6qqKrnd7vCWl5cXy+cBAABSQFTBpKOjQ5deeqnuv/9+jR8/XgsWLNB3vvMdPfroozEXUFlZqWAwGN6OHDkS87UAAEDfFlWPyYgRI3ThhRdG7Bs9erT+/Oc/S5K8Xq8kKRAIaMSIEeFzAoGALrnkktNe0+l0yul0RlMGUgTP3gGA5LKx9y+qGZPJkyervr4+Yt/bb7+tc889V9J/GmG9Xq+qq6vDx0OhkHbt2iWfz5eAcgEAQCqLasZk6dKluvzyy3X//ffr+uuv1+7du/XYY4/psccekyQ5HA4tWbJE9957r0aNGqWCggItW7ZMubm5mjlzZk/UDwAAUkhUwWTixInatGmTKisr9bOf/UwFBQV66KGHNGfOnPA5d9xxh44fP64FCxaoqalJU6ZM0fPPP6+BAwcmvHgAAJBaHMYYk+wiPisUCsntdmuqZijdkZHscnoNz/eInw33RnubjfeHk43vJSB28fwMOWHatV1bFAwG5XK5Yr4Oz8oBAADWIJgAAABrEEwAAIA1omp+Rfdxn7v30W/RP/G9BqQWZkwAAIA1CCYAAMAaBBMAAGANekwShPvc9rGh56S3/1909fH6Yt8N31tA/8KMCQAAsAbBBAAAWINgAgAArEGPCfqN7vScRNuXQv/DqeLtc2FMgeSxoTePGRMAAGANggkAALAGwQQAAFiDHpNu4r536unO15Sve9eiHSPGFLCHjWsbMWMCAACsQTABAADWIJgAAABr0GNyBtwHB06P7w2g77Kxp6QzZkwAAIA1CCYAAMAaBBMAAGANekyAfoT+EAC2Y8YEAABYg2ACAACsQTABAADWIJgAAABrEEwAAIA1CCYAAMAaUQWTn/zkJ3I4HBFbYWFh+HhLS4vKy8uVk5OjwYMHq6ysTIFAIOFFAwCA1BT1jMlFF12ko0ePhrcdO3aEjy1dulRbt27Vxo0bVVNTo8bGRs2aNSuhBQMAgNQV9QJr6enp8nq9p+wPBoN64okntH79ek2bNk2StG7dOo0ePVo7d+7UpEmT4q8WAADErPMiizY+1C/qGZNDhw4pNzdX559/vubMmaOGhgZJUl1dndrb21VSUhI+t7CwUPn5+aqtrT3j9VpbWxUKhSI2AADQP0UVTIqLi/Xkk0/q+eef19q1a3X48GF95StfUXNzs/x+vzIzM5WVlRXxHo/HI7/ff8ZrVlVVye12h7e8vLyYPhEAAND3RXUrZ/r06eF/jx07VsXFxTr33HP1pz/9SYMGDYqpgMrKSlVUVIRfh0IhwgkAAP1UXA/xy8rK0gUXXKC///3vuuqqq9TW1qampqaIWZNAIHDanpRPOZ1OOZ3OeMqIGw82AwD0Rzb2nMS1jslHH32kf/zjHxoxYoSKioqUkZGh6urq8PH6+no1NDTI5/PFXSgAAEh9Uc2YfP/739c111yjc889V42Njbrnnns0YMAAzZ49W263W/Pnz1dFRYWys7Plcrm0ePFi+Xw+/iIHAAB0S1TB5J///Kdmz56tDz/8UMOGDdOUKVO0c+dODRs2TJK0cuVKpaWlqaysTK2trSotLdWaNWt6pHAAAJB6HMYYk+wiPisUCsntdmuqZijdkSGJHhAAAJIhmp6TE6Zd27VFwWBQLpcr5o/Js3IAAIA1CCYAAMAaBBMAAGCNuNYx6Umb3n5driHkJgAAkiUZ65zwmx8AAFiDYAIAAKxBMAEAANawtsekv+l83y7etVtseN5BZ6xHAwDoCjMmAADAGgQTAABgDYIJAACwhrU9Jt+6YEy/flaOjT0iAAD0NGZMAACANQgmAADAGgQTAABgDWt7TD6rq36L/tiD0hfxdQSAvq03np3DjAkAALAGwQQAAFiDYAIAAKzRJ3pMupLo58wgOaL9OvJ1B4DUw4wJAACwBsEEAABYg2ACAACskRI9Jp3Re5AaeF4QANjts79fQ80dGnpB/NdkxgQAAFiDYAIAAKxBMAEAANZIyR4TAACQeJ/X+3fCtEt6J+6PwYwJAACwBsEEAABYI65gsnz5cjkcDi1ZsiS8r6WlReXl5crJydHgwYNVVlamQCAQb50AAKAfiLnHZM+ePfr1r3+tsWPHRuxfunSpnn32WW3cuFFut1uLFi3SrFmz9Oqrr8ZdbCqL9rkwAACkophmTD766CPNmTNHjz/+uIYOHRreHwwG9cQTT+jBBx/UtGnTVFRUpHXr1ukvf/mLdu7cmbCiAQBAaoopmJSXl+vqq69WSUlJxP66ujq1t7dH7C8sLFR+fr5qa2tPe63W1laFQqGIDQAA9E9R38rZsGGD9u3bpz179pxyzO/3KzMzU1lZWRH7PR6P/H7/aa9XVVWln/70p9GWAQAAUlBUweTIkSO6/fbbtW3bNg0cODAhBVRWVqqioiL8OhQKKS8vLyHX/hTPzukfeuLrzP8dAOhdUd3Kqaur07Fjx3TppZcqPT1d6enpqqmp0apVq5Seni6Px6O2tjY1NTVFvC8QCMjr9Z72mk6nUy6XK2IDAAD9U1QzJldeeaVef/31iH233nqrCgsLdeeddyovL08ZGRmqrq5WWVmZJKm+vl4NDQ3y+XyJqxoAAKSkqILJkCFDdPHFF0fsO/vss5WTkxPeP3/+fFVUVCg7O1sul0uLFy+Wz+fTpEmTElc1AABISQl/Vs7KlSuVlpamsrIytba2qrS0VGvWrEn0h4lLb/QNdLXuCL0KPa8n1n6h5wRAf9b5Z15P/JyNO5hs37494vXAgQO1evVqrV69Ot5LAwCAfoZn5QAAAGsQTAAAgDUS3mPSF8Vyjyze+2xd9SrwbJy+IdpeInpUAODzMWMCAACsQTABAADWIJgAAABr0GMSo0T3gNBTkpq6+romer2bnuhd6un+J/psAHwWMyYAAMAaBBMAAGANggkAALAGwQQAAFiD5lfAYvE2miaiUbWnG7NZdA7AZzFjAgAArEEwAQAA1iCYAAAAa9BjAsAq9JwA/RszJgAAwBoEEwAAYA2CCQAAsAbBBAAAWINgAgAArEEwAQAA1iCYAAAAa7COCQCrsa4J0Ht6+tlY3cGMCQAAsAbBBAAAWINgAgAArEEwAQAA1iCYAAAAa0QVTNauXauxY8fK5XLJ5XLJ5/PpueeeCx9vaWlReXm5cnJyNHjwYJWVlSkQCCS8aAAAkJqiCiYjR47U8uXLVVdXp71792ratGmaMWOG3njjDUnS0qVLtXXrVm3cuFE1NTVqbGzUrFmzeqRwAACQehzGGBPPBbKzs/WLX/xC1157rYYNG6b169fr2muvlSS99dZbGj16tGprazVp0qRuXS8UCsntdmuqZijdkRFPaQBSEOuYAImTyHVLTph2bdcWBYNBuVyumK8Tc4/JyZMntWHDBh0/flw+n091dXVqb29XSUlJ+JzCwkLl5+ertrb2jNdpbW1VKBSK2AAAQP8UdTB5/fXXNXjwYDmdTt12223atGmTLrzwQvn9fmVmZiorKyvifI/HI7/ff8brVVVVye12h7e8vLyoPwkAAJAaog4mX/7yl7V//37t2rVLCxcu1Ny5c3Xw4MGYC6isrFQwGAxvR44ciflaAACgb4v6WTmZmZn60pe+JEkqKirSnj179Ktf/Uo33HCD2tra1NTUFDFrEggE5PV6z3g9p9Mpp9MZfeUAACAqNjwLpytxr2PS0dGh1tZWFRUVKSMjQ9XV1eFj9fX1amhokM/ni/fDAACAfiCqGZPKykpNnz5d+fn5am5u1vr167V9+3a98MILcrvdmj9/vioqKpSdnS2Xy6XFixfL5/N1+y9yAABA/xZVMDl27JhuvvlmHT16VG63W2PHjtULL7ygq666SpK0cuVKpaWlqaysTK2trSotLdWaNWt6pHAAAJB64l7HJNFYxwTA52EdE6D7erOnJOnrmAAAACQawQQAAFiDYAIAAKwR9TomAJBM0d4zj6UnJd778vTBIBb8v/sPZkwAAIA1CCYAAMAaBBMAAGANekwApLRkPBuk88dMlXv/QG9gxgQAAFiDYAIAAKxBMAEAANagxwQAehg9J+gNyein6gnMmAAAAGsQTAAAgDUIJgAAwBoEEwAAYA2CCQAAsAbBBAAAWINgAgAArEEwAQAA1iCYAAAAaxBMAACANQgmAADAGjwrBwB6Gc/OgZQ6z7ZJNGZMAACANQgmAADAGgQTAABgDYIJAACwBsEEAABYg2ACAACsEVUwqaqq0sSJEzVkyBANHz5cM2fOVH19fcQ5LS0tKi8vV05OjgYPHqyysjIFAoGEFg0AAFJTVMGkpqZG5eXl2rlzp7Zt26b29nZ9/etf1/Hjx8PnLF26VFu3btXGjRtVU1OjxsZGzZo1K+GFAwCA1OMwxphY3/z+++9r+PDhqqmp0Ve/+lUFg0ENGzZM69ev17XXXitJeuuttzR69GjV1tZq0qRJXV4zFArJ7XZrqmYo3ZERa2kA0Gex4Fr/kGoLrJ0w7dquLQoGg3K5XDFfJ64ek2AwKEnKzs6WJNXV1am9vV0lJSXhcwoLC5Wfn6/a2trTXqO1tVWhUChiAwAA/VPMwaSjo0NLlizR5MmTdfHFF0uS/H6/MjMzlZWVFXGux+OR3+8/7XWqqqrkdrvDW15eXqwlAQCAPi7mYFJeXq4DBw5ow4YNcRVQWVmpYDAY3o4cORLX9QAAQN8V00P8Fi1apGeeeUavvPKKRo4cGd7v9XrV1tampqamiFmTQCAgr9d72ms5nU45nc5YygCAlMRD/vqHzl/XVOs5iVVUMybGGC1atEibNm3SSy+9pIKCgojjRUVFysjIUHV1dXhffX29Ghoa5PP5ElMxAABIWVHNmJSXl2v9+vXasmWLhgwZEu4bcbvdGjRokNxut+bPn6+KigplZ2fL5XJp8eLF8vl83fqLHAAA0L9FFUzWrl0rSZo6dWrE/nXr1umWW26RJK1cuVJpaWkqKytTa2urSktLtWbNmoQUCwAAUltc65j0BNYxAYBI9Jj0T32t58SKdUwAAAASiWACAACsQTABAADWiGkdEwBA72Fdk/6pv65zwowJAACwBsEEAABYg2ACAACsQY8JAPQxXfUa0IOSmvpLzwkzJgAAwBoEEwAAYA2CCQAAsAY9JgCQYrrTe0AfCmzFjAkAALAGwQQAAFiDYAIAAKxBjwkAAH1Qqq5rwowJAACwBsEEAABYg2ACAACsQY8JAPRDnfsRWNcEtmDGBAAAWINgAgAArEEwAQAA1qDHBADQ5RoY9KCgtzBjAgAArEEwAQAA1iCYAAAAaxBMAACANQgmAADAGgQTAABgjaiDySuvvKJrrrlGubm5cjgc2rx5c8RxY4zuvvtujRgxQoMGDVJJSYkOHTqUqHoBAEAKi3odk+PHj2vcuHGaN2+eZs2adcrxFStWaNWqVfrNb36jgoICLVu2TKWlpTp48KAGDhyYkKIBAL2LZ+vYp6u1Z/qqqIPJ9OnTNX369NMeM8booYce0o9//GPNmDFDkvTb3/5WHo9Hmzdv1o033hhftQAAIKUltMfk8OHD8vv9KikpCe9zu90qLi5WbW3tad/T2tqqUCgUsQEAgP4pocHE7/dLkjweT8R+j8cTPtZZVVWV3G53eMvLy0tkSQAAoA9J+rNyKisrVVFREX4dCoUIJwBgOXpOkq/zmKdKz0lCZ0y8Xq8kKRAIROwPBALhY505nU65XK6IDQAA9E8JDSYFBQXyer2qrq4O7wuFQtq1a5d8Pl8iPxQAAEhBUd/K+eijj/T3v/89/Prw4cPav3+/srOzlZ+fryVLlujee+/VqFGjwn8unJubq5kzZyaybgAAkIKiDiZ79+7V1772tfDrT/tD5s6dqyeffFJ33HGHjh8/rgULFqipqUlTpkzR888/zxomAJDC6DlBojiMMSbZRXxWKBSS2+3WVM1QuiMj2eUAAGJAMOl9yW5+PWHatV1bFAwG4+oX5Vk5AADAGgQTAABgjaSvYwIASD29fVuBW0enH4Nk396JBTMmAADAGgQTAABgDYIJAACwBsEEAABYg+ZXAECfxwJvqYMZEwAAYA2CCQAAsAbBBAAAWIMeEwBAyollYTH6UuzAjAkAALAGwQQAAFiDYAIAAKxBjwkAAOq6L6Uv9qB0rrkvPNSPGRMAAGANggkAALAGwQQAAFiDHhMAALohFXpQuqrRhh4UZkwAAIA1CCYAAMAaBBMAAGANekwAAEiAzv0ZfaHnxEbMmAAAAGsQTAAAgDUIJgAAwBr0mAAA0AP6Qs+JDeuWdMaMCQAAsEaPBZPVq1frvPPO08CBA1VcXKzdu3f31IcCAAApokeCyR//+EdVVFTonnvu0b59+zRu3DiVlpbq2LFjPfHhAABAinAYY0yiL1pcXKyJEyfqkUcekSR1dHQoLy9Pixcv1l133fW57w2FQnK73ZqqGUp3ZCS6NAAA+oTOPSk29oN81gnTru3aomAwKJfLFfN1Ej5j0tbWprq6OpWUlPzvg6SlqaSkRLW1taec39raqlAoFLEBAID+KeHB5IMPPtDJkyfl8Xgi9ns8Hvn9/lPOr6qqktvtDm95eXmJLgkAAPQRSf9z4crKSlVUVIRfB4NB5efn64TapYTfZAIAoG8INXdEvD5h2pNUSfec0H/qi7dDJOHB5JxzztGAAQMUCAQi9gcCAXm93lPOdzqdcjqd4def3srZof9LdGkAAPQZQy/ovOedZJQRtebmZrnd7pjfn/BgkpmZqaKiIlVXV2vmzJmS/tP8Wl1drUWLFnX5/tzcXB05ckTGGOXn5+vIkSNxNdH0d6FQSHl5eYxjHBjD+DGGicE4xo8xjN+ZxtAYo+bmZuXm5sZ1/R65lVNRUaG5c+dqwoQJuuyyy/TQQw/p+PHjuvXWW7t8b1pamkaOHBmeOXG5XPznSQDGMX6MYfwYw8RgHOPHGMbvdGMYz0zJp3okmNxwww16//33dffdd8vv9+uSSy7R888/f0pDLAAAwGf1WPProkWLunXrBgAA4FPWPivH6XTqnnvuiWiMRfQYx/gxhvFjDBODcYwfYxi/nh7DHln5FQAAIBbWzpgAAID+h2ACAACsQTABAADWIJgAAABrEEwAAIA1rA0mq1ev1nnnnaeBAwequLhYu3fvTnZJ1qqqqtLEiRM1ZMgQDR8+XDNnzlR9fX3EOS0tLSovL1dOTo4GDx6ssrKyU55nhP9Zvny5HA6HlixZEt7HGHbPe++9p5tuukk5OTkaNGiQxowZo71794aPG2N09913a8SIERo0aJBKSkp06NChJFZsl5MnT2rZsmUqKCjQoEGD9MUvflE///nPIx6MxhhGeuWVV3TNNdcoNzdXDodDmzdvjjjenfH617/+pTlz5sjlcikrK0vz58/XRx991IufRfJ93ji2t7frzjvv1JgxY3T22WcrNzdXN998sxobGyOukYhxtDKY/PGPf1RFRYXuuece7du3T+PGjVNpaamOHTuW7NKsVFNTo/Lycu3cuVPbtm1Te3u7vv71r+v48ePhc5YuXaqtW7dq48aNqqmpUWNjo2bNmpXEqu21Z88e/frXv9bYsWMj9jOGXfv3v/+tyZMnKyMjQ88995wOHjyoX/7ylxo6dGj4nBUrVmjVqlV69NFHtWvXLp199tkqLS1VS0tLEiu3xwMPPKC1a9fqkUce0ZtvvqkHHnhAK1as0MMPPxw+hzGMdPz4cY0bN06rV68+7fHujNecOXP0xhtvaNu2bXrmmWf0yiuvaMGCBb31KVjh88bx448/1r59+7Rs2TLt27dPTz31lOrr6/XNb34z4ryEjKOx0GWXXWbKy8vDr0+ePGlyc3NNVVVVEqvqO44dO2YkmZqaGmOMMU1NTSYjI8Ns3LgxfM6bb75pJJna2tpklWml5uZmM2rUKLNt2zZzxRVXmNtvv90Ywxh215133mmmTJlyxuMdHR3G6/WaX/ziF+F9TU1Nxul0mj/84Q+9UaL1rr76ajNv3ryIfbNmzTJz5swxxjCGXZFkNm3aFH7dnfE6ePCgkWT27NkTPue5554zDofDvPfee71Wu006j+Pp7N6920gy7777rjEmceNo3YxJW1ub6urqVFJSEt6XlpamkpIS1dbWJrGyviMYDEqSsrOzJUl1dXVqb2+PGNPCwkLl5+czpp2Ul5fr6quvjhgriTHsrqeffloTJkzQddddp+HDh2v8+PF6/PHHw8cPHz4sv98fMY5ut1vFxcWM439dfvnlqq6u1ttvvy1J+utf/6odO3Zo+vTpkhjDaHVnvGpra5WVlaUJEyaEzykpKVFaWpp27drV6zX3FcFgUA6HQ1lZWZISN4499qycWH3wwQc6efLkKQ/883g8euutt5JUVd/R0dGhJUuWaPLkybr44oslSX6/X5mZmeH/PJ/yeDzy+/1JqNJOGzZs0L59+7Rnz55TjjGG3fPOO+9o7dq1qqio0A9/+EPt2bNH3/ve95SZmam5c+eGx+p039+M43/cddddCoVCKiws1IABA3Ty5Endd999mjNnjiQxhlHqznj5/X4NHz484nh6erqys7MZ0zNoaWnRnXfeqdmzZ4efMJyocbQumCA+5eXlOnDggHbs2JHsUvqUI0eO6Pbbb9e2bds0cODAZJfTZ3V0dGjChAm6//77JUnjx4/XgQMH9Oijj2ru3LlJrq5v+NOf/qTf//73Wr9+vS666CLt379fS5YsUW5uLmMIK7S3t+v666+XMUZr165N+PWtu5VzzjnnaMCAAaf8tUMgEJDX601SVX3DokWL9Mwzz+jll1/WyJEjw/u9Xq/a2trU1NQUcT5j+j91dXU6duyYLr30UqWnpys9PV01NTVatWqV0tPT5fF4GMNuGDFihC688MKIfaNHj1ZDQ4MkhceK7+8z+8EPfqC77rpLN954o8aMGaNvf/vbWrp0qaqqqiQxhtHqznh5vd5T/rjixIkT+te//sWYdvJpKHn33Xe1bdu28GyJlLhxtC6YZGZmqqioSNXV1eF9HR0dqq6uls/nS2Jl9jLGaNGiRdq0aZNeeuklFRQURBwvKipSRkZGxJjW19eroaGBMf2vK6+8Uq+//rr2798f3iZMmKA5c+aE/80Ydm3y5Mmn/Kn622+/rXPPPVeSVFBQIK/XGzGOoVBIu3btYhz/6+OPP1ZaWuSP5gEDBqijo0MSYxit7oyXz+dTU1OT6urqwue89NJL6ujoUHFxca/XbKtPQ8mhQ4f04osvKicnJ+J4wsYxhmbdHrdhwwbjdDrNk08+aQ4ePGgWLFhgsrKyjN/vT3ZpVlq4cKFxu91m+/bt5ujRo+Ht448/Dp9z2223mfz8fPPSSy+ZvXv3Gp/PZ3w+XxKrtt9n/yrHGMawO3bv3m3S09PNfffdZw4dOmR+//vfm7POOsv87ne/C5+zfPlyk5WVZbZs2WL+9re/mRkzZpiCggLzySefJLFye8ydO9d84QtfMM8884w5fPiweeqpp8w555xj7rjjjvA5jGGk5uZm89prr5nXXnvNSDIPPvigee2118J/LdKd8frGN75hxo8fb3bt2mV27NhhRo0aZWbPnp2sTykpPm8c29razDe/+U0zcuRIs3///ojfNa2treFrJGIcrQwmxhjz8MMPm/z8fJOZmWkuu+wys3PnzmSXZC1Jp93WrVsXPueTTz4x3/3ud83QoUPNWWedZb71rW+Zo0ePJq/oPqBzMGEMu2fr1q3m4osvNk6n0xQWFprHHnss4nhHR4dZtmyZ8Xg8xul0miuvvNLU19cnqVr7hEIhc/vtt5v8/HwzcOBAc/7555sf/ehHET/8GcNIL7/88ml/Bs6dO9cY073x+vDDD83s2bPN4MGDjcvlMrfeeqtpbm5OwmeTPJ83jocPHz7j75qXX345fI1EjKPDmM8sJwgAAJBE1vWYAACA/otgAgAArEEwAQAA1iCYAAAAaxBMAACANQgmAADAGgQTAABgDYIJAACwBsEEAABYg2ACAACsQTABAADW+P8GYHap8vghYgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(mask, origin='lower')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Defining bounds for the upstream area of a collection of points\n",
    "The subsettools [`define_upstream_domain`](https://hydroframesubsettools.readthedocs.io/en/latest/autoapi/subsettools/subsettools/index.html#subsettools.subsettools.define_upstream_domain) function returns a tuple `(imin, jmin, imax, jmax)` of grid indices that define a bounding box containing the upstream area of a collection of lat-lon points (outlets), and a mask for that domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bounding box: (3874, 1883, 4054, 2180)\n"
     ]
    }
   ],
   "source": [
    "outlets = [[39.8195, -75.3820]]\n",
    "ij_upstream_area_bounds, mask = st.define_upstream_domain(outlets=outlets, grid=\"conus2\")\n",
    "print(f\"bounding box: {ij_upstream_area_bounds}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now visualize the upstream area mask and the grid point corresponding to the outlet. Note that to add the outlet to the figure, we need to convert the global grid coordinates of the outlet to local mask coordinates. To do this, we first compute the global grid coordinates using the `hf_hydrodata` function [to_ij](https://hf-hydrodata.readthedocs.io/en/latest/hf_hydrodata.grid.html#hf_hydrodata.grid.to_ij) and then subtract the lower-left corner of our mask `(imin, jmin)` from them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x152a326a2c10>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARgAAAGdCAYAAADJ82znAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABBBElEQVR4nO3dd3wUdf7H8dfsZrNpJCGkQwihd5AWgg0FKSKCxDtAEAvCgYCnWPhxVtQ7Tj0bSrEdiIKiInJERCkCIqFFmpRAIBBaQijpbcv8/kAXlpa2m8nufp6Px945s9+d/UwIb2a+853vKKqqqgghhBPotC5ACOG+JGCEEE4jASOEcBoJGCGE00jACCGcRgJGCOE0EjBCCKeRgBFCOI2X1gVUhdVq5eTJk9SpUwdFUbQuRwiPoqoq+fn5REdHo9Nd/xjFJQPm5MmTxMTEaF2GEB7t2LFjNGjQ4LptXDJg6tSpA8BN3IkXBo2rEdejGLz5765fCNb7MeFYd7YltSX67c1alyWqwYyJDSy3/T28HpcMmD9Pi7ww4KVIwNRmimKgTh0dgXodn7XewpR6Zex4R/7MXNofdy9WpHtCOnmF8ygKOn9fu1UGxYLOz0+jgkRNk4ARTqP26MA3v6+krv5ioEwL28lH+39CMRo1rEzUFJc8RRK1lE5Pq60KoYYCAMINP+Gn87Zrold0ROn9uG3bOX5+uDvq1t1aVCpqiASMcAiv2Bj2ToliScT7V4TK5fSKjin1DjL36e6EL+yG73dbKv193r4GAsMCZJiCg6mqSnF+KUU5RThipigJGFFpOn9/inu2tlt3Os5A+uBZwPXD5VL7b/qMpuYHaZzdEeXXHRX+XEzbaBKf7Y9foK8EjKOpYLFYSNt6hJ/mrCcvu6Bam5OAEZXXLJa1H33kkE2l9ZzH4y26kNrzwiVPa0EB1/un09vXQOKz/anfKBqj3schNQh7KlYC76hDZNNwZo+ej8VsrfK2JGCE5t6J2kbpfhMAA+99BCV55zXbBoYF4Bfoi1Hvgx59TZXoYfQE+NQhMLSQ4MhAzh7PqfKW5CqSqLBDC26g7+953PvFGodv26gYMCoGRs79nszHe1yznaIoclpUIy78nHX66kWEHMGI6/KKjGDvyw0BeK3rV/w1INep3zcq8Aw/D9vNFt8eNJi+0anfJZxPAkZclxoSRPpdjulvqai5DX/h1eHZrNx9Cz5Jlb/CJGoPOUUS16QYvLEEaDMg7rnQ/SyY9Rb6eiEgp0TX1aJrU1atXal1GVclASOu6eSkLiQtmafZ9zfwCmDprpUonVqX39gDvPfhuwy6b2C1t/N/Lz3Do0+Nc0BF5ZNTJMGBj7sw9cblV6xvaZyFQdH2So1B0fPwF0lMf3cE4TMd2CdjseC7YyteZ7Ixh4ZR3LEr6OWqlKNJwAhCI/MYG3RS6zKu6a8BuSSN/J2UwB6ELdxf7e0FrPmRsDdfwXA607bOFB5J9pPPU3B732pv/1rKykp5fcZrfP9TEgWFBbRt1Y6pTzxL+zbt+XbZYv711qts+3m7rf2qtSuZ8PR4Urem8e2yxbz/0XvAhVMigOkvvMaQgYlXfM+pzJP8+93p/LppAzqdjs4du/Dsk8/TILoB7334Lku+/9ZuO/PnfE585+5O2Wc5RfJwxYO60apeZvkNNTY/dj33DP2F4js6VGs7AWt+JGrKRLxO2++z1+ksoqZMJGDNj9Xa/vW8PuN1flzzI/9+8Q2WfLaU2AaxPPLYQ+Tk5pT72TvvGMDDI0bTrHEzNvyQzIYfkrnzjgFXtDOZTYx+7CH8/fxZ8NGXfPHxIvx8/XnksYcpM5Xx8MhH6N/7Tm5OuMW2nRvad3LC3l4gAePhXn37Q+bHrte6jAp5NXw306fOhXKmabwmi4WwN18BVC7vNlb+mOQk7K1XwWKpVp1XU1RcxJeLF/LMY1O49cZbadq4Ga8890+MRh+++d/X5X7ex8cHPz8/9HovwkLDCAsNw8fnypHMy3/6HqtV5Z/PTadF0xY0iWvK9Bf/zanMk2xJ2Yy/nz8+Rh+8vb1t2/E2VPz2jsqSUyThUhp6+WGNLgHVCKXmSn3Wd8dWu9OiyymoGLJO4btjK8UOPmXIOJ6ByWyiU4fOtnUGLwPt27TnUPohQm4Iccj37D+4n4zjR+l0q/2RXmlZKRnHMxzyHZUhASNcUkjDfPRGPaeKgjEeLqrQZ7zOZDu0nSPpdDrUy+7BMplNld5OUXERbVq25T+vvHnFeyF161W5vqqSUyQPpa8XwqEFN9DSUKh1KVUSqLNST2+mgd95SmL9K/QZc2iYQ9tVRsMGDTEYDPy2M8W2zmQ2sXvvLprGNaVucAiFRYUUFV8My/0H9tltw2AwYLVe//StTYvWHD12hHp16xEb08juVSegzsXtOOE08GokYDyUEuBP2m1zCddX7C9nbRWst9A8KAtLiH+5fTPFHbtiCo+8Sg/MBSoKpoioC5esHczP14/hiSN4fcZrrN+4jrTDB3n+1WcpKSnh3kF/oUPbDvj6+PLWzDfJOH6UZSv+x7dJ39pto35UA46fPM6+1L2cyzlHWVnpFd8zsP8g6gbXZfxT49i2fSvHThxjc8omXv3Py2Rmnbqwnej6pKalcvjIYc7lnKvSkVJFScB4IMVoxBISqHUZDmNUdLSMyUb1M4LuOqN+9Xqyn3weuLKb98/l7MnPOW08zFMTn6bv7X155sWnuOf+QRw9fpSPZ8wlKDCI4KBg3nj5TdZvXMvAYQP4/sdlTBozye7zfW/vy80JNzNq/EgS7uhG0o9JV3yHr48vn3/wBdGR0Ux85lHu/Gtfnn1lKqWlpQT4BwDw18FDiYuNI/GBe0i4o5vdUZWjKerlJ34uIC8vj6CgIHoySJ4qUAVnxySwbdpsrcuokhJTNEdynycuNhwfnyvDZO+pcAynrz9J0lXHwUREkT35OaeOg3ElFixkncnkw/ELyD56zu49s2piLUvJzc0lMPD6/1BJJ6+HOfBRVxb1ngFu+jypRhFZHPIPwyf92h2/Bbf3peDW3jKStwZIwLg56803kDHxYofetA5L6GZ0z3AB8NMpGLwq0IGp1zv8UrS4kgSMm8uL8yH1Ztc8HRKuTzp53Zg+OIiyOp411UGZasVilV/r2kKOYNxY2uxGHLh1ltZlOJj6x/+qcJXLzWlZkXhnVW8mfHGBqqpXDP6rLIl6N9Vsq5E1N87UugyH0+sKQTVjutbQDZe7Jlo7WbFgLrOQd6Z6AzHlCMbN6OrUIWNeQz6K+IAGXgFal+NwXrp8/Lz2kH0mGC8vX3SXzHZ38HwIhtwSLNTMKFV3ZcXCuZxzbF22k7KismptSwLGzSg+PuxJWAC4X7gAKIpKZOAijpxvyNGMYC7coqhyuswf3dl81BoaAu+uVFXFXGZh67KdrJ1b/Qm+JGCEy/HWn6dpvecxWULJs5r4obAhPyU65m5kT6eqKnlnCqt95PInCRjhknSKBaNXFvftG4iu1zHgXLmfETVPOnndSPHgbvzf5to5u7zwTHIE46q6tyf7H/Z307YN288tHvS45ubrRxH5hQ++HNO6FHENEjAuyHrrDRwaYuRwl/lal6Ip428B+C6Vpz/WZnKK5GK8GtQnfSwc/sscrUvRXFmgildkhNZliOuQgHExvVbsI+22uVqXUSukPjybsgXOm7BaVJ+cIrkInb8/U3Ylk+BTirtOtSDcjxzB1DK6Dq3I+6EJeT80Qd+mBQBK5zYUfhvOzT5mjDLBlp1nG31PxtfttC5DXIMcwdQSefd1pzhUR34jK4c6XOhfafrQOAIyepDXxMrhdguQfw+u1NPXyhddPmYK8VqXIq5CAqYW0DdrTL8p63kxbK/d+rT7pCNXuDYJGI0pBm8WrPmcuno/rUsRwuHkmFtDakIH/n3gFwmXamrnbeC19M3om8ZpXYq4jASMhlS9QkejUesyXJ5e0V34Ocqk3bWOBIyGDKfzab1xJCZVphgQ7kkCRkPWo8eJfdVKkeqYW+M9XVHjuuiDg7QuQ1xCAkZDlm6t+WH5QoJ0vlqX4hbWfvIRJ0e10boMcQkJGA0Y1kbxzpGNfPj5e1qX4ra8YhrwzpGNeMXGaF2KR5PL1DXg/IMJhD5w1Lb8WqPFtPKWK0fO8NdHVrPhnib4G4pp5e1HnYWF5JVFX7P9kbMhxNz7ew1W6FkkYJxAH1qPjEda0ODtFM7/tRP6oadZ3mL5JS3klMhZ/hGaCqGptuUv49Zct/2ushKGT51Mw7d/w1pS4uzyPE6lTpGmT59O165dqVOnDuHh4QwePJjU1FS7Nj179kRRFLvXuHHj7NpkZGQwYMAA/Pz8CA8P5+mnn8ZsNld/bzTm1bjRhXuJbm3KzknvY+nWmq6P/UZyh8Valyauob23D3smzUIJuv5D3EXVVOoIZt26dUyYMIGuXbtiNpv5xz/+QZ8+fdi7dy/+/v62dmPGjOHll1+2Lfv5XTwdsFgsDBgwgMjISDZu3MipU6cYNWoUBoOBf/3rXw7YJY0oCvq5JXzf7Ls/VuhYuUimVXAFFtWqdQluq1IBs2LFCrvlefPmER4eTkpKCrfccottvZ+fH5GRkVfdxk8//cTevXtZtWoVERERdOzYkVdeeYUpU6bw0ksv4e3tmvN7TDiQSm/fzYBr1u+pvi/yYWb3G7GcOa11KW6pWleRcnNzAQgJsX9kxIIFCwgNDaVt27ZMnTqVoqIi23vJycm0a9eOiIiLM5H17duXvLw89uzZU51yNNXMcAY/nYSLKxl8sC+vPz4Ky5mzWpfitqrcyWu1Wnn88ce58cYbadu2rW39fffdR2xsLNHR0ezatYspU6aQmprKt99+C0BmZqZduAC25czMzKt+V2lpKaWlFye4zsvLq2rZQgBwZ+qdnPqmEeFJMqevM1U5YCZMmMDvv//Ohg0b7NaPHTvW9t/t2rUjKiqKXr16cejQIZo0aVKl75o+fTrTpk2raqlOpRi8oWMLvJUN5TcWmiqylvF5fiMA8mbGEP6NhIuzVekUaeLEiSQlJfHzzz/ToEGD67aNj78wEVBaWhoAkZGRZGVl2bX5c/la/TZTp04lNzfX9jp2rPY8pkLXrBErln5GE4N7PqrVnfxcEsjiVuEsbhWO/zebtS7HI1QqYFRVZeLEiSxZsoQ1a9YQF1f+7fE7duwAICoqCoCEhAR2797N6dMXO9VWrlxJYGAgrVu3vuo2jEYjgYGBdi8hKqPb9r8ws1M3rcvwOJU6RZowYQILFy5k6dKl1KlTx9ZnEhQUhK+vL4cOHWLhwoXceeed1KtXj127dvHEE09wyy230L59ewD69OlD69atuf/++3n99dfJzMzkueeeY8KECRhdbOqCvPu603VyitZliAowmfVYpO+uxlXqCGb27Nnk5ubSs2dPoqKibK9FixYB4O3tzapVq+jTpw8tW7bkySefJDExkWXLltm2odfrSUpKQq/Xk5CQwMiRIxk1apTduBlXURilY0b0Vq3LEOXot38A/BRSfkPhcJU6glFV9brvx8TEsG7dunK3Exsby/Lly8ttV5vpWzSlKEoGaLmC0181JHyOdOhqQe5FqqKi98yktZVJuWuzUtUEgHL9fxeFE0nACLe0pdTEi+1vByCseCuSMdqQgBFuZ+SRnpx8qSmG/G1al+LxZMIp4XYO5dbD8JOES20gRzBVoPboQIOAo+U3FDXqvfOx5Fp8ycwIIZBDWpcjkICpkukLPqKzUW5srG1+GNIVS2oazZGhA7WFnCIJIZxGAqYKnhk9ntEZN2ldhvjDIVMBtz30COrR41qXIi4jp0hV4LUmhUN5jbQuQ/wh32rA+8dtyLDH2kcCRris7woD2FYYx+HCUOC81uWIq5BTJOGynv9oFFs76jl7o4RLbSUBI4RwGgkYIYTTSB9MFeV/E0VPBrO27Xdal+Kx6vc7yoEG8XbrmiwuQ7duu0YVictJwFRR6IfJ5BQnMO2p1rwYttfuvTOWQl49fSuvRSZjVAwaVej+VrT8Hlrar+uY9igR5c8YImqInCJVQ/BnyWwe1obzlouPZSlVTXye14Z9XSzsLAOTatGwQiG0JQFTTZZ9BxneuCf7yi6ETNt1Y1jRIRRUlRea38gDR3prXKEQ2pFTJAdQTWX8ffh4VL1Cs6w8LH88Z1s1lWFFnjYgPJcEjIMoyTtRgMtPiA581oI7/uLPylbLrvYxIdyanCI5WdicZM4tasCL2W20LsUj5DWzoHRtp3UZ4g8SMDUg9MNkNj/SiVPmAq1LcXuHEz+g6NUCdHXqaF2KQAKmxqhbd/NQs14SMjXg57aLmbLjV63LEEjA1Ci1tJSHEsfTd8go4r4fo3U5bkuv6EjwKeWuPefxalBf63I8mgRMDVO37oZNuzBkS/+6MxkVA5PqHiX19TCst96gdTkeSwJGIwEZ8HSm/OI7W1rPeZxt5aN1GR5LAkYjYXOS2TmxPRnmAtsr11qsdVlu57i5AJ1Z6yo8lxyna0jZuJMxsTfbltPejOfQMHlapCP9redI6h1O1roMjyUBo7VLnvctjzh1gnKepy6cS06RapGGP5iIWyZXl4T7kICpRQyrUoj7WqaudoRcazFD0u6A4hKtS/Focook3FKqyYuSIVYsZ89oXYpHkyMY4Za6GQ0s27kSr0YNtS7Fo0nACLelV3SM+Wk1o1KPceDjLlqX45EkYGqR3JHdMT99Tusy3Mpg/wJG1DnLtJu+48DsblqX43EkYGqRMzfA+nZLtC7DLY0KPMMHd8zVugyPIwEjPIa3YsErMgIURetSPIYEjPAYPX2tfP/bj3jVj9a6FI8hASM8zqS1q8gf2l3rMjyCBIzwOP38SjH5yWlSTZCBdrVE8eBuBDWXK0g15UxXC5BQobZ1DxSj/LrDqfW4KwmYWqLr89t4M+o3rcvwGOmDP4TBFWvbdO2DtDhyod/GfOKk84pyQxIwQpQjrec82HrhKZ2DWt2GJS9P65JchvTBCFEB8/LCGdTzL1jy87UuxaXIEYwQFZBv8cVy8LDWZbgcOYIRQjiNBEwtkZofwXF5ZlKtdMpcwK6CBlqX4ZIkYGoJU89T9Pz6Ka3LEFfRY+XjZMQXal2GS5KAEUI4jQSMEMJpJGCEuI770m8jaIe31mW4LLlMLcR1ZLzdnIhvNmpdhsuSIxghhNNUKmCmT59O165dqVOnDuHh4QwePJjU1FS7NiUlJUyYMIF69eoREBBAYmIiWVlZdm0yMjIYMGAAfn5+hIeH8/TTT2M2y/M9m0/bS48nxmldhgAsqpX+/YYRsGyH1qW4tEoFzLp165gwYQKbNm1i5cqVmEwm+vTpQ2HhxUt4TzzxBMuWLePrr79m3bp1nDx5kiFDhtjet1gsDBgwgLKyMjZu3Minn37KvHnzeOGFFxy3Vy7KkpeHz1mT1mWIPygnTqOWlmpdhktTVLXqz9bMzs4mPDycdevWccstt5Cbm0tYWBgLFy7k3nvvBWD//v20atWK5ORkunfvzg8//MBdd93FyZMniYiIAGDOnDlMmTKF7OxsvL3L71DLy8sjKCiIngzCSzFUtfxaydS7M2vmf6J1GR4t3VTAgK3jiH0gHWuhjH+5nFk1sZal5ObmEhgYeN221eqDyc3NBSAkJASAlJQUTCYTvXv3trVp2bIlDRs2JDn5wgPIk5OTadeunS1cAPr27UteXh579uy56veUlpaSl5dn93JX+hIL3xf5aF2GR/uluBEx9/4u4eIAVQ4Yq9XK448/zo033kjbtm0ByMzMxNvbm+DgYLu2ERERZGZm2tpcGi5/vv/ne1czffp0goKCbK+YmJiqll3r6Tbs4P22HSmwyiNPheurcsBMmDCB33//nS+//NKR9VzV1KlTyc3Ntb2OHTvm9O/UkrWkhKFdBvHGuSZalyJEtVQpYCZOnEhSUhI///wzDRpcvAksMjKSsrIycnJy7NpnZWURGRlpa3P5VaU/l/9sczmj0UhgYKDdy92ZM7NY9mwvEnYmal2KEFVWqYBRVZWJEyeyZMkS1qxZQ1xcnN37nTt3xmAwsHr1atu61NRUMjIySEi4MP9pQkICu3fv5vTp07Y2K1euJDAwkNatW1dnX9yO79ItZB0O1boMjzIzJ4aXVkqoO0qlRvJOmDCBhQsXsnTpUurUqWPrMwkKCsLX15egoCBGjx7N5MmTCQkJITAwkEmTJpGQkED37hceE9GnTx9at27N/fffz+uvv05mZibPPfccEyZMwGg0On4PXZzhvI61xTp6+lq1LsUjvL29N80mbda6DLdRqSOY2bNnk5ubS8+ePYmKirK9Fi1aZGvz9ttvc9ddd5GYmMgtt9xCZGQk3377re19vV5PUlISer2ehIQERo4cyahRo3j55Zcdt1dupNFzybw0abTWZQhRJZU6gqnIkBkfHx9mzpzJzJkzr9kmNjaW5cuXV+arhRAuSO5FEkI4jQSMEMJpJGCE+MPqYj1qtlxocCSZD0aIPzz3/BiaLtykdRluRY5ghBBOIwEjhHAaCRghhNNIwAghnEYCxhWoYFItWlchRKVJwLgA44+/MbjbQEpVmU5TuBYJGBegdm9L7HfnMLrZ9KC1SfenxxGyJl3rMtyOBIwLKAv2ZlZ9GZ/hDIdMBcT9byx1/7cHc2ZW+R8QlSID7YRH210WSfNxW5DJMJxDjmCEEE4jASOEcBo5RarlTvxfDx4euULrMtxS3313Yf5nBF6kaF2K25IjmFou6LCVT1ITtC7DLR09WxevNRIuziRHMLVcwFebMBR0BckY4YLkCEZ4pFLVhNUiv/7OJj9h4ZF6THuMxqP2aV2G25OAER5JZwLVVKZ1GW5PAkYI4TQSMEIIp5GAqeUULy+s3orWZQhRJXKZupY78FZn9iW+D8id1ML1yBFMbadHpmkQLksCRgjhNBIwQginkYARQjiNBIwQwmkkYGq5Vm+couO/H9W6DCGqRC5T13LmIxlEL9fTImI8AB8Nn80tPhoXJUQFyRGMC7CkpdPo2WQaPZvM+O0jWFFk1Lokl/bYya74ZctzpmqCBIyLaZC4h0d/fJAzlkKtS3FZhwaH45O0ResyPIIEjAtqNnELgx+frHUZQpRLAsYFpf+zO8/8+zOty3BZf1m1lcJ747UuwyNIwLigqGQLjycP07oMl/Vg4GlK68ivfk2Qn7IL8lm2hfpL5P4kUftJwLgonVklw1xAhrkAiyrPJawss6+Czs9P6zLcngSMi/JZtoUxsTczJvZm5uVFa12Oy9n27PscnttU6zLcngSMK1NVUFUWJ95Chy3Dta7GpegVHYqial2G25ORvG7AsvcAwZ90o82WR7H4qOx/ZLbWJQkBSMC4DZ9lW2iwDPSBgdx9az8+bryYcL2/1mUJDyenSG7GkpdH6a2ZzM3pSK61WOtyhIeTgHFTa9r50/WXcVqXITycnCK5sWZPn6F/8DBMoX6sWvhfrcsRHkiOYNyY+fgJrL/vx7BlP+3eeZQ9ZXLKdKnJ7VZzYFY3rctwaxIwbqBoSDz6Ni2u+b61qIjo1zdyxFy3Bquq/cYGneTtOxZqXYZbk1MkF6UYvNGHhgAwbvo3TFv8V5q9F3ndz/gou2uiNCFsJGBcVOFdN/DLzA9syyMenA0PalePEFcjAeOCDr2RwOeJ7wN6rUsR4roq3Qezfv16Bg4cSHR0NIqi8N1339m9/+CDD6Ioit2rX79+dm3OnTvHiBEjCAwMJDg4mNGjR1NQUFCtHXF3itFIxtftOL64DRPv/IHuPhIu1TXm2I28+tr9Wpfh1ip9BFNYWEiHDh14+OGHGTJkyFXb9OvXj7lz59qWjUb7OWRHjBjBqVOnWLlyJSaTiYceeoixY8eycKF0uF2NV1QkpwbHsavH+xgUCRZH+eVoY2I/Tta6DLdW6YDp378//fv3v24bo9FIZOTVOxz37dvHihUr2Lp1K126dAHgvffe48477+Q///kP0dHuf2ewPiIcxdvbtmw+dvyq7bwiI8Bg4OwtDfjt+dnIKZHjHDIVUFogk6c7m1P6YNauXUt4eDh169bl9ttv59VXX6VevXoAJCcnExwcbAsXgN69e6PT6di8eTP33HPPFdsrLS2ltLTUtpyXl+eMsmtMzLJ8Pmhw4V/OImsZic17Yi0quqJdtx+P8WLY3pouzyMMe+lpms+Voxdnc/g4mH79+jF//nxWr17Na6+9xrp16+jfvz8Wy4XHRGRmZhIeHm73GS8vL0JCQsjMzLzqNqdPn05QUJDtFRMT4+iyNeOn82bK7k1YenayrVOMRsYfTOPJer9pWJmolNUNOPJqgtZV1DoOD5hhw4Zx9913065dOwYPHkxSUhJbt25l7dq1Vd7m1KlTyc3Ntb2OHTvmuIJrkqKQ8XU7Roett1vd09eKxXjhj0LXviXZ38Qy0C+PAJ08Ya22U7y8OL64Da81XszfBv/IwXe7a11SreL0kbyNGzcmNDSUtLQ0ACIjIzl9+rRdG7PZzLlz567Zb2M0GgkMDLR7uRp9aD1Oj09gU8KHdDNeOZ/usd4GssclkDaiLimdv0KvyCBrZzqTYMbUp4vduuJB3TDf3rnC29BHhJP1t25sjZ9LR6ORySGH+Ve/RZwZmwCK4uiSXZLTx8EcP36cs2fPEhUVBUBCQgI5OTmkpKTQufOFP8w1a9ZgtVqJj3ffR0lYmtZn+3OzAN+rvp82QiaJqknpd31Ep8ihRO1rgPnYcbxiY2jwzEF2ZUUTmx5ra2c9lYW1pOSq2zA1i2b7s7OAix32d/mfwu/pT5n1aTvUS/oNPVWlA6agoMB2NAKQnp7Ojh07CAkJISQkhGnTppGYmEhkZCSHDh3imWeeoWnTpvTt2xeAVq1a0a9fP8aMGcOcOXMwmUxMnDiRYcOGecQVJFF7/NZlEd+tCWB2s6b8Z90iWnn7QRzw68U2Nz32N/y/2Vzhbb6S3Z0dNwBIuEAVAmbbtm3cdttttuXJky88YfCBBx5g9uzZ7Nq1i08//ZScnByio6Pp06cPr7zyit1YmAULFjBx4kR69eqFTqcjMTGRGTNmOGB3aqcTU3rw+bi3AbksWtsM8Msl7PBvNDdUoL9Lp2fs/oNEeuUA4K9sRv5Mr6/SAdOzZ09U9dqTJf/444/lbiMkJMT9B9UpCicWt6aObwn31V9NR6P8ItZGBkXPjRXIFn3zJuTOULnLfxNG5c8+NPs/04Sdieg/CcWfih/xuDu5F8mJlnX+gDhDgNZliGo4cbtKUFQPiqJVDrSfDVz7gXenU8No+s2mmivOBUjAOIHi5YWuSSP0yi9alyKqKX3whzBY6ypclwSMEyhtm7N8+UJAjl6EZ5OAEaIaiqxl3NtzKEpxKS3y92LRuqBaRkZzOVj+sO74v3e6/IbC5S0uCKTnc3/HcjgD84mTWFz8HjlnkCMYB8of2p38oXl802SV1qUIJ+q3fwBHzoRgOuFP03lyw+T1SMA4iL5FU5o+vpf5sevLbyxcWvFb9YlN2qJ1GS5BAsZB3vjxM9p4X/02ACE8lQSMgzwdP4iAxRZ2n4qm8aSLfTDh3xUxt6FcrhaeSQLGQSxZpznzUmdiisyYM7Ns64stIRpWJYS2JGAcyLAqResShKhV5DK1EMJpJGCczKrKxEPCc0nAOFnBnSaazx+vdRlCaEICxskseXnoS+QoRngmCRghKsiiWmn680P4ZcgtARUlV5GEqIBcazFf5DWl2dgDWAsLtS7HZUjACHEZk3rxnug/H9U7N7cVP7QJBiRcKkMCRohLtPhlFE3GXXzuVuhyC7+mNaHFhEOAnBpVlvTB1IC4hVl0eONRrcsQFWAx67GcP297nXqmCc3fLJWpGKpIAqYGlEUHkdfCrHUZogL6Nd/Lqck9bMu6X7Zj3SHPB68qOUVyAp2PD5YbWtiWDw31Iv3uDzWsSFTU+/U3M3/sIRa81UDrUtyCBIwztGzMT4s/1boKITQnp0hCCKeRgBFCOI0EjBCXmJbdmn9+8xety3AbEjBCXOKz3+Np9JxM5O0oEjCOpiioBr3WVQhRK0jAONjJpxNYtEQuSQsBcpnaoQ583IXXbllAkE6eLiAEyBGMQwXt8mb+yR7lNxS1UqdtQ6n3o4/WZbgVOYJxoIgZG0n37wHNta5EVIXxy7oELpQOXkeSIxghgAJrCYpV6yrcjxzBCAEMvfEv1MnYrHUZbkcCxkHCNgbTyO8sE/znaV2KqIT1JTDtkdEYTuwEVdW6HLcjAVNN+rAwUp9twoKYNwnV+2tdjqikc5YAvNakINHiHBIw1eAV04DsXjEc+utsQMJFiMtJwFSRYjRydERDfn9sltaliCoyqRZyLH5al+HWJGCqKH9pfba3ew+Q2wJcVYfkB2g4/ABQpnUpbksuU1dBRHIgc1t9ZptxXrimL7t8jNfKelqX4dYkYKrg75GraG6QPhdX9+X5bqT/GKd1GW5NAkZ4rG8PdqDB9I1al+HWpA9GeIRS1USu1b6vxWySX39nk5+w8AgJKSMJH3LIbl0T6y6NqvEcEjDCI6iqgmqWZ1PVNOmDEW6v/ZbhGL8K1roMjyQBI9ye14pgghZs0roMjySnSJWhKOiDAtHLnSsu45S5AJ2cGWlGAqYS9E0asXTdNxgUo9aliAoaffv91EuTSaS0IqdIFZRzfwIjvl8vo3ddxK6yEvrdPRLrkWNal+LR5AimgsoCFUbUOat1GaICXj3TkkXzbyd6mwyi01qlj2DWr1/PwIEDiY6ORlEUvvvuO7v3VVXlhRdeICoqCl9fX3r37s3Bgwft2pw7d44RI0YQGBhIcHAwo0ePpqCgoFo7IsSfvj/Rhuj/SLjUBpUOmMLCQjp06MDMmTOv+v7rr7/OjBkzmDNnDps3b8bf35++fftSUlJiazNixAj27NnDypUrSUpKYv369YwdO7bqe+Fkujp1MMtd/S7htKWQ/GJ5MkBtoahq1ecJVBSFJUuWMHjwYODC0Ut0dDRPPvkkTz31FAC5ublEREQwb948hg0bxr59+2jdujVbt26lS5cuAKxYsYI777yT48ePEx0dXe735uXlERQURE8G4aUYqlp+hZ1Z1pyUzl85/XtE9XV47VEi35WjF2cyqybWspTc3FwCAwOv29ahnbzp6elkZmbSu3dv27qgoCDi4+NJTr7Qk5+cnExwcLAtXAB69+6NTqdj8+arT7pcWlpKXl6e3aumtErxYkmH/9bY94mqu3XsWOrP36d1GeISDg2YzMxMACIiIuzWR0RE2N7LzMwkPDzc7n0vLy9CQkJsbS43ffp0goKCbK+YmBhHln1V+nohHFpwA/8IX0dDrwCnf5+oPv+D57CcP691GeISLnGZeurUqeTm5tpex44599KjV1wsp4a1JO22uYTLRN4u40T/cPTNGmtdhriEQwMmMjISgKysLLv1WVlZtvciIyM5ffq03ftms5lz587Z2lzOaDQSGBho93IWfWAgxwfVZ/uzMteuq9n5zCyO3huJPiLc9kJRtC7Lozk0YOLi4oiMjGT16tW2dXl5eWzevJmEhAQAEhISyMnJISUlxdZmzZo1WK1W4uPjHVlOlWTMi2HnMxIurmrPpFks3/6T7eXVqKHWJXm0SgdMQUEBO3bsYMeOHcCFjt0dO3aQkZGBoig8/vjjvPrqq/zvf/9j9+7djBo1iujoaNuVplatWtGvXz/GjBnDli1b+PXXX5k4cSLDhg2r0BUkZ1MUuc/InTy68kdy7k/QugyPVemRvNu2beO2226zLU+ePBmABx54gHnz5vHMM89QWFjI2LFjycnJ4aabbmLFihX4+Fwcm7BgwQImTpxIr1690Ol0JCYmMmPGDAfsjhD2BviV8Lx0o2mmWuNgtOLMcTAZL/XgnkEb+FeEzHbmLuKWjaHJlxb0P/+mdSluQbNxMO4g7t19JH1+k9ZlCAdKH/gRxVNy8IqMKL+xcCgJmMtkfBLNrielk9fd/Nr+W/6+YY3WZXgcuZv6Eh23w7x6HwAysM7dHVnUnnc7L7Itz+h+E5Yzcre8o0nAAPq6dcn4JJp59T4gSkbtuq2OxhxOLmkNwFvtvqKfX6ntvRl6mefHGTw6YMy3dya3iTdldRR+7z4LOXJxb+F6f3bHL9S6DI/ieQGjKHg1bADAqceK2NXtE40LEsJ9eVzA6END+T55mdZlCOERPCZgDr4fz3/7f4RBMSMXz4SoGR7xN+3Ax114qfdievpaudHHI3ZZVFLepwFYb71B6zLcjkf8bXvuxiRGBZ7RugxRi21o/y1pww2Ye3XWuhS34rYBo/PxQd80Dn3TOPx1peV/QHi89Ls/5OjDVq3LcCtu2wdT0L8Dv8z8QOsyhPBobhkwaW93J+metwB5FICouE6vjKf5l/uxaF2IG3G7gDnwcRdevOlbWnlLuIjKMeaqMqevg7ldH8yk7mt4MPB0+Q2FEE7nPgGjKOibNyFAX1J+WyGuoihMZzelg75ZY3R+FT8S1ofWkyk6L+M2AaMLCGDJz18yNuik1qUIF7VzyixS3/xj2lZFYeaq+eTf2a7Cnz8wtRkdlhxxTnEuyi1mtDP17sw7H8+kvbc8MlRUz3lLEUfNenSKShuDNxnmInKtBv579iZSu5iu+Tn/9WHMifuOIJ03m0sN/KtVPGrp9YdH6IODeGn7Kryxcv/MJ1zmedoeN6Od1Vsn4SIcoq7ej45GI+29fdArOuIMAXQ0Gnkq/GfKVsbanTLp27TAujoG6+oY/tlwKeF6f4yKgQ7exeV+j3pjR3Tf+dHNaKCj0YjF6Jz98YppgHV1DPrgIOd8QXnfr8m3OpC5V2eO3iXPvhHO1dArgFWtl9D+8Ynoyi6sK6pv5VCrPyetqlhfzdnRCZTWVShoZiK9+Q+29cZu58gb3p3ALzY5rGalcxsOJgZyoNVsWj32KHFfZmE5cMhh268Ilw4YfeNYDj1kJf12mXJBOJ9e0bFnYvnTqepRoHVTdKYrR9SMf2oJo4OufETy3PafMtLyMIFfOKRUAE72DOLAgxfq3TduFglp4wiUgKm4Lp8f5Ke45VqXIYSdAJ0PK75fUKnP3PfJE8S86hp9MJXh0gEjhLtY+7c3yBnj2G0G634FtH0olEsHzNCgFMB5z6kWoqaE6/0Jd8NpgV36KlKcQebQFaI2c+kjmC/yQ/BV7XdhUMAxgnS+GlUkhLiUSwfMV53rX/Ho2Kzfg3g65Oo95Rb14lwfesWlD96EcAlu97fs55sa0Pin0Vesb/HJeAa27217fVcop1dCOJtLH8FcjSUnl+bvlnDj8nF265vszMZy9pxtuUx1wx41IWoZtwsYAHX7HgK226+TSYSEqHlud4pUUd9md+aAqVDrMoSoMQUNdDU+nYTHBsz5G88x4NcJmFSL7SWEO9v9+CxS/xUCuprrHnCL6RqqSufnh+J98fOPb9tIH79r35IvhKsrspbx3vk2rGlX9RG+lZmuwS37YCrKWlQERReXTegBCRjhvvx03jwUvIP9m/rZ1m1Oaue0+6A8OmCE8EThen/mNvzFtjxmoJUthT2IfNvxIeOxfTBX82nmjewrKyq/oRBu5KOYX5k67gvUhA6gOHZuJQmYS+TedJa/bn9E6zKEqHHD6pzn+2/+i75OHYduVwJGCOE0EjBCCKeRgLmM9w9B3LL7Hq3LEMItSMBcJvTDZIq+idS6DCHcggSMEMJpJGCEEE4jASOEcBoJGCEEKaVltPxqAtZyHndbWXKrgBCCjUXNaPrEJhx957McwQjh4UyqhXyLc57tLgEjhIdrsXIs6zsHOWXbcop0mYPvduff/RdqXYYQNaLNe4/ScnEWFlOZU7YvRzCXOPLPBMb3WslfA3K1LkWIGuHd/RznuoY5bfsSMIBi8KasX1c+Hf7+NZ+pJIQ72t71SyzDz2G9+QanbN/hAfPSSy+hKIrdq2XLlrb3S0pKmDBhAvXq1SMgIIDExESysrIcXUal6END+Pqjd+juI48yEZ5na6evuHvOGqds2ylHMG3atOHUqVO214YNG2zvPfHEEyxbtoyvv/6adevWcfLkSYYMGeKMMirMfCqTkY1vY0WRUdM6hHA3Tunk9fLyIjLyyhsGc3Nz+eSTT1i4cCG33347AHPnzqVVq1Zs2rSJ7t27O6OcClFNZVhw7GxeQriKkYH7yN91A78khGItdNzjfJxyBHPw4EGio6Np3LgxI0aMICMjA4CUlBRMJhO9e/e2tW3ZsiUNGzYkOTn5mtsrLS0lLy/P7iWEcJy6ej+errcXRe/YbgKHB0x8fDzz5s1jxYoVzJ49m/T0dG6++Wby8/PJzMzE29ub4OBgu89ERESQmZl5zW1Onz6doKAg2ysmJsbRZQshnMDhp0j9+/e3/Xf79u2Jj48nNjaWr776Cl9f3yptc+rUqUyePNm2nJeX55CQ0fn4oHh725YN8oBZ4eGUwDpQUAhWx/xdcPpAu+DgYJo3b05aWhp33HEHZWVl5OTk2B3FZGVlXbXP5k9GoxGj0fEdsAf/25J9t35iWzYochVJeC6Doue7zf+j99/G45O0xSHbdPo4mIKCAg4dOkRUVBSdO3fGYDCwevVq2/upqalkZGSQkJDg7FIA0Nety117znPXnvMs6zELg6K3vYTwdAZFz2NvfcnRlx3z99HhRzBPPfUUAwcOJDY2lpMnT/Liiy+i1+sZPnw4QUFBjB49msmTJxMSEkJgYCCTJk0iISGhRq4gKV3akvq4geV1f/5jjZ/Tv1MIV5MYkEfKwF9ZrL+RRs9e++JLRTg8YI4fP87w4cM5e/YsYWFh3HTTTWzatImwsAvDkd9++210Oh2JiYmUlpbSt29fZs2a5egyriqvSQCHbp9TI98lhCv7V8Qugu4uZs2zVX+GNTghYL788svrvu/j48PMmTOZOXOmo7/6unR+fph9ZZyLEBVlUCzo69bFcv58lbfhMfcipc5uSfI/azbUhHBlk0MO88Xu5ej8qt6V4BEB478+jBU930OveMTuCuEwQTpfhv92oMo3Q3rE37gBYbtpbqjeuaQQnmpU4BkK/pHHwfmdOPLPyl1d8ogJp17f1YewTt9wt3+Rw7d9xlLII4evfrPmoqZJGBWDw79TiJqW3GExAFOyOrLtHxX/nEcETKOhu/j7ByPo2v9torwCqr29dFMBpj9ujEzKb0fxrVeZbkJRWHfQj0aGHEJ0EKq3P4LKtRaTa7XQ0AH1CFFbeUTAADT/21YShz3Fxreqf5l64k3DMB87fv1GqsqbTdsAkPFSD/aNtb8U3zHp7zRdYGLlornVrkeI2spjAgYg6Ps93HHioSr/pf61xMqrg+/DeuJgpT4XN/Mg/b8dZreudXYG1vM59Bs4gkVLPyZIV7X7tISozTwqYKz5+Rj2l3PkcQ1PZ97Aqk8SCN+1sdKftWRnQ3a2fS1//sdve7n57Sd57m8LZC5g4XY8KmAAKC2l7767+Kz5Ij7Nbc+qrFZXNOkems60sD225WnZrfnf8u40mln5cCmXqhL15kY239eEvwb85vjtC6EhjwsYS14e9Mpj9s6ufPvfnkS+e2Vo/PDALQybtpVW3n7sKysi6d1bafTf6t2TIYQn8ohxMFezsYP3VcMFoO6nyTx161AAnrptOCESLkJUiccGTHnMx05yZ6+/YD5atT6bytqf2IC478fUyHcJUVM87hSpwqwWLPsqd7WoOsxHMmg6vx4d9j5qW/fKo/OcMjhQiJoiRzC1iO6X7US+vdH22l7UqNLbKLKWccvuezhlLnB8gcLjtfE9zvn7u1W4vRzB1GIHC8M5HbKV8D9GAaeUlv/84INlEfj2TeeXA/XlsrdwiEt/79oYT/L2M7O5/bOKfVYCphbL7pFD95mTOXzPB5y2FPJss1tRzWatyxIepMhaxnNtemItuniqblZNQMX6JiVgarmWL6Zx5ztDQFVRzelalyM8RJPVD9Hi5VxQVaxFVf+9k4Cp5SxnzsKZs07/nokn4vl1XmcAvn3mdeIMchOmJ/mpyMCUty9exWz8ezGWg4ervV0JGA+Qay3mjp2jSGo/jzeyb2LdqaZXtCnYFEbMHyOV+978KO93+YI+fqaaLlXUsIcybmbP2UiyTwbT3Akj1SVg3NTWgsZEem0HYHdJc+oOOMjMnd34eXZ36n105cDBuly8JB83bBd//3oYSd1m00SOZNzawf+0pu43m6nrpO1LwLipXZ1UdtHRbt2mDgbqUbFRyQ3/spt7nnyGXU/WzBMfhHuSgBHCA1lUK3fdOYI6qTsv3tnvBBIwoloaf/M36u7RUVgf9j8yW+tyxCVGHunJvk+vnC0AABXC9qagmsofW1UdEjDimoLTzLTfMvy6bZp+UYKSvJOIpnF0aD+cLV3nyzzEGkvYmUhhqTfqr3WJ/uDaHbdqDdQiASOuyXfpFnyXVqytJS2dyCF6zh0pJcpLAqamnbYUsqM0GICQhwsJPHVI24L+IAEjhBvo+9towgft/2MpU9NaLiU3OwrhQtq8/ygDug2g7+D7tS6lQuQIRggXYigA8/ET6M6dp+tz4wHoOn67xlVdmwSMEC6i+457CU67MLraWlRkm2lxXUQPFGdea64GCRghaqldZSUcMYXYlutNNGM+vPWKdg2mO2EyegeRgBGilhoxazLRr18aHke0KqXKpJNXOI7VwiPdEnn8VBetKxG1hASMcChzZhbFFhkHU13t3nqUmGXZ5Tes5SRghKiFGqw4V6OTzjuLBIxwuO3ZDdhVVqJ1GaIWkIARDld3wEHuXfiE1mWIWkCuIglRC7T473iavJtmW7aeSdWwGseRgBFO0WThedrkP8qeSTJh1aXazniUwCNXjoprsvMMlmzX79S9nASMcArr7/uJDu4Ik7SupHYpamAhKrkM3Tr74f0WjepxNumDEU6jLzEzPy8Ui1pLx7Fr4PCQD0gbZkDfrLHWpdQICRjhNOq231nYrjFnrcVal1KrpA/6kND5Z7Quo0ZIwAinUk1lPNj+Lv6R1V7rUoQGJGCE01nOn2fzk11pvXFkpT532lLITY/9jRVFxqtsVIWNRbAk/8L/W2piAkhRWdLJK2qE15oUgsO608l7KL91WXTNdu23DCf/rP+FhTIdzRdv4cQrdbGbpe37ApTns1FOXewaVaP0qK+EwYAA5uWFM+3Xu8utqXPzI3zTZFVVd6lCJp6I5/td7a5Y73vYmxhq713QjiIBI2pMnUWbCDjWgTkf12dc8AkANpVY2Fp8scMzepoOdfs2u88tzuxEkXUPAE1WnmbA42lXzlidaUEZk8mumdFM87mb5qO3UZ4DT/aAJ50bMD+u7kTz/6vYs6jckQSMqFHKxp0s7RTDyLRDBOh8uP+riTS2+wu454rPWG47SRJ10akqn//xr75y+XZVUBUIf6aIlgVbsSqXtxBakD4Y4TLakk0YxVeEy58UFaIKcmiL+w1Yc1VyBCNqnLW0lAHjHkPVQ7P92RUeZFaPit1AWdF2MYuP0zl/PCkvyQPjnEUCRtQ8VcUnaQtQuRGsZ/FxaDvzkQwilpYQ12ksvw14h7p6v0pUU742ySOI3OLZgwzlFEm4jN8JIxvfaz5L2QqcxpffCavwNi1Zp2k+bguZThirH/WeEb9vNzt+wy5EAka4DKuiMIuOKHBFyFi50PE7m45V6uAtUr0cdktDkbWMImsZiipjczQNmJkzZ9KoUSN8fHyIj49ny5YtWpYjXMAGpT4vk8BZfO3Wn8GXl0lgg1K/Stt9ru1t/OVQ32rX99zpdiS2uI3EFrehW7+j2ttzdZr1wSxatIjJkyczZ84c4uPjeeedd+jbty+pqamEh4drVZZwARuU+mxUo2lLNvUo4Sw+/E5YtS5NWwsLKXqyMbfVG0NeQ68qd/yaVD3WwsIq1+FuFFXV5jguPj6erl278v777wNgtVqJiYlh0qRJ/N///d91P5uXl0dQUBA9GYSXIhNMC8fyiopk3z9ibcudOh6q0Ijfvx1PYMOSG2r1c4ocwayaWMtScnNzCQwMvG5bTY5gysrKSElJYerUqbZ1Op2O3r17k5zsuaMeRe1gPpVJs0kXb03Y81IPuE7ATMtujRXFI8KlsjQJmDNnzmCxWIiIiLBbHxERwf79+69oX1paSmlpqW05NzcXADOmK4eMC+FglpIS8vKv3gFcZC1jQ3c/rEXFRLIOcw3XpgUzFx5fW5GTH5cYBzN9+nSmTZt2xfoNLNegGuFxpi+l7vTrNfiypiqpVfLz8wkKCrpuG00CJjQ0FL1eT1ZWlt36rKwsIiMjr2g/depUJk+ebFu2Wq0cPXqUjh07cuzYsXLPA91JXl4eMTExst8eojbut6qq5OfnEx0dXW5bTQLG29ubzp07s3r1agYPHgxcCI3Vq1czceLEK9objUaMRvs5QXS6C1fYAwMDa80PvibJfnuW2rbf5R25/EmzU6TJkyfzwAMP0KVLF7p168Y777xDYWEhDz30kFYlCSEcTLOAGTp0KNnZ2bzwwgtkZmbSsWNHVqxYcUXHrxDCdWnayTtx4sSrnhJVhNFo5MUXX7zi1MndyX7LfrsSzQbaCSHcn9zsKIRwGgkYIYTTSMAIIZxGAkYI4TQuGTDuPo/MSy+9hKIodq+WLVva3i8pKWHChAnUq1ePgIAAEhMTrxgV7SrWr1/PwIEDiY6ORlEUvvvuO7v3VVXlhRdeICoqCl9fX3r37s3Bgwft2pw7d44RI0YQGBhIcHAwo0ePpqCgoAb3ovLK2+8HH3zwit+Bfv362bVxhf12uYD5cx6ZF198kd9++40OHTrQt29fTp8+rXVpDtWmTRtOnTple23YsMH23hNPPMGyZcv4+uuvWbduHSdPnmTIkCEaVlt1hYWFdOjQgZkzZ171/ddff50ZM2YwZ84cNm/ejL+/P3379qWk5OLE3iNGjGDPnj2sXLmSpKQk1q9fz9ixY2tqF6qkvP0G6Nevn93vwBdffGH3vkvst+piunXrpk6YMMG2bLFY1OjoaHX69OkaVuVYL774otqhQ4ervpeTk6MaDAb166+/tq3bt2+fCqjJyck1VKFzAOqSJUtsy1arVY2MjFTfeOMN27qcnBzVaDSqX3zxhaqqqrp3714VULdu3Wpr88MPP6iKoqgnTpyosdqr4/L9VlVVfeCBB9RBgwZd8zOust8udQTz5zwyvXv3tq1z13lkDh48SHR0NI0bN2bEiBFkZGQAkJKSgslksvsZtGzZkoYNG7rdzyA9PZ3MzEy7fQ0KCiI+Pt62r8nJyQQHB9OlSxdbm969e6PT6di82bUn3F67di3h4eG0aNGC8ePHc/bsWdt7rrLfLhUw15tHJjMz8xqfcj3x8fHMmzePFStWMHv2bNLT07n55pvJz88nMzMTb29vgoOD7T7jbj8DwLY/1/vzzszMvGKKVS8vL0JCQlz659GvXz/mz5/P6tWree2111i3bh39+/fHYrnw+ANX2W+XmA/G0/Tv39/23+3btyc+Pp7Y2Fi++uorfH19r/NJ4S6GDRtm++927drRvn17mjRpwtq1a+nVq5eGlVWOSx3BVHYeGXcRHBxM8+bNSUtLIzIykrKyMnJycuzauOPP4M/9ud6fd2Rk5BUd/GazmXPnzrnVz6Nx48aEhoaSlpYGuM5+u1TAXDqPzJ/+nEcmISFBw8qcq6CggEOHDhEVFUXnzp0xGAx2P4PU1FQyMjLc7mcQFxdHZGSk3b7m5eWxefNm274mJCSQk5NDSkqKrc2aNWuwWq3Ex8fXeM3Ocvz4cc6ePUtUVBTgQvutdS9zZX355Zeq0WhU582bp+7du1cdO3asGhwcrGZmZmpdmsM8+eST6tq1a9X09HT1119/VXv37q2Ghoaqp0+fVlVVVceNG6c2bNhQXbNmjbpt2zY1ISFBTUhI0LjqqsnPz1e3b9+ubt++XQXUt956S92+fbt69OhRVVVV9d///rcaHBysLl26VN21a5c6aNAgNS4uTi0uLrZto1+/fuoNN9ygbt68Wd2wYYParFkzdfjw4VrtUoVcb7/z8/PVp556Sk1OTlbT09PVVatWqZ06dVKbNWumlpSU2LbhCvvtcgGjqqr63nvvqQ0bNlS9vb3Vbt26qZs2bdK6JIcaOnSoGhUVpXp7e6v169dXhw4dqqalpdneLy4uVh999FG1bt26qp+fn3rPPfeop06d0rDiqvv5559VLkzdbvd64IEHVFW9cKn6+eefVyMiIlSj0aj26tVLTU1NtdvG2bNn1eHDh6sBAQFqYGCg+tBDD6n5+fka7E3FXW+/i4qK1D59+qhhYWGqwWBQY2Nj1TFjxlzxj6gr7LdM1yCEcBqX6oMRQrgWCRghhNNIwAghnEYCRgjhNBIwQginkYARQjiNBIwQwmkkYIQQTiMBI4RwGgkYIYTTSMAIIZxGAkYI4TT/D/xEyfdrWQyWAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(mask, origin='lower')\n",
    "imin, jmin, imax, jmax = ij_upstream_area_bounds\n",
    "outlet_global_i, outlet_global_j = hf.to_ij(\"conus2\", *outlets[0])\n",
    "outlet_local_i, outlet_local_j = outlet_global_i - imin, outlet_global_j - jmin\n",
    "plt.scatter(outlet_local_i, outlet_local_j, marker='o', c='r', label='outlet')\n",
    "plt.legend()"
   ]
  }
 ],
 "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}