{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculating cyclone density"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import all the necessary libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from functools import partial\n",
    "from pathlib import Path\n",
    "\n",
    "from octant import RUNTIME\n",
    "from octant.core import TrackRun\n",
    "from octant.misc import check_by_mask\n",
    "\n",
    "RUNTIME.enable_progress_bar = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the common data directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_dir = Path(\".\") / \"sample_data\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data are usually organised in hierarchical directory structure. Here, the relevant parameters are defined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset = \"era5\"\n",
    "period = \"test\"\n",
    "run_id = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct the full path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "track_res_dir = sample_dir / dataset / f\"run{run_id:03d}\" / period"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load land-sea mask array from ERA5 dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "lsm = xr.open_dataarray(sample_dir / dataset / \"lsm.nc\")\n",
    "lsm = lsm.squeeze()  # remove singular time dimension"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now load the cyclone tracks themselves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='671' class='' max='671', style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [671/671 00:03<00:00]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tr = TrackRun(track_res_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Classify the tracks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = [\n",
    "    (\n",
    "        \"aa\",\n",
    "        [\n",
    "            lambda ot: ot.lifetime_h >= 6,\n",
    "            partial(check_by_mask, trackrun=tr, lsm=lsm, lmask_thresh=0.5, dist=50.0),\n",
    "        ],\n",
    "    ),\n",
    "    (\n",
    "        \"bb\",\n",
    "        [\n",
    "            lambda ot: (ot.vortex_type != 0).sum() / ot.shape[0] < 0.2,\n",
    "            lambda ot: ot.gen_lys_dist_km > 300.0,\n",
    "        ],\n",
    "    ),\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "    <div>\n",
       "        <style>\n",
       "            /* Turns off some styling */\n",
       "            progress {\n",
       "                /* gets rid of default border in Firefox and Opera. */\n",
       "                border: none;\n",
       "                /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
       "                background-size: auto;\n",
       "            }\n",
       "            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
       "                background: #F44336;\n",
       "            }\n",
       "        </style>\n",
       "      <progress value='671' class='' max='671', style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
       "      100.00% [671/671 00:23<00:00]\n",
       "    </div>\n",
       "    "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tr.classify(c, True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "tr.rename_cats(**{\"aa\": \"cat1\", \"bb|aa\": \"cat2|cat1\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "  table.octant {\n",
       "      white-space: pre;\n",
       "      border: 1px solid;\n",
       "      border-color: #f9f9ef;\n",
       "      font-family: monaco, monospace;\n",
       "  }\n",
       "  th.octant {\n",
       "      background: #084469;\n",
       "      color: #fefefe;\n",
       "      border-left: 1px solid;\n",
       "      border-color: #0a507a;\n",
       "      font-size: 1.05em;\n",
       "      min-width: 50px;\n",
       "      max-width: 125px;\n",
       "  }\n",
       "  .octant-la {\n",
       "      text-align: left !important;\n",
       "      white-space: pre;\n",
       "  }\n",
       "  .octant-ra {\n",
       "      text-align: right !important;\n",
       "      white-space: pre;\n",
       "  }\n",
       "</style>\n",
       "<table class=\"octant\" id=\"140570351867832\">\n",
       "    <tr class=\"octant\">\n",
       "<th class=\"octant octant-la\" colspan=\"7\">Cyclone tracking results</th>\n",
       "</tr>\n",
       "    <tr class=\"octant octant-ra\">\n",
       "<td rowspan=\"4\">Categories</td>\n",
       "<tr class=\"octant\">\n",
       "<td class=\"octant octant-la\"></td><td colspan=\"4\">671</td><td class=\"octant octant-la\">in total</td>\n",
       "</tr>\n",
       "</tr>\n",
       "<tr class=\"octant\">\n",
       "<td class=\"octant octant-ra\">of which</td><td class=\"octant octant-ra\"colspan=\"4\">131</td><td class=\"octant octant-la\">cat1</td>\n",
       "</tr>\n",
       "<tr class=\"octant\">\n",
       "<td class=\"octant octant-ra\">of which</td><td class=\"octant octant-ra\"colspan=\"4\">44</td><td class=\"octant octant-la\">cat2|cat1</td>\n",
       "</tr>\n",
       "<tr class=\"octant\">\n",
       "<td class=\"octant octant-ra\">Data columns</td>\n",
       "<td class=\"octant octant-la\" colspan=\"6\">lon, lat, vo, time, area, vortex_type</td>\n",
       "</tr>\n",
       "<tr class=\"octant\">\n",
       "<td class=\"octant octant-ra\"                rowspan=\"2\">Sources</td>\n",
       "<tr class=\"octant\">\n",
       "<td class=\"octant octant-la\"                             colspan=\"6\">sample_data/era5/run000/test</td>\n",
       "</tr>\n",
       "</tr>\n",
       "</table>\n",
       "        "
      ],
      "text/plain": [
       "<octant.core.TrackRun> [671 tracks]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculate density"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, arrays of longitude and latitude are needed to define the grid (1 degree step in this case)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "lon_dens1d = np.arange(-20., 50.1, 1)\n",
    "lat_dens1d = np.arange(65., 85.1, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cell-based method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simplest method to calculate cyclone density is to loop over grid cells and sum up tracks' occurence in each of them. By default, the cell counts are divided by grid cell areas so the units are cyclone (tracks) per $km^{2}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "dens = tr.density(lon_dens1d, lat_dens1d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here it is assumed that both arrays represent grid cell centres. Otherwise, the `density()` method should be called with `grid_centres=False`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that if `subset=` keyword is not given, the computation is done for all existing categories:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'cat1': <xarray.DataArray 'point_density' (latitude: 21, longitude: 71)>\n",
       " array([[0.      , 0.      , 0.      , ..., 0.      , 0.      , 0.      ],\n",
       "        [0.      , 0.      , 0.      , ..., 0.      , 0.      , 0.      ],\n",
       "        [0.000621, 0.      , 0.      , ..., 0.      , 0.      , 0.      ],\n",
       "        ...,\n",
       "        [0.002655, 0.001327, 0.000664, ..., 0.      , 0.      , 0.      ],\n",
       "        [0.      , 0.      , 0.      , ..., 0.      , 0.      , 0.      ],\n",
       "        [0.      , 0.      , 0.      , ..., 0.      , 0.      , 0.      ]])\n",
       " Coordinates:\n",
       "   * longitude  (longitude) float64 -20.0 -19.0 -18.0 -17.0 ... 48.0 49.0 50.0\n",
       "   * latitude   (latitude) float64 65.0 66.0 67.0 68.0 ... 82.0 83.0 84.0 85.0\n",
       " Attributes:\n",
       "     units:    km-2\n",
       "     subset:   cat1\n",
       "     method:   cell,\n",
       " 'cat2|cat1': <xarray.DataArray 'point_density' (latitude: 21, longitude: 71)>\n",
       " array([[0.      , 0.      , 0.      , ..., 0.      , 0.      , 0.      ],\n",
       "        [0.      , 0.      , 0.      , ..., 0.      , 0.      , 0.      ],\n",
       "        [0.000621, 0.      , 0.      , ..., 0.      , 0.      , 0.      ],\n",
       "        ...,\n",
       "        [0.      , 0.      , 0.      , ..., 0.      , 0.      , 0.      ],\n",
       "        [0.      , 0.      , 0.      , ..., 0.      , 0.      , 0.      ],\n",
       "        [0.      , 0.      , 0.      , ..., 0.      , 0.      , 0.      ]])\n",
       " Coordinates:\n",
       "   * longitude  (longitude) float64 -20.0 -19.0 -18.0 -17.0 ... 48.0 49.0 50.0\n",
       "   * latitude   (latitude) float64 65.0 66.0 67.0 68.0 ... 82.0 83.0 84.0 85.0\n",
       " Attributes:\n",
       "     units:    km-2\n",
       "     subset:   cat2|cat1\n",
       "     method:   cell}"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dens"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output of `TrackRun.density()` is an `xarray.DataArray` with useful metadata and coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "arr = dens[\"cat1\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "OrderedDict([('units', 'km-2'), ('subset', 'cat1'), ('method', 'cell')])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arr.attrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAELCAYAAAD+9XA2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztvXmcXnV59//+zEz2DQigEMAECSBSQI2IYrUKslXB+mANVR9ULLUFl/rUFrRFpeX3aGu1WpeaImopshjkMVWURVRE2cJOgEjYQ8KahUDIMjOf3x/nOzPnnLmXM1tmcs/1zuu85nzXc51zn9zX/V2u65JtgiAIgmCotI22AEEQBEFrEAolCIIgGBZCoQRBEATDQiiUIAiCYFgIhRIEQRAMC6FQgiAIgmEhFEoQBEEwLIRCCYIgCIaFUChBEATBsNAx2gIMJxM1yZOZNtpiBEEBTZ5USHvT5lG5Lp1dRTk6O0dFjm11/83YwNpnbO8ylD6Ofss0P7umq3lF4JY7N19h+5ihXG+s01IKZTLTeJ2OGG0xgqBA+977FtJd9/5+VK6rZ9cW0p1PPT0qcmyr+2/G1V78yFD7eHZNFzddsVeluu273b/zUK831hnxKS9Jfy1pmaS7JV0oabKk70l6SNLt6TikTtuTJd2fjpNHWtYgCIKBYKC74r/xwIiOUCTNAT4GHGD7RUmXAAtT8adsL27Qdifgs8ACss/tFklLbK+t1yYIgmBbYsxWV5vyGg9si0X5DmCKpA5gKrCqYrujgatsr0lK5CqgpecfgyDY/ogRSh8jqlBsPw58CXgUWA2st31lKj5H0p2SviJpUo3mc4DHcumVKS8IgmBMYEyXqx3jgRFVKJJ2BE4A5gG7A9MkvQ84E9gfeC2wE/B3tZrXyOv3qUg6VdJSSUu3MjZ2jwRBMH7oxpWO8cBIT3kdCTxk+2nbW4EfAW+wvdoZm4HvAofWaLsS2DOX3oMa02W2F9leYHvBBGoNdIIgCEYGA1240jEeGGmF8ihwmKSpkgQcAdwraTeAlPdO4O4aba8AjpK0YxrpHJXygiAIxgwxQuljRHd52b5R0mLgVqATuA1YBPxM0i5k01q3Ax8BkLQA+IjtD9teI+kfgZtTd2fbXjOS8gZBEAwEA1vHyfpIFUbcsNH2Z8m2/+Z5a526S4EP59LnAeeNnHRBEASDx+NoOqsKLWUpvz3QvtOOhXTXmjCrGQsM9HPJ1y/X7Xj5vEK6c5Qsw8eKRbq2bC2kW+r/gKEr9EkvoVCCIAgGSWYpH/QQCiUIgmDQiK6aFg7jk1AoQRAEg8RAd0x59RIKJQiCYJAY2BJhpXoJhRIEQTAEuh1TXj2EQgmCIBgkmaV8KJQeQqEEQRAMEiO6Ysqrl1AoQRAEQyCmvPoI1RoEQTBIjNji9kpHFSQdI2m5pBWSzqhRPknSxan8RklzU/5sSb+U9Lykr9fpe4mkWn4Th41QKEEQBIMkM2xsq3Q0Q1I78A3gWOAA4CRJB5SqnQKstb0P8BXgiyl/E/APwN/U6ftdwPODuceB0FpTXtOmwMEH9yY7Hn+2UOz1G+o27Vq/vmHX7Qfu13u+afcZhbIJVy6tLGLZzUT7rFmV20JjOct9aVZRzs5HV1ZuW75Os/LRYKjPLt9+oO4/GtXvfOChhm0bPcux+Jyr0Oiz8DNFn67byz1VZRgX5Q8FVth+EEDSRWTxpO7J1TkB+Fw6Xwx8XZJsvwBcJ2mfcqeSpgOfBE4FLhkuYWvRWgolCIJgG2KLLg/bRE+tKLWvq1fHdqek9cBs4JkG/f4j8K/AxuEStB4x5RUEQTAEulGlA9i5J7psOk4tdVUlSm2lSLa9laVDgH1sXzagmxokMUIJgiAYJJkdSuXf5c/YXtCgvEqU2p46KyV1ALOARnGiXg+8RtLDZN/3u0r6le0/qir0QIgRShAEwSAxYqs7Kh0VuBmYL2mepInAQmBJqc4S4OR0fiJwjV0/wpftb9ne3fZc4I3A70dKmUCMUIIgCIZE1zDZoaQ1kdPJQp23A+fZXibpbGCp7SXAd4DzJa0gG5ks7GmfRiEzgYmS3gkcZfue8nVGkhFXKJL+miwKo4G7gA+SPZQFwFbgJuAvbG+t0bYrtQF41PbxIy1vEARBVYbbUt725cDlpbyzcuebgHfXaTu3Sd8PAwcOWcgGjOiUl6Q5wMeABbYPJNO6C4ELgP2BPwCmkAv7W+JF24ekI5RJEARjjm63VTrGA9tiyqsDmCJpKzAVWGX7yp5CSTeRLT4FQRBsVwxwUb7lGdEnYftx4EvAo8BqYH1JmUwA3g/8vE4Xk9P2uhvSnGAQBMGYwYguVzvGAyM6QpG0I5ll5zxgHfBDSe+z/d+pyjeBa23/pk4Xe9leJWlv4BpJd9l+oHSNU8ksQJnMVLj+jt6yziHIvv79ry+kd1z2XO/51hnFxzZ5RtEivRGaMb2Q7ly1ehDS1aafBfIALJKbWS/3szIfwD13bSh6KGjUdkB1h2hxPWoW2931o5CPFSvyZp9v+XMaK3Jva2yq7uAaF4z0WO1I4CHbT6dF9x8BbwCQ9FlgFzKXADWxvSr9fRD4FfCqGnUW2V5ge8EEJg3/HQRBENSlmlFj9ziJmTLSCuVR4DBJUyUJOAK4V9KHgaOBk2zX/LkmaUdJk9L5zsDhFH3aBEEQjCoGutxW6RgPjOhYzfaNkhYDt5LNQN0GLAJeAB4Brs/0DD+yfbakBcBHbH8YeAXwbUndZIrvC9t6T3UQBEEzYlG+jxGf/LP9WeCzVa5reylpC7Ht35FtKw6CIBiTGEWArRyxmhQEQTBITCzK54knEQRBMGg0nPFQtntCoQRBEAwSw7ixgq9CKJQgCIIhECOUPkKhBEEQDBJbMULJEQqlDrPOv76QzhvLTL21WLdrIB1vqB/XfqiUrftnX/lAscL0ab2nzeKel2mfXrTw19Qpved+oRhZtOv55xu2LVtZN2IgdZvR7x6mTe0973zyqW1ynVrXKtdv1LZMQ7kPO7iYvuGO2vVqMJLPvfx+bO+MFxuTKoRCCYIgGCRZgK320RZjzBAKJQiCYJBki/KxhtJDKJQgCIIhEJbyfYRCCYIgGCRhKV8kFEoQBMEQ6I4RSi/xJIIgCAaJzbAG2JJ0jKTlklZIOqNG+SRJF6fyGyXNTfmzJf1S0vOSvp6rP1XSTyXdJ2mZpC8M063XJBRKEATBIDGis7u90tEMSe3AN4BjgQOAkyQdUKp2CrDW9j7AV4AvpvxNwD8Af1Oj6y/Z3p8sntThko4d1M1WoOGUl6TnGpUDAlbb3nf4RAqCINh+GEZL+UOBFSmgIJIuIot4mw/bcQLwuXS+GPi6JNl+AbhO0j75Dm1vBH6ZzrdIuhXYY7gELtNshPKA7ZkNjhlksU2CIAjGHT3bhqscFZgDPJZLr0x5NevY7gTWA7OrdC5pB+AdwC+q1B8MzRbl/1eFPqrUCYIgaEEG5HplZ0lLc+lFthcVOuuP+12weZ1+SOoALgS+1jMCGgkaKpQqFx5J4YL+tE2t74pjx0uLrjXKsZW3HLRX7/nkDUX3F917vqR4nWfWN5Sj++ln+843bmxQEzR7x0L62ZMO6j2f9fCWQlnHVUtpRPsr+mZX7/tYsd+dbi3OU+9ywe1FOXLuYmB43a3k6edapImrkUauSDrKz+7NxdmKHRf3te33OZRcreSfHYAfWVm/bRP8xkOKfd/6+7p9dR5SmIVhwmPPUo/ORx6rWzZWGUC8+GdsL2hQvhLYM5feA1hVp87KpCRmAWsqXHsRcL/tf6sq7GCopFolvUvS/ZLWS3pO0oYK6ys9bf867S64W9KFkiZLmpd2KNyfdixMrNP2zLSbYbmkowdyY0EQBCONDVu72ysdFbgZmJ++HycCC4ElpTpLgJPT+YnANbYbjlAk/ROZ4vnEgG5uEFQdq/0zcLztWT1rJ7ZnNmskaQ7wMWCB7QOBdrKH9EXgK7bnA2vJdi6U2x6Q6r4SOAb4ZtoFEQRBMCboMWwcjjWUtCZyOnAFcC9wie1lks6WdHyq9h1gtqQVwCeB3q3Fkh4Gvgx8QNJKSQdI2gP4DNmusVsl3S7pw8P4CApUNWx80va9Q7jGFElbganAauCtwJ+l8u+T7Vr4VqndCcBFtjcDD6UHeChwPUEQBGOEAUx5NcX25cDlpbyzcuebgHfXaTu3TrfbzJS/2bbhd6XTpZIuBv4fsLmn3PaPGrW3/bikLwGPAi8CVwK3AOuSNobaOxlIeTfk0jXrSToVOBVgMo1dfQdBEAwn4RyySLMRyjty5xuBo3JpAw0ViqQdyUYa84B1wA/JjHbK1JoDrLSbIe2SWAQwUzs13e0QBEEwnESArT6a7fL6IICkw23/Nl8m6fAK/R8JPGT76dTmR8AbgB0kdaRRSq2dDFBtx0MQBMHoUd3GZEwjqbz4X4s1tj/QqELVNZR/B15dIa/Mo8BhkqaSTXkdASwls9w8EbiIbMfCj2u0XQL8QNKXgd2B+cBNFeUNgiAYcQx0tsYI5RVAo8V6kbmFaUizNZTXk40odpH0yVzRTLIdWw2xfaOkxcCtQCdwG9n01E+Bi9J2ttvIdi6QdjIssH1W2t1wCZnbgU7gNNsDirYbBEEwkrTQGspnbP+6UQVJn2/WSbMRykRgeqo3I5f/HNkIoym2Pwt8tpT9INmOrXLdJeT2Xds+BzinynWCIAhGg1ZQKLYvGY46zdZQfi3pOuAPbDfVTtsTHbOL7m86n61vvTualOXccuDLes8n3v1IsfKUyYVk58rHi33lrNA7KdIxsWhb2v1s0fi2+8UX68rVNqVogZ6XEYCSnJtyt7T0898plB3w7b8spOdetraQ7rqzb/f6/GLVfpQ9BQzUGnwsULYcn/3TkhV6g3sqfy5MKE0q5OzhOvYobqAsvzvlvrqvK3ohKD/rPCrVLb972zOtEmAr2fh9mGyt+uf5NXNJf2/7n6r003TyL00z7TRYQYMgCFqZblTpGON8G3gz8CzwtbR23cO7ajfpT9VF+dvSLoAfkvMu3MwOJQiCoKVxa0x5AYfaPgggBej6ZtqVexIDMIysqlB2ItNcb83lNbVDCYIgaGUMdHa3xC6v3jnvZM5xqqSzgGvI1tErUUmh9NijBEEQBH20yhoKmTeUY2z/vCfD9tmSVtHfLVZdqnob3kPSZZKekvSkpEuT07EgCIJxja1Kx1jG9vvyyiSXf67tCVX7qTpW+y7Zdt7dyfxp/U/KC4IgGNe0yKJ8PyQtal6rSFWFsovt79ruTMf3gF0GerEgCIJWwh7WEMBjjUbBwGpSVaE8I+l9ktrT8T6yRfogCIJxjOjqbqt0bIcMOKRp1bv8EPCnwBNk8UxOTHlBEATjmlZYQ6mF7WMG2qbqLq9HgeObVhxjtL+yGEO7a1n9uNcdc3YvpDsfr+/YuGw1rJKFul/cVEjnrczbdyrGCC+jHXcoyvHAQ8Vr/7pvYDgUi+PyPXjWtEK6e6/ijGbbbcuL5bnnV7ain7hiUiFd9kIw5wu/6z0/9quHFcp2Prroru3xtxaf1x7P9n1O3S8t2tt233I32yMdL5/Xe9696olCWfnZlp9l22sO7D3XQ0Xrds8rWr93NXg+bZOKnhLK/x+61xQ9FjT6v1Wm3NfmfXcrpCc81xtiibZ1LxTKyu//WKOFfHkBIGkBWYTHl5HpBwHusVFpRiWFImkX4M+Bufk2tmOUEgTB+MUFDzatwAXAp4C7aOxRpyZVp7x+TBbk/moyT8E9RxAEwbhmOHd5STpG0nJJKySdUaN8kqSLU/mNkuam/NmSfinp+WTpnm/zGkl3pTZfk9RImKdtL7H9kO1Heo6qz6KqpfxU239XtdMgCILxgGHY1keSg8ZvAG8jCzB4s6Qltu/JVTsFWGt7H0kLgS8C7wE2Af8AHJiOPN8iC5N+A1m8+mOAn9UR47OSzgV+wQDCvfdQdYTyE0nHVawbBEEwThBd3dWOChwKrLD9oO0tZAEITyjVOQH4fjpfDBwhSbZfsH0dmWLpk07aDZhp+3rbBv4LeGcDGT4IHEKmdN6RjrdXER6qj1A+Dnxa0mZgK30LNTOrXigIgqAVGcYdXHOAfLyClcDr6tWx3SlpPTAbeKZBnytLfc6pUxfgYNt/MBCh81Td5TWjUbmkV9peNlghgiAItkfsASmUnSUtzaUX2c5bo9fqqLzkX6XOUOrfIOmA0jRbZaqOUJpxPjXiy0vaD7g4l7U3cBbwemC/lLcDsM72ITXaPwxsALqATtsDttwMgiAYSQawbfiZJt9hK4E9c+k9gLL9Qk+dlZI6yDZLraE+K1M/jfrM80bgZEkPka2hDP+24QrUfKK2l5PNx/UsOD0OXGb733obSv8KrG/Q91ts1xvOBUEQjCrDuG34ZmC+pHlk35ULgT8r1VkCnAxcT2Zgfk1aG6kjm1dL2iDpMOBG4H8D/95AhgEbM+YZLoVS5ZEeATyQ34KWtq/9KcU4K0EQBNsFRnQPk1uVtCZyOnAF0A6cZ3uZpLOBpbaXAN8Bzpe0gmxksrCnfZrRmQlMlPRO4Kg0dfWXwPeAKWS7u+rt8AI40nYhLrekLwD9tjDXYrgUShUWAheW8v4QeNL2/XXaGLhSkoFvl+YbAZB0KtmWOCYzdRjFDYIgaM5w2jXavpxsa28+76zc+Sbg3XXazq2Tv5T+W4nrcaKkTbYvAJD0TWBSkza9DJdC2dKoUNJEMtctZ5aKTqK/kslzuO1VknYFrpJ0n+1r8xWSklkEMFM7FT7bRu4gyi4tuh8vpjvmvqxYvvrJum3buhu/Um2Tcq5ZSnXLblo0qfFnl+/LhxTdX7SvL8rFpuLH0vWSWX3nE9sLZROefr6QXv2GonKevWNx48fklX312+57sFDmF0pubUrPksl9bj467yv+lpj2szuKVQ99RSG98ZC+6eCuicWZ1hdf9YZCetefluQqPeuudevqyvjUkY02wsCuV/e5Ocm/GwBtOxXd53SvWUcj/GzOrUnp/ejY7aXFyrOK+2M6c+5UCu8Z0PbshkJaOxTlyrsM2rh/0dXO1GXFeyq/l17xKPUoy8GMolufF19SDK8xafVzfYkJxa+ktoMPKKbXF12zsHlzIdm116695+2PFn0bbjis9B7+aDFDZmCL8tsD7wKWSOoGjgXW2P6rqo2rBtg6XNK0dP4+SV+W1Pvp2D6sfmtIgt1qu/ctTQtK76K4aF/A9qr09yngMrJ92kEQBGMHVzzGMJJ2krQT2bTYh4G/BZ4Dzk75lag6+fctYKOkg9OFHiEzkKlKrZHIkcB9tlfWqI+kaZJm9JwDRwHbp/e/IAhalhbxNnwLsDT9/SXZ7ts/zuVXouqUV6dtSzoB+Krt70g6uUpDSVPJXAn8Ramo35qKpN2Bc20fB7wEuCy5nekAflArRGUQBMFo0grOIW3Pa16rOVUVygZJZwLvB/4wbQGuFGfY9kYyS85y/gdq5K0CjkvnDwIHV5QvCIJgm2ODt8/gWQUkvdr2rUOtU/VJvIfMyOVDtp8gM93/l4ptgyAIWha72jHG+a6kHXvWUmodZFuWG1LV9coTki4F5qesZ8gWyYMgCMY3Y19ZVGEW2XpJQ9f2zTqpGmDrz8lsPXYCXk42QvkPMmPFIAiCccp2seDelHo2LAOl6pTXacDhZNvISIaIuzZsEQRBMB5ogW3Dw0XVRfnNtrf0BPpKNiTj5BEFQRDUofUMG4dEVYXya0mfBqZIehvwV8D/jJxYw0M/i908B+5TTN+9opDsfLh+1MuOvfYopL2+aJHcdcDcQrr9nodzlRvr4RcW7FVIT31oViGtzZ19iYefaNiXNxWtiJ27p44D9yuWrS5aFe/541I46eeK9+hdduw7Lz1LrSz25bVFS3HtunOfHE2e5Qtzihbam3bsG1TvetNzhbK2zqJ1/7o3zy2kd/hd0eSpY1KfxX5Zxl1+WF/mcn1NLsrYvXuxLk0s5bVL3yZIlZ/zxqL3g42l92Pa2voTBY3eYYC2F/v+f0y5rngd71my0J9essJ/em0xPSl3z5tLjjNK9zDr7mJbuvreNU8qfiXp4ceLdUvPunPv3Ypy3Lq8r69S3en3N/4cBk0olF6qTnmdQbYgcxeZPcnlwN+PlFBBEATbDS005SXpUkl/LGlQe6Gr7vLqlvTfwLXJJX0QBEEA242yqMi3yMIAf03SD4Hv2b6vauOqvryOB24Hfp7Sh0haMghhgyAIWgeTTXlVObYDbF9t+71kARMfJnPK+ztJH5TU1Ji96rDms2SOGdeli94OzB2UxEEQBC1Eixg29iJpNvABMieRtwFfJVMwVzVrOxBfXut7dnkFQRAEie7W+V6U9CNgf7Kw7u+wvToVXSypqZPIqgrlbkl/BrRLmg98DPjdYAQOgiBoJbQdjT4qcG4K8tWLpEm2N9te0Kxx1SmvjwKvJPPn9QOyGPCfGKikQRAELUXVHV7bj9L5pxp511dt3HSEkjwLf972p4DPDECwIAiCFmf7WXBvhKSXkrnUmiLpVfT59JoJ1WOrNx2h2O4CXjMYIYMgCFqeYRyhSDpG0nJJKySdUaN8kqSLU/mNkubmys5M+cslHZ3L/2tJyyTdLelCSbUsvo8GvgTsAXwZ+Nd0fBL4dDXpq6+h3Ja2Cf8Q6A3qbPtHVS80GrSV4m93ze2z/m3bULQiZ++iBXLHhmLs6s6VuRjiTz7T+MLXF+OiMytn7b77S4plDxZjc3dNLOp431+0ds6/l+X7c8lCuWwp356Xo6v4hnfvV7z/9meLMebL21T01Jre8859i/HX25ptacnFds8/V+jv3WDWz+4tpHfYuX400olTiq/ztBuLfXeXrPDzlK3dlbOizxp31y0v99v2WDEeezkOfJmuHft+ALZ1tBcLVxX7mnZP0eFr97o+bwHle+jYtRgnvixn9+a+z4H8OcD69YVkP68T5eeV/8wnFp/d5vnFd779F7cU06/Yty/xWNH7Q1dJDkpJnix6ZVDuHS9/hu4aoXmnYeo2zQZ9gywg4UrgZklLbN+Tq3YKsNb2PpIWAl8E3iPpALKgha8EdgeulrQv8FKyNe8DbL8o6ZJU73uFW7C/D3xf0v+yfelg76GqQtkJeBZ4a14GoKFCkbQfxZjxewNnkYWX/HP63CF/urwQlNofQ7ZlrZ1ssegLFeUNgiAYecxw7vI6FFiRggsi6SLgBCCvUE4APpfOFwNfV7b99gTgItubgYckrUj9PUr2PT9F0lay6atV5QtLep/t/wbmSvpkudz2l6vcQFVL+Q9WqVej3XLgEOjVvo+TxVH5IPAV21+q17aitg6CIBhVhnGX1xzgsVx6JfC6enVsd0paTxYRdw5wQ6ntHNvXS/oSmWJ5EbjS9pU1rj0t/Z0+lBuoGg/lazWy1wNLbf+44rWOAB6w/UhFe5Yq2joIgmB0qa5Qdi7ZciyyvSiXrvXFWO69Xp2a+ZJ2JPvenEdmmP7D3Gikr6L97fT3803uoSFVtw1PJhtp3J+Og8imwU6R9G8V+1gIXJhLny7pTknnpZsuU0tbz6lRLwiCYHvgGdsLcseiUvlKYM9ceg/6T0/11klhRGYBaxq0PRJ4yPbTtreSLVO8oZ6Akv5Z0kxJEyT9QtIzkt5X9QarKpR9gLfa/nfb/56EfAXwJ8BRzRpLmggcT7aoD5kDspeTKanVZLsJ+jWrkdfvt4CkUyUtlbR0K5trNAmCIBg55GpHBW4G5kual74zFwJln4lLgJPT+YnANbad8hemXWDzyMK130Q21XWYpKlpreUI4F7qc5Tt54C3kympfYFPVZKe6ovyc8jm2Hr2WEwDdrfdJanKt/ixwK22nwTo+Qsg6T+Bn9RoU0Vbk7T8IoCZ2mn7MR8KgqA1GCY7lLQmcjpwBdlGpPNsL5N0NtnywhLgO8D5adF9DZnSIdW7hGxJoBM4LZl83ChpMXBryr+N9H1Zhx4HkMcBF9peMxCXW1UVyj8Dt0v6FdnI4U3A/ydpGnB1hfYnkZvukrRbzkfMnwB312jTq63JFvMXAn9WUd4gCIKRx0B301rVu8t2u15eyjsrd74JeHedtucA59TI/yyZg98q/I+k+8gW8P9K0i7ApiZteqm6y+s7ki4nWygX2TbfntFCw+GQpKlkO7X+Ipf9z5IOIfs4Hu4pk7Q72fbg4+pp66o3FgRBsC1oJV9ets+Q9EXguTQD9QLZon4lqu7y6pl729v22ZL2knSo7ZsqCLiRbFtbPu/9dequIhtq9aT7aesgCIIxRQsplMQryOxR8vrhv6o0rDrl9U2ygd1bgbOBDcClwGsHIGQQBEHr0UIKRdL5ZBumbge6UrYZZoXyOtuvlnQbgO21aRfCmEJtbbRP73Nz4ReLU3/tjz/bV/Zc0Q2FNxf3FpTfkba8OwkXJ03bZs0spqeV3Kvk2LLLtEJ64ovFujNvKLpi6Xb9CdrudSX3GCU5yu5C8umue39fKOqYs3ux7oRicLbys8y7+eiaXHyNJpRdr5TlyJX3u26J5xfsWUhPv6PPNUf37JLrmbbi4mGz57Pp4Lm955MfW1e68MZiX4+V3Lhs6XNz0zH3ZUU5ZkwppicWn4+2dBbTt+eiak8qujQp46eKbn80se9zKrspaZvY+L9o4dlPL/r/ayZz507F97hjTc5V0eathbJJd5bcB80ofW4PPNx7nn+ug6Gfq5Y8Tz1dv2yQDGAH1/bCAjI3LYO6q6oKZWuyXDdAWqgZxqWoIAiC7ZQWCrBFtkHqpWTmHAOmqkL5GpnLlF0lnUO2//nvB3PBIAiCVqLFRig7A/dIugn6DPtsH1+lcdVdXhdIuoVsYV7AO203Mo4JgiAYH7SWQvncUBo3VCiS8r7Cn6JoS7KT7TX9WwVBEIwTWmwNxfavJb0MmG/76mT20d6sXQ/NRii30Od4bC9gbTrfgcykf96gpA6CIGgVWkihSPpz4FQyX40vJ/OS8h9ks1NNaejLy/Y823uTGRe+w/bOtmeT+XkZ08G1giAItgmtFVP+NOBw4DkA2/cDu1ZtXNU55GvzAbBs/wx48wCEDIIgaEmG0TnkWGCz7d6928m4sbL0VRXKM5L+XtJcSS+T9BmyCI5BEATjm9YaofyQ2FokAAAgAElEQVRa0qfJIjy+jcxD/P9UbVxVoZwE7EK2dfiydH7SAAUNgiBoLSqOTrajEcoZZKHZ7yLzsXg5AzARqbpteA3w8cFIt01pE2pgHdz5eD/v931NS+26Dj2gkJ5wb1+sr62vKFpv644Him2fLlozd+y+W991rr2tUNZdum7ZUri9bFWcs+hXyaq6s3TdMvm+8jIBsKnoKeDFVxYt2Cc/XrzWi3P6rM4nritFMOgu/u/p95nkjHC94fliUcljwfTrSh4Mcs+nrfSsJm4upjtL5Srd48Rr76xbt/zcG30ufra42dGritdpZv1d+CxK98Ck0rObUPovm7Po75hdcJmHX1pMdy0rekfozv1/KLftZ6pXkqN91ZOFdP4T79pQ9EJRpvx/Lf98ys+9/O50l/ou/x/Ivz9DtbqvzPajLJpiuxv4z3QMmGbbhj9n+3NDrRMEQdCytIBCkXQXDe7E9kFV+mk2QvmwpOcayUEWp+RzVS4WBEHQSojtajqrEW9Pf09Lf89Pf98LbOxfvTbNFMp/AjMq1AmCIBh/GNQCXg1tPwIg6XDbh+eKzpD0WzIv801pqFBsf37wIgZBEIwDWmOE0sM0SW+0fR2ApDeQhXyvRFXnkEEQBEEtWkuhnAKcJ2kW2Z2tBz5UtXHVbcODQtJ+km7PHc9J+oSkf5F0n6Q7JV0maYc67R+WdFdqu3QkZQ2CIBgMw7ltWNIxkpZLWiHpjBrlkyRdnMpvlDQ3V3Zmyl8u6ehc/g6SFqfv3Hslvb7e9W3fYvtg4CDgENuH2L4119fJjeQfUYVie3kS6BDgNWSLO5cBVwEHpp0DvwfObNDNW1IfC0ZS1iAIgkExTIaNKebUN4BjgQOAkyQdUKp2CrDW9j7AV4AvprYHkG2QeiVwDPDN1B/AV4Gf294fOBho6ine9nO2a0Ura2g+UkmhSNpX0i8k3Z3SB0kaaDyUI4AHbD9i+0rbPSHgbgD2GGBfQRAEo09alK9yVOBQYIXtB5P7k4uAE0p1TgC+n84XA0dIUsq/yPZm2w8BK4BDJc0E3gR8B8D2Ftul8KQDomE0saojlP8kG0VsTULdSaYNB8JCcu7vc3wI+FmdNgaulHSLpFMHeL0gCIKRZ/hcr8wBHsulV6a8mnXSj/L1wOwGbfcms3z/rqTbJJ0rqfIiew0a3klVhTLV9k2lvM6aNWuQ4s8fT+YXJp//mdTPBXWaHm771WRDwNMkvalG36dKWipp6ZbuTf17CIIgGEEGsIayc893VTrKP5Jr/fovf4HXq1MvvwN4NfAt268CXiBzrzJYGo5Qqu7yekbSy+mLKX8iA4s5fCxwq+1efw1pceftwBG2a2o926vS36ckXUY2JLy2VGcRsAhg1uSXmh36XIJ0P/Z4ob+yy4dGtN90T1GWXNtyWVczFw+52+vY7aXFopLrkY6XFj1Fdz/xVCGdd0VRdlNSpny/+fr93GMcdnAhOfmJF4rlq58uJKds2tqX2Lq1UNbPPcb00g+i5/rKXXaP0kDmMt3PNPZPWn7WZfJ9t++0Y7Hs+eL9t0+fXrdt2cVHs/es3Ff+PgbaV/7ZlmWm9Dl0vLx6+KLy/53uZxvLVXCf0uRZlj/jwq/a8lfBxAkN5Wzm5mWbUH2X1zNN1oJXAnm/TnsAZX9RPXVWJk/As4A1DdquBFbavjHlL6aBQpE0L02Z1cv7bQP5K49QTgO+Dewv6XHgE8BfVmwLmSPJfLTHY4C/A463XdMKU9I0STN6zoGjgLsHcM0gCIKRpep0VzWlczMwX9K8NKuzEFhSqrME6NlpdSJwTfpBvgRYmHaBzQPmAzfZfgJ4TNJ+qc0RwD3U59IaeYt7b9c+vdENVHUO+SBwZPpib7Nd+WdBCiH5NjLPlT18HZgEXJWtJ3GD7Y9I2h041/ZxwEuAy1J5B/AD2z+vet0gCIKRRjSZAxoAtjslnU4W0LAdOM/2MklnA0ttLyFbXD9f0gqykcnC1HaZpEvIlEUncJrtrtT1R4ELkpJ6EPhgv/uQ9ifbITZL0rtyRTOByVXvoZlzyE/WySfdxJebXSCNQGaX8vapU3cVcFw6f5Bsi1sQBMGYZThdr6RAhpeX8s7KnW8C3l2n7TnAOTXybweamV3sR7YEsQPwjlz+BuDPq8gOzUcoPX689gNeS9/w6x2U1jKCIAjGJS1gKW/7x8CPJb3e9vWD7aeSLy9JVwKv7pnqkvQ5Sju2giAIxiUtoFByrEgRG+eS0w+2K7lfqbrLay8gv81jS7pgEATB+GX7isZYhR8DvwGuBrqa1O1HVYVyPnBT2rpr4E+A/xroxYIgCFqO1lIoU23/3WAbV93ldY6knwF/mLI+aPu2Rm2CIAjGA60QDyXHTyQdlzYHDJhKCkXSXsAzZI4de/NsPzqYiwZBELQKLTbl9XHg05I2k7naEmDbMxs3y6g65fVT+gZ2U4B5wHKyfctjBwk62vvSryzuTm57us95ptcVHWlqTtGqunP5ikK6PWfd27Zjydv+2qKvtbZddymku3bfqa+fR54slPWzBJ5S3PLdz8p8S59VuvYtWj63/b5g4Ern6w8slm/t+ynV8URRZt/9AI0oW7QX5NhhVvE6pfvvd48v5GxZ5xZdFbl0D2Wr6kJZg2cDwMyiRTovvFi3b00oyVjue/eiB4O2dX0eDlR6lzx/r2LbKaV7uO+Rkhx9124vxTLqJ1cJ555l+Xlo6tRi3SeL3g7yn+mWNxVDhk+YXXx2HeuLz451pcjguf8DzSzj+71LDehc/UQhXbbCb99S7Lvr+aLniTz9vA40djRRjepGi9sFtptF6G1I1SmvP8inJb2aoqFiEATB+KQFFIqk/W3fl77b+5GPidKIQUVstH2rpNcOpm0QBEGrIFpmyuuTwKnAv9YoM/DWKp1UXUPJW8y3kXmvfLpO9SAIgvFDCygU26emv28ZSj9VRyj5ebVOsjWVWk7EgiAIxg8GdbeARklImkDm+LcnVMivgG/b3lq3UY6qCuUe2+VYJu8mrOWDIBjntMiUVw/fAiYA30zp96e8D1dpXFWhnEl/5VErLwiCYHzRWgrltbbzTnmvkXRH1cbNvA0fS+b9d46kr+WKZjKAiI1BEAStSouNULokvdz2AwCS9mYALliajVBWAUvJwvfeksvfAPz1AAUNgiBoPVpLoXwK+KWkB8k2sb2MGvFT6tHM2/AdwB2SLrAdI5IgCII8LeYc0vYvJM0nC1ki4D7blU1Am015XWL7T4HbpP6PzfZBNZqNGt0T29m4d58l7dQ7i3Gx81bFZYvajrUly/mOooVy96a+Z9rWUXxs5TjgKlkKt23os2DufLIYI76tbM38WDGEdFvJ+rkgx9pS7PbddyukJ64oWhmTs7rufqq467t8HUpydT1ctO4md8/aVHzfyn2VLaPznga67l5OQxpYVZc/o/Z5RQv17ulFrwPdJe8HXW95Te95x2/uLJS5s7SpZc3axnLmua1RhNVBuHCtSPl5eAAyb5nZXkh3XH1XId3V0cRiv/y8Cp2XPv/Su9U2p++97by/iceGTY2/2zpyXhq61pQ8WJQ8OlByWjEYRGv58pI0Gfgr4I1kY6/fSPqPFNirKc2mvD6e/r59kMLtB1ycy9obOIvMU/HFZC7wHwb+1Ha/t1/SycDfp+Q/2f7+YOQIgiAYMdxCQ5Tsu3kD8O8pfRKZt/maUSLLtDUqtL06nf6V7UfyB5kWa4jt5bYPsX0I8BpgI5mDyTOAX9ieD/wipQtI2gn4LPA64FDgs5J2LNcLgiAYTeRqx3bCfrZPsf3LdJwK7Fu1cUOFkuNtNfKOrXqRxBHAA0kZnQD0jDa+D7yzRv2jgatsr0mjl6uAYwZ4zSAIgpHDAzgqIOkYScslrZBU64f2JEkXp/IbJc3NlZ2Z8pdLOrrUrl3SbZJ+0kSE2yQdlmv3OuC31aRvvobyl2Qjkb0l5SeZZwzkIomFwIXp/CU9ox/bqyXtWqP+HOCxXHplyguCIBgzDNcaiqR24BtkP+BXAjdLWmI7vyh3CrDW9j6SFgJfBN4j6QCy79hXArsDV0va13bPkt3HgXvJTD4a8Trgf0vqCU2yF3CvpLvI3Ng3XDdvtobyA+BnwP+lOC21wfaaJm17kTSRbOvxmVXbkK13lemn5yWdSubUjEmTd+jXIAiCYCQZxkX5Q4EVth8EkHQR2WxOXqGcAHwunS8Gvi5JKf+itCPrIUkrUn/XS9oD+GPgHDInkI0Y0ixQs23D64H1ZAszpJHEZGC6pOkDCLB1LHCr7Z59FU9K2i2NTnYDnqrRZiXwR7n0HmR+ZcoyLgIWAcyYtcf2M1MZBMH2jxnIovzOkpbm0ovS91cPtWZlXlfqo7eO7U5J64HZKf+GUtueGZ1/A/6Wok/GmqQliUFTaQ1F0jsk3Q88BPyabGfWzwZwnZPom+4CWAKcnM5PBn5co80VwFGSdkyL8UelvCAIgjHDABbln7G9IHcsKndVo/uytqpXp2a+pLcDT9m+pUb5sFN1Uf6fgMOA39ueR7bAXmkNRdJUsjnBH+WyvwC8LSmpt6U0khZIOhcgTan9I3BzOs4eyDRbEATBNmH4FuVXAnvm0nuQeSupWUdSBzALWNOg7eHA8ZIeBi4C3irpv6ve2kCpqlC22n4WaJPUZvuXwCFVGtreaHt2mj7ryXvW9hG256e/a1L+UtsfztU7z/Y+6fjuAO4rCIJgxOkJsDVM24ZvBuZLmpfWnReSzebkyc/unAhcY9spf2HaBTYPmA/cZPtM23vYnpv6u8b2+4Z63/Wo6m14naTpwLXABZKeIpxDBkEw3rGHzbAxrYmcTja13w6cZ3uZpLOBpbaXAN8Bzk+L7mvIlASp3iVkC/idwGm5HV7bDLnCw5A0DdhEppDfSzbMuiCNWsYMs6bs5jfM6/Nj5ifqB5XsWld0y9C+Q3GHWD93IXvldiw/WzLqnzqlmN74YjGdc3nC1pKLipmldbLnNjQufzHnAaHUV9e8ousV3V10Y6GJE/vOJxR/S3SVXM+0z3lp8bqle+rOubFpRtsus+v25a3F3yXl5969sXidvHuRsruP7je9qnjda28rpMufsWb1PdvORx4j6KPsxqXsTqf8/6fR51LuSxPru3HpPqRoQ9dx/8pihdL/h+7VRf8pbdP63Lp0Pv1MQzmu2nrhLbYX1BWmAjN22MOvetPHm1cEfvM/fzvk6411Ko1QbOedU4X7kyAIgsR2ZAU/4jQzbNxA7eUkkRm5NDOSCYIgaF0MtFAI4KHSzA6l6b7lIAiCcU3ok16qLsoHQRAENYgprz5CoQRBEAyF1nJfPyRCoQRBEAwWt1aAraESCiUIgmCQZIaNMULpIRRKEATBUIgRSi+hUIIgCIZAjFD6aCmF0jWlg+de2WeVPWN10St+1/o+a/B+1rszpxfTJQtu1vdZsHtLyRJ4RrGtd92pWN7Z9xOm+7GirzeV+vI+exbS3PtgsX7O6r7rFXMLZR1PFa3du+cXy5+f17cLfMatjxf7LVm+dz3+RFGukvVznvZZs4p9zW4SqTn3vDxjaqGo7ami/8/8/ZbZ/Nr5hfSqj2wupPe+v2jt37m6eE+UrL2DPsqfd9kyfiC071h8Pyj9fyHntcE3LSsUdZbkaCu9p22ziqZwXU/Xd97R6B0eNAOIxjgeaCmFEgRBsG0xCsPGXkKhBEEQDIWY8uolFEoQBMFgiW3DBUKhBEEQDIUYofQSCiUIgmAohD7pZcQViqQdgHOBA8ke/YeATwD7pSo7AOts94sAmcJWbgC6gM5WjyUQBMH2R2wb7mNbjFC+Cvzc9okprOVU2+/pKZT0r8D6uq3hLbafaVAeBEEwOhjoCoXSw4gqFEkzgTcBHwCwvQXYkisX8KfAW0dSjiAIgpFAOEYoOdpGuP+9gaeB70q6TdK5KZxwD38IPGn7/jrtDVwp6RZJp46wrEEQBAOnJ658s6MCko6RtFzSCkln1CifJOniVH6jpLm5sjNT/nJJR6e8PSX9UtK9kpZJqhaveJCM9JRXB/Bq4KO2b5T0VeAM4B9S+UnAhQ3aH257laRdgask3Wf72nyFpGhOBZisacy4+t7esrxlfJmy1Wznoyvr1EzXyVnWt5ctwUvWu2VDp+4n+yz2u18sxZsv0XbPlobledpfKFqGdz1atH4v3+PUZX33ULZALlu7l/8D5GN1Q+nZ7v6SYtNVxTjf2nXnQrp7ZZ+3AD9a9EjQVZSi3z2sOeUNvecTXyju13zZwjsK6fI9BiNHIyv0cmx3yukB0L1xY8P0qDBMIxRJ7cA3gLcBK4GbJS2xfU+u2inAWtv7SFoIfBF4j6QDgIXAK4Hdgasl7Qt0Av/H9q2SZgC3SLqq1OewMdIjlJXASts3pvRiMgWDpA7gXcDF9RrbXpX+PgVcBhxao84i2wtsL5jYNmWYxQ+CIGiAyZxDVjmacyiwwvaDaXngIuCEUp0TgO+n88XAEWnp4ATgItubbT8ErAAOtb3a9q0AtjcA9wJzBnezzRlRhWL7CeAxST07uo4AejTjkcB9tmsODSRNSxqVNE12FHD3SMobBEEwUGRXOiowB3gsl15J/y//3jq2O8k2NM2u0jZNj70KuJERYlvs8voocEHa4fUg8MGUv5DSdJek3YFzbR8HvAS4LFO+dAA/sP3zbSBvEARBRQzdlU3ld5a0NJdeZHtRLq3aFyhQr07DtpKmA5cCn7D9XEV5B8yIKxTbtwP97Edsf6BG3irguHT+IHDwSMsXBEEwaMxA1lCeaWJLtxLIuxvfA1hVp87KtGwwC1jTqK2kCWTK5ALbP6oq7GAY6TWUIAiC1mb41lBuBuZLmpdmdBYCS0p1lgAnp/MTgWtsO+UvTLvA5gHzgZvS+sp3gHttf3nQ91iRcL0SBEEwBIbLDsV2p6TTgSuAduA828sknQ0stb2ETDmcL2kF2chkYWq7TNIlZGvUncBptrskvRF4P3CXpNvTpT5t+/JhEbpEKJQgCIKhMIyGjemL/vJS3lm5803Au+u0PQc4p5R3HbXXV0aEUChBEASDxYau8F/fQyiUIAiCoRCuV3oJhRIEQTAUQqH00lIKxV1dDd2tDIS8qxUoulvZsv8exbLf3lVMz5lU7GvK5L6y3DkAW4uuR9xZTJfxlj7XLL73gYZ1+7XNucco3x/txQ1/bS8r3mP3I0X707YpfV4J/PBjNKL7oUcry9g2c3oh7T13K6R3/dXq3vPOBx4q1q18lSAYJgxETPleWkqhBEEQbFsMjjWUHkKhBEEQDIWY8uolFEoQBMFgMbHLK0colCAIgqEQI5ReQqEEQRAMmurBs8YDoVCCIAgGixmIt+GWJxRKEATBUIgRSi+hUIIgCIZCKJReQqEEQRAMFht3dY22FGOGUCh1yFuVA3Q++VTv+YSNLxbKukp1eXFTIaldZveedz/6eKGs+8ViX/0s2JvIlSdvvQ6gCcWPVx196a4NGwplXWvWFvsq3UPbjKIFe9eadZVkqiVX2159kUn99LMN5aCUjtnqYMwRlvK9jHiALUk7SFos6T5J90p6vaTPSXpc0u3pOK5O22MkLZe0QtIZIy1rEATBgLGrHeOAbTFC+Srwc9snpihkU4Gjga/Y/lK9RpLagW8AbyMLb3mzpCW279kGMgdBEDTHA4op3/KM6AhF0kzgTWRRxrC9xfa6xq16ORRYYftB21uAi4ATRkbSIAiCQRIjlF5Gesprb+Bp4LuSbpN0rqRpqex0SXdKOk/SjjXazgHybmxXprwgCIIxQrYoX+UYD4y0QukAXg18y/argBeAM4BvAS8HDgFWA/9ao22tsJX91LykUyUtlbR0K5uHTfAgCIKm9Livr3KMA0ZaoawEVtq+MaUXA6+2/aTtLtvdwH+STW/VartnLr0HsKpcyfYi2wtsL5jApHJxEATByOLuakcFmm1EkjRJ0sWp/EZJc3NlZ6b85ZKOrtrncDKiCsX2E8BjkvZLWUcA90jKR036E+DuGs1vBuZLmpcW8xcCS0ZS3iAIgoFgwN2udDQjtxHpWOAA4CRJB5SqnQKstb0P8BXgi6ntAWTfka8EjgG+Kam9Yp/DxrbY5fVR4IKkFB4EPgh8TdIhZJ/Hw8BfAEjaHTjX9nG2OyWdDlwBtAPn2V62DeQNgiCohoc1wFbvRiQAST0bkfI7W08APpfOFwNfl6SUf5HtzcBDklbQN/PTrM9hY8QViu3bgQWl7PfXqbsKOC6Xvhy4fOSkC4IgGBpVRh8VqbUR6XX16qQf3euB2Sn/hlLbnk1MzfocNlrKUn4Da5+52osfAXYGnhmxCz3XpPzJuunmcjU2Om/MxkG37C9Xua/B992/7X2VW47s5zg4xqJMEHINhB6ZXjbUjjaw9oqruy/ZuWL1yZKW5tKLbC/KpatsRKpXp15+rWWNEdsh0FIKxfYuAJKW2i6PikadkGtgjEW5xqJMEHINhOGUyfYxw9FPospGpJ46KyV1ALOANU3aNt3cNFyMuOuVIAiCoBJVNiItAU5O5ycC19h2yl+YdoHNA+YDN1Xsc9hoqRFKEATB9kq9jUiSzgaW2l5C5nXk/LTovoZMQZDqXUK22N4JnGa7C2Bbbm5qVYWyqHmVUSHkGhhjUa6xKBOEXANhLMoE1N6IZPus3Pkm4N112p4DnFOlz5FCHic+ZoIgCIKRJdZQgiAIgmGhpRSKpH9JcVfulHSZpB1yZTXdEmwjud4taZmkbkkLSmWjKdeYiDeTHIQ+JenuXN5Okq6SdH/6W8uB6EjLtaekX6Y4PsskfXy0ZZM0WdJNku5IMn0+5c9LrjjuT645Jm4rmUrytSdHsD8ZK3JJeljSXcpiLy1NeaP+frUiLaVQgKuAA20fBPweOBPquyXYhnLdDbwLuDafOZpybWuXDE34Htn95zkD+IXt+cAvUnpb0wn8H9uvAA4DTkvPaDRl2wy81fbBZM5Vj5F0GJkLjq8kmdaSuegYDT4O3JtLjxW53mL7kNx24bHwfrUcLaVQbF9puzMlbyDbcw05twS2HwLybgm2hVz32l5eo2g05Roz8WZsX0u2YyXPCcD30/n3gXduU6EA26tt35rON5B9Uc4ZTdmc8XxKTkiHgbeSueLY5jL1IGkP4I+Bc1NaY0GuOoz6+9WKtJRCKfEh4GfpfKzGVhlNucbqM+nhJbZXQ/bFDuw6msIo8+r6KuBGRlm2NK10O/AU2aj8AWBd7sfUaH2W/wb8LdDj3Gr2GJHLwJWSbpF0asobU+9Xq7DdbRuWdDXw0hpFn7H941TnM2TTFRf0NKtRf1i3t1WRq1azGnnbatvdaF57u0LSdOBS4BO2n8t+eI8eyb7gkLRGeBnwilrVtqVMkt4OPGX7Fkl/1JNdo+povGOH214laVfgKknVnf8EA2K7Uyi2j2xULulk4O3AEe7bE10ptspIylWHEZdrjF67Ck9K2s32amXhDp4aDSEkTSBTJhfY/tFYks32Okm/Ilvf2UFSRxoNjMZneThwvKTjgMnATLIRy2jL1eN0FttPSbqMbLp3THyGrUZLTXlJOgb4O+B423mXhPXcEow2oynXWI83k3cxcTJQb5Q3YqQ1gO8A99r+8liQTdIuPbsXJU0BjiRb2/klmSuObS4TgO0zbe9hey7Zu3SN7feOtlySpkma0XMOHEW2SWbU36+WxHbLHGSL2o8Bt6fjP3JlnyGba14OHLuN5foTshHBZjLfw1eMEbmOI9sN9wDZ1NxofW4XkoWC3pqe0ylk8++/AO5Pf3caBbneSDZFc2funTpuNGUDDgJuSzLdDZyV8vcm+zGyAvghMGkUP88/An4yFuRK178jHct63vOx8H614hGW8kEQBMGw0FJTXkEQBMHoEQolCIIgGBZCoQRBEATDQiiUIAiCYFgIhRIEQRAMC6FQgiAIgmEhFMo4Q9LzzWsNuM/je9zfS3rnYLwWS/pV2bV/hfrLJR1fo2xu3hV+qyPpA5J2z6UvkLRG0omN2gXBcBMKJRgytpfY/kJKvpPMHf624L3O4myPGNs4zMFg+QDQq1CcWaiPJa8HwTghFMo4RRn/IunuFHzoPSn/j9Kv/8XKgpVdkFyQIOm4lHedpK/lgih9QNLXJb0BOB74lxTM6OX5kYeknSU9nM6nSLpIWTC0i4EpOdmOknS9pFsl/TA5Z2x2P69RFnTqeuC0XH57us+b07X+IuW3SfqmsiBVP5F0ec8vemUBmc6SdB3w7nQfP0/ean8jaf9UbxdJl6a+b5Z0eMp/c7r/25UFm5rRQO5P5WT7fC7//6XrLVPykJvu5Xu5z+yvk8wLgAvS9abUu1YQjDTbnXPIYNh4F1mApoOBnYGbJfUEAHsVWdCvVcBvgcOVRbr7NvAm2w9JurDcoe3fSVpC5nZjMYDqe+b9S2Cj7YMkHQTcmurvDPw9cKTtFyT9HfBJ4Owm9/Nd4KO2fy3pX3L5pwDrbb9W0iTgt5KuBF4DzAX+gMx1+b3Aebl2m2y/Mcn0C+Ajtu+X9Drgm2RxPr5KFjzqOkl7AVeQef79G+A0279NynBTLYElHUXmv+1QMs+8SyS9yVl8mA/ZXpMUxM2SLk3yzrF9YGq/gzMHkacDf2N7aZNnFAQjSiiU8csbgQuduUJ/UtKvgdcCzwE32V4JoCzuxlzgeeBBZ4HAIPO/dWq/XqvzJuBrALbvlHRnyj+MbMrst0kZTQSub9SRpFnADrZ/nbLOJ4tECZkzwINy6wmzyL7E3wj80HY38ISkX5a6vTj1PR14A/DDnHKclP4eCRyQy5+ZRiO/Bb4s6QLgRz3PsgZHpeO2lJ6eZLsW+JikP0n5e6b85cDekv4d+ClwZaPnEgTbmlAo45dGQT025867yN6TwQYB6aRvanVyqayWIzkBV9k+aQDXUJ2+eso+avuKQqb0x036fCH9bSMLEnVIjTptwOttv1jK/4Kkn5I5krxB0pG2a8caEvEAAAI6SURBVMXgEPB/bX+7JNsfkSmr19veqMxF/WTbayUdDBxNNq33p2SB5IJgTBBrKOOXa4H3pHn5XchGDI1c599H9ut4bkq/p069DUB+zeBhsukl6HNj3nP99wJIOpDMiy5koZsPl7RPKpsqad9GN2J7HbBe0htT1ntzxVcAf6ksrgmS9lXmxvw64H+ltZSXkHnIrdX3c8BDkt6d2it9qUM2Qji9p66kQ9Lfl9u+y/YXgaXA/nVEvwL4UM8akaQ5yoJAzQLWJmWyP9morWc6sM32pcA/AK9O/ZSfeRCMCqFQxi+XkblAvwO4Bvhb20/Uq5x+hf8V8PO0WP0ksL5G1YuAT6XF6JcDXyL7Qv8d2VpND98Cpqeprr8lKTPbT5PtWrowld1A/S/kPB8EvpEW5fMjhnOBe4BblW0l/jbZiOtSMlf5PXk31rkfyBTUKZJ6XKCfkPI/BixIC+r3AB9J+Z9IC+d3JFl+1q/H7F6vBH4AXC/pLrLY6zOAnwMd6f7/MT0DyMLn/ipNQ34PODPlfw/4j1iUD0abcF8fVEbSdNvPK1s0+AZwv+2vjJIsv2KIC9G5+5lNptAOb6RUtyckfY/c5ogg2BbECCUYCH+efh0vI5uW+XaT+iPJGuB7qmHYOAB+ku7nN8A/tpAyuQB4M3V2lwXBSBEjlCAYYST9AdnOszybbb9uNOQJgpEiFEoQBEEwLMSUVxAEQTAshEIJgiAIhoVQKEEQBMGwEAolCIIgGBZCoQRBEATDwv8PmBfOdk0McHYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "arr.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, these data are not very useful without weighting it by grid cell areas, but it can be retrieved by calling the function again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "dens_nonweighted = tr.density(lon_dens1d, lat_dens1d, weight_by_area=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAELCAYAAADZW/HeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztvXmcHWWV///+dHf2hRASMIQlYRlQWQKERcM4bCKgBhdQGHEQcaLjhjqjgs5PFPX3xRWXUSQCwmhkkWXki7KJLKJsCQQChAxLAoSELEIWSEjS3ef7x1M3XfX0Xep239t9+/Z551Wv1LPWqbq376lnOefIzHAcx3Gclv4WwHEcx2kMXCE4juM4gCsEx3EcJ8EVguM4jgO4QnAcx3ESXCE4juM4gCsEx3GchkFSq6SHJd2YpC+TtFjS/OSYVs/rt9Wzc8dxHKcqzgIWAmNTeV80s2v64uI+QnAcx2kAJO0EvBO4uL9kaKoRwlANs+GM6m8xHCeDWlszaevo6Jfr0tmZlaOPvBT01/1XYj2vrDazib3p4x1HjrK/v5zvfuY9uulx4PVU1mwzm51K/wj4EjAmavptSV8DbgfONrNNvRC5LE2lEIYzikN1dH+L4TgZWseOy6Q71qzpl+vaxtcz6c5N2XRfydFX91+JP9k1z/W2j7+/3MEDt+ySq27rpKdeN7PpxcokvQtYaWbzJB2RKjoHeAkYCswGvgyc1yuhy1D3KSNJn5f0uKTHJF0haXjehRJJp0t6KjlOr7esjuM41WBAZ85/FZgBzJS0BLgSOErSb8xsuQU2Ab8CDqnn/dR1hCBpMvBZ4E1mtlHS1cApSXHZhRJJ44FzgemE5z5P0g1m9ko9ZXYcx8mLYWyx3k+Bmdk5hNEAyQjhP8zsNEmTzGy5JAHvAR7r9cXK0BdTRm3ACElbgJHAspzt3gHcZmYvA0i6DTgOuKIuUjqO4/SAHG//vWGOpImAgPnAJ+p5sbpOGZnZi8D3geeB5cBaM7s1Kf62pEclXSBpWJHmk4EXUumlSZ7jOE5DYBgdlu/I3afZnWb2ruT8KDPb18z2MbPTzOzVut0MdVYIkrYFTgSmAjsCoySdRhga7Q0cDIwnLJR0a14kr9tTlTRL0lxJc7dQt8V3x3GconRiuY6BQL0XlY8BFpvZKjPbAlwHvDXnQslSYOdUeieKTDeZ2Wwzm25m04dQbKDhOI5THwzowHIdA4F6K4TngcMkjUwWRY4GFkqaBFBhoeQW4FhJ2yYjjWOTPMdxnIahmUYIdV1UNrP7JV0DPAS0Aw8T9tLeVGyhRNJ04BNm9jEze1nSN4EHk+7OKywwO47jNAIGbGmiMMR132VkZucSto+mOapE3bnAx1LpS4FL6yed4zhOz7EBNB2Uh6ayVB4ItAwbnkn3lbWoU55qP5d0/bhu67jGsMxtFIvgmKb6GzDoaB594ArBcRynpwRL5ebBFYLjOE6PER1Fd8gPTFwhOI7j9BADOn3KyHEcxzFgcxOFlXGF4DiO0ws6zaeMHMdxBj3BUtkVguM4zqDHEB0+ZeQ4juOATxk5juM4hBHCZmutXHGA0DxjHcdxnD4mGKa15DryIKlV0sOSbkzSUyXdn4QRvkrS0HreT1ONEDR0CG07dnnMtrXrM+VxkPE0lczn2ya9oSsxdnSmrH3R07lljK8Tm/FX275cXxqRTZdzZVDJnUAjuhvo7bMr536i2r7SVHIZUe5ZNuJzzkO5zyL+uxso95SXGi8qnwUsBMYm6e8AF5jZlZJ+AZwJXFjLC6bxEYLjOE4PMRMd1pLrqISknYB3AhcnaREcgRZiz19OCBdQN5pqhOA4jtPXdNZuhPAj4EvAmCS9HbDGzNqTdN3DCPsIwXEcp4cEO4SWXAcwoRDuNzlmFfqR9C5gpZnNS3WfK4xwLfERguM4Tg8xxBbL/TO62symlyibAcyUdAIwnLCG8CNgnKS2ZJRQNIxwLfERguM4Ti/oMOU6ymFm55jZTmY2BTgF+LOZfQi4AzgpqXY68Pt63kvdFYKkz0t6XNJjkq6QNFzSHEmLkrxLJQ0p0bZD0vzkuKHesjqO41RDwVI555RRT/gy8AVJTxPWFC6pmfBFqOuUkaTJwGeBN5nZRklXE7TfHOC0pNpvCWEzi22l2mhm0+opo+M4Tm/ozLGDqBrM7E7gzuT8WeCQml6gDH2xhtAGjJC0BRgJLDOzWwuFkh4gzI05juMMKAqLys1CXe/EzF4Evg88DywH1kbKYAjwYeDmEl0MT1bj75NU1/23juM41WLkWz+otIbQKNR7ymhb4ERgKrAG+J2k08zsN0mVnwN3m9lfSnSxi5ktk7Qb8GdJC8zsmegas4BZAMMZSftzL9REdjs8mql64e9bTzvGjMgUVWMxW431cLV0swCtwiK0kvVobyysq2nbm7rV0ogWs40iU6XPN5azUeTua8yoZpdRw1Pvsc4xwGIzW2VmW4DrgLcCSDoXmAh8oVRjM1uW/P8sYU7tgCJ1ZpvZdDObPoRhtb8Dx3GckojOnMdAoN4K4XngMEkjEzPso4GFkj4GvAM41cw6izWUtK2kYcn5BMI+3SfqLK/jOE5uDGrmuqIRqOtYx8zul3QN8BDQDjwMzAZeA54D7g16guvM7DxJ04FPmNnHgDcCF0nqJCiu883MFYLjOA1FMy0q133yy8zOBc7Nc10zm0vYgoqZ/Q3Yt77SOY7j9BxDHiDHcRzHCVNGzbSo3Dx34jiO0+eo1vEQ+hVXCI7jOD3EqL2lcn/iCsFxHKcX+AjBcRzHwUw+QhgM6J75mXR7OhFZQ1cVsaKOFp0dRx+USQ+dlzHqzlhJty9/qaq+y8VrrhQztzdxgmtpAVvuHmppMV5tbOtyVsFx25hycrfttUcm3ZvY371hoMaJzstAsTHIgysEx3GcHhIC5LT2txg1wxWC4zhODwmLyr6G4DiO49BclsrNcyeO4zh9TMFSOc9RiSSa5AOSHkmiTH4jyb9M0uJU9Mi6BQ3zEYLjOE4v6Kzde/Um4CgzezWJFXOPpJuSsi+a2TW1ulApXCE4juP0EDNqFvzGzAx4NUkOSY6qNjH2Fp8ychzH6SGGaO9szXUAE5IIkIVjVtyfpFZJ84GVwG1mdn9S9G1Jj0q6oBAWoB74CMFxHKcXVGGpvNrMpperYGYdwDRJ44DrJe0DnAO8BAwlhA/4MnBezyUuTdkRgqR1FY71kv63HoI5juM0OoVtp7VYVM70a7aGECXyODNbboFNwK+AQ2p+IwmVpoyeMbOxZY4xhGA3juM4g5DguiLPUbEnaWIyMkDSCEII4iclTUryBLwHeKxed1Npyuj9OfrIU8epEeXcHLTc83gmHa9Gde6x89bz1sjdhHaYkK287lXK0flyl8uESq4IWsaPy6Q3HLLr1vMRL27IlNncBWX7atu16x6eP2XnTNn4Jzsy6ZE3P5JJV3IhUSu6PY8Kz6fc82uLn93+2XsefkfX84r7iV1VpJ8dQOdLq3LJUIzW/d6YSduixSX70l5Ts3KseqVkv9W6VGkEahgveRJwuaRWwsv61WZ2o6Q/S5oICJgPfKJWF4wpqxCS4PZlqVRH0ucJUdAMWACcQbjxK4HxhPCaHzazzUXangOcCXQAnzWzWyrJ4ziO01eYwZbO2riuMLNHgQOK5B9VkwvkINcuI0nvk/SUpLWptYN1OdpNBj4LTDezfYBW4BTgO8AFZrYn8ArhRz9u+6ak7puB44CfJ5rTcRynIailYVojkHfb6XeBmWa2TWHtwMzG5mzbBoyQ1AaMBJYDRwEFI4vLCfNiMScCV5rZJjNbDDxNHRdTHMdxekInynUMBPIqhBVmtrDazs3sReD7wPMERbAWmAesMbOCR+mlwOQizScDaT/TRetJmlXY17uFTdWK6DiO02Pqtcuovyi7hiDpfcnpXElXAf8DXb+6ZnZdhfbbEt70pwJrgN8BxxepWswar9gT7FbPzGYT9uYyVuP71KrPcRxnMAXIeXfqfANwbCptQFmFQNg2tdjMVgFIug54KzBOUlsyStgJWFak7VIgvS2iVD3HcZz+YQC9/eeh0i6jMwAkzTCzv6bLJM3I0f/zwGGSRgIbgaOBucAdwEmEnUanA78v0vYG4LeSfgjsCOwJPJDjmo7jOH2CAe1NNELIeyc/zZmXIfHDcQ1ha+mC5HoF0+svSHoa2A64BEDSTEnnJW0fB64GngBuBj6VmHU7juM0BINtDeEthCmeiZK+kCoaS9hCWhEzOxc4N8p+liI7hszsBsLIoJD+NvDtPNdxHMfpDwbKj30eKq0hDAVGJ/XGpPLXEaZ8Biyt47IWoPWyWu0tsZxpa+OWp1/IlMWWuN2sPlNWwPFQqy1qm7ZEhu7Wp2m5WuO2e2QtYi2S87U3dL1LvPvb92bK5lx0bCa9460rM+m09e2O38n2G9MZZwzA4O7xZzgisjDvKHNP3azah5d2ktk26Q1lrxv31fFo/k2H1dQdaBTsEJqFSmsId0m6B9jXzL7RRzI5juMMGAaKjUEeKrq/NrMOSeP7QhjHcZwBhQ2uKaMCD0u6gWBHsNW7aSU7BMdxnGbGgPbO5tlllFchjAf+TnA5USCPHYLjOE7TMqjWEAoU7BEcx3GcLDbYFIKknQh2BzMII4N7gLPMbGkdZXMcx2l4GmVROZnWr8TLZvaRUoV5p4x+BfwWODlJn5bkvT1ne8dxnKbDarioLGk4cDcwjPDbfI2ZnStpKjnixwBvJMSeKXkJ4GflZMirECaa2a9S6cskfS5nW8dxnCZFdNRuUXkTcJSZvSppCHCPpJuALxDix1wp6ReE+DEXFmn/VTO7q6y0Ulnzgbx3slrSaZJak+M0wiKz4zjOoMZMuY7K/ZiZWSF27ZDkMPLFj8HMrs5xjbJ18o4QPgr8F3BBIuDfkryGpm2vPTLptJWrRRaflSw108RWm7GFcNx32so3tjyO0TZjMun25yJr3LldFsTdHDtVYWzdzYp17OhM0iZPzNZf8FS2PHWPsRVz24uRxXRkBb7dRX/ben77ZdlYzuOO3JJJrzhi+0x6h1SsZ5u4baZsoFrEpmMdp+McQ/dnG1vUZ2IbP7882/EukzLJ9jLPp2XokKxM0d9DbLle7m8rJu5ry+7ZdNu6rjgmWpuN5d3t+99gFHwZ5WSCpLmp9OzEff9WkqiQ84A9CNM7z5AvfkxZJM02s1mV6uXdZfQ8MLNaIRzHcZoaC+sIOVltZtPLdhcceE6TNA64nrAuUOSq3SljQCzghDwC5t1lNBH4V2BKuo2ZNfwowXEcp57UY5eRma2RdCdwGPnixwCsAp4jG1zMkvT2RVtE5J0y+j3wF+BPFJmpcBzHGYwYtbNDSF68tyTKYAQhwNh3yBc/BoIX6aOTGZ2471xzb3kVwkgz+3LOuo7jOIME0dFZsxHCJODyZB2hBbjazG6U9ARwpaRvAQ+TxI8pwo+AbQmByWK+m0eAvArhRkknmNkfc9Z3HMcZFNRqhGBmjwIHFMkvGj+mSL2SNgZmVjGgGeTfdnoWQSlslLRO0npJ63K2dRzHaUrMarfttLdIOrC3dfLuMhpTrlzSm5OQl3H+XsBVqazdgK8BbwH2SvLGEbZVTSvSfgmwnrBu0V5phd5xHKevaSDndr+SdASUXeW+hCKjkAJ5p4wq8Wugm+Yxs0XANNi6v/ZF4Hoz+1GhjqQfAGvL9H2kma2ukZyO4zg1pYptp/VmG4INQzmFsKpMWc0UQh4VeTTwjJk9t7WRJOADZN1qO47jDAgM0dkg8RDMbEpv+6jVneTRkacAV0R5/wisMLOnitQv9HurpHmSilrZSZolaa6kuVvYVKyK4zhO3bCcx0CgViOEskgaSrB0PicqOpXuSiLNDDNbJml74DZJT5rZ3ekKien3bICxGp957uXM6WOXAJ2Rq4q0OwHIuhSI21bSqt3cRKSI3VzEbjDK9WX77pkpa12/MVv59ayC7JjY5TajY3hrpqxlVdZlwPJ/zC4bjZ+wbyY9YllX/ZZFizNl3dyCRM8yHew9/oyG37Egkx5y8N6Z9Ib9U24ehmWf/MaD3pJJT7y5vFxpNxCxjKuO2olyTPxzl+f32N1Ey/ise5LY7UOMrV1fsix2+xC7GEm7o4i/Zy2RG4jYbUr6u7bhjdnrjFz4Usm6AJ1LSnu+r+QWZcOkbPk26e9e6rsBkWsOuru2YHPW1UlHyuVK64vZz+XVQ3bNtv2fa+g11lzxEGo1QijmijXN8cBDZraikCGpDXgf2UXnDGa2LPl/JcGMu+LWK8dxnD6liYYIuRSCpBmSRiXnp0n6oaSt6tbMDqvQRbGRwDHAk6WC7EgaJWlM4Rw4Fngsj7yO4zh9RaNsOy0g6VpJ75RU9Qt/3gYXAhsk7Q98ieAv479zCjeSEEgnjr/cbU1B0o6SCsZvOxD8gT8CPAD8wcxuzimv4zhOn2CW7+hDLgT+GXhK0vmS9q7UoEDeNYR2MzNJJwI/NrNLJJ2ep6GZbQC2K5L/kSJ5y0i88iXWefvnlM9xHKfPMQNrkF1GBczsT8CfJG1DmJ25LfFl9EvgN2a2pVTbvHeyXtI5wIeBPyQ2BUMqtHEcx2l6GnCEgKTtgI8QQmo+DPyYYCt2W7l2eRXCBwnh3T5qZi8RAjR8r6fCOo7jNA0Ntqgs6TqCd+qRwLvNbKaZXWVmnwFGl2ub13XFS5KuBQr7HFcTdv04juMMYvp2wTgnF8eOSCUNM7NNldz/5N1l9K+EmJ4XJVmTgf/piaSO4zhNRYONEIBvFcm7N0/DvIvKnyLYANwPYGZPJcZijuM4g5cGMkyT9AbCy/oISQfQ5VJoLGH6qCJ5FcImM9scXA9tNSpreFOLshbCkZWvoiDy5YJ7x1atsaVp+z5TsvUfW1JGyiyxNeWoxdlA8tqU2iDwYlk/Vd0scy11T22RBaityPoP3CmOfLEuayGaDnAfP8uWSK74+ShljVrpWa7fdUQm/fr4rj++He6NLXyz3/mXj5qaSY+/J2vy0payvo2vO+H6bEB67TChpJyxFW/aWhZAFSyV030rfs7RZ/ha9P0YHdVPUylAfcvGLrlH3Bd9V3aZlK28TXbquWXVK9n00NQek8h6OLaYH7uotC9LGxbtVXl+eTYdPestu2ctrFsffLKrr6juqMV18thfu4hpOxO2878B6ARmm9mPJX2dEMK48If1lRKxad5BWEjeCfhhKn898JU8MuRVCHdJ+gpB87wd+CTwf3O2dRzHaV5q92rcDvy7mT2UGOXOk1TYFXSBmX2/rBhmlxMirr3fzK7tiQB5FcLZwJnAAuDjwB+Bi3tyQcdxnKaiRgrBzJYDy5Pz9ZIWEqaAciHpNDP7DTBF0heK9P/DIs0y5N1l1CnpN8DdSYwDx3Ecx6hmymiCpLmp9OzEOWc3JE0hBLK5H5gBfFrSvwBzCaOIV4o0G5X8X3ZraTlyKQRJMwl2B0OBqZKmAeeZ2cyeXthxHKcZqMLobHWeqI+SRgPXAp8zs3WSLgS+SVA/3wR+AHy0uxx2UfL/N3JLFJHXMO1cwi6jNckF5wNTenpRx3GcpqFT+Y4cSBpCUAZzzOw6ADNbYWYdZtZJcD9R1uuzpO9KGitpiKTbJa2WdFqe6+dVCO1mVi7MpeM4zqBElu+o2E/YxnkJsDA93y8pveXrvVT2+nysma0D3gUsBf4B+GKee8m7qPyYpH8GWiXtCXwW+FvOto7jOM1JbY3OZhD8xS2QND/J+wpwajJNb8ASwsaechT27p4AXGFmLxdMBiqRVyF8BvgqwZ/Rb4FbKG4N5ziOM4hQzewQzOweisenL2ZzUI7/K+lJYCPwSUkTgdcrtAFyKITEs+k3zOyLBKXgOI7jFGgwE10zO1vSd4B1ZtYh6TXgxDxtKyqEpMODeitkfxDHtk1bNbaty1pPslfWqrUtit2atvqMY+jG6J752Yx0LNvYAjSKR7x5dHZZZ2Qcrzh1Ht9fN8vkKB3H1E0TW1cPWVXaAhZAKUvV9jdlt0q3vli2acZyNbamja3Lt/1DeYvhNENeHZpJj34gK0i52MaxtXGl2Nbp8rjfOJYv40s/d4AtE7t2CbZVsNQd/eTfM+n0tWOZ43jMsZyZ2OBRnHDWZOt2s/qPn1c6MTR7D3G85qE3P5hJZ+ImR/fbEclBlFQUC12p73j8POr2u91gCiHhjQR7hPRvfMWgZnmnjB6WdAPwO+C1QmZhFdxxHGdQYuTeQdRXSPo1sDswH+hIso0aKoTxwN+Bo1J5RvewmLFgewFXpbJ2A74GjCOHbw5JxxECO7QSXLqen1Nex3GcPiHPDqI+ZjrwJrPqw/LktVQ+o2qRQrtFwDTYuhbxIiGOwhlU8M2R1P8ZIR7zUuBBSTeY2RM9kcVxHKcuNJ5CeIzgIG95pYoxeS2Vf1Ikey0w18x+n/NaRwPPmNlzObdAHQI8ncRWRtKVhIURVwiO4zilmQA8IekBws5QAPJ4lsg7ZTQc2JuwhgDwfuBx4ExJR5rZ53L0cQpwRSpdyTfHZCC94rgUODTuVNIsYFYQMpfLb8dxnJrRgFNGX+9pw7yWynsAR5nZT83sp8AxhFXs9wLHVmosaSgwky6FciFh0WMaYVjzg2LNiuR1e/RmNtvMppvZ9CEMK9LEcRynjpjyHX0ljtldBAO2Icn5g8BDedrmVQiT6fKkR3K+o5l1kBqSlOF44CEzW5EInMc3x1IgHT1lJ2BZTnkdx3HqjxFC2eQ5+ojehDzOO2X0XWC+pDsJb+5vA/5/SaOAP+Vofyqp6SJJkxLf31DaN8eDwJ6SphIWo08B/jmnvI7jOH1CA04Z9Tjkcd5dRpdI+mNyERG2iRbe1ss6TZI0krBTKO1/47vFfHNI2pGwvfQEM2uX9GmCm4xW4FIzezyPvI7jOH1G4ymEHoc8zrvLSIRdQruZ2XmSdpF0iJk9UKmtmW0AtovyPlyi7jKCQ6ZC+o9U78fDcRyn72g8hdDjkMd5p4x+TpgFOwo4jxC0+Vrg4OplrR9qbaV1bJfpeuy6YcgLXWb/cVD1bm4for67me6nyyLXBC1jSwcs2rBjtmzkpp0y6TgQfLmpx9gVQSxHOToezbqEGBIFu2d4doE+fj5ptwBbRme/RtUs7bfF1414+fDo+cztcguRdvkA0DEkuyRW6fm8mgpY3y0Ae+S6pHNJ9LmkXD207bVHVo4xI7J1h7dm5Xi9I5NOB4aPXULE2IrVmXT6c4jdPLRsLN9X+tnbNtln2Tk0+5lqc3smvXlidkff0FUbuupu2pIpG/lI1j2JRS5ULOWepTN2oVElmWcQeyqJ3FzUgryurfuYHoc8zqsQDjWzAyU9DGBmryQ7hxzHcQY3Dea6IrVZ55fVts2rELYklsMGkLhT7cN1c8dxnMakUUYIkhZQZgLLzPar1EdehfATgsuJ7SV9GzgJ+M+cbR3HcZqXGikESTsTHNC9gfDCPdvMfixpPMEn3BTCJpwPFDHkhRAhDcIuI4BfJ/9/CNjQvXp38u4ymiNpHmFhWcB7zGxhhWaO4zjNTW3XENoJXhsekjQGmCfpNuAjwO1mdr6kswlrBF/uJorZcwCSZpjZjFTR2ZL+Slj/LUtZhZBopgIrydoSjDezlytdwHEcp6mpkUJIbLOWJ+frJS0kGJWdCByRVLscuJMiCiHFKEmHJxHYkPRWsobFJak0QphHuF0BuwCvJOfjgOeBqaWbOo7jDALyK4QJkuam0rPNbHaxipKmAAcQjMt2KBjymtnyHEZmZwKXStomkW4t8NE8ApZVCGY2NRHuF8ANhZgFko4n+DNyHMcZ1FQxZbTazKZX7E8aTdjW/zkzW5fTO/RWzGwesL+ksYDMbG3U/+lmdnmxtnl9GR2cDmBjZjcB/1SVlI7jOM2I5TxyIGkIQRnMSUWkXCFpUlI+iTB9X1kss3WxMkg4q1SbvAphtaT/lDRF0q6SvkqIoOY4jjN4sS7jtEpHJRKPEJcAC83sh6miG4DTk/PTgbwxaEpeqlRB3m2npwLnEraeGnB3ktdYtKhscPQ4oHumaWSJ/PqR+2bSaWvLDftnrWtH3PdUJt0RWUSmLULjAOOd0XVjS83W2KozZTEc32t7BUvMdF/dLIQ3Z61L1x4yOZMevThrffvq1K41quEvZ9vGlPtMKlmMj//zlpLlQ17POtodEt1De/QsFfU98uZHtp53VHju5T6XbtbDS6LrbipvBd+a/iyie4gD1scW5KzrsqhuG5ENZt8xeWL2unMXZNKdqb+HtknZtlnb6u5yDH2+dDCu9jWxiXCW+G8t/Wzj5x5/d2Lr87g8/f3ordVzbmq3y2gG8GFggaT5Sd5XgPOBqyWdSVi7PbmX1ykpcd5tpy9TZpjhOI4zaKndLqN7KP32fnRtrgJlrlF+ykjS1yv2nKOO4zhOMyJqN2VUM5lCyIByeX8t1bbSCOFjktaVKRchTsHXK/TjOI7TfBio8Zz4XAscGOVdAxwEYGafLtWwkkL4JTAmRx3HcZzBSeP4MtobeDOwjaT3pYrGAuVd3yZUskP4Rs/FcxzHGQQ0iEIA9iL4MxoHvDuVvx741zwd5N1l5DiO4xShUbydmtnvgd9LeouZ3duTPuqqECTtRfDSV2A34GsE/xzvBjYDzwBnmFm3vWqSlhC0WwfQnsfKz3Ecp09pEIWQ4ukkYtoUUr/xZlbRfUVdFYKZLQKmASTxFF4k2DLsBZyTxE3+DnAOpZ01HWlmq0uUOY7j9B+Nuaj8e+AvwJ8IL9O5yRtT+R+ACwlOlvaRtB8w08y+VcW1jgaeSVy0PpfKv48QX8FxHGfg0XgjhJFmVs4baknyuq74JeEtfguAmT1K2G5aDaeQcp+d4qPATSXaGHCrpHmSZhWrIGmWpLmS5m7u3FilSI7jOL2j0ewQgBslndCThnmnjEaa2QOR1732UpVjkvjLMwlKJZ3/1aSfOSWazjCzZYm719skPWlmd6crJO5jZwNsM3JHs4nbdpWlgndDd5P5cgy/I2vmbykT+bgsdntQjthlROy6oW1KNqh8HNw9baofu3mIie83XT8OyL75uIMz6ZHLI1cNL67KpMemgq7HQdUrBbdP33N8D+VcEcR0vrSqZBkUcc9WDAAuAAAaIklEQVQRke47dt0Q30M5FyKxi4RK37NubjFS91FtX+ln282tQ/w57PfGkv3Ev1fx304ludLllZ5l/BmXfSuNXXdExN/jfqHxRghnAV+RtInwEi/AzGxspYZ5FcJqSbvTFVP5JJJADjk5HnjIzFYUMiSdTtgidbSZFX2kZrYs+X+lpOuBQwh+lBzHcfqfKjyZ9hVmVsl2rCR5FcKnCG/he0t6EVgMnFbFdU4lG23tOMIi8j+ZWdFYn5JGAS1J5KBRwLHkCAHnOI7TV4gyjoH6GEl7m9mTkmIrZQDM7KFKfeR1bvcscEz6R7oKIUcCbwc+nsr+L2AYYRoI4D4z+4SkHYGLzewEYAfg+qS8Dfitmd2c97qO4zh9QQPtMvoCMAv4QZEyA46q1EGlmMpfKJEfrpD12V2UZASwXZS3R4m6y4ATkvNngf0r9e84jtOvNMiUkZnNSv4/sqd9VBohFOai9gIOJgRqgGBU5nP5juM4DaIQCiRR1/4NeFuSdSdwkZmVD1xCTl9Gkm4FDixMFSUur3/Xc5Edx3GagBpuKZV0KWGjzUoz2yfJ+zrBD1FhK9pX0uGMS3AhMAT4eZL+cJL3sUoy5F1U3oXgZqLAZoJZtOM4zuCmdiOEywjrq/8d5V9gZt+vop+DzSw93f5nSY+UrJ0ir0L4NfBAsvXTgPfSXWjHcZxBR60Wlc3sbklTatBVh6TdzewZAEm7kdOFRd5dRt+WdBPwj0nWGWb2cI9EdRzHaSKqmDKaIGluKj07MaytxKcl/QswF/h3M3ulQv0vAndIepawK3ZX4Iw8Aub1ZbQLsJrgmG5rnpk9n6d9X2ESnUO7bqn98DdnyoctXdtVNwqM3rlHZEEcBSRvTVlXtrwhG7ycyGK2JbI23jC1y3p65MKXMmWx1WbHmGww+9bYyjdtbXzw3tm6Dz6ZSa+fOS1bvrnrmztqcTYQ3oj7nqIcnbFFcTq9w4RMWXz/Nixrbaq1XYHh2/eZkpUxuof4+WSuU+bZAHRMzJa3rs+6Nsn0HVnExn1v2T1rfTtkVdc9KPoubT5o92zb0dk/s9EPPJdJp+XoZrVbwVI3bfUdy8zY0dm6z2dtSdPPa81JB2TKhu+4byY9bHVkqbwqshBOW1tXsEyuZGGfpv25FzLp2Apa8WdexnK5m9V3fjFKU51h2uoeeGy+EPhmcpVvEraUlvVaama3S9qTsBlIwJNmtinPxfJOGf2BrtseAUwFFhGi8ziO4wxe6rjLKPLu8EvgxkptJA0HPgkcnkj3F0m/MLOKKjDvlFHmdSGxhPt4ieqO4ziDAlFfx3WSJplZYWj3XuCxHM3+mxBH5qdJ+lTCOvDJlRr2KB6CmT0k6eDKNR3HcZqc2m07vQI4grDWsBQ4FzhC0rTkKkvI9yK+V7TL6I6a7jKKLJZbgAPp2hfrOI4zODFQZ200gpmdWiT7kh509bCkw8zsPgBJhwJ/zdMw7wgh7T2vnbCmcG1VIjqO4zQhjRJTOcWhwL9IKmz62QVYKGkBwQ32fqUa5lUIT5hZxjJZ0sm4tbLjOIOdxlMIx/W0YV6FcA7df/yL5TmO4wwqGm2EkIQp7hGVvJ0eT/A+OlnST1JFY6kiYprjOE7T0mAKoTdUGiEsI1jHzQTmpfLXA5+vl1CO4zgDgr6Pl1xXKnk7fQR4RNIcM2v4EUH7qBZWH9S1/r3DnSuzFdZ1WZfGFo1tq7JR5zrLxCOOLS+7xZtNXQdgyKtdFqOx5WUcX7fl6Ww520TR8FJypK1lAWyvqdmmj2YtaDMWw5HVqqLr2DZZK9fORxdm5Ujdc0tkLdrNYja2TE1ZNuue+dnrEFEmXnXLxugz2nfPbNMJ2fKhkfX5ug+9Zev5uGuynljiz1TLsxbm5f4YWm+fl01H5fX6Q4qfR2ckczk2bJ+N+zV2zoOZtOIYylH7+HlliMq6xZROeQnojD6jmEpWzuk42nHM7W4eBpaU7SoXoqEC5PSaSlNGV5vZBwjbmLrpwXKr1Y7jOIOC4iHhBySVpozOSv5/V086l7QXcFUqazfgawRLuqsILrSXAB8o5rBJ0unAfybJb5nZ5T2Rw3Ecp14005RRN19aaVIm0580s+fSB8FXRlnMbJGZTTOzacBBwAaCg7yzgdvNbE/g9iSdQdJ4gqXeocAhwLmSto3rOY7j9BtWxTEAKKsQUry9SN7xVV7raOCZRJmcCBTe9i8H3lOk/juA28zs5WT0cBu92F/rOI5TD9SZ7xgIVFpD+DfCSGA3SY+misaQ0xQ6xSnAFcn5DoXRh5ktl7R9kfqTgfQK69IkL5ZxFjALYMhoH0A4jtO3DJQf+zxUWkP4LXAT8H/ITuusN7OX815E0lDC1tVzqpBNRfKKLWzPBmYDjNx+5wEyMHMcpykwmmpRudIawlozW2JmpyZTPRsJj2B0EjQnL8cDD6V8e6+QNAmCe1dgZZE2S4F01JqdCHYRjuM4DYMs3zEQyLWGIOndkp4CFgN3EXYG3VTFdU6la7oI4Abg9OT8dOD3RdrcAhwradtkMfnYJM9xHKdxGISLyt8CDgP+18ymEhaIc60hSBpJWJS+LpV9PvD2RMm8PUkjabqkiwGSKalvAg8mx3nVTFM5juPUm0KAnGYZIeR1brfFzP4uqUVSi5ndIek7eRqa2QZguyjv7wSlEtedC3wslb4UuDSnjI7jOH2LWc3WECRdSrD5Wmlm+yR548lhs1Ur8iqENZJGA3cDcyStpAGd27Vt7GT84xu6MiIXEmXZvCWTjN1TpIPBt7zw90xZ635vzPa1NnJdkXYxsevOmbJugeCj4OXdyocP60q8no2b/dreGb3LqLsWZdKZe4pcVcRm/opcV7RFcqeDu8d0cy+wy6RsOvV84qDpcdvYxUg6UHrsLmHdXlmZx865N5OOr7Xtgq6+O8q5XhgglHUfUYE3XPC3TDoOSB+7I2mP3GKU+1y6BbePSLtrWXfyYZmybR54MZOO/x66uXpJ0e15vFSfmF413GV0GfBfBMPdAgWbrfMlnZ2kv1yzK0bknTI6kbCg/HngZuAZ4N31EspxHGegUKspIzO7G4inxfPYbNWMXCMEM3stlXT3EY7jOBAWi2sUQrMEeWy2akYlw7T1FF8fFyEU29i6SOU4jjNQyK8PJkiam0rPTuyoGoZK7q/HlCt3HMcZ7FSxg2i1mU2vsvsVkiYlo4NSNls1I+8aguM4jlOMwk6jSkfPyGOzVTPy7jJyHMdxYqx2u4wkXQEcQZhaWkrw9nw+cLWkM4HngZNrc7XiuEJwHMfpIcEwrTaLymZ2aomibjZb9cIVguM4Tm8YRN5OHcdxnDLUaoTQCDSVQtg8poWlR47cmp76TLY8bV0ZW0/axGwsBW3KWi5nrI0jq+bOodnH2L5Hdqtwy6aOredtj2UD37dGfb16yK6ZdGxtTMraeO3hUzJFoxe/lklvPCwbdH71/kO2nu98U9b6XZHlsS1anEmXs+SNLYA7Jk8sUTMhZW26aUL2cxi5OCtX24jSVq6rZu6WSc8998JM+p1z359Jty96OttBFUHoBxuxlW9nL55VHNw+tjbuHN669XzMDfMzZe2RHK3R91Q7TMj2tWRpSTl6Y8ldkgHkuC4PTaUQHMdx+hZD9TVM61NcITiO4/QGnzJyHMdxarnttBFwheA4jtMbfITgOI7jAL6o7DiO4wR822kVSBoHXAzsQ9ClHwU+B+yVVBkHrDGzaUXaLgHWAx1Aew8cQzmO49QPAzpcIVTDj4GbzewkSUOBkWb2wUKhpB8Aa8u0P9LMVpcpdxzH6ReE+QghL5LGAm8DPgJgZpuBzalyAR8AjqqnHI7jOHXDFUJudgNWAb+StD8wDzgrFYHtH4EVZvZUifYG3CrJgIuKBZOQNAuYBTC8dTRTL352a1kc9zVNN6vFRxeWvZG0ZXPLlJ0yZa3rN5Zvm4r7GscIJkqOuiu/NeWQ9R2ZtBZkH+PQ6B53uqPrHmLL49jaOKatTEzdTW+anCkb9kQ2Du6W3bN9tz22ZOv58AWRRWx03fhzeurCQ7v6iZ7d8VMPzaQ7N0WWyU7dKGcF3P5cFPc4SivdT4XrdP/7WVO8Yl/SRAqh3vEQ2oADgQvN7ADgNUKQ6AKnAleUaT/DzA4Ejgc+JeltcQUzm21m081s+tCWETUU3XEcpwJG0GJ5jgFAvRXCUmCpmd2fpK8hKAgktQHvA64q1djMliX/rwSuBw6pq7SO4zhVIrNcx0CgrgrBzF4CXpBU2FF0NPBEcn4M8KSZFfVGJWmUpDGFc+BY4LF6yus4jlMdBp2d+Y4cSFoiaYGk+VH85T6hL3YZfQaYk+wwehY4I8k/hWi6SNKOwMVmdgKwA3B9WHemDfitmd3cB/I6juPkw6jHGkK/7aysu0Iws/lAN/sBM/tIkbxlwAnJ+bPA/vWWz3Ecp1cMkPWBPNR7DcFxHKepqfEaQmFn5bxkB2Wf4q4rHMdxekP+H/sJ0brA7CJb6WeY2TJJ2wO3SXrSzO6uiZw5cIXgOI7TU8ygI/ec0epK7nfSOyslFXZW9plC8Ckjx3Gc3mCW76hAI+ys9BGC4zhOb6jdLqN+31nZVArBtrSXdVdRDWlXFZB1V7HiiO0zZRMveyjbeN9scHttM2breVvqHIDNWzJJ21jedUW6fPgdC8rWjUm7F4jvj6FDsmJF7iiGznsmk24d1+XKojUqi/88Wh98smx5mjgg+4ZpO2bSe8zZtPVc92QDsjfRZg9noGBAjWIqN8LOyqZSCI7jOH2LgTXPq4grBMdxnN4wQNxS5MEVguM4Tk8xqtll1PC4QnAcx+kNPkJwHMdxwhqCKwTHcRzHyO3JdCDgCsFxHKc3+AjBcRzHAVwhOI7jOIAZ1tFRud4AwRVCCeKg4Z2LugK277Du1UxZexywfv3GTHrTHl2WzbHFbxw0vJsFcQW50qSthwFaR2T7UiodW3R3RoHPh65dn227w4Rs/SVdge7KyVRMrs0H7b71fNjTKzNlcUD24XGAdsdpNGpkqdwI1N25naRxkq6R9KSkhZLeIunrkl5MwsTNl3RCibbHSVok6WlJZ9dbVsdxnKqpkXO7RqAvRgg/Bm42s5OSMJojgXcAF5jZ90s1ktQK/Ax4O7AUeFDSDWb2RKk2juM4fYpZU+0yqusIQdJY4G3AJQBmttnM1pRvtZVDgKfN7Fkz2wxcCZxYH0kdx3F6SBONEOo9ZbQbsAr4laSHJV2c+PkG+LSkRyVdKmnbIm0nA+kJ5KVJnuM4ToMQFpXzHAOBeiuENuBA4EIzOwB4DTgbuBDYHZgGLAd+UKStiuR1U7OSZkmaK2nuFjYVaeI4jlMnCu6v8xwDgHorhKXAUjO7P0lfAxxoZivMrMPMOoFfEqaHirXdOZXeCVgWVzKz2WY23cymD2FYjcV3HMepgHXmOwYAdVUIZvYS8IKkvZKso4EnJE1KVXsvxcPEPQjsKWlqshh9CnBDPeV1HMepBgOs03IdeejvnZV9scvoM8Cc5Ef9WeAM4CeSphGe5xLg4wCSdgQuNrMTzKxd0qeBW4BW4FIze7wP5HUcx8mH1S5ATiPsrKy7QjCz+cD0KPvDJeouA05Ipf8I/LF+0jmO4/SOvG//Odi6sxJAUmFnZZ8pBNkA2Q6VB0mrgOeACcDqfhanGC5XdTSiXI0oE7hc1VCQaVczm1ipcjkk3Zz0l4fhQNqsf7aZzU71dRJwnJl9LEl/GDjUzD7dGxmroalcVxQ+XElzzSwelfQ7Lld1NKJcjSgTuFzVUEuZzOy4WvSTkGtnZT2pu+sKx3EcJxe5dlbWE1cIjuM4jUG/76xsqimjFLMrV+kXXK7qaES5GlEmcLmqoRFlohF2VjbVorLjOI7Tc3zKyHEcxwFcITiO4zgJTaUQJH0vCcTzqKTrJY1LlZ2TmIMvkvSOPpbrZEmPS+qUND0q60+5GiIAUeLxdqWkx1J54yXdJump5P9iHnHrLdfOku5IAjs9Lums/pZN0nBJD0h6JJHpG0n+VEn3JzJdlSxK9jmSWhPPxjc2ilySlkhakATjmpvk9fv3qxFpKoUA3AbsY2b7Af8LnAMg6U2EFfs3A8cBP0/MxPuKx4D3AXenM/tTrpSZ/PHAm4BTE3n6g8sI95/mbOB2M9sTuD1J9zXtwL+b2RuBw4BPJc+oP2XbBBxlZvsTvAUfJ+kw4DuEoFN7Aq8AZ/ahTGnOAham0o0i15FmNi1lf9AI36+Go6kUgpndambtSfI+wj5eCObfV5rZJjNbDDxNcQ+r9ZJroZktKlLUn3I1TAAiM7sbeDnKPhG4PDm/HHhPnwoFmNlyM3soOV9P+KGb3J+yWaAQ1HtIchhwFMGbcJ/LVEDSTsA7gYuTtBpBrhL0+/erEWkqhRDxUeCm5LxRg+30p1yN+kwK7GBmyyH8MAPb96cwkqYABwD308+yJdMy84GVhFHxM8Ca1MtQf32WPwK+BBS8vW3XIHIZcKukeZJmJXkN9f1qFAacHYKkPwFvKFL0VTP7fVLnq4Th/pxCsyL1a7rfNo9cxZoVyeurfcD9biY/UJA0GrgW+JyZrQsvvv2HmXUA05I1suuBNxar1pcySXoXsNLM5kk6opBdpGp/fMdmmNkySdsDt0l6sh9kGBAMOIVgZseUK5d0OvAu4GjrMrKou0l4JblK0J+m6v1uJl+BFZImmdnyJH7Gyv4QQtIQgjKYY2bXNZJsZrZG0p2E9Y1xktqSt/H++CxnADMlnUBw4jaWMGLob7kKXpQxs5WSridMlzbEZ9hoNNWUkaTjgC8DM81sQ6roBuAUScMkTQX2BB7oDxkj+lOufjeTr8ANwOnJ+elAqVFW3UjmwC8BFprZDxtBNkkTC7vnJI0AjiGsbdwBnNQfMgGY2TlmtpOZTSF8l/5sZh/qb7kkjZI0pnAOHEvY5NHv36+GxMya5iAsyr4AzE+OX6TKvkqYa10EHN/Hcr2X8Ea+CVgB3NIgcp1A2I31DGFqq78+tysIsbW3JM/pTML88+3AU8n/4/tBrsMJUxyPpr5TJ/SnbMB+wMOJTI8BX0vydyO8TDwN/A4Y1o+f5xHAjY0gV3L9R5Lj8cL3vBG+X414uOsKx3EcB2iyKSPHcRyn57hCcBzHcQBXCI7jOE6CKwTHcRwHcIXgOI7jJLhCcBzHcQBXCIMOSa9WrlV1nzML7rMlvacnXlMl3Rm7Bs9Rf5GkmUXKpqRdaTc7kj4iacdUeo6klyWdVK6d48S4QnB6jZndYGbnJ8n3ENxp9wUfMrO6Wlf3sZv0nvIRYKtCsGAh3EhW584AwRXCIEWB70l6LAke8sEk/4jk7fsahWBDcxIXDkg6Icm7R9JPUkFQPiLpvyS9FZgJfC8JRrJ7+s1f0gRJS5LzEZKuVAhmdBUwIiXbsZLulfSQpN8lzuUq3c9BCkFj7gU+lcpvTe7zweRaH0/yWyT9XCHIzI2S/lh4o1YIqPI1SfcAJyf3cXPiLfMvkvZO6k2UdG3S94OSZiT5/5Tc/3yFYDFjysj9xZRs30jl/09yvceVeOhM7uWy1Gf2+UTm6cCc5HojSl3LcSox4JzbOTXjfYQAK/sDE4AHJRUC+BxACNqzDPgrMEMh0tRFwNvMbLGkK+IOzexvkm4guC24BkClPYP+G7DBzPaTtB/wUFJ/AvCfwDFm9pqkLwNfAM6rcD+/Aj5jZndJ+l4q/0xgrZkdLGkY8FdJtwIHAVOAfQmujxcCl6bavW5mhycy3Q58wsyeknQo8HOCn/8fE4K/3CNpF+AWgufR/wA+ZWZ/TZTZ68UElnQswX/VIQTPoDdIepuF+BAfNbOXkx/4ByVdm8g72cz2SdqPs+Dg7tPAf5jZ3ArPyHHK4gph8HI4cIUFV8orJN0FHAysAx4ws6UACn73pwCvAs9aCOQDwf/QrG695udtwE8AzOxRSY8m+YcRppz+miiTocC95TqStA0wzszuSrJ+TYgEB8GZ2X6p+fRtCD/ChwO/M7NO4CVJd0TdXpX0PRp4K/C7lHIblvx/DPCmVP7YZDTwV+CHkuYA1xWeZRGOTY6Hk/ToRLa7gc9Kem+Sv3OSvwjYTdJPgT8At5Z7Lo5TLa4QBi/lnPpvSp13EL4nPQ0C0E7X1OTwqKyYIy0Bt5nZqVVcQyX6KpR9xsxuyWRK76zQ52vJ/y2EIC/TitRpAd5iZhuj/PMl/YHgCO8+SceYWTEf/AL+j5ldFMl2BEHZvMXMNii4uB5uZq9I2h94B2Fa7AOEQFCOUxN8DWHwcjfwwWReeiLhjb2c6+0nCW+nU5L0B0vUWw+k58yXEKZnoMsNcuH6HwKQtA/BiyeE0KczJO2RlI2U9A/lbsTM1gBrJR2eZH0oVXwL8G8KcQ2Q9A8KbpDvAd6frCXsQPDQWazvdcBiSScn7ZX8KEN4Q/90oa6kacn/u5vZAjP7DjAX2LuE6LcAHy2skUiarBDEZRvglUQZ7E0YNRWm01rM7Frg/wMOTPqJn7nj9AhXCIOX6wkulB8B/gx8ycxeKlU5eQv+JHBzsti6AlhbpOqVwBeTxdTdge8TfpD/RlirKHAhMDqZKvoSiTIys1WEXTNXJGX3UfoHNc0ZwM+SReX0G/vFwBPAQwpbUS8ijHiuJbjaLuTdX+J+ICiYMyUVXCgXYk9/FpieLAg/AXwiyf9csvD7SCLLTd16DPd6K/Bb4F5JCwixh8cANwNtyf1/M3kGEMJP3plM410GnJPkXwb8wheVnd7i7q+d3EgabWavKkya/wx4yswu6CdZ7qSXC6mp+9mOoJBmlFOKAwlJl5Fa3HecPPgIwamGf03eTh8nTGtcVKF+PXkZuExFDNOq4Mbkfv4CfLOJlMEc4J8osbvJcUrhIwTHqTOS9iXsfEqzycwO7Q95HKcUrhAcx3EcwKeMHMdxnARXCI7jOA7gCsFxHMdJcIXgOI7jAPD/AF4ImptnG4wpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dens_nonweighted[\"cat1\"].plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Area weights are calculated by `octant.grid.grid_cell_areas()` function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Grid centres vs bounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is what coordinates of the density array look like, e.g. longitude:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray 'longitude' (longitude: 71)>\n",
       "array([-20., -19., -18., -17., -16., -15., -14., -13., -12., -11., -10.,  -9.,\n",
       "        -8.,  -7.,  -6.,  -5.,  -4.,  -3.,  -2.,  -1.,   0.,   1.,   2.,   3.,\n",
       "         4.,   5.,   6.,   7.,   8.,   9.,  10.,  11.,  12.,  13.,  14.,  15.,\n",
       "        16.,  17.,  18.,  19.,  20.,  21.,  22.,  23.,  24.,  25.,  26.,  27.,\n",
       "        28.,  29.,  30.,  31.,  32.,  33.,  34.,  35.,  36.,  37.,  38.,  39.,\n",
       "        40.,  41.,  42.,  43.,  44.,  45.,  46.,  47.,  48.,  49.,  50.])\n",
       "Coordinates:\n",
       "  * longitude  (longitude) float64 -20.0 -19.0 -18.0 -17.0 ... 48.0 49.0 50.0\n",
       "Attributes:\n",
       "    units:    degrees_east"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arr.longitude"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Density is calculated assuming these are grid cell **centres**.\n",
    "\n",
    "It can be also assumed that `lon_dens1d` and `lat_dens1d` are grid cell **boundaries**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "dens_b = tr.density(lon_dens1d, lat_dens1d, grid_centres=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how the shape and values of longitude (and latitude) change:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray 'longitude' (longitude: 70)>\n",
       "array([-19.5, -18.5, -17.5, -16.5, -15.5, -14.5, -13.5, -12.5, -11.5, -10.5,\n",
       "        -9.5,  -8.5,  -7.5,  -6.5,  -5.5,  -4.5,  -3.5,  -2.5,  -1.5,  -0.5,\n",
       "         0.5,   1.5,   2.5,   3.5,   4.5,   5.5,   6.5,   7.5,   8.5,   9.5,\n",
       "        10.5,  11.5,  12.5,  13.5,  14.5,  15.5,  16.5,  17.5,  18.5,  19.5,\n",
       "        20.5,  21.5,  22.5,  23.5,  24.5,  25.5,  26.5,  27.5,  28.5,  29.5,\n",
       "        30.5,  31.5,  32.5,  33.5,  34.5,  35.5,  36.5,  37.5,  38.5,  39.5,\n",
       "        40.5,  41.5,  42.5,  43.5,  44.5,  45.5,  46.5,  47.5,  48.5,  49.5])\n",
       "Coordinates:\n",
       "  * longitude  (longitude) float64 -19.5 -18.5 -17.5 -16.5 ... 47.5 48.5 49.5\n",
       "Attributes:\n",
       "    units:    degrees_east"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dens_b[\"cat1\"].longitude"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:test_env]",
   "language": "python",
   "name": "conda-env-test_env-py"
  },
  "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}