Week8_EDS220_MultiDataset-interponly.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8e88be16-a155-4c09-b8cc-0f82867d0f82",
   "metadata": {},
   "source": [
    "# Ocean Remote Sensing: Combining Datasets\n",
    "## Marine Heat Wave Focus Continued\n",
    "## EDS 220 Fall 2021\n",
    "\n",
    "This week, we will build on the marine heat wave exercise we began last week, and expand on our Python analysis skills while becoming familiar with new sources of environmental data.\n",
    "\n",
    "The major topics/learning objectives for this week are:\n",
    "- Additional familiarity with probability distribution functions\n",
    "- Interpolation of spatial data fields to a common grid\n",
    "- Combining multiple datasets to gain insight into an environmental system\n",
    "- Visualizing multiple datasets on a single plot (color/quiver plots)\n",
    "\n",
    "These exercises build on the marine heat wave diagnostics discussed last week and presented in the \"marineHeatWave\" package developed by Eric Oliver for the Hobday et al. (2016) publication (see References) and available on Github:\n",
    "https://github.com/ecjoliver/marineHeatWaves\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f08e2001-c2eb-40d3-b800-03354782e371",
   "metadata": {},
   "source": [
    "### Datasets\n",
    "\n",
    "Several datasets will be used during this exercise:\n",
    "\n",
    "- SST: NOAA's 1/4° Daily Optimum Interpolation Sea Surface Temperature (OISST) version 2, or OISSTv2\n",
    "- Surface currents: NCEP Global Ocean Data Assimilation System (GODAS)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f8912e4-8032-44e1-ba4f-6a6dcb5bf9c2",
   "metadata": {},
   "source": [
    "### Science Goals\n",
    "\n",
    "Last week, we learned about the importance of marine heat wave events for human and ecological systems. This week, we'll apply our new knowledge of environmental data analysis to answer more complex questions about these events:\n",
    "\n",
    "- How does ocean circulation change during sustained marine heat waves?\n",
    "\n",
    "### Step Zero: Configure Environment\n",
    "\n",
    "First, load in the required packages. Note that we'll be working with daily OISSTv2 data stored in a \"pickle\" file, so the pickle package is included in addition to other packages from last week."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "79113a11-1b9d-412a-a5f2-5b4b4f42664c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure environment\n",
    "import numpy as np\n",
    "import scipy\n",
    "import scipy.ndimage as ndimage\n",
    "import scipy.stats as stats\n",
    "import matplotlib.pyplot as plt\n",
    "import xarray as xr\n",
    "import cartopy\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.feature as cfeature\n",
    "from cartopy.feature import NaturalEarthFeature\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e8b3359-c1dd-482c-a57c-fca886e4f1b3",
   "metadata": {},
   "source": [
    "### Review from Week 7: Marine Heat Wave ID, properties\n",
    "\n",
    "The commands for extracting the Northeast Pacific portion of the OISSTv2 data from the THREDDS server are included below, but commented out since they have already been executed and stored in a file: 'oisstv2_nepac.pickle'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d82b8e12-aea8-4663-94ed-9a77a7d0c896",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "range(2018, 2021)\n"
     ]
    }
   ],
   "source": [
    "## NOAA THREDDS server: base URL for all files \n",
    "#base_url=\"https://psl.noaa.gov/thredds/dodsC/Datasets/noaa.oisst.v2.highres/\"\n",
    "\n",
    "## Naming convention for files: \"sst.day.mean.[year].v2.nc\"\n",
    "yrs=range(2018,2021)\n",
    "print(yrs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "382fdf15-5c8d-4383-b6a9-b382557db180",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2018\n"
     ]
    },
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] No such file or directory: 'https://psl.noaa.gov/thredds/dodsC/Datasets/noaa.oisst.v2.highres/sst.day.anom.2018.v2.nc'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[0;32m/opt/anaconda3/envs/eds223/lib/python3.8/site-packages/xarray/backends/file_manager.py\u001b[0m in \u001b[0;36m_acquire_with_cache_info\u001b[0;34m(self, needs_lock)\u001b[0m\n\u001b[1;32m    198\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 199\u001b[0;31m                 \u001b[0mfile\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_key\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    200\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/envs/eds223/lib/python3.8/site-packages/xarray/backends/lru_cache.py\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m     52\u001b[0m         \u001b[0;32mwith\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_lock\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 53\u001b[0;31m             \u001b[0mvalue\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cache\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     54\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cache\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmove_to_end\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyError\u001b[0m: [<function _open_scipy_netcdf at 0x7fbc93dbfc10>, ('https://psl.noaa.gov/thredds/dodsC/Datasets/noaa.oisst.v2.highres/sst.day.anom.2018.v2.nc',), 'r', (('mmap', None), ('version', 2))]",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m/var/folders/8s/qv6gw27d6nsbh7wpqhj27p600000gn/T/ipykernel_98518/1578530640.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m     \u001b[0;31m# Open the dataset into an xarray object\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m     \u001b[0mmysst\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mxr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mopen_dataset\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"https://psl.noaa.gov/thredds/dodsC/Datasets/noaa.oisst.v2.highres/sst.day.anom.\"\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m\".v2.nc\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      8\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m     \u001b[0;31m# Define logical mask to be true in the Northeast Pacific (30-50N, 180-240E)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/envs/eds223/lib/python3.8/site-packages/xarray/backends/api.py\u001b[0m in \u001b[0;36mopen_dataset\u001b[0;34m(filename_or_obj, engine, chunks, cache, decode_cf, mask_and_scale, decode_times, decode_timedelta, use_cftime, concat_characters, decode_coords, drop_variables, backend_kwargs, *args, **kwargs)\u001b[0m\n\u001b[1;32m    495\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    496\u001b[0m     \u001b[0moverwrite_encoded_chunks\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpop\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"overwrite_encoded_chunks\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 497\u001b[0;31m     backend_ds = backend.open_dataset(\n\u001b[0m\u001b[1;32m    498\u001b[0m         \u001b[0mfilename_or_obj\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    499\u001b[0m         \u001b[0mdrop_variables\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdrop_variables\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/envs/eds223/lib/python3.8/site-packages/xarray/backends/scipy_.py\u001b[0m in \u001b[0;36mopen_dataset\u001b[0;34m(self, filename_or_obj, mask_and_scale, decode_times, concat_characters, decode_coords, drop_variables, use_cftime, decode_timedelta, mode, format, group, mmap, lock)\u001b[0m\n\u001b[1;32m    280\u001b[0m         \u001b[0mstore_entrypoint\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mStoreBackendEntrypoint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    281\u001b[0m         \u001b[0;32mwith\u001b[0m \u001b[0mclose_on_error\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstore\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 282\u001b[0;31m             ds = store_entrypoint.open_dataset(\n\u001b[0m\u001b[1;32m    283\u001b[0m                 \u001b[0mstore\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    284\u001b[0m                 \u001b[0mmask_and_scale\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmask_and_scale\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/envs/eds223/lib/python3.8/site-packages/xarray/backends/store.py\u001b[0m in \u001b[0;36mopen_dataset\u001b[0;34m(self, store, mask_and_scale, decode_times, concat_characters, decode_coords, drop_variables, use_cftime, decode_timedelta)\u001b[0m\n\u001b[1;32m     22\u001b[0m         \u001b[0mdecode_timedelta\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     23\u001b[0m     ):\n\u001b[0;32m---> 24\u001b[0;31m         \u001b[0mvars\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mattrs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mstore\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mload\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     25\u001b[0m         \u001b[0mencoding\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mstore\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_encoding\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     26\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/envs/eds223/lib/python3.8/site-packages/xarray/backends/common.py\u001b[0m in \u001b[0;36mload\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    122\u001b[0m         \"\"\"\n\u001b[1;32m    123\u001b[0m         variables = FrozenDict(\n\u001b[0;32m--> 124\u001b[0;31m             \u001b[0;34m(\u001b[0m\u001b[0m_decode_variable_name\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_variables\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    125\u001b[0m         )\n\u001b[1;32m    126\u001b[0m         \u001b[0mattributes\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mFrozenDict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_attrs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/envs/eds223/lib/python3.8/site-packages/xarray/backends/scipy_.py\u001b[0m in \u001b[0;36mget_variables\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    169\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mget_variables\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    170\u001b[0m         return FrozenDict(\n\u001b[0;32m--> 171\u001b[0;31m             \u001b[0;34m(\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mopen_store_variable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mk\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mv\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mds\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvariables\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    172\u001b[0m         )\n\u001b[1;32m    173\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/envs/eds223/lib/python3.8/site-packages/xarray/backends/scipy_.py\u001b[0m in \u001b[0;36mds\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    158\u001b[0m     \u001b[0;34m@\u001b[0m\u001b[0mproperty\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    159\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mds\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 160\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_manager\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0macquire\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    161\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    162\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mopen_store_variable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvar\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/envs/eds223/lib/python3.8/site-packages/xarray/backends/file_manager.py\u001b[0m in \u001b[0;36macquire\u001b[0;34m(self, needs_lock)\u001b[0m\n\u001b[1;32m    179\u001b[0m             \u001b[0mAn\u001b[0m \u001b[0mopen\u001b[0m \u001b[0mfile\u001b[0m \u001b[0mobject\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mreturned\u001b[0m \u001b[0mby\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m`\u001b[0m\u001b[0;31m`\u001b[0m\u001b[0mopener\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;31m`\u001b[0m\u001b[0;31m`\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    180\u001b[0m         \"\"\"\n\u001b[0;32m--> 181\u001b[0;31m         \u001b[0mfile\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_acquire_with_cache_info\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mneeds_lock\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    182\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mfile\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    183\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/envs/eds223/lib/python3.8/site-packages/xarray/backends/file_manager.py\u001b[0m in \u001b[0;36m_acquire_with_cache_info\u001b[0;34m(self, needs_lock)\u001b[0m\n\u001b[1;32m    203\u001b[0m                     \u001b[0mkwargs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    204\u001b[0m                     \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"mode\"\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_mode\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 205\u001b[0;31m                 \u001b[0mfile\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_opener\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_args\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    206\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_mode\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m\"w\"\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    207\u001b[0m                     \u001b[0;31m# ensure file doesn't get overriden when opened again\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/envs/eds223/lib/python3.8/site-packages/xarray/backends/scipy_.py\u001b[0m in \u001b[0;36m_open_scipy_netcdf\u001b[0;34m(filename, mode, mmap, version)\u001b[0m\n\u001b[1;32m     97\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     98\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 99\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mio\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnetcdf_file\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmmap\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmmap\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mversion\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mversion\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    100\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0mTypeError\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m  \u001b[0;31m# netcdf3 message is obscure in this case\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    101\u001b[0m         \u001b[0merrmsg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/anaconda3/envs/eds223/lib/python3.8/site-packages/scipy/io/netcdf.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, filename, mode, mmap, version, maskandscale)\u001b[0m\n\u001b[1;32m    245\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfilename\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfilename\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    246\u001b[0m             \u001b[0momode\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'r+'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mmode\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'a'\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 247\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'%sb'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0momode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    248\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mmmap\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    249\u001b[0m                 \u001b[0;31m# Mmapped files on PyPy cannot be usually closed\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'https://psl.noaa.gov/thredds/dodsC/Datasets/noaa.oisst.v2.highres/sst.day.anom.2018.v2.nc'"
     ]
    }
   ],
   "source": [
    "nepac_sst=[]     # Blank variable to store output data \n",
    "\n",
    "for y in yrs:\n",
    "    print(y)\n",
    "    \n",
    "    # Open the dataset into an xarray object\n",
    "    mysst=xr.open_dataset(\"https://psl.noaa.gov/thredds/dodsC/Datasets/noaa.oisst.v2.highres/sst.day.anom.\"+str(y)+\".v2.nc\")\n",
    "    \n",
    "    # Define logical mask to be true in the Northeast Pacific (30-50N, 180-240E)\n",
    "    mask_lon = (mysst.lon >= 180) & (mysst.lon <= 240)\n",
    "    mask_lat = (mysst.lat >= 30) & (mysst.lat <= 50)\n",
    "    # Read in data \n",
    "    nepac_sst.append(mysst.where(mask_lon & mask_lat, drop=True).squeeze())\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3974239d-ef82-46ec-9015-e013ef1078a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "#with open('oisstv2_nepac.pickle','wb') as f:\n",
    "#    pickle.dump(nepac_sst,f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed157ec4-56dd-408f-b55b-334d5bfce3b5",
   "metadata": {},
   "source": [
    "Now we can load in the previously stored information from the pickle file with `pickle.load`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cabfc1e-5337-47b4-a0ee-b04308fbef2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load in NE Pacific OISSTv2 data\n",
    "#with open('oisstv2_nepac.pickle','rb') as f:\n",
    "#    nepac_sst=pickle.load(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "281ea5c6-b7de-414a-8ff2-be5f720c6c44",
   "metadata": {},
   "source": [
    "The data are stored in a list at the moment, so they need to be translated into an xarray Dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9566f1df-8aef-42d6-92c5-d5569a29342a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Concatenate list of xarray Datasets into a single Dataset, along dimension \"time\"\n",
    "nepac_fldts=xr.concat(nepac_sst,dim='time')\n",
    "nepac_fldts\n",
    "\n",
    "# Calculate a regionally averaged time series for the NE Pacific\n",
    "nepac_ts=nepac_fldts.mean('lon').mean('lat')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05943744-1921-4039-adf6-6892996b721f",
   "metadata": {},
   "source": [
    "Now we determine the SST threshold appropriate for a marine heat wave using `np.percentile`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35738e49-cae0-4406-be50-860ac6087967",
   "metadata": {},
   "outputs": [],
   "source": [
    "thr=nepac_ts.reduce(np.percentile,q=90)\n",
    "print(thr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dff7042a-c915-4de1-a23f-f047a770517a",
   "metadata": {},
   "source": [
    "and apply this as a logical mask:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "928889a8-09ca-4e53-bc96-105a50768dbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mask for MHW\n",
    "mhw_msk = (nepac_ts > thr)"
   ]
  },
  {
   "attachments": {
    "b26ea220-4576-47ef-b3e9-8c9e9fa81b1a.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAy0AAAEmCAYAAABvSH/sAAABSGlDQ1BJQ0MgUHJvZmlsZQAAKJFjYGASSSwoyGFhYGDIzSspCnJ3UoiIjFJgf8rAzcDHIMwgyqCSmFxc4BgQ4ANUwgCjUcG3awyMIPqyLsisDuWI9w/KTnKdD/0Ve8ZYZDWmehTAlZJanAyk/wBxenJBUQkDA2MKkK1cXlIAYncA2SJFQEcB2XNA7HQIewOInQRhHwGrCQlyBrJvANkCyRmJQDMYXwDZOklI4ulIbKi9IMDrlJoXqBDgYmxkYBBKwL0kg5LUihIQ7ZxfUFmUmZ5RouAIDKVUBc+8ZD0dBSMDI0MGBlCYQ1R/DgKHJaPYPoRY/hIGBotvDAzMExFiSVMYGLa3MTBI3EKIqcxjYOBvYWDYdqggsSgR7gDGbyzFacZGEDaPPQMD693//z9rMDCwT2Rg+Dvx///fi////7sYaP5tBoYDlQAUPWGbpe5waAAAAFxlWElmTU0AKgAAAAgABAEGAAMAAAABAAIAAAESAAMAAAABAAEAAAEoAAMAAAABAAIAAIdpAAQAAAABAAAAPgAAAAAAAqACAAQAAAABAAADLaADAAQAAAABAAABJgAAAADCzJUvAAACC2lUWHRYTUw6Y29tLmFkb2JlLnhtcAAAAAAAPHg6eG1wbWV0YSB4bWxuczp4PSJhZG9iZTpuczptZXRhLyIgeDp4bXB0az0iWE1QIENvcmUgNS40LjAiPgogICA8cmRmOlJERiB4bWxuczpyZGY9Imh0dHA6Ly93d3cudzMub3JnLzE5OTkvMDIvMjItcmRmLXN5bnRheC1ucyMiPgogICAgICA8cmRmOkRlc2NyaXB0aW9uIHJkZjphYm91dD0iIgogICAgICAgICAgICB4bWxuczp0aWZmPSJodHRwOi8vbnMuYWRvYmUuY29tL3RpZmYvMS4wLyI+CiAgICAgICAgIDx0aWZmOk9yaWVudGF0aW9uPjE8L3RpZmY6T3JpZW50YXRpb24+CiAgICAgICAgIDx0aWZmOlBob3RvbWV0cmljSW50ZXJwcmV0YXRpb24+MjwvdGlmZjpQaG90b21ldHJpY0ludGVycHJldGF0aW9uPgogICAgICAgICA8dGlmZjpSZXNvbHV0aW9uVW5pdD4yPC90aWZmOlJlc29sdXRpb25Vbml0PgogICAgICAgICA8dGlmZjpDb21wcmVzc2lvbj4xPC90aWZmOkNvbXByZXNzaW9uPgogICAgICA8L3JkZjpEZXNjcmlwdGlvbj4KICAgPC9yZGY6UkRGPgo8L3g6eG1wbWV0YT4Kpl32MAAAQABJREFUeAHsfQd8FNX2/0nbhFSSQCAk9CA1FIEAUsxDQUBFURRFffoQFRWxPfvzPfX/3kMs71mxdx7F3kUBRYpI5wcooICUJEBCEtJ78p9zlzs7O5nZnd2d3Z2ZPffz2Z07t5x7zvdMuWfOLWEtQgAKhAAhQAgQAoQAIUAIEAKEACFACBgUgXCD8kVsEQKEACFACBAChAAhQAgQAoQAIcAQIKOFLgRCgBAgBAgBQoAQIAQIAUKAEDA0AmS0GFo9xBwhQAiEAgJr166F3r17axL17bffhjFjxqiWzc3Nhddff101nzIIAUKAECAECAEzIkBGixm1RjwTAoSApRAYO3Ys7Nu3z1IykTCEACFACBAChICeCJDRoieaRIsQIAQIAQ8RaGxs9LAGFScECAFCgBAgBEIPATJaQk/nJDEhQAj4iMDjjz8O06dPd6Jy++23w7x581jaW2+9BX379oWEhATo0aMHvPLKK2LZ1atXQ2ZmJixYsAA6duwIf/nLX4Cn8UJIv2fPnqx+v3794JNPPuFZ7IiLPt52222QlJQEffr0gVWrVjnlS0/efPNNxktycjKcd955cPjwYWm2z/HS0lK44IILoH379oBtYDwvL0+ki8PVHn74YRg9ejSTZ+LEiXDy5Ekx//PPP4f+/ftD27ZtAcvu2bNHzOvWrRs8+eSTMHDgQIiLi4Prr78eTpw4AZMnT2a0zj33XMD2KRAChAAhQAhYHwEyWqyvY5KQECAEdEbgyiuvhK+//hrKy8sZ5aamJnj//fdh5syZ7DwtLQ2+/PJLlo8GzJ133gnbtm0TuTh+/DiUlJQwA+LVV18V03kEDRac51JWVgb/+Mc/4Oqrr4Zjx47xbNi4cSMzhrDz/+ijj8Ill1zC6IkFTkc+/fRT+Pe//w0ff/wxFBUVAQ5DQ97VAhoOaj80pJRCc3MzM7zQGDpy5Ai0adMG5s6d61R08eLFgDgUFhZCfX09PPXUUyz/t99+Y/w888wzjL8pU6bAhRdeyMpwAh999BGsWLECsOwXX3zBDBaUCWXHtp977jlelI6EACFACBACFkaAjBYLK5dEIwQIAf8g0LVrVzjzzDMBjQIM33//PcTGxsLIkSPZ+fnnn888JWFhYXD22WcDehfQCOEhPDycGRvR0dGsk8/T+fGyyy6DTp06AZabMWMG9OrVCzZt2sSzAY2iO+64A6Kiolg+TuL/6quvxHweQQ/PAw88wDwtkZGR8OCDD8KOHTtUvS2nTp0Ctd/999/PyTodU1NT4dJLL2Xyo2fpoYcegh9//NGpDHqTzjjjDCbr5ZdfznjAAsuWLQPEasKECUyWv/71r1BTUwM//fSTWB89Sh06dICMjAxmdI0YMQKGDBkCiN20adNg+/btYlmKEAKEACFACLhGID8/H6qqqlwXMmguGS0GVQyxRQgQAsZGAL0qS5YsYUyiJ4F7WTDhm2++YQZMSkoK81ygV0Y6JAqHUsXExKgK+O6778LgwYNFr8fu3bud6mMHHg0iHtCIKigo4KfiEb0fOGyNe0+QHxxahi8tvUJ1dTXcdNNNgDwkJibCuHHjmOGD3icecBgcD2jcVVZWslPkGevxgEZa586dnfhDg4UH9OLIzzktXoaOhAAhQAgQAsoINDQ0wPDMHjAjwfHcVS5pzFQyWoypF+KKECAEDI4AekNwLgrO38A5J9xoqaurY54H9Brg/Av0XOCwJzQWeJAaHDyNH9HQuOGGG+CFF16A4uJiVn/AgAFO9dHokNLDYVnomZEHNADQ2yL1nqAn46yzzpIXZefx8fGg9sMhWUrh6aefZiuf4ZA1HC63Zs0aVkzKn1I9TEOepXNssM7Ro0eZV0WtDqUTAoQAIUAIeIcALoc/KCwedrdUs/eTd1SCV4uMluBhTy0TAoSAiRFAb0lubi6bz9G9e3c2BAvFwTkbaLhgPg7JQq/Ld999p1lSdNujUYP1MeBcEPS0SAPODcG5HPjV7IMPPmCT19Ewkoc5c+bA/Pnz4ZdffmFZOEcGy6sF9Fqo/XBomVKoqKhgw77Qm4PzdHCOjdaAQ8VwWBsuJICyoAGEw77UjCqtdKkcIUAIEAKEgDMC+IxdcMvdcGlEIpwXngjXp/d3LmCCMzJaTKAkYpEQIASMiQB6V1auXCl6WZBLnNeBBgV2yHE1LRw6NnXqVM0C4Gphd999N4waNYoNhdq1axdbeUtKAOd1/P7779CuXTs2h+TDDz8EnFsiDzjn47777oMrrriCDd1Cjw0aUXoGnFuD3hvkBef0TJo0STN5nIuzaNEithIa1seJ9viz2WyaaVBBQoAQIAQIAfcIoJdlQFgstA0Lh/Hhsab0toQJ7njHmAX3MlMJQoAQIAQIAUKAECAECAFCgBAwCQLoZellS4K/RXZgRguyvaKpCvJbGuHLZscS9EYXhzwtRtcQ8UcIEAKEACFACBAChAAhQAh4iYDUy8JJmNHbQp4Wrj06EgKEACFACBAChAAhQAgQAhZCQMnLwsUzm7eFPC1cc3QkBAgBQoAQIAQIAUKAECAELISAkpeFi2c2bwt5Wrjm6EgIEAKEACFACBAChAAhQAhYBAFXXhYuopm8LeRp4VqjIyFACBAChAAhQAgQAoQAIWARBFx5WbiInnhbcB+tP/3pT2yJ//79+8Ozzz7LyYhHXN9r3rx5kJWVBQMHDoRt27aJeb5GyGjxFUGqD7NmzYK0tDTA5VR5mDFjBtvRG3f17tatG4vzPNw3Ai9mXO7022+/5cnwww8/wLBhw9gSrZj42WefAS7ZygPW69WrFz9lS6NedNFF4jlFPEMAHyxcd9nZ2U6Vn3/+eaYffCjhkrk8cN316dOHdMdBCeKR60967z3yyCNsc0a89/D39ddfixxy/cnvPdwkE++9e++9l5XFe+/iiy9uVY8n4LLEnizjzOvR0RkBJf3xEk899RTbr+fkScfKPqQ/jk7wj0q6e/jhh1knDe+7iRMnQkFBgcgo6U6EgiIBQqCxsZHty3JxRILLFiOEfcG07tuCe4/hflp79uyBn3/+GV588UX49ddfnejjsvq4JD/+Xn31Vbj55pud8n05IaPFF/SoLkPguuuug+XLlzuhsWzZMtixYwf7XXrppXDJJZewfLy4ly5dyja7wzq33HILNDU1sd29cedu3E0bb7S9e/eyfSo2bNjA6mEHG+OJiYmAG+vh+U8//USb0Dmh7tkJbmCIupPv24HG4+effw47d+5kesI9QzDgBoWoVzxiR/jWW28l3XkGue6lle49bOTOO+8U7z++6aTavYflX3rpJVi7di3TJ957uLkjv/cwX3rv4Tnee6NHj8YoBR8QUNMffs1csWIFdOnSRaRO+hOhMERESXf33HMPe27iu++CCy6Axx57jPFKujOEykKOif/973/ivizuhOfeFulHEqU66enpcOaZZ7Is3JOsb9++kJ+f71QUP3r9+c9/Zh9dcO+uU6dOwbFjx5zKeHtCRou3yFE9EYFx48ZBSkqKeC6NoHHx/vvvw5VXXsmS8WLGje5w12vcRRw9Lps2bWJ5zc3N7CIPDw9nRgnuCJ6UlAT79+9n+XhjoPGDHSYM2JGinbMZFF7/yXWH+sIOLHpXUEcY0IuGAQ0Z9KCR7hgchviT688VU77ee/jxgd979MHAFdLa89T0h0bnE088wZ6HnBrpjyNhjKOS7vCjGg9VVVWi/kh3HBU6BhKB9957D+6/JBuGXnmG21/OzN5wyZld4a677mJed/S84w89JWrh0KFDsH37dsDNjqUB+2qdO3cWkzIzM1sZNmKmhxEyWjwEjIp7hgB+ve3QoYM4rEvtYsav/tdffz37eosdZ7TeMeCu4NhB2rdvH6OBVjsaK+iNQU/A8OHDPWOISrtF4LfffmNf3fFBlJubC5s3b2ZGpFx3GRkZ7EFEunMLacALvPDCC2yYCg5hKS0tZe3L9Sd9kcyePZt9AMAPB/zeww8C8nsPz+ne86868eMA3luDBg1yaoj05wSHYU8eeugh1mHDr9zc00K6M6y6LM0YPqsToiI1yxgfFcHeA1u2bAH+u/HGGxXrV1ZWAn7IeuaZZ9gIGGkh7MPJA/YT9AhktOiBItFQRWDJkiWilwULubqYzzvvPDZhC8dyY8CLHIegYEcJDRU0YNBI2bhxI7PucVx+TEwMK0t/+iGADzp05+J41QULFjDvCuqNdKcfxv6khOOHDxw4wIaHoSufD+/Toj8cq8wDv/fw/sN7Lycnh+49Do6fjtXV1fCvf/1L7OxKmyH9SdEwbhz1h8P7rrrqKsCPBxhId8bVF3HmOQK4IhkaLHiN86H/Uir4QQzvAR7y8vKgU6dO/NSnIxktPsFHlV0hgJ3fjz/+mHV6eTlPL2buaeEdJxxDWVNTAz/++CONqeeg6nxEHeECCGg0YkcVh+sVFxeDXHf49dDVg4h0p7NiNJJDz2ZERATT2w033CAOv5Trz92LhHtapPdebW0t4KR9ms+iURkeFkNj848//mBelm7CAiaoIxw/fvz48Vb3H+nPQ3ADXHzmzJnw0UcfsVbp3gsw+NSc3xBAAxxHxaBHHoeSKQVcpOXdd99lxjp+/MRh/vgBTY9ARoseKBINRQRWrlwJuMoUPrB5wIsZJ+LX1dWxlzOuLoEdY7XQr18/tgLLunXrYMiQIawYrsyCk/axU0xBXwTQUMFVo3AyPgYcKlZfXw/t2rWDCy+8kE3EJ93pi7ne1KQTHj/55BNxVT9v7z0c4im9915++WWaS6a30k7Tw1X8cKERHCuOP3x24nKhHTt2ZKu1efPsJP35SVkKZPF9xgMO88P3Hwa69zgqdDQ7AuvXrwecK/P999+LK8Tiwjz4XsAfBlz8pUePHmzOMn44W7hwoW5iax/spluTRMhqCOAke/z6iqtO4Ev20UcfZZY4vmD5BHwuMy6he/nllwMaI7h0Hi6Xh1+F1QJ2onFuRXl5OURFRbFiOK/ltddeo46TGmga0/GLCX4N5LrDiXO4XO5f/vIXpj9cRhcn3b/99tvM64Lnl112GdMd6gKHPpDuNILtp2JK9x7qE1cvwnunm/C1Hg18DN7ee2VlZeK9hx8KcGImLYChj0KV9IdfMZUC6U8JleClKekOO284/xK90127dhU7caS74OmJWtYXgTFjxigOd5S2gu8e7Nv5I4QJHZfWM2b80RLRJAQIAUKAECAECAFCgBAgBAgBvyOAC+k8mx4FiTZt/om39h2FTtfNhTlz5vidN28boOFh3iJH9QgBQoAQIAQIAUKAECAECAFCICAIkNESEJipEUKAECAECAFCgBAgBAgBQoAQ8BYBMlq8RY7qEQKEACFACBAChAAhQAgQAoRAQBAgoyUgMFMjhAAhQAgQAoQAIUAIEAKEACHgLQJktHiLHNUjBAgBQoAQIAQIAUKAECAECIGAIKBtSYGAsOJ9I2effTaUlpZ6T4BqBg0B3IASf7SzfdBU4HXDpDuvoTNERdKfIdTgNROkP6+hC3rFYOgO2ywpKYETJ06w5dBxc0Bcsh6XaMa9uDIyMiA5OZnleQNQU1MT29cLNxJMTExke7HFxcUB37sG28NtDuShqKiIbZ6K6T179oTY2Fh5EUOe48bKy5cvNyRvVmaq9RVkQmnRYFn15tMm5JxYRgQ2HTkFOV3aEhgyBOJ3rZOlGO90ta0b5NYfMh5jxJEmBEh/mmAybCGr6q8ye4xhMdeLMb3fe3ttw1yytvVwMssvPn4Q3vx/EyElNQ2ibRGwa9culp6fnw/4O3LkCOCeXd6E999/H26+4UYYnl8NlzTHQkRYExyL7ALv26ph95Gj8J/nn2P7gHHahw8fZntJ8fP9+/fzqOGPw4a5xtvwApiUQRoeZlLFEdvWRQCNFTMYLNbVAElGCBACwUSAnn+eo9+nfgvgTy0M7WofjZLasQfc8u8N0GngFXAwr5wVRyPlwQcfhLZt28IPP/zAdjO/55571EippuPG0T9v2QyH+naCv8eWwdGWBkgPi4Qwocapqkq46frZMHbsWBg4cCD79erVy2vPjioTlGFpBMhosbR6jS/cseIyqKisND6jAeDQbMZKXm0zlFfXBgAZasIfCJD+/IFq4GhaXX9WNlz8+d7TYrjEJbWHUefdBGdNuZVdsOPHj4d//vOfkJqaCjfccAP88ccfsHTpUrbJYPv27Zk35ODBg5oubjREFn/4AfxWXgL/jChmddLBPqinoaUZ1q1bB+Xl5WxIeENDg9Pu6tOnT4f58+dDfHw8DB8+XGzvlltuYUPX6urqxDSKhCYCZLSEpt4NI/X+giIoKjllGH4CzQg3VMz4gt5b3gRFZfYvdYHGjdrzHQHSn+8YBpNCKOjPjM9FLdeEv997WgwX5HPYn66FeY+vgdXrt0JGZhcoPVUGmd0HMBHy8vLglVdegZMnT7I5KjgnRmtIapsEccLclCbBKPk4vArG1kbB6+Fd4dzwFAgPCwMcFvb7jp0wPCoZ2kRGwbJly5iR8tFHHzGPT1VVFTOYcD7M/fffDy+99BKbdxMdHa2VBSpnUQTIaLGoYkks4yJgZkPFuKgSZ4QAIWBFBPjz0oqy+VMmrYZL935j4fand8LV930Bk655Ag7u2wZvvvkmoGcFJ9fjb+fOnYAeFK2hXbt2kCfMj2nTLgU+bC6B25qPwuzmw3BTRFu4OCwVbDYbnGqog80NpTCrOQXmzLpenIC/cOFC5n356aef4IwzzoAnn3ySNbtt2zax+bKyMmhubhbPKRI6CJDREjq6JkmDjAC9fIOsAGqeECAETIuAVb0u/lSIVsMlTPB+dOo+GIaMncnYmTVrFpvXgnNPcIUx9HCMGzcOsrKy4L333hNX+3LFO86P+W7lSugnGB7odRkXlQK/tdRDD4hiXhOse0ZkIuSEx8KDdcnQtGMvM2YeeOABuPTSS2HFihWMPK4siiuKcS9LS0sLm3vzxBNPuGqe8iyKgCVWD7OobkgsCyFAL1wLKZNEIQQIgaAgwJ+jobC6mF4Ao+GitrIYTs7nq4phexGRNphwyV1QfGI/pHUeBJtXvc7YQKMB56JguPfee6GmpgZOnXI/rDs7ewDs2L0bDuw/AA/ddy+8vXkrnCyrgA6QBAMabZBQ3wL45bxbWBTcXh0lxOLgZFUTrP90HfSOASiLsgnzJquxWejevTszXNAYwnDO+HPYkf5CCwEyWkJL3yRtEBDgL9ogNE1NEgKEACFgOQSkz1QyYNyrV6vhgh6XKdc5to8YOOpSeOL2wU4GyvHjx1mDOHG/f//+MG3aNFUGXn75ZbY3DK5MtlSYr8LDzz//DE899RTgHJZP4SSMGDEC0HOy+5xrIE1YbeyiiHiABoCpLTFgE6yaOc1HWFWciI/7zMyZMweG5zgm6nO6dLQ+AjQ8zPo6JgmDiID05RpENqhpQoAQIAQsiQA9Y7Wp1d1QMb4kspRaevdB8N/PW+Bf/7OvAsaHaGGZhx9+GD799FNp8VZxnED/0EMPtUofOXIkfPDBB3D06FFYv34925Ry6tSpcGvjEZgXdgLqm1pYHVwuOTUiEqb0tO8bg96WHTt2wAvPv9CKJiWEBgJktISGnknKICBAL9MggE5NEgKEQMghgM9a/gs54XUUWMlwQfKxCSnwt9cOwMy7/ie2dt5558Fbb73FznHSPg4b27t3r7iEMe50jxP4ced4pYBenczMTDjrrLPY3jDXXHMNK3asvgruDz/B6Ay98gzA3/W97UYLLsWMc2sioyLZymb19fVKpCnNwgiQ0WJh5ZpBtLHZWdCjS6YZWNXMY6i8PM9Ni4Ke6WmacaGCxkKA9GcsfXjKDelPGTEzPH+D9d5z5W3haKLhomS8pHboAThc7IGX9sDlt74EO3f/Bpdccgl8/PHHbHgXrvLVt29fSEhIYBtGpqXZ3w1ffvklJ616HDRoEAwYYF9qGb0wRxur4fLGP+CzQ8ehprEJhrRLhNdusA9Du+6sbPjHP/7BhogNHjxYlSZlWBMBmtNiTb2SVAFGAF+UFAgBQoAQIASMgYD0mUzzXhw64YaL2uR8XhINF+kkfZ6eltEH8Dcg52JY89lTcMtt90Bm/2ugd1JXqK85BS3QDFtWPAY43Av3eenYsSOv6vKIm1oWFhbCc889B5GRkYD7wty5Zger8/W918CFXTvB6r9dD2d2S4fyviPg0UdB0ypmLhulTNMhQEaL6VRGDBsBAekL0Qj8EA+EACFACBACyghIn9dkwNgxcjU5n6OoZrhgfkJyRzj/uqegY//7eXHYt3EhbPnuMbbPy8yZM8VlisUCLiLh4eFsnszo0aPhnHMcK4NlpqXClCfeE2uuWvQyZNsiYcOGDTBq1CjIycmB8ePHs00pExMTxXIUsSYCNDzMmno1jVTb9x+FfGE1EDMEfPHxnxn49TePm0oaIe9kib+bIfp+QoD05ydgA0SW9Ocd0EZ4hhvlvce9Lq6QVBoqpla+a/YM6NprJNz/wN/gur/Mhg8//JAtj6xWXikdPSw84CaUOWPOZqcZnbvCxPOnw4tfbYG04RPh8ccfZ+lbt26FBQsWwGuvvcar0dHCCJDRYmHlmkG0ypo6qKsT1jY0eJB+qTM4qwFjr7yxBeoaHC+YgDVMDemCAOlPFxiDRoT05xv0wXymG+m956vhMnR4O1ERMbGpMGbG55Bz8WLYdzwd7r7/ceiZ1RueeeYZNikfd7LH1cJwg0i1MHHiRFi7di3Lxon2n3zyCdx+++2Q3CELjuYXw8dL7HvHfPbZZ2zvFtxDBsOLL74Izc3NLE5/1kWAjBbr6pYk0wmBYL7cdBKByBAChAAhQAjIEKBnuwwQF6daPS64KljbtP7Qf8y9MG7mcuh79uPw5H8/hD+dM4XtZD9mzBi46qqr2Epjas116NCBZc2dOxdSUlJg4/4M2L1lFezZ8YNY5bLLLgNcTSw9PZ2lYVzqpRELUsRSCJDRYil1kjB6I0AvNb0RJXqEACFACBgHAXrGA2jxtrjSmNTbIi/XvvMoSEq7GLr0ewbOzP0UsgY+CFFRUZCbm8s2p7zzzjthxYoVIF2+uFevXrB69WrYvn07W4ns6C/LYfiEB5xI4z4vWO73339n6bfddhvgcDIK1kaAjBZr65ek8wEBepn5AB5VJQQIAULAJAjgsz7Un/daDBdX3hY1w+Xg//0P9my5B3asuQyO7X9MMEKi4NVXX4X8vHz473//y4Z03XXXXcyjgsskT5o0CaZPn86MmoN5iRCdMAQS218ALc1nw9DxX8HgcR9Cj/53wZIlS9gSy3iJoRGEE/O3bNlikiuO2PQWATJavEWO6hEChAAhQAgQAoSAZRAgw8V9p9+d4SI3Xrr0s++v8n//939w/z3XQljtx5CRkQEjRo5gc1f+9re/wa5du9jyxWh0oAfm+x93Qfe+syEj63ZI73E3xCX2Y9dYWFiEsLFkPKR0nMj2g/n8889ZekNDAzNYpkyZYplrkQRRRoCWPFbGhVIJAUKAECAE/IjA8dWbWlHvmJvTKo0SCIFAIiA1XGh5ZGXk0XBR2sOFl0bDZevmk+w0b6/dsBg4cCDgb968eXD48GHYsWMHvPPOO8xb0qNHD+ZpadOmDVRVVcHVV54Hn359HH7fdit07D5LmLjfBL9t/zsktu0OZwx9CepqCuCjj+yT9bER9LSg4VJUVMRZoKNFESCjxaKKNYtYHZMTobwpwnDsSl9chmPOIAxlxIRDqS1G4KbaIBwRG54gECz9KRkrnG9pHhkwHBXlY7D0p8yNNVOV3gN6GDJGfe+hFnGYmLuNJ7Ec97ioGS/c43IybzDEJyZhFTF07doV8HfRRRfBkSNH4JdffoETwtYH1dXVEBMTAz///DP07l4KKw8cgHD4LySntIdBgwYBemu2rJrE6Oz6SSTHDBZuuPz1r3+Fa6+9DrKzBzgKUMwyCJDRYhlVmlOQXplpUNp8yjDMK72kDMOcwRjpmxgBJ2zCZl71tFeLwVSjiR25/qQGgxoBXw0JLW3wtrGsr+1xWlY8yvVnRRmNKBN/R/hivBjtvSfHmc9v0cN4mXDRebDxsyjYt28f9O7dW94UdOnShf2kGbNmzZKeivF3330X1qxZA2+88YaY1q1bNzh06BAzXJKTk+Hpp59m81twaWUK1kOA5rRYT6ckkRcI4IuIv4y8qE5VCAFTI6DVmMByWstKAfGlnpQOxQkBoyAQCu8LbrxowRw9L9z7Ii2PO90PGDkNZt9wIxw7dkya5XH8z3/+M7z++uuQJ0ziLz5ZzOa/oMHCQ2lpKVxyySWwaNEinkRHiyHgsdHS1NQEQ4YMgQsuuIBBMWPGDBg8eDD7ocWLcaWAednZ2Sx/2LBhYpGCggIYP348cxNWVlay9EceeQRww6DCwkKxXHx8vBiniHUQqKiuhdr6uoALxI0Ufgw4AxZosKyhBWqFccQUzInA/jXboKrkVECNEHMiZUyu6f4Lvl68NVyC9d7zBjFPDBekr2S4TLpmAdSFdYD+A7Jh86bN3rDhVOeOO26H1HapbFlknoHzYbp3786MGjwGM1A/2X/oezw87Nlnn2UTp8rLyxlXy5YtE7m7++67ISnJeeyimClEfvjhB2jXzrF7KuY999xz8Pzzz8PBgweZdTxnzhxWBcuhm2/BggXsnP6sicCOA3lQUFYLkGXfTMqfUnr7gvEnT2amvbm0EfIjSgFSzSxF6PHOPSW/xKVBYUUlDPUBAqRFQ7h8ANCHqnT/+QCejlXxveLpULFAvvf0EJUbLlqGi2F7aLhI57rExiXDzL++Dz9++gScJyxpfPz4Ma/3VGlpaYEPP/zQSazo6GioqamBfv36AQ4RC3agfrL/NOCRpyUvLw+++uormD17diuO8EJ6//334corr2yV5yoBLVJ0H+IPafCAYxrRICopofHyHBM6AhvC5anxgeU9rUNYEwJWQ8DbIVrucOBGkLty3ub7m763fFE9QoAjECrvFzReuAHDZVc7Kg0XO/vieyG1Y3fWt1Or5yq9qbEJcnJarzBYV1cHHTp0YP1TPmLHFR1/5lm9n4x9c9xPZ8AA5YUOcFNQdF7wEViPPfaYrnB75Gm544474IknnoCKiopWTKxdu5ZdNLhDqVIICwuDiRMnst1Nb7rpJrjxxhtZsblz58I111zDhFy8eLFYFYeDIThosT766KNiulIEl8hbu2u/UhZkdWoP6al278+x4jLYX6C+JN7Y7CyRxvb9R6GyRnnYEq78gRPpMKCbF7+aqIXBPTMhIRZXWAL4Pa8QjpfaPVTy8vFtomFIVmcxWU0eLBAMmaIKj4q8DU+OhKSoMHa+p7wJ8mubxbyGNIcMWmQ6UlgClXVNgLpxpSdp++emRYlGyNo2XVX11K08H3CyKgYcSoFfJtWCK5mkdRIjwyAnxXHbrCxUHx7VR2g7U1hhC0OegNFeASu1gDLxsKmkEcobHQY8T8cjrhhkFJkqBB5P1VTDyiZlDMwokxX1hDJt3XUEQPCsSEO+LQ6qZc+4XW2SoTrccS1Ky7dvqIYe9fYhvJXhkfBLmxRH9uZDLB7XNYMd5ffTb7K2ecXY5gbIrhG8dafDRpVyccJ9ZvX7ydNr71B1E5SFtb7/jPSM8FQmMzz3VGVa9QPg+09LPwLfe/Utjm/G5utH2PtbxREd+a3Ljkmp6ZDe2T7Zvqa6HA7/tg3Ky9o4lfnTxGmwYsUKmDp1qjgyZ9eu3ZCfr9yPSkhIgNGjRzMaY8aOgcjISLjvvvtgxIgRbHNKXOoYQ0pKCmzbtg2Ki4sB+4+4Ktmvv/7K8pT+cBNLHnDSvlK/FvMzMjLFlcjKysrYJH9eT+lo1H6yEq/epF133XWA/XacX6QWxo4dC19++aVatk/pjt6XGzLIAFpXQ4cOdRpHyKvh7qSuvCx4UXTq1InNU5kwYQL06dMHxo0bx5a9w9UglAKu543WGg47cxWam5shH4cYKYT6iAo4WmXvAFYIc2aKVMph1U1HHKtY5ZdWQ12dcmcMl+jlK17hfAw2vEmhbUxqKSiDGJudt6KSKqioVOYzurYJGqTtu+AzkDKFNdbbJYtwDPtrtCVDjLAuOoYiWzmUN0hkKnYYZba4OLcyocHSJHjYfj2hrifGg6T91TZH5yvvRAnUNdiNkZZIG+OJ81wZ196+upWQWhvWYB/KxEq0/nMpk6R4tPDArLY5OmuHIxzzriTFWLQ6MhH2syWBQTBCaqEowoGNvKyTTJGCTC3KBhYuMcxW7BIIBFum4+EF0CTgeliiG6lcZpTJinoqKSuG4swsqWpYHA2W5ohI2DVqsph3rKwS6gTvt1KoibZBVZy9A1Lf0MSGlrUqJ1QNt9lAfj+dyoxuVRQToiOEjwpJjvmKhSVliuXCI2xg9fvJ02uvLKxcuP/CWt1/RnpGeCqTGZ57LmUS3n9a+hH43gsXPoDxsubtRxxi92tTmL0rGdZQBhH1xSytrr4WqkoLQFhfEhok2xr0yIiDmqTubLoALm+MAecvq3lIbMLzBA2K48ePw+TJk9meLI2NjWxKwaWXXgoYxwn4OFIHDRk0VPbs2QM4heHkyZOMvtLf8uXLxWT0jtTXn+7riKn2CPLGDara2lrAedhqwcj9ZDWePU3Hfrt08QNP6/taXrPRgkYH7j769ddfAyoOL4irr76azUPBi+bjjz+GrVu3qvKDBgsGNHymTZsGmzZtYkaLagUho23btjBz5kxYuHChq2JsZ9Qrxii7qpwrthVOM52T1M66YFmNQet8DA9o5mgu6z+ZVF3eLdXCMrenscFnjv25owCWUK6sSBzvqyQTepTQ4Dynn1QvzjK14qP+kKMtfCKKQWjPKSCfkuGFruZeaJZJaKDeYXzkuqKJ+5dwnPBOc1VWVSYngYQT48jUGC4YLHFJcH2qHPfTPJtQJjChnrwdPoVejULBmMne8I14kWWLMfeRoe6L2Oe7CM+H4/tabySpVF2tffu8GWvfT55eeytr7R8MWt9/xnlGeCqT0jVhTzORTGW7xXceqLzH+XvP6Z1o8n6E0nyXzn0dGuVzXCKEjyIPTIuEp556yu0HaUdtYKN0pOcYx/kraLBIA64cdtVVV0mT/BJ/8803W9E1cj+5FbN+TNiwYQPbVwf7/ajn/v3769aawz/phuT8+fOFZebymIW1dOlStuIXX1Zu5cqVzHOSmSnteDoI4vAt7nrD+Hfffac6Hs5Ryx6766674JVXXmHWtDyPzv2LgKeTC9W4aWV0qBWkdEKAENCEAJ+f4q3BoqkRHQpx/nyZrO9LXR1EIBKEgMcIhOI7T+tcl3DBu5pzztWAm0CWFEs+KrpBGec4/+Uvf4EHH3xQLMkNFpx+gAE9LRdeeKGYH+iIEfvJHcYOZR+P8Dnq7pfQszOsW7cOcIVf/nv11Vc9gvHMM8+Ew4cPs41Ab7vtNrj44os9qu+usGajxRUhNGLkQ8PQhTZlyhRWDXc6HTNmDLO8cBLV+eefD9LxhK5o4ypi6JnBiVYUzIsAPsT5z7xSEOeEgHsEpAaFNO6+prYS3BDQVjr4pTi/ZHwEXxfEASFgBASuvP09GDzqfHhv0Xua2bn88sshThhy/u9//1usg+fp6elsGgEOzcIP4omJTsMvxLLBjpiln4x99S1btog/Pv9cK36IP9+iBG0AnHPkapieVrq8XJgwDlB5xi8vYYLjwIEDYdWbT5uAU/Ox6K8vRtyLgxMQtwvzfsapuMb91b75NGE8jnFxgw22TJjUkm885oLAEe+cu2ra1467ljZctS/Nw8n0+84cD0O3fCdN9lucy+6JDLyO35gyMWG6/4yvPP6ek3Pq7r0nL2+mc6VhYsg/Hx7GZVn31XOQv/N9WPTuO9CjZw/F4V+8LB5x7nIEzoGTBJzv0rdvX/jpp5/Y3n6SLL9H0ROBnXsjh9zcXFh6+UhIEhZ60hJeWLkJEsZeBHzrEbU6OKcF92rcvXt3qyI49whXckPvF04DmT59OvO8cG9YqwoeJujiafGwTSpOCIirf+HKajE2bTcUwWYsBHAFOb4gg7E4Cyw32AnX2hH3pKxcCq1tyOupncc3N4ItyrkToFZWj3TOPxoiWowRLWX04MusNOj+M77m1D66heJ7T77pZEqH7vDzhvWQ1SuLbXmBw4i2b9+uqlTcFuOXX35hq8qOHDmSGTG4iBOuGIabkVMIDAI4qmrUqFGwb98+wCkhb7zxBrz88svshxzgHjq4HPKgQYMAF9NCD5NeBgvS1zwRHwtTCD0E8EuR2oPXVzREukmtF1EQ83xthOoTAn5EgHfE/dgEIx2odgIhBzdG+BHb5PJJ0/zNC9EnBAKBAL7L1DwugWjfqG20yFYp/Oyzz6BHjx4wZMgQVZZx80jsJPOAc1goBBYBXCnYVcDlkPHnr0CeFn8hS3Q1IYD7vBSfcCwhiA94Mlg0QRf0Qqi7ImEVwVANvKPtjfye1PWkrCe8HLTFQ3FVjSdVdCmrJA8aK2SweAZvqN9/nqFlrNK4Z1tRifOqV8bi0HtutE7I75dzEaSmdWINnXPOOWzJ4f/85z/eN0w1QwIBMlpCQs2+CenPr0S4MWV5VS1N0vdNRUGpzXUXlMaD2Ch2upU63p6ypIWGljKetsvLF0XFQkVdPT+lo8kQCNX7z2RqYuzKP8ThJtMVlVVmFMUnnqVDxHDI0NgL7mCT61etWgU4PwWNFgtMs/YJI6rsGgEyWlzjQ7mEACFACIgI6G1EID0lmmrpIiMUIQQIAVMhIDdcTMW8h8xq9bYM+9N1gNtg8IAbiePGkBQIATUEaE6LGjKUTggQAoTAaQSUDAs9wfE3fT151YMWykvDwfRAkmgQAsZEAA0XpZXE0NvCVxKLiLI5MV9TUwMxMaq7VTuVpZPQRICMltDUu8dS+3NCvsfMUAVCwI8IhJoB4UcoXZImw8UlPJRpQQTQ2+LP4dZmgYwbLhGRzkYL7ulBRotZtBgcPsloCQ7u1CohQAgYCAEyVIKjDDnu5H0Jjh6oVULAHwioeVuwLTRcNvxm3+4AlzPGfVi++OILmDlzpj9YIZoWQYDmtFhEkSQGIUAIEAJmRwCNGOnP7PIQ/4SAHIFQmtuCsrua35IWto3BgwbLI488As899xykpKTA7NmzYejQoXDbbbfJ4aPzEEeAPC0hfgEEW/zEyDCIjqTLMNh68KZ90p03qBmnTmxzA0TLdpg2Dnd2TqSeGPLCOGuH7j9nPMx0llqSBydj083EMuNVyeDyZbhbp8xuwsaSkYKXpZEZLdhIUlIS24sFVxfDSfnLli2DCRMmAG4oSUaM6S4Z3RkmT4vukFqToNLDSg9Jc1IiIbNdih6kiEaAEbCK7qQd4wBDGNTmsmtKIT0pPqg8eNI498B4UsfKZa1y/1lZR2qyoe4yOnRQyzZcOr7/1foAaulSIdS8LW3axMF3G49Ki8IVV1zBDJi6ujrIysqCoqIiWLx4MeA5BUKAPnHTNeAWAS0PJbdEqAAhQAgQAjogIDUyyfuiA6BEghAIIgI2WzScO2kanJWTDTNmzADc9Z6HnTt38igdCQGGAHlaTHIh8C8d/BgItgPZViDkoTYIAULAWghIDRhrSUbSWB2B8JpKy4io5cOmmrfl6NEDsHL5JzBu3Dgng8Uy4JAguiJAnhZd4dSfmNrDgKf7Mp5Uzi2nKU/35/nKwgY4HFEIuan+bIVo+wMB0p0/UA0czY1xaVBYUgbZgWvSLy1xw8VbrwuvL2fOW3pyOv46p/vPX8j6ny7X3bDemf5vzMcW/N0vyMzsAenp6ZCZaXwsfISSquuAABktOoDoLxJaHhZYxhvDRQttf8lFdAkBQoAQ0BsBbnxoMTZ4WVc8YBkttFzRoDxCwOwIYP9Cr/6C0hLIlZVlcOzYMVi3bh307t3b7HAR/35GgIwWPwPsLXlPHhK8rNx44ene8kD1CAGjIMA7mdJOpFKaUfglPoKHgNp1wdM94YzXkV53ntSnsoRAKCAg73t4InNCQltWHJc5njVrFuCqYRQIATUEyGhRQ8aE6WSkmFBpxLIqArzDKC3gLo06l1K0QjuudK2ENiIkPSHgHQKuvC2+GCxSbkaMGEEGixQQiisiQBPxFWEJbiIZH8HFn1oPLgLY2fS2w+lL3eBKTa0bGQFvr0cjy0S8BR8BM73r9TJO5Kg3NzexpI0bN8qz6JwQaIUAeVpaQUIJhAAhECwE9OocIh3yugRLi9Zsl64pa+qVpNKOgD8Ml7raGsbAeeedp50RKhmyCPjF09LU1ARDhgyBCy64gAH7yCOPQEZGBgwePJj9vv76axHwe+65B4YNGwY//vgjSzt06BBzET7//PNimblz58Lbb78tnlOEECAErIeAXgYLR0ZvepwuHQkBQoAQ0AsBM3lb9JJZSsdWe4Cdvvvuu9JkU8epD+w/9fnFaHn22Wehb9++TlzfeeedsGPHDvabMmUKy9u7dy87rlmzBl588UWxfFpaGiCN+vp6MS2UIv74mmFU/PokRkD7pESjskd8uUBAT935y8BwR9ddvgvxTZ/Vra4CUmPbmF6OUBVAz/svVDEMltykOwfyo0aNYidJSUmORJPHqA/sPwXqbrTk5eXBV199BbgShLuA1mh4eDjzrLS0tIjF27dvD+eccw688847YlqoRULFcMmMCYfE2JhQU68l5NVLd/42HNToq6VbQjkahOjQWAMJMTYNJakIR8BI14xe9x+XjY6BQ0Cuu1Dxtuy1DWsFMi53jAaLzWaNZxH1gVupWNcE3ee03HHHHfDEE09ARUWFE6MvvPACoPsPh4I9/fTTkJycDP3794fq6moYM2YMPPnkk07l77//fpg8eTJbAs8pQ+GkqqoK1u7ar5ADkNWpPaSn2i34Y8VlsL+gSLEcJo7NzhLztu8/CpU1deK5NNIxORF6ZaaxpIrqWthxIE+a7RQf3DMTEk53yn/PK4TjpeVO+fwkvk00DMnqzE9Py9MRogqPimkYaUjr7CTTqc1rYG+5fSKbU8HTJ+emRYnJm0oaobzRYRyKGUIkQzAe+gpeDwxlDS2wubSRxZX+hidHQlKUfVnCPULb+bXNSsUgMTIMclIclxhuqKUUcHPJrLhmwAc5hjyBntllQjnwa5rVZTpQXAiNTcp61XLtVR3Oh/a2eOhRb98dujI8En5pk4LwKYb+NSUQ32y/Ng8K9YqiYhXLxTY3QHZNqZj32eZDYlyMCJsrYkCPA3bgMZyIbAOHohNYXOlvRFWhmLyrTTJUhzvuLzFDiLRvqPa7TLg5pFrQKpN8c0kryISY+FtPejz3kE9fnxH47JTff8F8lushE9KQByvKVJoo9AVOf69j79xVP7D3u1x2PPetH6FEMTh9o+II+3MeOep6xpnwa1FX+Hbj04DzWb799ltFRhMSEmD06NFi3vLly8W4PNKvXz/o0qULSz5y5Aj8+uuv8iLi+aRJk8T4+vXrW/VZeWZGRiZkZw9gp2VlZbBhwwaepXgMRh9YkRGLJjp6lDoI+OWXXwIO7Ro6dCisXr1apHjzzTfDww8/zDwqeLz77rvhzTffZPnSuStiBSHSvXt3yMnJgcWLF0uTFePNzc2QX1armFcfUQFHq+wd9YrKSihSKYeVNx05JdLIL62Gujrlzlh5UwSUNtvL1tbXQYELmi0FZRBjs/NWVFIFFZXKfEbXNkGDtH1OM7q9yBOLCOlSmSojO0FRhLIhhOVX2xwdm7zIEqhrUTZGSm0xcMJmH6ZVG9YA+RGODp8zAwCNtmSIibJ31ops5VDeoCJTZCRU2xwdUHzBykN9QyM0tDTDNkGO/QIPGMoba00tE5exOjLR0jIdqzgFVYKuDkcqu/W1XHvNmW2gJtoGVXH2YUr1DU1QWOF4sXEs+TEyIR5sUXbjuriqBirqlIeQRkcIZZLieTW287t4Ios0CUOkCk97HCpq66G42m7AyIqx010pQ8XkY2WVUCd4i5WCGWSqb2yCprAI2DVqsiiC2WXigvhTT3ttCeDrc4/z6cszoqK+lj075fdfsJ7leshktPeTv2TC916E8HFudWI31oT4zi0uh5bI1h4Hn/sRXBDJUdqP8HffqCmMdzXtH3gbhD5UeNh2SIneA206d4bU1FQ4ePCghDtHFD0w0o/gauWwRnl5uWioYPzkyZMOQrKY1PhB74jadITCwkLIz7fzXVtbCwUFBTJKjtNg9YEdHFg/FiYMy1L+9O6F7A888AC89957ECl0VlG5eNFccsklsGjRIpEaTrTHCfq7d+8W06QRaT7OeZk+fTqMGzeOGTDXXXedtKgYHzhwIKx682nxPNQiZnYto/flcEQSXJ9aHWpqM728vujOSMNsTK8ILwVAT01hZhZcuO8nLymEZjWjrErny/0XmpozjtTudGel4eHyIWFbDyeLisAPzgvvPgMWLFgAl156qZhuhgiOGtqyZYsTq8HqAzsxITnJzc2FpZePhCRhJI+W8MLKTZAw9iKYM2eOluJBKaPrnJb58+cDWqxoeCxduhTGjx/PDBYcs8jDJ598AgMG2F1tPE3t2KdPH0B3H1qvFAgBQsAaCJDBYg09khSEACHgHwTM/CFSiogrgwXLrf9hF/Nc4MdtKwTqA/tfi9xn59eW7r33XrZqWFhYGHTr1g1eeeUVze099NBDbPlkzRWoICFACBgOATJUDKcSYogQIAQMjAAaLmb2uLgzWLZuFoZuCX3CmpoaKCoqYlMLDKwOn1ijPrBP8DlV9pvRgm4p/GHAIWNaAxo10qFjgwYNAnQhUiAECAFzIUCGirn0RdwSAoSAsRAwq+GiyWARoG45vagKDqt64403jAW+j9xQH9hHAFWq6zo8TKUNSvYzAmb+GuNnaIh8kBAggyVIwFOzfkeArm2/Q0wNSBAw21AxucEiEYVFmYfldGJF8UFISUmBl19+WV6MzgkBRQTIaFGEhRIJAULAWwSoU+ctclSPECAECIHWCJjFcFEyWKQT76UGC0p5YNuPUFJSIq4m21pySiEEnBEgo8UZD9OekbfFtKqzFONksFhKnSSMCgJ0nasAQ8mEgAQBqcEiSWbRvT9thcaGQpg2bRpce+218mw6JwQUESCjRREWcyaa0XDBDQh7pjv2kjEn8qHJNenO3HrHDRi7pSjvsWNuyQLDfbANF7r/AqNnf7Tije7M4m2R4jW0q/p+b1iupakERo0aBTExp3fZlFamOCGggAAZLQqgUBIhQAgQAoQAIeAOgWAbLu74o3xCwNAItFRDu3btDM0iMWcsBMhoMZY+fObGjN4Wn4UmAoZAgDpwhlADMRFgBOi6DzDg1JxhEehT77zZoitG+5w1FOprCqFLly6uilEeIeCEABktTnAE7wRdv/Kft9yYyXDZVNIIeSdLvBWV6gURAanuqOMWREV42fSuNslwrKzSy9pUTYpAMK5/6f0n5YXixkcgVHU3dLizV6WivBDi4uKMrzDi0DAIkNESRFVIjRQlNtzlK9XhaWYxXMobW6CuoZGzTUeTIICdtGMH8qCmqhqC0WEzCUyGZrM6PArqmpoMzaOZmAv0fUDPTjNdHc68eqM7s7zT5d4WtXktDfVV0NTcADk5Oc7g0JmhEZg1axbbCHTAgAGKfLa0tMC8efMgKysLBg4cCNu2bVMs522i3zaX9JahUKmHBokngZc3y4PLE9morHERCHRHzLhIEGeEgHsE8H7pmEudMPdIUYlQR6AkfwvU19XAiRMnID09PdThMI381113HcydOxf+/Oc/K/L8zTffwO+//85+GzduhJtvvhnwqFcgT4teSGqkw70nGou3KuZJfTJwWsFHCRoRwM4XGSwawaJihIAEAbpvJGBQlBCQIcCHiG1dcR/LGTRokKwEnRoZgXHjxrENQdV4/Oyzz5hBExYWBiNHjoRTp07BsWPH1Ip7nE6eFo8h874C95Z4T8FRE2lpMUqwjJ7tOjigmFURoE6XVTVLchEChIAZEdDyrjebXDERNXDvvffCY489ZjbWTcVvVb+REBGvbd5Q/S8nYd26dfD666+LMt54442AP60hPz8fOnfuLBbPzMwETNPLm0ZGiwit+SLcGHH3QCPDxXy6JY4JAULAvAig4U/DxMyrP+Lcfwhs3XwS7p9xApY9XgQPP/wwREdH+68xouwxAmPGjIFFixZ5XI9XwDkt8oBeF70CDQ/TC8kg0kHjhRswQWSDmiYECAFCgBAgBAgBHRFw91FSx6b8Qmrr4eRWdHESd0NDA0yYMKFVHiWYGwH0rBw9elQUIi8vDzp16iSe+xoho8VXBDXWD4RR4aoNoz74MmLCITGOdsPVeBkZqlj7hmpIiLYZiidiRjsCpD/tWHlT0t/DLOnZ6Y1WjFHHyrrbaxumCjJ6WXav/ruY//PPPzPjRUygiOkRmDp1Krz77ruAHhfUb1JSkm5DwxAcGh4WoEskUEO0uOFiVCNFDnffxAg4YUsEqKe9WuTYGP28R30lVMW1MTqbxJ8KAqQ/FWBMkkzPTpMoSoFNrbozy3tcQUSWJPWyoMGCoez4NjZBe8WKFRAfH8/S6M88CFx55ZWwevVqOHnyJKBX5dFHHxUNzzlz5sCUKVPg66+/Zksex8bGwltvvaWrcGS06AqncYjJjRd+bhwOiRNCgBAgBAgBQoAQUELAjAaLmpeFGywtzc1QkLcbfly1hwwWJaWbIG3JkiUuucT5Ky+++KLLMr5kktHiC3oe1g2Ut0XKltGNlbKGFqgNa5CyTPEgIuDJkJbK8Eiob6DNCYOoLp+aJv35BF/QK9OzM+gq8JqBUNCd1MvCgSo9/n+QmtrOaXUpnkdHQkALAjSnRQtKOpYx49cTHcVvRWpzaSPknyxtlU4JgUfAE4MFufulTQoUVFQGnlFqURcESH+6wBg0IvTsDBr0PjfsTndm7CdIvSxSg4V7WRC0Xd8/AxdPPd9n/IhA6CJAnpYg6D4YHpcgiElNmgQBT40Vd2JJl3rVm7a7timfECAECAFCILAISA0WactSg2Xbt0uhvGQ93HffS9IiFCcEPELAY09LU1MTDBkyBC644ALW0D333AN9+vSBgQMHwrRp09jul0ocdOvWDbKzs2Hw4MEwbJhjdYmCggIYP348XHTRRVBZaf9q+8gjjwBO4CksLBRJWW3CFhouZvyaIiqEIqZDAA0IpZ+egkgNFqSL5/I0PdsjWoQAIUAIWAkBs/UL5AaL1MvC9bL3p61QWfozlJWVhcTQMOonc83rf/TYaHn22Wehb9++Iie4zvbu3bth586dcMYZZ8D8+fPFPHnkhx9+gB07dsCWLVvErOeeew6ef/55mD17ttOGNu3atYOnn35aLGfVCDdezPagsqo+rCaX1Ejxt2yujBNXef7mi+gTAoQAIUAI6I+A3GCRtiD1smB6SqdLIDGpPcybd7u0mCXj1E/2n1o9Mlpwk5ivvvqKGRicpYkTJ0JkpH2U2ciRIwHLeBLQIg0PD2c/6U6as2bNgmXLlkFJSegshUsGjCdXDpV1hQA3VlyV0TOPjBI90SRaVkCA7gkraDGwMpjp46WSwaLmZUEUIyPjoGv/x+HtdxfBr7/+GlhgA9ga9ZP9C7ZHc1ruuOMOeOKJJ6CiokKRqzfffBNmzJihmIfLoKGBg8ebbroJbrzxRlZu7ty5cM0117ANaBYvXizWxeFgaLigxYrrQLsKVVVVsHbXfsUiWZ3aQ3pqEss7VlwG+wuKFMth4tjsLJaHK25tKmmE8sYWp7INaZ3ZecfkROiVmcbiFdW1sOOAuqE2uGcmJMTaN0/8Pa8QjpeWO9HkJ/FtomFIlp0+pn0DHVlWVKFjZ1Feto+wt0mmsCkjhrzaZthbrr6C07lpUbyaokw8Eze7wrXjMeDKJjhRUC0MT46EpKgwlr1HaDtf4EEpJEaGQU6K4xJbWXwAWqcAAEAASURBVNh6lbBD1U1QFlYNeXHNlpGJYxEsPVUdzgeIS4P+NSUQ32zX40FbPBRFxXLWnI6xzQ2QXeNYDGGjUFctdKurgA6NNSy7Slg9rLy+EXb3GwK7FXTb6toTyjHeZMRxk0PcMwQDrmiFE8TVgr9lOhHZBg5FJ6g1DyOqHENWd7VJhupwx/0lrWQGmfJtcVBdUydlG8wuExfGSHpSeu5xPn15RvBn58om5+dqsJ7lesgkf+dymlaTCXVXY6vn4onv3AaVfoy3/Qi1fhE27E3fCOt98wdAXQ0+r1fjqVOobe4FMakDWFp4SwXENm6HLh2PO8oIXaoUYfPJDRs2QEZGBuv3YeauXbshP1+5H5WQkACjR48WaSxfvlyMyyP9+vWDLl26sOQjR464NI4mTZokVl+/fr1qvzYjI1OY2mCXCYe3Ie+uglH7ya54NlOeo0fphusvv/wS0tLSYOjQoWxjGXnxf/3rX8zjctVVV8mz2DleFJ06dWLzVHBIGc6DGTduHHTt2hXWrFmjWGfevHlsDszdd9+tmM8Tm4W1v/PLavmp07E+ogKOVtmNjwphzkyRpFxYo+OhgZW27Dt909i6QV5kCdS1yDruxXaDo7wpAkqbT7F2auvroEBC06lx4aSloAxibHbeikqqoKJSmc/o2iZoOGKniTREeaLbg5zP6shE2G+zG0LljbVQFKFsCCGd1TZHB1RRJiwkhFKBHtvkUYjjEsT5EY5OLCsg+Wu0JUNMlL2zVmQrh/IGFZkED1y1zdEBPRzh6PBxcmiwNAmG7LbITpaRicsWLD01Z9o3fIxMiAdblN0QLa6qgYo65+ud8xkdIZRJiuenUFhSJsblkabYNlAYY2PJZaUV0FTfAIcj7B8F5GWVrj3Om7RsTbRN3KQSl1AudLEimb9lqqith+Jqu1Em5ZHHd6UM5VE4VlYJdYKnWCmYQSY0WJojImHXqMmiCGaXiQsSbD1FZabD3tPMKD33OJ++PCPKwsrZs1N+/wXrWa6HTK3euaeJWk0mfO+Fh4UL7+duTEJ85+ZFVwlfDJXfpT73I7hyJEdXfSNJMRbdJPRNasPsH70qT+ZDvdDvUQrxMRWQWGPvRw3rLPR3SgugY6yjH9XSIvRzug2Huro6wOkCMTH2fgzOX+ZzmuV0bTabk0Fx8OBBeRHxvLy8XDRUMI4bIKoFqfGD3pH6euX3I/LGDara2lrAedhqwcj9ZDWezZau2WhBo+Pzzz9nO12i4vCCuPrqq9k8lHfeeQdQWatWrWKeFCUQ0GDBgIYPTtjftGkTM1qUyvK0tm3bwsyZM2HhwoU8SfGIlvjsJEeH36lQk5Au7YMp96/sVeoPOaoKm7Srh2qo7NLdkZ3VwRF3FevS1lWuU16OrKzzfivVwg7yp4ujBlOdqjqfeCCT0670rmi2SNrHZ479uePcLj+rdxhUuQo0yxJbYIMtAya1CN4Bi8jERT++ejWPOh0Vh43opCccFuZryHZDgPOfk4q6y7TrTqmOgkxa+HOYBUpEvUtzJ5M3VD2haUSZ0Ku178zxkL3hG1F8s8skCiKJBFomdn9Irn2l556DPcmz1MNneVkbybPTQVCIIU3JsGqF565YXMdnuUiTtX/6zEOZHDTkMWvJZH/vtYfc09eJR8PCZH0DOVLSc3k/QprnHMe+SaZz0ukzPgyMd4nSeioWAxwahuYJv/K27rYbDDgJn4fa6qOQt2e+MOz/pGpfkZc1+hFHFsmDkfvJcl7Neq55TgtOsEdr9NChQ7B06VK24teiRYsArdUFCxYwgwZX/FIKOHyLDynD+HfffQcDBtjdbUrlpWl33XUXvPLKK9DY6LDWpfnBijsbEcHiwvzt4jAz7rUxvzQOCVx1zl3lOSgYM8YNFuTOG91J6xtTwtDhCocOck9c6EjtX0kDeX17c//5V3qirhUBqe48Mli0NuBDOTRSpD8fSLGqfc5yfLKxCSNHyspKYePGjb6SNWR96if7Xy2ajRY1VnBOChokOOQLlzOeM2cOK4outClTprD4iRMnYMyYMTBo0CDIycmB888/H6TjCdVoYzquIoaeGXQnUiAECIHgIaBXh0wvOsFDglomBJwRwGuarmtnTOjMPQJGMFikBgr3qrjn3LsS4RExkN5lIrz11lveETBpLeon66e4MGHFLvuED/1oBpwS7hGz8c5pAW8XGwzkQ8eK3h2cyL/b1h4ui+FO5aCoUfdG3XlT/NXBcdeuN4Kq8eqL7vzBpzeyhXIdXJzhSJ9hkLtzdSjD4LPsaveHz4TdEPDl/nNDmrL9jADX3eTh/fzckp28v40RLoR89TDpssd8mFh5ySY4efgZWLJkCVucidc12xH3G5Ru32FE/nNzc+GNR+6AxPg4Tey9uvRTSM3KFp0PmioFuJDPnpYA80vNWQwBXHmsvEp58qHFRA2IOIHsQPmiu0DyGRDgTdgIriantjiDCcUJCsvBvI59uf+CAhY1KiKAuiurU151UyykUyRQBos7dvkwscSUHPjPf/4DuGgTThWgQAh4ggAZLZ6gRWUJgRBDwJ+dMqTtT/p6qspMvOopN9FSR8As1666BJRjdQSCbbAMHd7OCWJuuFx77bUwYsQIOHr0qFM+nRAC7hDQvHqYO0KU718ErDg0zL+IEXWzIOCu8xfsoWRS/jAebH7MolfikxAgBJQRYHu+qSxvrFzD89RgGCxDu5ayVcSk3KLhIh0mhobL9Lv3sI3KH3vsMWlRihMCbhEgT4tbiKgAIeA5AtKOrrQ2pqvlScv5EteLvl50fJEF63LMpEdfafpS3yi4+CID1SUECAFCwB8IoOEiD3KPS0tTPbRN6cgWa2pocN4cVV6XzgkBKQJktEjRMGicvCwGVYwbtqSdbB53U0W3bKt3rDme/pZTjb5aum4KJEKGR4A8boZXUcgyGAwvixRsd4ZLbGIGXHDrLgiLSoIvvvhCWpXihIBLBCxjtOAqXoFcycslqjpmksGiI5ghRsqXjrUvdQMNM/LKf3q27Q4Dd/l68kK0jIkAGi7SnzG5JK6MhkBUYWjO5ZB7XJIzxsI///lPo6mH+DEwApab06JkuJit4282fn25vhMjwyA60nKXoS+Q6FoXO9aefBH2pCNuRN1J+fdEbjnoUjryPOm5p/hK63oS5/z4IpO8vdjmBoiOiJAn07kPCEj1w3XmAzmXVY14/7lkmDJFBJjuoqPEc70jfeq3sA0i9abrCT2l+S1YnxsuOM8la8i1sOHDzzwhS2VDHIGQ6C1KDRkjGwRG5s1f90lOSiRU21IA6sv91UTI0+WdJ2mHSg4KLyNPd3VudN1xmVzJrSQfr6eUF6g0NR6k6Z7KJec9u0YYe54UL0+mc50QQP1I9aUTWZGM0e8/kVGKtELArjtte2e0qmyiBDXDBUVA42Xllw0QH0/PIBOpNOishoTRIkVZbwMGDQ0pTWlbWuKhaKhowYXK6I+APztQ+nOrH0Uut5ZOPi/rSetYRwttLTQ9aZ+X1attLfxRGc8Q8Lfh4hk3VJoQMBYC0bALBg8ebCymiBtDI2CZOS3eoIzGhvTnDQ1v66CxQgaLt+hRPULAcwR4J1+tprt8tXqY7ktdTlcPGpwWHQkBQiC0EcAhYkYISpPyka+WlhbYuXYRW0HMCHwSD+ZAIOQ8La7UIveY+MOo8AdNVzIZPW9lYQMcjiiE3FSjc0r8yREwo+64YcC9E/xcLlugz33hA+tyeTzhe2NcGhSWlEG2J5WorGEQMOP9ZxjwgswI192w3plB5iR4zRcf3w/792yBCy74PHhMUMumQyCkPS3utCX1wvC4vI7c0JHn83PyrHAk6EgIGAcB7Ox70+FXksAbwwPreFNP3r4eNOQ06ZwQIAT8h0BLpM1/xA1EeevhZEVuUjtmQafOWbBs2TLFfEokBJQQIE+LEiou0rQaKUiCvCougKQsQiDEESBDI8QvABI/ZBFoSOsMUFZrOPn5kLJA7PPS2FAHBUf3w9ChQw2HAzFkXATI0+In3ZDB4idgiSwh4AMCap6VQHpb9PKu+AADVSUECIEgIeDJh09fWeRGiDs6WE5aVn7urr5avpqXBcv/a7ZguAlh/vz57Eh/hIAWBMho0YISlSEECAHTI+DOMHGXrxUAVx4UV3la6auV8ydttTYp3TUCpBPX+IRKLh9eHkiDRSu2UmNFXsdVnrysJ+dff/AxlJ06yaosXLjQk6pUNsQRoOFhfrgAyMviB1CJJCHgAwJaDRIsp0dHE2nwNvWg54PoVDUICJDOgwC6wZo0ooEih0iLUYJlvBkupuZl2bKpCDZ9eQtjZefOXdCtWzc5W3ROCKgiQEaLKjSUQQgQAoSA9wgEo+MqNZa855xqIgJq+uPGqBJKanWUylKaMRFAY8ObD49mMFKkiGsxWHh5aVlvDBhOZ+vmk3B49zJI75AMJ44dhujoaJ5FR0JAEwJktGiCiQr5C4E+iRFQHZkokK/2VxNE108ImEV3rjqZcmis0OnUarh0q6uAptg2cgi8OlfCzRPcvWo0SJWUZA0GK2a5/4KBjTdtSo0OaZzTUjJklMrx8q6OWZ3aQ31EhasiuuV56ylRY4AbMN4aL3m//g/uvWMeGSxqAFO6SwTIaHEJD2X6G4HMmHDYb4sBqPd3S0RfbwRId3ojqh89LYbLoDHZUGpLdfIoeGJouOu883xPaOqHgG+UOO++UfFvbbr/9MNXi/GhpYxWjtJTk+BoVYvW4rqX44aHL4Q5DTXjBTeVVBoiFh7RBg4dPuJL01Q3hBHwy0T8pqYmGDJkiLBp0AUM2pKSEpgwYQL06tWLHUtLS0XI77nnHhg2bBj8+OOPLO3QoUMQFhYGzz//vFhm7ty58Pbbb4vnFCEECAFCgBBwjQB2vF39lGorlZeX42Xk6WrnWN5MwWz8mglbo/GKhoiexojR5OP8oIEh/fF0PY7ceNFKK637JFi+/DutxU1ZjvrA/lObX4yWZ599Fvr27Sty/fjjj8M555wDv//+OzviOYa9e/ey45o1a+DFF19kcfxLS0sDpFFfT5/fRVAsGsmrbYbyauOtV29RuHUVy2q6C7XOqlb9cSOFH3W9iIiY1who1Z/XDVi8YjCNlWPFZVBRWWkZhLlBJBcIvS3SMHR4O+jQbRzk5R+HCy+6BEqKS6TZlolbvQ+8fPly6N27N2RlZQHvz0uVt3r1akhKSoLBgwez32OPPSbN9imuu9GSl5cHX331FcyePVtk7LPPPoNrr72WnePx008/ZXG0RsPDw5lnpaXF4Spt3749M27eeecdkYa7CI435T93ZSnfOAjsLW+CorJy4zBEnGhGwCy6CzVjRKsCA6k/s+jALHyijgOpP63XlBnKGcG7sr+gCIpKTpkBLo94VPK6yA2XhNQsmDBrHRzIj4IeWb3g9dff8KgNoxcOVh84ULhgv/3WW2+Fb775Bn799VdYsmQJO8rbHzt2LOzYsYP9/v73v8uzvT7XfU7LHXfcAU888QRUVDgmmZ04cQLS09MZk3gsLCxk8f79+0N1dTWMGTMGnnzySSch7r//fpg8eTLMmjXLKV3ppKqqClYWNjiyVv3A4rjrLE54w/GjGPDrBj4s1MLY7Cwxa/v+o1BZUyeeSyMdkxOhV2YaS6oQvAQ7DuRJsyFKwsvw5EhIigpj+XuEDnq+4FlQComRYZCT4lCHkzyyCjgBE8czY8CvbfjyUgvnpkWJWZtKGqG80WEcihlCJEOg11egi6GsoQU2lzayuNKfnjIdqm6CsrBqyItrtoxMHDMr6UlJpgrhWjpVUw0rmyT3Hi8oHI107e3uNwSqDudLuHNEh2Z3Cclrj997kOrAwq/PCEEHGKR6iG1ugOwax9fYjXH256qDI0cMFw7o0FjDEk5ECuPioxMcmbLYiCr7OwaTd7VJhupwx3NQWrR9QzX0qLd/8W4zerjze0RaUIjr+dzjpH15RnD9ye+/YD3L9ZDJ3+8nthP9rv2M1fg20TAky77BISasPZ3O5ZAe9e5HHCksgfoWxzdjpX6EtP3BPTMhIVaY+ymE3/MK4Xip8oe+YMqEvGHfqE+mfYnkmupyOPzbNkyG8jLHgh+92gt9sYJqGD75ccjb/xPcfPNfICIiXOwjsgqSv4SEBBg9erSYgl/51UK/fv2gS5cuLPvIkSOKnWled9KkSTwK69evd+qzihlCJCMjE7KzB7CksrIy2LBhgzS7VTwYfeBWTPgxYdOmTczD0qNHD9bKFVdcAeiYQOwDERy9ZB1a+/LLL9nQrqFDhwK6h7QE6dwVafnu3btDTk4OLF68WJqsGG9ubobDEXbDxKlAcTk0CH3v/JN2A6pM6NsXlakPRdp0xPHlI7+0GurqlDtj5U0RUNpsL1tbXwcFMpphEl4abckQE2V/YRbZyqG8Qbn96MhIqLaliOwfjnC8dMXE0xFcbYtNXhfOyxtroShC+QGGxVfbHJ2AvMgSqGtRNkZKhcnwJ2y4ihdAbVgD5Ec4OhIsUfKnp0xosDQJc5i2RXayjEwcKivpSUmm4+EFgu4alO89oYLhrr0eCs8Igc9tFrqflPSk9ozg995qW1deDQLyjJDoIUr4agdJ8WL7hSVlYlwewZXOCmNsLLmith6Kq+0GjLRcuM2evzc9V0w+cVJ47uGLQBKaTw89rhc6WfWJp597DYF77nFWfHlGlIWVs2en/N0XrGe5HjL58/3UEilcG5J3dXRtEzRI3/mSPC4LP+JKX3ziPA7r8rUfUVknjDIRPlTyPodSP4K3jceWgjKIsdn7DkUlVcLQMpV+RBBlQj4dfaOVwsfPcKgqLcBkwDusQeg38dBDMFxi4rfBqO5tYeq4t4WJ+b/DqVOnFFcUswn3tPQj+MGDBzmZVsfy8nLRUMH4yZP2DSxbFRQSpMYPekfUpiPgR/b8fPuH6draWigosMukRDNYfWAlXnjab7bBEG9TfvfxMvxYGLEF9qxbJ3i/XudJcOONN7IfT8jPz4fOnR3GfmZmJmzcuJFni0c07gYNGgSdOnWCp556CtBJoUfQ1WhBa/Xzzz+Hr7/+GlC5eNFcffXV0KFDBzh27BizpPGIc1a0hAcffBCmT58O48aNc1kcLfHrU6tVygjp0qkxEt25HNPapa0KPYXkrA5OiU5LI7ZI2scPJfaPJU7lxZN6h/GRK/n6KeaLEQlN1KCrsvWHxFrsyeE4k8WQZokjzRVNHWXCr4T40p0aedyhJ5PL5ADROnpSkqkxXNBdXJL6vWfwa09JJgiha4/fe7nB1hPYPdq4yli2QymaY85Duk4bJ5plEp55QXjuOYTz/hmxstb+7Gz97gvOs1wPmRw05DHfZHL5rj/dVI7mdz72DTLlDCqfq9BEr06+YCQ5tSnrRygTFFJVaCqVd6KvVEBM810mkZQksrfveMkZOK0m1iTkbBL2bcFwaFcx/Lr2n/Ddt18DfvQ2Q3jzzTdbsRmsPnArRnxIwJFPixYtUqUgncrBC+HiWdJw5plnwuHDhyE+Pp7ZAxdffDGb0y4t423c4Z/0loKk3vz58wEtVlwBbOnSpTB+/Hgm/NSpU4HPT8HjRRddJKmlHu3Tpw9zOaH16o/A58A4GRn+aIhoEgKEACFACLhEAI0P+U9eQZ7vbLDIS9M5IUAIBBMB+RwX+fwWnJiPoVv2ldChxwR45ZVXgsmuz22brQ/sjcDoWTl69KhYFfv86E2RhkTBe40GC4YpU6ZAg+DFduX1ktZ1F9fVaFFrDOenrFixgi15jEc81xoeeughZghpLe9tOW7AeFuf6hEChAAhQAjoi4DcSNGXOlEjBAgBfyOg1XDpOWQ2mw6AHVyrBTP0gbViPnz4cOY1+eOPP9iQOnRQoGNCGo4fPw7cI4NzYHAKR2qqq+E70tqu4zgYwi8hNzcX8IcBmV21ahWLu/vr1q0b7N69WyyGY+JQ4EAFNF60uJJd8YP1yXvjCiHKIwQIAUKAECAECIFQQAANF+kmlPKNJ9HjsnXzAIiNjWUd4kBN6vYn9mbtA7vDJFKYf/3CCy/AeeedB7iSGC6WhfNVXn75ZVZ1zpw58OGHH8JLL70EWLZNmzZs5JV8CJm7dtTy/Wa0qDXor3Q9DQU9DBd/yWk1urjCFJuwLR2DbjUhLSoP6c7ciiX9kf7MjYB5uceVSvkkfPNK4RnncsNFXrtH10pdv8jL6dO5fgjgkC/8SQMaKzzghvD480cIyPAwfzCuRNNXD4mUJnlKpGhQnBAgBAgBQoAQIAQIAX0QkM9vOfDLajYJHxduokAIqCFgGU+LmoC+pPvicdHT8+OLDFSXECAECAFCgBAIZQR8eZcbETfpUCsl/uTzSJTKGC0tJjYJdu7caTS2iB+DIUBGixuF+PKwC4bh4s7bZDQPEm5mh3tDuF6O2Y2SKDsoCJDuggK7bo2S/nSDMiiESH9BgV2XRnHzatwLzpPli7Fhd8YKZ05ezgxGTL/hF8Fnr9wIP//8M4wcOZKLQkdCwAkBSw0PQ8ncddqdpNd44ktHH/nxB09y1rW2o7WcnL6/znEHZPnGb/5qi+jqiwDpTl88A02N9BdoxPVtj/TnGZ6+vMc9a8l96cqaOtXNq5VqoxEiN0SUyqml+VJXjabe6eHh4TBv3jyQzo3Quw2iZ34ELGe0GFUl/jIWvKXrbT2j4kt8EQKEACFACBACVkLAV2NFioUZDJc2/e+F/fv3A+5CT4EQUELAkkYLdsj1DHrS48aCrzQ5HV/l9JUPX9un+oQAIUAIEAKEgL8RMNu7zh9Ghj9o6qm37WsXQ1VVFSQnJ+tJlmhZCAHLzmnBB5SR3MFK14zSQ9QVz0rlleh6msbpumrbU5pUnhAgBAgBQoAQIAQ8Q8DfhgXSD8YcFyW55Pu1DBgxDTp1fQYum3EFLHr3HXFXdc8QpNJWRsCSnhauMN4Z5+feHPWg4Um72J7azxM63pQNtKze8Eh1CAFCgBAgBAgBQsB7BJQMCO+pua+ptb1ff22E0Zd/DLv2lsHAQWfCrl2Ojcbdt0IlQgEBSxstqEBuAHijzFDsxPuClzcYUx1CgBAgBAgBQsCfCITiu9wdnmhIaDUm3NFSyuf0PW0jOiYJRl26BFK6T4ex486GX375RYk8pYUoApYdHibXp/yh5WoolLysnFYonHMMXOGkBw4ZMeFQaosRSAnLP1IwFQKkO1OpqxWzpL9WkJgqgfRnKnU5MdsxORHKmyKc0oJ1omRUeDN8TImOpzINHd4Otm4+CWFhYdD3rLsgMioWxo49G77/fiUMHjzYU3JU3oIIhIzRItcd75TL0+ncGQFXOOlh0PRNjIATtkSAemGvFgqmQoB0Zyp1tWKW9NcKElMlkP60qcvVO0wbBf1L9cpMg9LmU/oT1omiHgaIHqz0Gj4HCg8XwNizx8P5kyfCwhcXQkpqih6kiYZJEbD88DCT6sUUbOPLgP9MwTAxSQgQAoQAIRAyCBjRYHEHvjdeDnc0zZw/evpj0HPgQvjxp3y4cuaVZhaFeNcBATJadACRSHi/qWdZQwvUNjQQhCZEgHRnQqVJWCb9ScAwYZT0p640o39Mq6iuhdr6OlUBQslwwRXEpAGHiMlD9p/OhdyrFsN3330Hd955J7S0tMiL0HmIIBCyw8NCRL8BEdOXYWKbSxshP0J4aKUGhFVqREcESHc6ghkEUqS/IICuY5OkP2cwzeRV2XEgDwrKagGyOjgLQWcMAT63RQpHZFQbeO211+CGG26ALVu2wNq1a6XZFA8RBMjTEiKK9peYvhgs/uKJ6BIChAAhQAhYGwHuTeFHq0kbSt4WJd0peVx+PHIxK7pu3Tp48sknxWrffvstZGRkQJcuXWDp0qXME3PTTTdBXFwcWzb5+PHjkJaWxib4v/7662I9ipgPATJazKcz4pgQIAQIAUKAEAhZBMzkVfFFSaFiuMiHiHHM0HCRGi/dOpfxLPj3v/8txpcsWQIFBQWQlZUFV155JYSHh8Orr74K1dXVMHBgNqSnp0NRURErj54aXJ2sgYali/iZKULDw8ykLYPxSl4WgymE2CEECAFCwMIIhIqxIlUhGi5GWc1LypfecTRcth5OViTLDZdft2xg+X3PnABFR7bDww8/DDfeeCMsX76cpf/www/w2GOPQV1dHfTt2xeGDRsGq1atgtWrV8MHH3wAV111FWRnZ8Onn34KZWVl0K5d6/kzigxQomEQIKPFMKogRgiB1gi4ekmT0dgaL0ohBAgB6yHg6jloPWlbSxRKhgtKr2a8XHPpGLjmUvv2CNs3r4X5f78FXnj+eTglGCA8cCOFn/fu3RtuueUWfsqO9913n9M5nZgHAc3Dw2prayEnJwcGDRoE/fv3h3/84x9MyhkzZrBNf3Djn27duqluAIR5aOFiObR+eUCX3vjx4+Giiy6CyspKlvzII49AbGwsFBYW8mIQHx8vxilCCIQCAu5e1Jjvrkwo4EQyEgKEgDURoGecQ6+hMlQMJUavi9LPgQbAkOFjYednb8FXr/9Hmgw//vgjfPTRR05pgTyhvrJ/0dbsaYmOjhZ2Jf2eGQ84FnDMmDEwefJkWLZsmcjh3XffDUlJSeK5PIKuO7k77rnnnoPnBUv54MGDsGjRIpgzZw6rhuWefvppWLBggZwMnRMClkbAU0OElyfPi6UvCxKOELAEAvx5ZQlhAixEqHhctMDKjbhe3TrD1KlT4ZdffoEDBw6wqtOnT2fDvxIThY2rAxyor+xfwDV7WnDiEvd2oNGCP0zjAdfNfv/999kkKJ6m5djU1MQmTeHEKena27NmzWIGUUkJ7ZSuBUezlhmeHAkZ7ZTHsZpVJl/49uWF7ktdb3gm3XmDmnHqkP6MowtvODGD/vCZJP95I6vV6gzumQmdOrb3SizeWfeqskUqyTHIP3QA0FjAgKuG9ezZE7744ougSEt9Zf/CrtnTgmyggTF06FDYv38/3HrrrTBixAiRO1wzu0OHDtCrVy8xTRpBRU6cOJEZOnhR4eQpDHPnzoVrrrmGeWgWL14sVkEDCQ2XZ599Fh599FExXSlSVVUFa3ftV8qCrE7tIT3V7v05VlwG+wvsK0goFR6bnSUmb99/FCprlDd/6picCL0y01hZ3CQK11xXC/hwSoiNYdm/5xXC8dJyxaLxbaJhSFZnMU9NHixgFJnwZXR0/Y+QX9ss8i2NJEaGQU6K4xJbWai8iSTu05IX1wyZMXYbOk+gt7e8SUrKKX5uWpR4vqmkEcoblTeayhDo9U2MYGVxIzbc10AtYAcgKcpuhO8R2vZVJmynj9C2JzJxo8O3a68jEzGq8Cj4W6aKphbIP1kKK5uU9WpVPfFryArXnnyPJCvIhPoJlWsP9Se//7Q89xrS7O8a6ftJj2c5Yo/PvbbDx2EU6J2r3o+oaLK/mxAnT/sRRaVLoTjC/qzH+jxEt4mH7r0dw+/37ljNs1odO2SeAcntOrH00pMFcCLvt1ZleEKfwbk8Cn/s2wJ1Nfah/GLi6UhSajqkd+7Nzmqqy+Hwb9vkRcTzrmecCW1i7Z6QY0f3QVnxMTFPGpHKhMYK9o3kvbiU9h0hOTUVcLQP7uWCHpe//vWvzNvSo0cPKTkxPmnSJDG+fv16qKioEM+lkYyMTGFqwwCWhJP3N2ywLwYgLSOPG7WvLOfTjOeOHqUG7iMiImDHjh1w6tQpmDZtGuzevRsGDLArE5ecw6Xm1AJeFJ06dWLzVCZMmAB9+vSBcePGQdeuXWHNmjWK1ebNm8fmwOCF6Co0NzdDPm7UpBDqIyrgaJW9U1shzJkpUimHVTcdOSVSyC+tFlagUO6MlQsPm9Jme1nc1ZZtEiXWdI60FJRBjM3OW1FJFVRUKvMZXdsEDdL2XfBpJJmKbe2hvEFFpshIqLaliIAcjnDMURIThUh9WARsi+wE+2124668sRaKIpSNO6y32mY3GDGeF1kCdS3KxkipQO+Ezf5QrA1rsG9iiZUUQqMtGWKi7MZQka3cZ5mwierIRM0yfZ8kPFhP61+Xay+6PdQntofY0waFkkwtkTaGRExYsyY9KclUH1YFhyOUh4RaUU9WkwnvvdW2buw6wD+r3E9W0xPqRkkm1J/8/lN77vH7HWnB6feL0/sprgtUtFRBWGM9KyL9i9bwLOf0q+LaQsLpZxm9c9X7EbY2bcQ+h3f9iEPQFObchWuoiobCSIf+qkrVP6iWhtVBQ9khpubKygqoKpWbAo4roPDAT+JJ5cl8qBf6PUohrKEMIuqLWVZdfa1As0CpGEs7ebgZovk7v6QIaiuETaYVQkt0FJRH2Ptbm4R8pb7e5MmT4GRJKdhsNjbfGkfo7Nu3j82TxqkHSoGvOIZ5eXl5gkwO3KTlcW51fr4dR5yvgvOw3QWj9pXd8W2GfOcrXiPHbdu2hdzcXLbMHBotjY2N8PHHH8PWrVtVKaDBggE3+EGDZ9OmTcxoUa0gZGA7M2fOhIULF7oqBgkJCXDFGLvx5LIgtBWyM10X4bldsKzGoHVXWw9o5mguG3yZXM6lqHcYH7kKu96jV2O3YPhMjTwuWC+n8carUqGsqI36Q2IU7DaJ49wpVi3QLHGkuKLZgmVPF0XbyW4/OepKYxKZhuWOYTnKGEhoupGpUqpvaVzarlJc67WnVFchDeVQ0pOjqEOmPdVNUBfbHi6LqXBkS2MG0pNWmcCNnsBCMvF7L9dCMomXXwjIxPXX+v5r/dzjXlwRH6WI5Lmj+DxTeO4pkXFOC/77yZkf9bNAvnNx1MVvJ6vAqU2tz3KJntwthZzW8yx1gSU5aZK4u2haT3clHPmd+zrirmL/v70zga+ivP7+CQlJQCAESNiCgARERPbFWuTFiqhArSh1X5EqiLV/xIXW1xba16Va3s/f4gLqn0VR0dZa1FeoUgW0ikgFrUWQUEGyQBII2TdC3jnPZS5zJzOTmXtn5s7c+T18wp2ZZzvP98xz7znzLDOxz3YpurdRknBcBLPwVaJZC/+PcFpuuOEGxVXnDlesWNFq4V6zlVsV2AcJ2piVkV/MwyMsHGpra2njxo1itITP5eOcHG2HgKdvyUNvfPzee++FR2g4v1G45557aPny5cIxMkqHuPgRMPWDqCMeT8OqqNYeqdHJ4onL3GZlu5XHnhAwSiGstMOvuosSTcJlg/78rVKz+rPSp5mIpsNyEhWXZbU8f1N2RnqeJl5ZVR1z4eq1HTEXGIcCbGuDYo11HJoRrhK2chiFIwemnZaioiK64IILpLeLDqOxY8cST/GaPn26EGrt2rUtpobxENrUqVNF/OHDh8VuY7xdMm+bPG3aNFLOJzRqGe8ixiMz/LIgBO8SCNIPmV5b9a63prVo87VWbrTxXpMn2nYgHwgEnYDVvmzksASdJdrvbQIXjB9FTzzxBK1fvz6ugsJWdhZ/krRjl/YqZmfrtbV0dqT+vmKJrWWisOgJWPnh48X5PCf7tq7SlAYfhNaMACttl5vbWplyOrc/W2uL33TnNj+v1wf9eV1DxvKZ0Z+V75bW+ruVsowlRywvJue1GeamtRvzam2KmHHu+MbaNspyshlPrfkzPf3Sn2j1iy+Jh+NOto7fN7h9O09r827gZRwPPLJK2vlXe92pWvI/v7yMzuqfGX71iDreC+dRrWnxguCQwbsElD9urf0QercVLSVTtqtlbOgKp0mUNidSW/T0hesgEHQCZr+vOJ2Z78Cg80T740fgrhtmUv/ePemWW26lkSNHUGlpqdixduXKlWKWUPwkQ812EYDTYhdJlKNJoLUfuUbpiVOz9MSp6pxRIr/ZH1DNyhy82Fo7lFVbMfatlKusA8cgAAIgECsBr37fxtou5Pc+AbtHWeQWT7vgh9Rh+BW0+f23KFd6TUdVZbnY0GnHjh3U9uQOoXJafPqPAJwW/+ksoSVWGvFe+EFVypPQ4HUaZ8UB0ykCl0EABOJEAKMjcQKPauNKoHuPPnTVjfOEDLwC4vI3VtEbb7xB11xzTVzlQuWxE4DTEjtDlBADAX6pJr+jRitoOQxuODJa9WrJZ3TNjLFvRz1GMjgdxy8P5fc3IPiTAPTnT73JUpvVn57jYvW71O/fVzI3L3wa/e5Zkc+v61mcGmVhdmom9fV1VFZWRu+++y6cFis3l0fTwuLwqGKCItbI3D4RL9Vsrd1aP5xGP75yeq00clxrdUYbz+Vr1cvlOV13tDJr5dNrx7guJ18eqnh/g1Z+XPMmAejPm3oxK5UV/fH3kJ++c8wy8Gs6q797fm1nvOV+8L+up48+XE9Dhgyh1atXx1sc1G8DATgtNkBEEfElYObH2EwaJ1qhNPjjJYMT7TJbptxmPefNbDlIBwIgEBuBWPqg3I9jkwC5QYDIyVEWNV92WDhs2rRJLMhXx+M8OgIbNmygX/ziF9QkrRmaPXs2LVy4MKIgnpLH8Ty61b59e1q1ahWNGhVatxyRMIoT0+9piaJsZAEBEJAI8A9+EH/0g9hm3PAgkEgEgvrd5QcdqqdB+UFmt2V8YPGT1L17D8rKynK76oStjx2VefPmiffh7Nq1i1599VXiT2Xgd+Xs3btX/D333HM0d+5cZXRMx3BaYsKHzLES4P3q//N9fqzFIH8cCPB7IvYVFbeoGYZOCySevKCnP08KC6FaEHBKf3L/xUOHFshtu4DfPdtQGhY0Zfo1dPjwIdqzZ49hOkSaJ7Bt2zbKzc2lM844g1JTU8U6oXXr1kUUwOc33XSTGN0699xz6dixY8Qv3bQjYHqYHRRRBgiAQCBHk6B2EPArATglftUc5NYjoB59apvSliZfMoMmXzSF/ueF52nKlCl6WRP2+pcHMym9vbmXS+Yfa0dHPv6YXnjhhTCP22+/nfhPDgUFBdSnTx/5lHJycuizzz4Ln/OBVhq+1rNnz4h00ZzAaYmGGvKAAAiECcD4CaPAAQiAAAiAgIcIPPjIcvrvpSvpxptm0cyZV9BTS5/E+hYD/UyYMIHWrFmjm4LXq6hDUlJSxCUzaSIyWDgJ7PQw9siVfxaYISkIgMBJAs0pqWABAiAAAiDgMgE3F7S73DRbq0tOTqFJM+6nu/+wg954cx09++yzpGVU21ppAhfGIysHDx4MtzA/P5969eoVPucDM2kiMlg4CazTYoERkoJA4AlojaY0Zp8aIg48IAAAAR8R0OrPPhIfop4kwI6L35wX9RQut5R5WkYWzbxzBS367SM0Z27oxZNu1Z1I9YwdO1YssP/uu++ooaGB1q5dS5dddllEE/n8xRdfFM7h1q1bKSMjw5apYVwJpoedRM0dyW+dP+IuwQkIOEyghaEjbaKAAAIg4C8CLfqxv8SHtBoElLZLvJwCDbFcvWSm3bnDLqTbfv03evZXE+j+exfQgNwBrsqYCJWlSC+Ufuqpp+jiiy8WWx7PmjWLzj77bFq2bJlo3pw5c2jq1Kliu2NesM9bHq9cudK2pgfSaTFzc9tGGAWBAAiIRfqxvCcCCEEABEAABFonoHRgOLUX7R2WSS1n6y2zJ0X3PmfT+Ivn0nkTzqd33lpHY8eNtafgAJXCTgn/KQM7K3LgNS5PP/20fGrrZ+CcFqMOLMfFqzPZqlmfFJbbK4sakit9Ii3EVBKA7pQ0/Hc8uFMy1aR0kgSv8Z/wAZRYPULS90g5VR+upKohOWTmgYA6fwAReqbJbn53Ku0Z2cbxAghZFqV80coll2U2/6U3PEJJbZLpzrvuoo8/2kJpaWlmsyJdnAlgTYuGAqx2AI0icMkkgZ5dM6hjhw4mUyOZlwhAd17ShnVZctLbUKf26dYzIofrBLQcDmX/43itNCyoUZzrDUGFgoBSd24iscNBsFtetrfkP7vLlssb3bdMPgx/TvzxfDpUUkHp6enEL0BE8AeBwIy0WHVEOL0XO7g/bitICQIgAAIgYAcBPWdEq2wrabXy41riE2C7xqo95BYVWS6ztpecPhr52nfsQgv++A3t/WojLfzltTRgwAC68MILoykKeVwkkDBOSyw3rx5vZZlmO5FeWbiuTaBImuJQWVUlRXbWToCrniUA3XlWNaYEy687QRXH67Adiyla8Ulk5ISg/8VHJ3bUCt0ZU5RtL7XdJV83zq0dy6Mt/zyQ2SJy4LDJdOXc5+nyK66UdrjqRcueeYp+9KMftUiHC94ggOlhJvXAnSWWDmOymsAlyyssoZKjxwLX7kRosFXdGRlgicDDb23YXdFEJeUVfhMb8p4kYLX/AZx3CEB35nQh213yp7lc+qm0polx6rPHX07/+/mD1JzSBaMt+vg8EWO701JXV0fjxo2j4cOHi23QfvOb34iGLlq0iHr37k0jRowQf++++24YwH333UdjxoyhzZs3i2v79+8XbyxdunRpOM1d0oKpVatWhc/jdQDHJV7kUS8IgAAIBIcAnPzg6BotdY+AnuOS1q4jpaTEbhInug3snqa0a4pdQ6pyeReGDz74gL788kvauXMnbdiwgfjlMhzmz58vrvF1ebu03bt3i7gtW7ZEbJGWnZ1NTz75pHh5jUjgof/guHhIGRDFVwRgiPlKXRAWBEAABBKOgJ7jMvJ/3UTtT8uQ3t9yJv3qVw/SiRMnLLc9CDawZSg2ZrDdaeH9mTuc3A2qsbGR+I+v6YWmpiZq06aNSNPc3BxOlpWVJYbpVq9eHb7mpQMrjos8tKn36aV2QRYQcJoAHBenCaN8vxNAH/G7Br0rvxXbxbutcEaycZNn0yU/20JZg26hZ5Y9H37gbqW2oNjAVpjYmdaRhfjsiIwePZry8vJo3rx5NH78eFq/fr14i+aLL74opoItWbKEMjMzxRSympoamjBhAj3xxBMRbVu4cCFdeumlxG/cNArV1dW0e+cmzSTdcwZRZrdeIq6stJAO53+rmY4vDh4xKRz33Z7tVF/LC8RbhoyuPalnnzPFGpfexz6mnfvyWyY6eaXDWVOpXWropOjgHio/UtQiLY819U0to5G5fcJxHxm8bZz3eOctEznwgj6eH6sXzj8nNxy1I+8gVdXWh8+VBz0yO9HAnGxxqbKmzrBNIwbkUMeTW6XuzS+mQ2Xa8+I7tEtrtU3fFx+lqvom0Y5EaZPMNZH0pNWmKmkqaHlFJendq8b3Xg9qW3xQFNtb2nr3LOmdIRzKG5vp87Lj4ljrv7GZKZTRNvQQ5BtpTUaBtJhcK3RKSaJxXU59vW0sbtRKJq7x+0p4+18OvDid13rohcnZbcNR244elxayn3rQEo6QDvzQpv01TVSeJL2jpespyf3eJrklftZTY7b0O6Dx/a/uT/J3p7r/xeu7XGaf6N97dvzmsu4amk89M3bzN/dIcsiucdI24nuhtqaCDnz7hXxbtPjsO2gUtWvP74mS7Bgd24jj0tp1oP5njuFDEfRsPY600qbRkr0nL8yvLvoHNdVXUlFhDWVLu8BnDz+TTmscRXv27KGOHTvROecMFXWXl5fTp59+Ko6N/nPbBjaSJdHiTv2q29iy5ORkMQ3s2LFjNGPGDPr6669p7ty59NBDD4kRFf5csGABrVixQtSqXLuiFKN///5ifcwrr7yivNzimIfwqsu0HYeypHpqLN8v8lRVVUrp9A384n2fhMuuKi2QpqZpG/hJjeWU3HBEpD3Y0EBHyqUdeBShKekU1uwD2yktNfQuhIqjJVQtyaAV9qcmUeP3pxakF6jKVObhlzEerA4ZS7zzVolB2m3KMstqqL5e23iraEqmshOh+uukdhcalNlcWE7pqaE2lxytlnb/imy/LGtaXVOrbWKHpUkaYdslvSQtUdoktz+R9KTVppLKRqE7vXu11XsvLUsUe+S00+hwasgZqEtqpILklnvqy/UfT82k9LahtCWpFVTRqHPvpaRQTWoXORsdSC4OH6sP+AWLeXIflXbTKknWdsI536bUkGPPx/kpR6m+WdvBKpPKO5wa+kH2apvYYWmSRsE3pfbl5ojg9zbJ7fCznpp1vnvV/Un+7lT3v3h9l8vsE/17z47fXNZdG+nBiqxTt35zQ7ZJyFZy2jaqb6iT7K1C+bZo8Vl64IQp26ixOo2KUxrC+fVsPU5gtU1dGkMPy5rr8+l4YwP16Rx6YHXiRBN1GTeaqiT7aseOL6igIMSM16sUFuq3SRbSbRtYrjcIn0nSlCztR4U2tX7x4sV0mmSU3HvvveESeaH99OnThTMTvqg4UMbzmpeZM2fSxIkThQNzyy23KFKGDocNG0ZPrwkt4m8R6aEL6u37PCRa3EThp4T8o3vNhNCTjLgJgootE3BCd2be7G1ZUGTQJMCjTweSM+i2rtJoC4InCFiZFuZE//MEhAAIES/dYWpY5M0lj7TIV//5eak4ZKdl00tT6IEFt9Hdd98tR0d88uZR27dvj7imdeKGDaxVL1+bNGkSXTJ3HaW3D83M0UsnX//onf+mi0am05w5c+RLnvs8NT5pk2glJSXEIywcamtraePGjTR48GAqKjo1LerNN9+koUPNGamcd8iQIfTOO+/YJKG7xbCjIv+5W7M/auMpD2ecnuMPYSFlBAEndMdGmxXDLUIgnFgiwFOoBvQ8NXJkKTMS207A6n3vRP+zvVEoUJNAPHQHh6WlKtQL8keP7SYStWmTTGnpHYQN2zKX8RXYwMZ8Yo09NY8p1pJO5mfn5Oabbyae08fTtq666ioxqnLjjTeKKWO8SKlfv360fPly0zU++OCDNHLkSNPpvZAQoype0AJk8CsB2YDDyItfNQi5zRKQ73Wz6ZEOBKwSgMNilZj0zt12XaJyWmADW2dtJYftTgtP1dqxY0cLGV566aUW1/QusFPD62DkwO98iWbrOTm/W59wVNwijXqCQoANOjguQdF2sNoJZyVY+o5Ha4PmrLANZleby4t30Suv7qPzzjuPpkyZYlp9QbaBTUOKIaHtTksMsvg6KxyW6NTHO5oVSBsE0OmdoysAueJGALqLG3pbKuadwnjhPYX2C7ClTBRiTMBORwX9z5i1l2Pd0J1dxruXOSplk20wuxyX8699h9rUfihmCvH6lUOHDon3BraVNoLp3Lkz8WJ7BPcJwGmxgbncWWwoKnBF8BbMejuaBQ6Gzxrslu4w2uLMjcHbNevtfuZMjfEvVc9psGM0T69sp1rtVv9zSv4gl+u07oLqsERzT/G6FuWCfF7Xwgvy252WRaMn3UmduvSkJOnfHaOyKC0tXewA2y29nH73u99FUx3yxEgATkuMAOGwxAgQ2UEABEDAQQJmnAllGjMOjDK9g6KjaBCwTAAOi/TOvRinickL8hn+0PEzhA7OkJwbOZw/LBNOiwzD5U84LS4DR3UgAALWCJgxIq2ViNRBIBCtYxFtviAwRRu9TQAOi7f1A+liJ2D7lsexi4QSQAAEQAAEQCA6Aux0wPGIjh1y+ZdA0ByW1jRlZRaMeutjZdn1tZW09b3nqODgf8RlfrVhg/RScYT4EMBISwzcrXSKGKpBVhAILAGMsgRW9ZYbDkfFMjJkSBACQXRY3LK//rJsDm378BVa/1IG1dc1SGtwaykjI4Nyc3MT5O7xVzPgtPhLX5AWBEAABEBAIgAnBbcBCJBtW/z6iaVbDgszufruF+meexfQ0dJiuvPmSykzM5PKyk6tb/ETt0SQFU5LImjRx23okdmJKpqwdaAfVei07jDK4uxd0Tu9DZWlpkuVSFuOezTAMdFXjNP9T79mxMRKALqzTtBpR0W9i5gsYRtpa+NeOf3F39r/909668+raM6tP6X58+fLSfDpIgE4LS7CRlUtCQzMyaayE8daRuCK5wk4qTs4LM6r/6xOyXQ4VXpJS4P0rhYPBDgo1pTgZP+zJglSWyVgl+4SeVqYHU4Kl2GFkZ7jIuuXnZc5/7WYePcwhPgQgNMSH+6oFQRAAARAQCIAZwW3AQhYJ2DFGLdeevxy2OGs2Cm90SJ9O+tBWeYIwGkxx6lFKq91rBYC+uRCZU0d1TXU+0RaiKkkAN0pafjvuLyxmeqSGuMmOJyV2NCj/8XGL565obt40jeuW3ZS+IWT8rFxDsS6SQBbHrtJG3W1ILBzXz4VHippcR0XvE8AuvO+jowk/LzsOBWUxmdBKRwWI82Yi0P/M8fJi6li0R2PsGCUxXmtwmFxnnE0NWCkJQpqGGWJAhqygAAIBJ4AnJXA3wIJDaA1ZyIa26G1MhMFaDRsEqXtaId5AnBazLMSKdGxLAJDchCIggAbt1iMHwU4D2eBw+Jh5UC0mAiYdSzkdEeSq6T69mvWKafRjMRFEAg4ATgtFm4AOCwWYCEpCMRIQDZy4bzECNID2WVdekAUiAACthKI1sloSkpJ2GleVgHDtrJKLLjp4bSY1D06lUlQSAYCNhMwMnjh0NgM2+bijHRnc1UoDgRcJRCts+KqkD6ozG+2FS/QP3+YD8AmqIhwWkwo1m+dykSTkAQEEoKA2iiGE+MNtar14g2pIAUI2EMADos9HP1WCjssCPElAKelFf5wWFoBhGgQ8BAB2VhOROdFbhvjjkf7lPV7SOUQBQRcJQCHxT7csK/sYxmUkkw7LXV1dTRx4kSqr6+n48eP08yZM2nx4sV033330dtvv02pqak0YMAAWrlyJXXu3LkFv379+lHHjh0pOTmZUlJSaPv27SJNYWEh3XDDDSLu5Zdfpg4dOtCiRYvo8ccfp/3791N2drZIx9erqnjxmjsBnckdziMG5FBzYbk7laEWWwl4WXdsYMfDsLcDsBnnQJ0mmraOzUyh46nSk8PmGkOx1XUZJkakawS83P9cg+BiRXY6K30HjaLSAydclD6YVdmpM7MEg2YrK7kcPXqUrr76amG7s83/+uuvU2Zmy9EpPX9AWZbesen3tKSlpdEHH3xAX375Je3cuZM2bNhAW7dupYsuuoi+/vpr+uqrr2jQoEH06KOP6tVFH374ocgrOyyc8I9//CMtXbqUZs+eTWvWrAnn7datGy1ZsiR8joPEJNCxfTqlp6YlZuMSvFVe1128jW2uX+9PeWuo0yjjzB5H09aMtkmU3ratYRXRlGtYICJtI+D1/mdbQxOwoHbtO1FaanoCtiyxm2Tm3S1BtpUfe+wxuvDCC2nv3r3ik8/1gpY/oJdWed2005KUlCRGQThzY2Oj+ONrU6ZMESMnfP3cc8+l/Px8PjQdmpqaqE2bNuKvubk5nG/WrFn02muvEXtubgeMsrhNHPWBgDME3Da6lQ6IUYvMpjMqQx1nd1vtLk8tL85BwC8E4vHE3i9sgiZna45LkGxlte7XrVtHN998s7jMn3/961/VSWI+Nz09jGtiB2P06NGUl5dH8+bNo/Hjx0cIsGLFCjE0FHHx5Ins4PDnHXfcQbfffruIueuuu+jGG2+kjIwMeuWVV8JZeToYOy5PPvmkmIYWjsBBQhHYm19MJUeriU5vOaUwoRqagI3xi+5k4zuaKVRm1CaXbyat02lYFrPt/KaiiUpSK4g0Hvh6qU1OM/Nr+X7pf37l66TcRQf3UMXREsoe4GQtwS47no6ml2zlHV8codS0RlM3Q/731fRx9XZ64YUXwunZVpft9fBFnYPDhw9Tz549RSx/FhcXa6bU8wc0E6suWnJaeD0KTw07duwYzZgxQ0wLGzp0qCjy4YcfFiMu119/vaqK0Ok//vEP6tWrl2gETykbPHiwWCPTt29f2rJli2aeu+++m0aMGEELFizQjJcvVldX0+6dm+TTiM/uOYMos1svca2stJAO538bEa88GTxiUvh0R95BqqqtD58rD3pkdqKBOaG1NpU1dbRzn/7oEs875mF8Dvwjc6hMMhI0Qod2aTQyt0845qN/5YWP1Qe5vbKoZ9cMcbnoSDnlFZaok4TPzz8nN3zsxTZ9X3yUquqbqKhPl4Rpkww8kfSk1aa8omLKL6mkj5Kb5OiIT+/dez2EfG2LD0bIObhTMuWkhwad8+tO0G7JmNcLk7PbiilfHB/gN4VIAAAT0ElEQVTuTxp9Nb7fET3onLrCVtu0v6aJypOkkey+p1q77ehxqjjeTI2eaxNRovcnq9/luw4Uie9Odf+L772XmHpK6dsrKjviuz3bqb625Vrc0sPfU1JTPZ15suvV1lTQgW+/ONURVUe8BoanlHFgh6f8SJEqReg0rV0H6n/mmHCcnl3ECaK1jfTaxGVmdO1JPfuEWmWmTXTSAnXCNrJi71ltU3VhO26ubvCqrawrsCJiwoQJEUs1FFHicPLkyXTo0CH1ZWIfwGzQ8wfM5LfktMgF8kL7SZMmiXUt7LSsXr2a3nnnHfr73/9O7EFpBXZYOPDCenZ4tm3bJpwWrbTyNa7nuuuuo2eeeUa+pPl54sQJqi7TdhzKkuqpsXy/yFdVVSml0zfwi/d9ItJVSAtTC8pqpE0HtL3TiqZkKjtxTKSta6inwvI6caz1Hy8yT08NxfOIQmWVdtq0uiZq/D5UJpdTYFBmQ3IlHawOTaWrlDYnKDFIu01ZpgfbxA5LkzQtcNfhxGmTfB8kkp602lRS2Sh0p3evevbeS8sSzUk63iA+a1I6Ud7J+eUVx+uoJDnywUJzSqrcfPogI4foZJ/y9HdElzMoLzkktlabOKY8qYaapO/rTamnvJb8lKNUlyw5cBrfKfje89Z3ufzdqe5/0JP9eupU8C/LdgT3sarSAmqQbAR1aKyrpBSpn8k2R31DnWSbFKqThc950b68BoZHaKolW0YrNFanUXFK6HuN4/XsIo6LxjbifHpt4rikxnJKbjjCh2SmTVuk7910yd6y2zaqS2ovbdz0L1P2HstqtU3/qTr55cqZDYLXbGUDUU1Hbdy4UTdt9+7dqaioSIy28Ke8kZY6QzT+gFxGkrSOJGT9yld0PktKSqittGiTlVBbWyvWsjzwwANidOWee+6hzZs3U1ZWyBhQF8EjIexY8O5hfMwjLb/+9a/pkksuUScV57x7GE8Pu/fee6m0tJTGjh0rQPCuDFph2LBh9PSazVpRUV3DmpaosEWViUeU+Ef3mgmhEbuoCkGmuBBINN3J06oSZWqU3B69m2NjcSMdSM6g27rWRCRJlPZHNCoBTxKt/3lZRXZPNeIREHYoRl9wg5eb7YpsdttbdutKC8L5wzJpzJgx4V1wlWm8ZCvz4EL38Suk6WGhUTqlnFrH336+nG79SVeaM2eOVnSr13g34a5du9LChQuJF+HzmnTeCVgZrPoDyrx8HJoTob6qcc5e0wUXXEDsILATwY7H9OnTidekVFZWinOeyiU3lrcynjp1qiiJ57nxkNPw4cNp3LhxNG3aNF2HRV017yLGIzO81bJbwY2b3q22oB4QAAFzBNhYTySDXW6P2TbJ6c3RQioQCA4Buw3r4JBrvaV22lt2ltW65NopgmQrqwmws/L+++/TwIEDxSefc7DLH+CyTI+0cGKvBrtHWrid+JJyR9t4WugOZydqge6coOp8mfIIDI+07JemymGU03nmTtSA/ucEVeMy7TKKMdISyTlWe8suvURKpX9mNNKin8v9GLdHWtxoYVRrWtwQLN51cCeItSPFuw2oHwRAAATUBOSRF7HYXmPtijo9zkEABEDAiwTcdla8yCBoMsFpMdC4Vx0XdUf1s3PFu6bxJgQI/iMA3flPZ0qJoT8lDf8dQ3/+05ksMe/yxYvmEUIErNpaahsIHINDAE5LK7qWO4cXHANZFrXI8nUvyKiWrbVz3uZZuWtaa+kR7x0C0J13dBGNJNBfNNS8kwf6844urErC2xIrd/mymj+o6WVbJ6jtR7stLMQPOizuLPHqMGbrjpd8Qb830H4QAAEQAAEQAAHnCMC+cY6tn0o2vXuYnxrlpKxmHQi7ZLDaUa2mt0tOlAMCIAACIAACIAACdhOAXWM3Uf+Wh+lhUepO3YnsnpqlLj9KMT2fTd4BZ9zpnT0vKwSMJADdRfLw2xn05zeNRcoL/UXycPrMzt9kefew7AHnOS22r8u3k7mvQUD4MAE4LWEUsR0YdS52aIziY6s5MrfdzlNk6TgDARAAARAAARAAARAAAfcJwGlxgbkbDgucFRcUiSpAAARAAARAAAQcJ+CG3eR4I1CB7QQS2mmRDXk/3fyyzGY0ze2ykt5MmUgDAiAAAiAAAiAAAiAAAl4jkHBOi5YRL1/zuvMiy2n2JrGa3my5SAcCIAACIAACIAAC8SLA9o3XbbZ4sQlyvQm1e1hrRnxr8fG8EbwsWzy5oG4QAAEQAAEQAAEQAAEQSBinxazRbzadm7eGF2Vys/2oCwRAAARAAARAAASUBNg2gn2kJILjhJse5ieVojMS5fbKoobkSj+pDbKeJADd+ftWgP6gP38T8K/03XMGUVlSvX8b4IDkRmt02VbCVDEHoPuwyIQZabHC3gvOghdksMLMqbQ9u2ZQxw4dnCoe5TpIALpzEK4LRUN/LkB2sAroz0G4Dhed2a0XdejQ0eFaEqt42EyJpc9oWxPYkZZoPHernUb9ZMBq/miVinwgAAIgAAIgAAIg4CcCRqMt3A4jG0ptb/mp3ZDVPIHAOi2MyIrjYtRZ9HBHk0evrES9XnSknCqrqqTmdU7UJiZsu6A7f6sW+oP+/E3APentNojLSgupqqqSst1rQsLXJNtbdusq4cH5rIGBnB6m1BHf6PLNrrwuH7cWL6fDZ3QE8gpLqOTosegyI1dcCUB3ccUfc+XQX8wI41oA9BdX/DFVfjj/Wyo7WhJTGcisTcDIntPOgat+IhDokRalonCjK2ngGARAAARAAARAAE/u3b0HmHes9pgyv1J/yuvqVinTqeNw7h0CcFq8owtIAgIgAAIgAAIg4BECMGQ9oogYxDByVJTFyumgcyUV7x0HfnqY91QCiUAABEAABEAABOJJAMZr/OjHk73svMSv9ajZiACcFiM6iAMBEAABEAABEAgMATaY42k0Bwa0QUPj7TjEu34DNIGPgtMS+FsAAEAABEAABEAg2ATgrARb/+rWw3FRE/HGOda0eEMPkAIEQAAEQAAEQMBlAhhVcRk4qgOBGAgkjNOSNeaiGDAga7wIXCHpbe/evZQ1cGC8REC9URKA7qIE55Fs0J9HFBGlGNBflOBU2bJU526cnj9shvjdGzgw043qfFiHN+y5eNwbPlSWqyInhNOSmppKY8aMcRUcKgMBEAABEAABEAABEAgegW7dugWv0R5ocUI4Ldu3b/cASogAAiAAAiAAAiAAAiAAAiDgBAEsxHeCKsoEARAAARAAARAAARAAARCwjQCcFttQoiAQAAEQAAEQAAEQAAEQAAEnCMBpcYIqygQBEAABEAABEAABEAABELCNAJwW21AGt6BZs2ZRdnY2DR06NAzh6quvphEjRoi/fv36iU858tFHH6Xc3FwaPHgw/e1vf5Mv04cffig2VHjggQfEtXXr1tGMGTPC8ZxvoGKXsbfffpt+8pOfhONxYI1Ac3Mzybo755xzIjIvXbqUzjzzTDr77LNJ1gcngO4iMMX9RNafsu8tWrSIevfuHe5/7777blhOWX+sW2Xf27Rpk+h7999/v0jLfe/yyy9vkU++wH3vsssuk0/xGSUBLf3JRf3hD3+gpKQkKi0tlS+F+x/0F0YStwMt3T300EM0bNgw0femTJlChYWFYfnQ98IocAACUROA0xI1OmSUCdxyyy20YcMG+VR8vvbaa7Rz507xd+WVV9IVV1whru/atYvWrl1L//73v4mNqXnz5lFTUxOxAb18+XLasmULHT9+nHbv3k0/+MEP6NNPPxX5OJ6PO3XqRMXFxSL9J598Quedd15EvTgxT4ANItbd+vXrIzKx8/jWW2/RV199JfS0YMECEc86Y71CdxG44nqi1fdYoPnz54f739SpU4WMyr7H/fXOO+8UfY8jn332Wfroo4/EOfc97ldy3+N4Zd/jc+57P/zhD/kQIQYCevo7ePAgvf/++3T66aeHS4f+wig8caClu/vuu098b/Jv3/Tp0+m3v/2tkBW684TKIEQCEIDTkgBKjHcTJk6cSF26dNEUg52N119/na699loRz09wr7nmGkpLS6P+/fuLEZdt27aJuBMnTogni23atBFOSVZWFmVkZFBeXp6ILygoEM4PG0wc2JCC0yJQRP2fWnesLzZgeXSFdcSBR9E4sCPDI2jQncDhif/U+jMSKta+xw8f5L6HBwZGpM3H6emPnc7HH39cfB/KpUF/MglvfGrpjh+qyaG6ujqsP+hOpoJPEIiNAJyW2PghdysE+Olt9+7dw9O62PHo06dPOBdPY+Fr/NT/tttuE09v2XA+66yzRBoebWEDac+ePaKMc889VzgrPBrDIwFjx44Nl4UDewh8++234qn7+PHjadKkSfT5558LJxK6s4evG6U89dRTYpoKT2EpKysTVar1l5OTI/oeR86ePVs8AOAHB3Lf4wcC6r7H5+h7zmqQHw7w9+Lw4cMjKoL+InB49uTBBx8Uv3Evv/xyeKQFuvOsuiCYjQT+9Kc/iSnl/ODZ6FUkPNLPU1x5mcBjjz1mSQI4LZZwIbFVAq+++mp4lIXzskOiDuywcLj44ovpiy++IJ7LzYGv8xQUNpR4VIUdGHZSPvvsM9qxY4e46dPT00Va/GcfATZKjx07Rlu3bqXf//73YnSF9Qbd2cfYyZLmzp1L+/btE9PDevbsSfL0PjP6W7JkSVg0ue9x/+O+N27cOPS9MB1nDmpqaujhhx8OG7vKWqA/JQ3vHrP+eHrf9ddfT/zwgAN05119QTL7CPDayr/85S/Eo5B6gZcD8LIAnpbO0ybZRuRPswFOi1lSSGeZABu/fAPzlCI58NNd/kKXAz+B6tWrl3za4lMeaZENp44dO1JtbS1t3rwZc+pb0LLnAuuIN0Bgp5ENVX5qcuTIEYLu7OHrdCk8spmcnCz09rOf/Yzk6Zdq/eXn5xv2PXmkRdn36urqiBftYz2LM1pkZ/O7774Toyz9pA1MWEejRo2iQ4cOteh/0J8zOrCr1Ouuu47eeOMNURz6nl1UUY6XCfAoPY+gGAX+PeIRljPOOINSU1PFcgGePmk2wGkxSwrpLBPYuHGj2CGMv7DlwDsO8UL8+vp68eO8d+9eYRjL8erPIUOGiB1YPv74Yxo5cqSI5l3JeNE+OzQI9hJgR4V3jeLF+Bx4qlhDQwN169aNfvzjH4uF+NCdvcztLq2oqChc5Jtvvhne1S/avsdTPJV9b9myZVhLFiZs7wHv4scbjezfv1/88Xcnjz736NFD7NYWzXcn9GevjoxK498zOfA0P94hkwP6nkwFn24S4Ie8/97yHu3+5J+m/gr27BTTusaMGSN2k+TP5557zlaRjaZKmqkoxUwipAEBIwK8yJ6fvvLWnPwju3jxYrE+hX9g5QX4cn7eQveqq64idkbatm0rhs/5qbBeYCOa11ZUVFSI9JyO17U8//zzMJz0oJm8zlMW+GmgrDtea8Tb5d56661CfzzUy4vuV61aJUZd+PynP/0pdGeSrxvJtPoe65N3L+K+0096Ws8OPgdl30tJSaGnn35ajMjoySn3vfLy8nDf4wcF/COGDTD0qFm7rqU/XtunFaA/LSrxu6alO94Rk9df8uh03759iR18DtBd/PQU5Jp//vOf0423Pkynn/WrVjEcP15N5Ue3SzbZN2Q07X7y5Mli5FddIE+LNPMKCqOpkuoytc6TpAJaLjLQSolrIAACIAACIAACIAACIAACnifA5n2XrLMpZ9B8Sm9/avt0LcEPH3hRcrbb0f49L2hFW7rGG/jw2mQeqVEHXp/MD0fl94Tx+4s4/PKXv1Qn1TzH9DBNLLgIAiAAAiAAAiAAAiAAAv4kwKPlq/7nUSo+sMawATzKUlbyCX2zc6lhOjsieTMlnkbJa/d46jnPyLHyomI4LXZoAWWAAAiAAAiAAAiAAAiAgIcIsEPQUF9CdTXf60p1pOAN6tr9QmrXrp1uGjMRvIaSlwjwaMq0adPEjrCcr7CwkOSXHPPUZN5Vj3eL5YX7vFyAp0+aDZgeZpYU0oEACIAACIAACIAACICAjwjw7ly33/l/Nde28ChL3s4FVFy0K2anxQ0kGGlxgzLqAAEQAAEQAAEQAAEQAAGXCRiNttg1yuJWk+C0uEUa9YAACIAACIAACIAACICAiwT01ra4uZbFrubCabGLJMoBARAAARAAARAAARAAAY8R0Bpt8dsoCyOF0+KxGwvigAAIgAAIgAAIgAAIgIBdBNSjLX4cZWEWWIhv1x2BckAABEAABEAABEAABEDAgwSU720pL9lk23tZ3GwqnBY3aaMuEAABEAABEAABEAABEIgDAd5JbPYdj1BdXbFvdgxTYoLToqSBYxAAARAAARAAARAAARBIQAKh0ZYhlNH1h7R/zwu+ayGcFt+pDAKDAAiAAAiAAAiAAAiAQLAIYCF+sPSN1oIACIAACIAACIAACICA7wjAafGdyiAwCIAACIAACIAACIAACASLAJyWYOkbrQUBEAABEAABEAABEAAB3xGA0+I7lUFgEAABEAABEAABEAABEAgWATgtwdI3WgsCIAACIAACIAACIAACviMAp8V3KoPAIAACIAACIAACIAACIBAsAnBagqVvtBYEQAAEQAAEQAAEQAAEfEcATovvVAaBQQAEQAAEQAAEQAAEQCBYBOC0BEvfaC0IgAAIgAAIgAAIgAAI+I4AnBbfqQwCgwAIgAAIgAAIgAAIgECwCMBpCZa+0VoQAAEQAAEQAAEQAAEQ8B0BOC2+UxkEBgEQAAEQAAEQAAEQAIFgEfj/aL6UBm/8jm0AAAAASUVORK5CYII="
    }
   },
   "cell_type": "markdown",
   "id": "e5b89348-77b5-4a9d-ac56-5971becc0392",
   "metadata": {},
   "source": [
    "Recall that the composite SST anomaly pattern associated with these events resembles the Northeast Pacific \"blob\" of 2014-16:\n",
    "\n",
    "![image.png](attachment:b26ea220-4576-47ef-b3e9-8c9e9fa81b1a.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3d5a5d2-2918-47ef-84c5-45388b33c803",
   "metadata": {},
   "source": [
    "Now let's take things slightly further. Often, we care not only about the _composite_ patterns across all heat wave events, but about individual properties of those events. Some relevant examples:\n",
    "- Temporal persistence: how long does a heat wave last?\n",
    "- Maximum temperature: how hot does it get?\n",
    "- Cumulative thermal stress: what is the \"effective\" combination of heat and time during the heat wave event?\n",
    "\n",
    "To accomplish this, we must be able to identify individual marine heat wave events within a time series, and perform appropriate operations on the portion of the time series within those events. \n",
    "\n",
    "The `scipy` package in Python has some excellent built-in tools to identify \"events\" within a time series: here we'll use the `scipy.ndimage.label` functionality. This treats the input time series as an \"image\", which can be n-dimensional (hence the name), and \"labels\" parts of the input data which are considered to be connected to one another in some way. There are various options for specifying that connectivity, but here we'll just use the simplest one: a data point is connected to the preceding point if it is equal to the value at that preceding point.\n",
    "\n",
    "Using the previously defined mask for marine heat waves, `mhw_msk`, this command identifies MHW events:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf91ffe1-0f94-4581-baa0-0264605460df",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Determine number of events\n",
    "events, n_events = ndimage.label(mhw_msk.to_array())\n",
    "n_events"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3643c97e-b8e3-4136-9785-8b0907f7ae3b",
   "metadata": {},
   "source": [
    "Note that we had to convert the `mhw_msk` object to a \"DataArray\" first using the `to_array()` method in order for this to work correctly.\n",
    "\n",
    "Using this set of events, we can now do things like collect all the temperature data within the events. Before we start doing that, though, let's specify an output variable to hold the information: the Python \"dictionary\" data type is useful here since it will allow us to store many different aspects of MHW behavior."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55085afc-bf96-4bfb-8051-1bc60e038a38",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define empty variable to put MHW data into (type \"dictionary\")\n",
    "mhw = {}\n",
    "\n",
    "# Specify start, end time variables within the dataset\n",
    "mhw['time_start'] = [] \n",
    "mhw['time_end'] = [] \n",
    "mhw['duration'] = [] \n",
    "mhw['dates'] = []\n",
    "\n",
    "# Specify variables to hold MHW parameters\n",
    "mhw['temp_mhw'] = []\n",
    "mhw['time_peak'] = []\n",
    "mhw['intensity_max'] = []\n",
    "mhw['intensity_mean'] = []\n",
    "mhw['intensity_cumulative'] = []\n",
    "mhw['duration'] = []"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef1c2265-4f74-48c1-a636-dd7fb3d97613",
   "metadata": {},
   "source": [
    "Now we can go ahead and compile the relevant information into `mhw`. First, let's simply identify and store the starting and ending times associated with each MHW event, with the requirement that there is some minimum duration it must exceed in order to be considered. Here I'll use the 5-day threshold suggested by Hobday et al. (2016):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d83cc8b-5451-4c0a-aad7-ad6a799d2857",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find all MHWs, collect temperature data within them into a list\n",
    "minDuration=5    # set some minimum duration for the MHWs in order to count them\n",
    "\n",
    "# Loop over set of all events\n",
    "for ev in range(1,n_events+1):\n",
    "    # Grab the part of the time series within the event, figure out how long it is\n",
    "    event_duration = (events == ev).sum()\n",
    "    # Figure out whether the event is long enough to count\n",
    "    if event_duration < minDuration:\n",
    "        continue\n",
    "    #print('Marine heat wave length = ' + str(event_duration) + ' days')\n",
    "    # If event is long enough, add data to an array\n",
    "    mhw['time_start'].append(nepac_ts.time[np.where(events == ev)[1][0]])\n",
    "    mhw['time_end'].append(nepac_ts.time[np.where(events == ev)[1][-1]])\n",
    "    mhw['duration'].append(event_duration)\n",
    "\n",
    "#print(mhw['time_start'][0])\n",
    "#print(mhw['time_start'])\n",
    "\n",
    "# Add number of events to MHW dictionary\n",
    "mhw['n_events'] = len(mhw['time_start'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0119050d-2e12-46d5-8187-a5ca1190a7c6",
   "metadata": {},
   "source": [
    "Next, we can apply this time information to extract specific portions of the temperature time series. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9a09371-4119-429c-8421-d50c728c434f",
   "metadata": {},
   "outputs": [],
   "source": [
    "for ev in range(mhw['n_events']):\n",
    "    # Get starting, ending point within the time series for the MHW event in question\n",
    "    tt_start = np.where(nepac_ts.time==mhw['time_start'][ev])[0][0]\n",
    "    tt_end = np.where(nepac_ts.time==mhw['time_end'][ev])[0][0]\n",
    "    \n",
    "    # Extract time series information\n",
    "    temp_mhw = nepac_ts.anom[tt_start:tt_end+1]\n",
    "    mhw['temp_mhw'].append(temp_mhw)\n",
    "    mhw['dates'].append(nepac_ts.time[tt_start:tt_end+1])\n",
    "    \n",
    "    # Get peak value = date of maximum temp during MHW\n",
    "    tt_peak = temp_mhw.argmax('time')\n",
    "    \n",
    "    # Time/date of maximum value    \n",
    "    mhw['time_peak'].append(mhw['time_start'][ev] + tt_peak)\n",
    "    \n",
    "    # MHW Duration\n",
    "    mhw['duration'].append(len(temp_mhw))\n",
    "    \n",
    "    # MHW intensity \n",
    "    mhw['intensity_max'].append(temp_mhw[tt_peak])\n",
    "    mhw['intensity_mean'].append(temp_mhw.mean())\n",
    "    mhw['intensity_cumulative'].append(temp_mhw.sum())\n",
    "  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dabdc16-ea2f-4428-90b5-279ad900c240",
   "metadata": {},
   "source": [
    "### Summarizing MHW Properties: tables, PDFs, filled line plots\n",
    "\n",
    "OK, now let's take a deeper dive into the properties of these events. The below code selects the MHW event with the longest duration, then displays the parameters associated with that event.\n",
    "\n",
    "Once we've put in this code, try playing around with it to return the properties of:\n",
    "- the event with the largest maximum intensity\n",
    "- the event with the largest cumulative intensity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05a2dd4f-ea6a-4352-8461-0cbd57afd644",
   "metadata": {},
   "outputs": [],
   "source": [
    "ev = np.argmax(mhw['duration']) # Find longest event\n",
    "dates = mhw['dates'][ev].values\n",
    "print('Maximum intensity:', mhw['intensity_max'][ev].values, 'deg. C')\n",
    "print('Average intensity:', mhw['intensity_mean'][ev].values, 'deg. C')\n",
    "print('Cumulative intensity:', mhw['intensity_cumulative'][ev].values, 'deg. C-days')\n",
    "print('Duration:', mhw['duration'][ev], 'days')\n",
    "print('Start date:', mhw['time_start'][ev].values)\n",
    "print('End date:', mhw['time_end'][ev].values)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73d6d19e-8c8f-4bfd-90a4-1cd5c42645af",
   "metadata": {},
   "source": [
    "Oftentimes, though, a more effective way to present this information is via a plot. The below code generates a line plot of the SST anomaly associated with the longest-duration MHW, and shades in the part of the time series above the MHW threshold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76452fd8-9493-4c58-84ab-36f6cc244588",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a blank figure\n",
    "plt.figure(figsize=(14,10))\n",
    "# Plot the SST anomaly time series\n",
    "plt.plot(nepac_ts.time, nepac_ts.anom, 'k-', linewidth=2)\n",
    "# Plot a horizontal line showing the MHW anomaly threshold\n",
    "plt.plot(nepac_ts.time, np.repeat(thr.anom.values,len(nepac_ts.time),axis=0), 'g-', linewidth=2)\n",
    "# Set the x and y limits of the plot\n",
    "plt.xlim(('2019-06-01','2020-01-31'))\n",
    "plt.ylim(0,2.2)\n",
    "\n",
    "# Find indices for MHW of interest and shade accordingly\n",
    "# index of MHW start\n",
    "t1 = np.where(nepac_ts.time==mhw['time_start'][ev])[0][0]\n",
    "# index of MHW end\n",
    "t2 = np.where(nepac_ts.time==mhw['time_end'][ev])[0][0]\n",
    "# Fill in parts of time series above the threshold\n",
    "plt.fill_between(nepac_ts.time[t1:t2+1].values, nepac_ts.anom[t1:t2+1], np.repeat(thr.anom.values,len(nepac_ts.time[t1:t2+1]),axis=0), color='r')\n",
    "\n",
    "# Set plot axis labels etc\n",
    "plt.ylabel('SST [$^\\circ$C]')\n",
    "plt.title('SST (black), threshold (green), detected MHW events (shading)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "582bea76-45a7-40bf-9fea-6ac964d5d90a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform kernel density estimation on mean duration values\n",
    "mhw_dur=np.array(mhw['duration'])   # duration is a list - convert to numpy array\n",
    "\n",
    "kde = stats.gaussian_kde(mhw_dur)\n",
    "\n",
    "# Display as PDF\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "x_eval = np.linspace(0.6, 120, num=200)\n",
    "plt.plot(x_eval, kde(x_eval), 'k-')\n",
    "plt.xlabel('Duration (days)')\n",
    "plt.ylabel('Probability density')\n",
    "plt.title('PDF of MHW durations')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4ca2e9f-4d83-4044-bd5c-3cd8b1daca4d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "89e3e026-dc06-41f7-ae47-81f0bc4476a0",
   "metadata": {},
   "source": [
    "### Identifying marine heat wave impacts\n",
    "\n",
    "Now that we've taken a look at the properties of MHWs themselves, let's combine this with additional data to determine the impacts of those events. \n",
    "\n",
    "To make things easy for ourselves, we'll use the monthly-resolution current data from the NCEP Global Ocean Data Assimilation System (GODAS), which is available on the same NOAA-hosted THREDDS server as the OISSTv2 information. This will let us visualize the differences in ocean circulation associated with the MHW!\n",
    "\n",
    "Important things we'll need to consider:\n",
    "- _Temporal resolution_. Here we'll be comparing the *daily* OISSTv2 information with the *monthly* average currents; so we'll want to pick longer-duration MHW events!\n",
    "- _Spatial resolution_. To directly compare SST and currents, it's useful to *regrid* data to the same spatial resolution!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36ed3401-8982-47ca-8d1d-dfc2fabc9521",
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOAA THREDDS server: base URL for all files \n",
    "base_url=\"https://psl.noaa.gov/thredds/dodsC/Datasets/godas/\"\n",
    "\n",
    "# Naming convention for files: \"sst.day.mean.[year].v2.nc\"\n",
    "yrs=range(2017,2021)\n",
    "print(yrs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dfc9a4e-2a2c-41af-a46b-5c7e57e83764",
   "metadata": {},
   "source": [
    "Read in the GODAS current data; this uses the same syntax that we applied for OISSTv2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccd7ed2f-feb9-4662-b532-7f4f0bd7584d",
   "metadata": {},
   "outputs": [],
   "source": [
    "nepac_uvel=[]     # Blank variable to store output data \n",
    "nepac_vvel=[]\n",
    "\n",
    "for y in yrs:\n",
    "    print(y)\n",
    "    \n",
    "    # Open the dataset into an xarray object\n",
    "    myuvel=xr.open_dataset(\"https://psl.noaa.gov/thredds/dodsC/Datasets/godas/ucur.\"+str(y)+\".nc\")\n",
    "    myvvel=xr.open_dataset(\"https://psl.noaa.gov/thredds/dodsC/Datasets/godas/vcur.\"+str(y)+\".nc\")\n",
    "    \n",
    "    # Define logical mask to be true in the Northeast Pacific (30-50N, 180-240E)\n",
    "    # NOTE: masks the same for U and V since they're on the same grid\n",
    "    mask_lon = (myuvel.lon >= 180) & (myuvel.lon <= 240)\n",
    "    mask_lat = (myuvel.lat >= 30) & (myuvel.lat <= 50)\n",
    "    # Read in data \n",
    "    nepac_uvel.append(myuvel.where(mask_lon & mask_lat, drop=True).squeeze())\n",
    "    nepac_vvel.append(myvvel.where(mask_lon & mask_lat, drop=True).squeeze())\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95274931-7ac7-4357-824d-fbfbe258ee21",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Concatenate data into a single xarray Dataset for u, v velocities\n",
    "uvel_fldts=xr.concat(nepac_uvel,dim='time')\n",
    "vvel_fldts=xr.concat(nepac_vvel,dim='time')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4f83415-ee7b-4804-a649-34b4d1fe504f",
   "metadata": {},
   "source": [
    "Now let's take a quick look at the SST and velocity datasets. It becomes pretty obvious quickly that the coordinates are quite different between the two datasets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f75dc4f-5dc0-4604-93a3-ead4f5054525",
   "metadata": {},
   "outputs": [],
   "source": [
    "# NE Pacific SST: OISSTv2\n",
    "nepac_fldts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40c845bc-49f5-4a5a-9a3a-418d707e0513",
   "metadata": {},
   "outputs": [],
   "source": [
    "# NE Pacific zonal (east-west) current: GODAS\n",
    "uvel_fldts"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e795da2b-2de5-4d45-ad43-bffe22bcbf55",
   "metadata": {},
   "source": [
    "There are several different ways of interpolating data in Python. Here, let's use the xarray interpolation functionality: this is a wrapper to the interpolate method within scipy, so functions essentially the same way.\n",
    "\n",
    "Here, we'll interpolate the OISSTv2 data to the GODAS grid:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b19104c8-a199-4536-9766-a0ead278f42d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Interpolate SST to GODAS grid\n",
    "ssta_int = nepac_fldts.interp(lat=uvel_fldts.lat, lon=uvel_fldts.lon)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af0248df-a7da-461c-ad1a-849d380f5607",
   "metadata": {},
   "source": [
    "You can easily see the difference between the two datasets by making maps:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "916eac7d-ddd8-4fa7-bd33-eb881d36322d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a blank figure\n",
    "plt.figure(figsize=(14,10))\n",
    "# Plot raw data\n",
    "plt.pcolor(nepac_fldts['anom'].mean(dim='time'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47d8a8b2-33e1-461a-900d-5b14fd9f3ab5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a blank figure\n",
    "plt.figure(figsize=(14,10))\n",
    "# Plot interpolated data\n",
    "plt.pcolor(ssta_int['anom'].mean(dim='time'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c22f6809-473f-4573-8162-76f3aecd14ae",
   "metadata": {},
   "source": [
    "But what we *really* want to do is to: first look at how velocity changes during a marine heat wave relative to the average (\"climatology\"), then make a plot of the total anomaly.\n",
    "\n",
    "Let's get started with that, by taking our interpolated current fields and converting them to anomalies using the procedure we applied to the ERSSTv5 last week. We can apply groupby to accomplish this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "619b72f7-d99a-41af-bd65-0edd0458e43d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the climatology and remove it from the current time series\n",
    "# Zonal (east-west) current\n",
    "\n",
    "\n",
    "# Meridional (north-south) current\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "604e0d70-9e8a-4835-959d-bbe15efba835",
   "metadata": {},
   "source": [
    "Now we need to pull out the portion of the current anomaly fields that correspond to the MHW event of interest. Here, we'll use the 'longest-duration' event identified above, since that will contain multiple months of data and therefore we can get a reasonable sampling from GODAS."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82497cb0-f67b-4e55-a7be-c468f82c96ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "#### Currents: GODAS\n",
    "# index of MHW start\n",
    "   # pick first time exceeding the MHW start time\n",
    "# index of MHW end\n",
    "   # pick last time below the MHW end time\n",
    "\n",
    "#### SST: OISSTv2\n",
    "# index of MHW start\n",
    "   # pick first time exceeding the MHW start time\n",
    "# index of MHW end\n",
    "    # pick last time below the MHW end time\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c8a9e5b-9ae4-4d09-9338-989fa52fe7c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate composites over the MHW\n",
    "# zonal current\n",
    "\n",
    "\n",
    "# meridional current\n",
    "\n",
    "\n",
    "# SST\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8708d6de-7602-4be7-b9ea-11c320bd124f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "c6d06bc0-6eb6-4aaa-9257-1d292edf7689",
   "metadata": {},
   "source": [
    "Finally, we make a plot! Let's include both SST and current anomalies here, and plot them on the same map:\n",
    "- SST in colors\n",
    "- currents as vectors (arrows)\n",
    "\n",
    "For the current information, we'll use the `quiver` method in matplotlib, which draws arrows having an x-component equal to the east-west current, and a y-component equal to the north-south current."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c130d61-a556-4bd0-95f6-d214024ffcf4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an object containing the Plate-Carree map projection, centered at 180E\n",
    "\n",
    "# Create a figure object\n",
    "\n",
    "# Create an axis object using the specified projection\n",
    "\n",
    "# Make a filled contour plot of the MHW SSTA using the specified map projection\n",
    "\n",
    "\n",
    "# Add arrows for currents\n",
    "\n",
    "\n",
    "# Set plot parameters\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1451694b-b382-4489-80de-9871e81dda44",
   "metadata": {},
   "source": [
    "In the plot above, it is hard to read the arrows! There are too many of them, so Python scales them automatically to a size that's difficult to read. To get around this problem, we can plot only a subset of data. Let's do an example plotting only every third point in the current arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1764759a-2584-48ea-b54f-31a16d6407cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an object containing the Plate-Carree map projection, centered at 180E\n",
    "\n",
    "# Create a figure object\n",
    "\n",
    "# Create an axis object using the specified projection\n",
    "\n",
    "# Make a filled contour plot of the MHW SSTA using the specified map projection\n",
    "\n",
    "\n",
    "# Add arrows for currents\n",
    "\n",
    "# Set plot parameters\n",
    "ax.gridlines()\n",
    "ax.set_extent([180, 240, 30, 50], crs=ccrs.PlateCarree())\n",
    "coast_10m = cfeature.NaturalEarthFeature('physical', 'land', '10m',\n",
    "                                        edgecolor='k', facecolor='white')\n",
    "ax.add_feature(coast_10m)\n",
    "ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,\n",
    "                  linewidth=2, color='gray', alpha=0.5, linestyle='--')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf5d1b08-7bd3-4c4e-a2f9-f9cf02e8bf38",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d62ffbec-7613-4588-9486-127806cc3c47",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}