Fetch_SNOTEL_SouthPlatte.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fetch and export SNOTEL sites and daily time series data\n",
    "2/29/2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import OrderedDict\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "from shapely import geometry as sgeom\n",
    "import ulmo\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "##USER INPUTS## - most of this can be read directly from the .ctl file or the .par file \n",
    "#-> need to read in text file\n",
    "#start date\n",
    "st = \"2000-10-01\"\n",
    "#end date\n",
    "ed = \"2020-09-29\"\n",
    "\n",
    "\n",
    "#Snotel bounding box\n",
    "Bbox = {\n",
    "    'latmax' : 39.825325,\n",
    "    'latmin' : 38.507589,\n",
    "    'lonmax': -104.39097, \n",
    "    'lonmin': -106.503086,}\n",
    "\n",
    "# Snotel projection\n",
    "stn_proj = 'epsg:4326' #WGS84\n",
    "#model projection\n",
    "mod_proj = 'epsg:32612' #UTM zone 12N"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Use ulmo to get SNOTEL station metadata within bounding box\n",
    "## Get SNOTEL metadata over modeling domain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#function to extract data\n",
    "def sites_asgdf(ulmo_getsites, crs=stn_proj):\n",
    "    \"\"\" Convert ulmo.cuahsi.wof.get_sites response into a point GeoDataframe\n",
    "    \"\"\"\n",
    "    \n",
    "    # Note: Found one SNOTEL site that was missing the location key\n",
    "    sites_df = pd.DataFrame.from_records([\n",
    "        OrderedDict(code=s['code'], \n",
    "        longitude=float(s['location']['longitude']), \n",
    "        latitude=float(s['location']['latitude']), \n",
    "        name=s['name'], \n",
    "        elevation_m=s['elevation_m'])\n",
    "        for _,s in ulmo_getsites.items()\n",
    "        if 'location' in s\n",
    "    ])\n",
    "\n",
    "    sites_gdf = gpd.GeoDataFrame(\n",
    "        sites_df, \n",
    "        geometry=gpd.points_from_xy(sites_df['longitude'], sites_df['latitude']),\n",
    "        crs=crs\n",
    "    )\n",
    "    return sites_gdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert the bounding box dictionary to a shapely Polygon geometry using sgeom.box\n",
    "box_sgeom = sgeom.box(Bbox['lonmin'], Bbox['latmin'], Bbox['lonmax'], Bbox['latmax'])\n",
    "box_gdf = gpd.GeoDataFrame(geometry=[box_sgeom], crs=stn_proj)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# WaterML/WOF WSDL endpoint url \n",
    "wsdlurl = \"http://hydroportal.cuahsi.org/Snotel/cuahsi_1_1.asmx?WSDL\"\n",
    "\n",
    "# get dictionary of snotel sites \n",
    "sites = ulmo.cuahsi.wof.get_sites(wsdlurl)\n",
    "\n",
    "#turn sites as geodataframe \n",
    "snotel_gdf = sites_asgdf(sites)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>code</th>\n",
       "      <th>longitude</th>\n",
       "      <th>latitude</th>\n",
       "      <th>name</th>\n",
       "      <th>elevation_m</th>\n",
       "      <th>geometry</th>\n",
       "      <th>easting</th>\n",
       "      <th>northing</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>335_CO_SNTL</td>\n",
       "      <td>-105.777893</td>\n",
       "      <td>39.803921</td>\n",
       "      <td>Berthoud Summit</td>\n",
       "      <td>3444.239990234375</td>\n",
       "      <td>POINT (-105.77789 39.80392)</td>\n",
       "      <td>947142.300804</td>\n",
       "      <td>4.419059e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>938_CO_SNTL</td>\n",
       "      <td>-106.113068</td>\n",
       "      <td>39.303501</td>\n",
       "      <td>Buckskin Joe</td>\n",
       "      <td>3398.52001953125</td>\n",
       "      <td>POINT (-106.11307 39.30350)</td>\n",
       "      <td>921455.974142</td>\n",
       "      <td>4.361858e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>415_CO_SNTL</td>\n",
       "      <td>-106.170952</td>\n",
       "      <td>39.489540</td>\n",
       "      <td>Copper Mountain</td>\n",
       "      <td>3215.639892578125</td>\n",
       "      <td>POINT (-106.17095 39.48954)</td>\n",
       "      <td>915353.758295</td>\n",
       "      <td>4.382249e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>936_CO_SNTL</td>\n",
       "      <td>-105.593452</td>\n",
       "      <td>39.656269</td>\n",
       "      <td>Echo Lake</td>\n",
       "      <td>3230.8798828125</td>\n",
       "      <td>POINT (-105.59345 39.65627)</td>\n",
       "      <td>963935.465583</td>\n",
       "      <td>4.403600e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>485_CO_SNTL</td>\n",
       "      <td>-106.196808</td>\n",
       "      <td>39.379910</td>\n",
       "      <td>Fremont Pass</td>\n",
       "      <td>3474.719970703125</td>\n",
       "      <td>POINT (-106.19681 39.37991)</td>\n",
       "      <td>913778.691339</td>\n",
       "      <td>4.369955e+06</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          code   longitude   latitude             name        elevation_m  \\\n",
       "0  335_CO_SNTL -105.777893  39.803921  Berthoud Summit  3444.239990234375   \n",
       "1  938_CO_SNTL -106.113068  39.303501     Buckskin Joe   3398.52001953125   \n",
       "2  415_CO_SNTL -106.170952  39.489540  Copper Mountain  3215.639892578125   \n",
       "3  936_CO_SNTL -105.593452  39.656269        Echo Lake    3230.8798828125   \n",
       "4  485_CO_SNTL -106.196808  39.379910     Fremont Pass  3474.719970703125   \n",
       "\n",
       "                      geometry        easting      northing  \n",
       "0  POINT (-105.77789 39.80392)  947142.300804  4.419059e+06  \n",
       "1  POINT (-106.11307 39.30350)  921455.974142  4.361858e+06  \n",
       "2  POINT (-106.17095 39.48954)  915353.758295  4.382249e+06  \n",
       "3  POINT (-105.59345 39.65627)  963935.465583  4.403600e+06  \n",
       "4  POINT (-106.19681 39.37991)  913778.691339  4.369955e+06  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CSO_gdf = gpd.sjoin(snotel_gdf, box_gdf, how=\"inner\")\n",
    "CSO_gdf.drop(columns='index_right', inplace=True)\n",
    "CSO_gdf.reset_index(drop=True, inplace=True)\n",
    "#CSO_gdf.index = CSO_gdf['code']\n",
    "\n",
    "#add columns with projected coordinates \n",
    "CSO_proj = CSO_gdf.to_crs(mod_proj)\n",
    "CSO_gdf['easting'] = CSO_proj.geometry.x\n",
    "CSO_gdf['northing'] = CSO_proj.geometry.y\n",
    "\n",
    "CSO_gdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAD4CAYAAADl7fPiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUV0lEQVR4nO3dfYwcd33H8ffHlwNf0riHyKXNnR1sVHpAbdkmK9dgQaljYpfSxI0VNRYPrVIpVZXykBYjW00qua2aUCPUQFtFVqIIVHAKwXFJCDWmtomIgp27nB/xOcVQEp/T+qh0BZer45hv/9g5c3Oevdu924fZ289LWnl3dn+z3xnvfm7mN7/ZUURgZjZmTqMLMLN8cSiYWYpDwcxSHApmluJQMLOUKxpdQJZrrrkmFi5c2OgyzGat/v7+H0VEV9ZzuQyFhQsX0tfX1+gyzGYtST8s9Zx3H8wsxaFgZikOBTNLcSiYWYpDwcxScnn0wZrLroEhtu0+yZmRUbo7O9i0tpf1y3saXZZN05RbCpLmSjoo6bCk45K2JtOXSnpW0lFJT0iaV6L93Um7Y5J2SJpb7YWwxtk1MMSWnUcZGhklgKGRUbbsPMqugaFGl2bTVM7uw3lgdUQsBZYB6yStBB4CNkfEEuBxYNPEhpJ6gI8AhYhYDLQBt1ep9pa2a2CIVffvZdHmr7Hq/r0N+xJu232S0QsXU9NGL1xk2+6TDanHZm7KUIiic8nD9uQWQC/wdDJ9D7ChxCyuADokXQFcCZyZUcWWq7/OZ0ZGK5pu+VdWR6OkNkmHgLPAnog4ABwDbk5echuwYGK7iBgCPgW8CLwM/E9EfKPEe9wpqU9S3/DwcMUL0kry9Ne5u7OjoumWf2WFQkRcjIhlwHxghaTFwB3AXZL6gauBVya2k/Q64BZgEdANXCXpAyXeY3tEFCKi0NWVOSTbEnn667xpbS8d7W2paR3tbWxa21v3Wqw6KjokGREjwH5gXUQMRsRNEXEDsAM4ldFkDfCDiBiOiAvATuAdMyvZ8vTXef3yHu67dQk9nR0I6Ons4L5bl/joQxOb8pCkpC7gQkSMSOqg+EX/pKRrI+KspDnAPcCDGc1fBFZKuhIYBW4EfKbTDG1a28uWnUdTuxDT+etcrUOJ65f3OARmkXK2FK4D9kk6AjxHsU/hSWCjpBeAQYqdh48ASOqW9BRA0vfwGPA8cDR5v+1VX4oWU42/znnqrLR8UR5/zblQKIRPna6tVffvZSijD6Kns4NnNq9uQEVWT5L6I6KQ9ZyHObeoPHVWWr44FFpUnjorLV8cCi3KhxKtFJ8Q1aLGOiV9IpNN5FBoYT6UaFm8+2BmKd5SqDP/9oDlnUOhjsYGDI2NRBwbMAQ4GCw3vPtQR3k6u9GsFIdCHXnAkDUDh0IdecCQNQOHQh15wJA1A3c01pEHDFkzcCjUWZ4GDPnwqGVxKLQoHx61Utyn0KJ8eNRKcSi0KB8etVK8+9Ciujs7Mn95Ka+HR93/UT/eUmhRpQ6P/uabu3Jx5anx/HuS9eVQaFFZP/664YYevtI/lLsvn/s/6su7Dy1s4uHRVffvLfnla+Smuvs/6stbCnZJXr98Hh5eXw4FuySvXz4PD68vh4Jdktcvny9NV18t2afgw1vZ8nxuRp6Gh892LRcKHt47OX/5rOV2H3x4y2xyU4aCpLmSDko6LOm4pK3J9KWSnpV0VNITkuaVaN8p6TFJg5JOSHp7tReiEnntYTfLi3K2FM4DqyNiKbAMWCdpJfAQsDkilgCPA5tKtH8A+NeIeDOwFDgx46pnIK897GZ5MWUoRNG55GF7cgugF3g6mb4H2DCxbbL18C7g4WRer0TEyMzLnr689rCb5UVZfQqS2iQdAs4CeyLiAHAMuDl5yW3AgoymbwSGgUckDUh6SNJVJd7jTkl9kvqGh4crXY6y+fCW2eQUEeW/WOqkuKvwYeBV4DPA64GvAh+JiNdPeH0B+A6wKiIOSHoA+HFE3DvZ+xQKhejr66tkOcysApL6I6KQ9VxFRx+STf/9wLqIGIyImyLiBmAHcCqjyWngdLJlAfAY8LZK3tPM6mvKcQqSuoALETEiqQNYA3xS0rURcVbSHOAe4MGJbSPiPyW9JKk3Ik4CNwLfrfIyWAN4ANjsVc6WwnXAPklHgOco9ik8CWyU9AIwCJwBHgGQ1C3pqXHtPwx8IWm/DPibKtZvDeDfN5jdKupTqBf3KeTbqvv3Zv5qU09nB89sXt2AiqxSVetTMAMPAJvtHApWMQ8Am90cClYxDwCb3VruLEmbuTyfYm0z51CwafEp1rOXdx/MLMWhYGYpDgUzS3EomFmKQ8HMUhwKZpbiQ5J14rMKrVk4FOrAPytvzcS7D3Xgn5W3ZuJQqAOfVWjNxKFQBz6r0JqJQ6EOfFahNRN3NNaBzyq0ZuJQqBOfVWjNwrsPZpbiLYUW5cFUVopDoQV5MJVNxrsPLciDqWwyDoUW5MFUNhmHQgvyYCqbjEOhBXkwlU1mylCQNFfSQUmHJR2XtDWZvlTSs5KOSnpC0rxJ5tEmaUDSk9Us3qZn/fIe7rt1CT2dHYji5d7uu3WJOxkNKO/ow3lgdUSck9QOfFvS14HPAh+PiG9JugPYBNxbYh4fBU4AJYPD6suDqayUKbcUouhc8rA9uQXQCzydTN8DbMhqL2k+8NvAQzOu1sxqrqw+hWTz/xBwluKl6A8Ax4Cbk5fcBiwo0fzvgE8AP5tRpWZWF2WFQkRcjIhlwHxghaTFwB3AXZL6gauBVya2k/Q+4GxE9E/1HpLulNQnqW94eLiSZcilXQNDrLp/L4s2f41V9+9l18BQo0syK0tFRx8iYgTYD6yLiMGIuCkibgB2AKcymqwCbpb0H8CjwGpJ/1Ri3tsjohARha6urkrKyp2xEYNDI6MEPx8x6GCwZlDO0YcuSZ3J/Q5gDTAo6dpk2hzgHuDBiW0jYktEzI+IhcDtwN6I+ED1ys8njxi0ZlbOlsJ1wD5JR4DnKPYpPAlslPQCMAicAR4BkNQt6alaFdwMqj1i0LsiVk9THpKMiCPA8ozpDwAPZEw/A7w3Y/p+irses153ZwdDGQEwnRGDPnnJ6s0jGmugmiMGvSti9eZTp2ugmj+/5pOXrN4cCjVSrRGD1dwVMSuHdx9yzicvWb15SyHn/EvQVm8OhSbgk5esnrz7YGYpDgUzS3EomFmK+xRywtdhsLxwKOSAhzJbnnj3IQc8lNnyxKGQAx7KbHniUMgBX4fB8sShkAMeymx54o7GHPBQZssTh0JOeCiz5YV3H8wspem2FDzIx6y2mioUPMjHrPaaavfBg3zMaq+pQsGDfMxqr6lCwYN8zGqvqULBg3zMaq+pOho9yMes9poqFMCDfMxqral2H8ys9sq56vRcSQclHZZ0XNLWZPpSSc9KOirpCUnzMtoukLRP0omk7UdrsRBmVj3lbCmcB1ZHxFJgGbBO0krgIWBzRCwBHgc2ZbR9FfiziHgLsBK4S9Jbq1K5mdXElKEQReeSh+3JLYBe4Olk+h5gQ0bblyPi+eT+T4ATgDsEzHKsrD4FSW2SDgFngT0RcQA4BtycvOQ2YMEU81hI8ZL2B0o8f6ekPkl9w8PD5VVvZlVXVihExMWIWAbMB1ZIWgzcQXF3oB+4GnilVHtJvwB8BfhYRPy4xHtsj4hCRBS6uroqXAwzq5aKjj5ExAiwH1gXEYMRcVNE3ADsAE5ltZHUTjEQvhARO2dWrpnVWjlHH7okdSb3O4A1wKCka5Npc4B7gAcz2gp4GDgREZ+uYt1mViPlbClcB+yTdAR4jmKfwpPARkkvAIPAGeARAEndkp5K2q4CPgislnQoub236kthZlWjiGh0DZcpFArR19fX6DLMZi1J/RFRyHrOIxrNLMWhYGYpDgUzS3EomFmKQ8HMUhwKZpbiUDCzFIeCmaU4FMwsxaFgZikOBTNLcSiYWYpDwcxSHApmluJQMLMUh4KZpTgUzCzFoWBmKQ4FM0txKJhZikPBzFKuaHQBlm3XwBDbdp/kzMgo3Z0dbFrby/rlvgyn1Z5DIYd2DQyxZedRRi9cBGBoZJQtO48COBis5rz7kEPbdp+8FAhjRi9cZNvukw2qyFqJQyGHzoyMVjTdrJocCjnU3dlR0XSzanIo5NCmtb10tLelpnW0t7FpbW+DKrJWUs5Vp+dKOijpsKTjkrYm05dKelbSUUlPSJpXov06SSclfU/S5movwGy0fnkP9926hJ7ODgT0dHZw361L3MlodTHlBWaTy8lfFRHnJLUD3wY+CnwW+HhEfEvSHcCiiLh3Qts24AXgPcBpilet3hgR353sPX2BWbPamtEFZqPoXPKwPbkF0As8nUzfA2zIaL4C+F5EfD8iXgEeBW6psH4zq6Oy+hQktUk6BJwF9kTEAeAYcHPyktuABRlNe4CXxj0+nUzLeo87JfVJ6hseHi6zfDOrtrJCISIuRsQyYD6wQtJi4A7gLkn9wNXAKxlNlTW7Eu+xPSIKEVHo6uoqq3gzq76Kjj5ExAiwH1gXEYMRcVNE3ADsAE5lNDlNegtiPnBmeqWaWT2Uc/ShS1Jncr8DWAMMSro2mTYHuAd4MKP5c8CbJC2S9BrgduCrVardzGqgnC2F64B9ko5Q/JLviYgngY2SXgAGKf71fwRAUrekpwAi4lXgT4DdwAngSxFxvPqLYWbVMuUhyUbwIUmz2prRIUkzay0OBTNLcSiYWYpDwcxSHApmluJQMLMUh4KZpTgUzCzFoWBmKQ4FM0txKJhZikPBzFIcCmaW4svGmc0C1bz2qEPBrMlV+9qj3n0wa3LVvvaoQ8GsyVX72qMOBbMmV+1rjzoUzJpcta896o5GsyY31pnoow9mdsn65T1VuwCxdx/MLMWhYGYpDgUzS3EomFmKQ8HMUhwKZpZSzlWn50o6KOmwpOOStibTl0n6jqRDkvokrSjR/u6k3TFJOyTNrfZCmFn1lLOlcB5YHRFLgWXAOkkrgb8FtkbEMuAvkscpknqAjwCFiFgMtFG8HL2Z5dSUg5eieFnqc8nD9uQWyW1eMv0XKV6OvtR7dEi6AFw5yevMLAfKGtEoqQ3oB34F+IeIOCDpY8BuSZ+iuMXxjontImIoef5FYBT4RkR8o8R73AncCXD99ddPY1HMrBrK6miMiIvJbsJ8YIWkxcAfA3dHxALgbuDhie0kvQ64BVgEdANXSfpAiffYHhGFiCh0dXVNa2HMbOYqOvoQESPAfmAd8PvAzuSpLwNZHY1rgB9ExHBEXEhef9kWhZnlRzlHH7okdSb3Oyh+0Qcp9g38RvKy1cC/ZzR/EVgp6UpJAm4ETlShbjOrkXL6FK4DPpf0K8wBvhQRT0oaAR6QdAXwfyT9AZK6gYci4r1J38NjwPPAq8AAsL0Gy2FmVaLiwYV8KRQK0dfX1+gyzGYtSf0RUch6ziMazSzFoWBmKQ4FM0txKJhZikPBzFIcCmaW4lAwsxSHgpmlOBTMLKWlLwaza2CoalfVMZstWjYUdg0MsWXn0UuX8B4aGWXLzqMADgZraS27+7Bt98lLgTBm9MJFtu0+2aCKzPKhZUPhzMhoRdPNWkXLhkJ3Z0dF081aRcuGwqa1vXS0t6WmdbS3sWltb4MqMsuHlu1oHOtM9NEHs7SWDQUoBoNDwCytZXcfzCybQ8HMUhwKZpbiUDCzFIeCmaXk8ifeJQ0DP5ziZdcAP6pDOeVwLZfLSx3gWrK8ISIyr8+Yy1Aoh6S+Ur9bX2+uJb91gGuplHcfzCzFoWBmKc0cCnm6JqVruVxe6gDXUpGm7VMws9po5i0FM6sBh4KZpUVEw2/AbcBx4GdAYcJzW4DvASeBteOmv4bi/tkLwCCwIWO+C4FR4FBye3DcczcAR5N5fwZQDet4D9CfvF8/sHrcc/uTeY7VeG0t18kU7S9bJzOoJXO5JrR9/7jnDyXzX1aqfQ3rqOhzUuN1UvFnperfx0YHQrKwbwF6k4UujJv+VuAw8FpgEXAKaEue2wr8dXJ/DnBNif/sYyXe8yDwdoph8HXgt2pYx3KgO7m/GBia8B9dqOM6maz9ZetkBrVkLtckn4ElwPcnWy+1qqPSz0mNa6n4s1LtWy52HyLiRERk/WLqLcCjEXE+In5AMX1XJM/dAdyXtP9ZRJQ9SkzSdcC8iHg2imv788D6WtUREQMRcSZ5eByYK+m1k9VYw3WS2b7UOplBLZXaCOyY7AV1quOSRqyT6XxWqi0XoTCJHuClcY9PAz2SOpPHfyXpeUlflvRLJeaxSNKApG9Jeue4+Z6eON8a1zFmAzAQEefHTXtE0iFJ90rSFO1nWktmeypfJ5PNa0wly/V7XB4K5bavRh3V+JxUq5YxM/2sTEvdQkHSNyUdy7jdMlmzjGlB8Rej5gPPRMTbgGeBT2W89mXg+ohYDvwp8EVJ+4B/BlaO1QB8juI+a63qKM5E+jXgk8AfJY+/CfxyMv/XAJ8AXqrxOinV/h+5fJ28a5q1ALw/IpYA70xuHyw5E+nXgZ9GxLGxzwmXr5dP17COrM/JPBq7TlKflUrbz0Tdfo4tItZMo9lpYMG4x/OBM8B/Az8FHk+mfxn4w4z3PA+cT+73SzoFbAKGgH0RsRhA0kbg3bWqI3mP+cnrPhQRp5Ka1kx4zR9Q3Gf8lxrWUqr9LWSsk2nWQkQMJf/+RNIXKW5Cf77EfG4n2UrI+pyMrZda1VHic/KrNGidZH1WKmk/U3nfffgqcLuk10paBLwJOJjs3z3Bz7/INwLfndhYUpektuT+G5P234+Il4GfSFqZbIJ9CJjsP3qmdXQCXwO2RMQz46ZfIema5H478D7gWC3XySTtK10nJedVyXJJmkOxJ//RcdMqXS8zqqOKn5Nq1NJJ9T4r01PrnsxybsDvUkzY88B/AbvHPffnFHtwT5L0/CbT3wA8DRwB/o3i5h/AzcBfJvc3UOysOQw8D/zOuPaFZKWeAv6e4mZfreq4B/hf0offrgWuonjY6UhS5wP8vKe6JrVM0f6ydTKdWqZYrom1vBv4zoTPQ2b7WtVBhZ+TWq4TpvFZqfbNw5zNLCXvuw9mVmcOBTNLcSiYWYpDwcxSHApmluJQMLMUh4KZpfw/zeMtINV5B1cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "CSO_gdf.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get SNOTEL SWE data in domain over time period of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Emilio code\n",
    "\n",
    "def fetch(sitecode, variablecode, start_date, end_date):\n",
    "    print(sitecode, variablecode, start_date, end_date)\n",
    "    values_df = None\n",
    "    try:\n",
    "        #Request data from the server\n",
    "        site_values = ulmo.cuahsi.wof.get_values(\n",
    "            wsdlurl, 'SNOTEL:'+sitecode, variablecode, start=start_date, end=end_date\n",
    "        )\n",
    "        #Convert to a Pandas DataFrame   \n",
    "        values_df = pd.DataFrame.from_dict(site_values['values'])\n",
    "        #Parse the datetime values to Pandas Timestamp objects\n",
    "        values_df['datetime'] = pd.to_datetime(values_df['datetime'])\n",
    "        #Set the DataFrame index to the Timestamps\n",
    "        values_df.set_index('datetime', inplace=True)\n",
    "        #Convert values to float and replace -9999 nodata values with NaN\n",
    "        values_df['value'] = pd.to_numeric(values_df['value']).replace(-9999, np.nan)\n",
    "        #Remove any records flagged with lower quality\n",
    "        values_df = values_df[values_df['quality_control_level_code'] == '1']\n",
    "    except:\n",
    "        print(\"Unable to fetch %s\" % variablecode)\n",
    "    \n",
    "    return values_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "335_CO_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "938_CO_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "415_CO_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "936_CO_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "485_CO_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "1057_CO_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "505_CO_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "248_CO_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<suds.sax.document.Document object at 0x7fb5589ff278>\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Unable to fetch SNOTEL:WTEQ_D\n",
      "248_CO_SNTL has no data\n",
      "531_CO_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "935_CO_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "970_CO_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "602_CO_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "937_CO_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "1014_CO_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "939_CO_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "1100_CO_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "802_CO_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "842_CO_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n"
     ]
    }
   ],
   "source": [
    "stn_swe = pd.DataFrame(index=pd.date_range(start=st, end=ed))\n",
    "\n",
    "for sitecode in CSO_gdf.code:\n",
    "    try:\n",
    "        swe = fetch(sitecode, variablecode='SNOTEL:WTEQ_D', start_date=st, end_date=ed)\n",
    "        stn_swe[sitecode] = swe.value\n",
    "    except:\n",
    "        print(sitecode, 'has no data')\n",
    "        stn_swe[sitecode] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>335_CO_SNTL</th>\n",
       "      <th>938_CO_SNTL</th>\n",
       "      <th>415_CO_SNTL</th>\n",
       "      <th>936_CO_SNTL</th>\n",
       "      <th>485_CO_SNTL</th>\n",
       "      <th>1057_CO_SNTL</th>\n",
       "      <th>505_CO_SNTL</th>\n",
       "      <th>248_CO_SNTL</th>\n",
       "      <th>531_CO_SNTL</th>\n",
       "      <th>935_CO_SNTL</th>\n",
       "      <th>970_CO_SNTL</th>\n",
       "      <th>602_CO_SNTL</th>\n",
       "      <th>937_CO_SNTL</th>\n",
       "      <th>1014_CO_SNTL</th>\n",
       "      <th>939_CO_SNTL</th>\n",
       "      <th>1100_CO_SNTL</th>\n",
       "      <th>802_CO_SNTL</th>\n",
       "      <th>842_CO_SNTL</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2000-10-01</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000-10-02</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000-10-03</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000-10-04</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000-10-05</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.3</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            335_CO_SNTL  938_CO_SNTL  415_CO_SNTL  936_CO_SNTL  485_CO_SNTL  \\\n",
       "2000-10-01          0.0          0.0          0.0          0.0          0.0   \n",
       "2000-10-02          0.0          0.0          0.0          0.0          0.0   \n",
       "2000-10-03          0.0          0.0          0.0          0.0          0.0   \n",
       "2000-10-04          0.0          0.0          0.0          0.0          0.0   \n",
       "2000-10-05          0.0          0.0          0.0          0.0          0.3   \n",
       "\n",
       "            1057_CO_SNTL  505_CO_SNTL  248_CO_SNTL  531_CO_SNTL  935_CO_SNTL  \\\n",
       "2000-10-01           NaN          0.0          NaN          0.0          0.0   \n",
       "2000-10-02           NaN          0.0          NaN          0.0          0.0   \n",
       "2000-10-03           NaN          0.0          NaN          0.0          0.0   \n",
       "2000-10-04           NaN          0.0          NaN          0.0          0.0   \n",
       "2000-10-05           NaN          0.0          NaN          0.0          0.0   \n",
       "\n",
       "            970_CO_SNTL  602_CO_SNTL  937_CO_SNTL  1014_CO_SNTL  939_CO_SNTL  \\\n",
       "2000-10-01          0.0          0.0          0.0           NaN          0.0   \n",
       "2000-10-02          0.0          0.0          0.0           NaN          0.0   \n",
       "2000-10-03          0.0          0.0          0.0           NaN          0.0   \n",
       "2000-10-04          0.0          0.0          0.0           NaN          0.0   \n",
       "2000-10-05          0.0          0.0          0.0           NaN          0.0   \n",
       "\n",
       "            1100_CO_SNTL  802_CO_SNTL  842_CO_SNTL  \n",
       "2000-10-01           NaN          0.0          0.0  \n",
       "2000-10-02           NaN          0.0          0.0  \n",
       "2000-10-03           NaN          0.0          0.0  \n",
       "2000-10-04           NaN          0.0          0.0  \n",
       "2000-10-05           NaN          0.0          0.0  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stn_swe.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>335_CO_SNTL</th>\n",
       "      <th>938_CO_SNTL</th>\n",
       "      <th>415_CO_SNTL</th>\n",
       "      <th>936_CO_SNTL</th>\n",
       "      <th>485_CO_SNTL</th>\n",
       "      <th>1057_CO_SNTL</th>\n",
       "      <th>505_CO_SNTL</th>\n",
       "      <th>248_CO_SNTL</th>\n",
       "      <th>531_CO_SNTL</th>\n",
       "      <th>935_CO_SNTL</th>\n",
       "      <th>970_CO_SNTL</th>\n",
       "      <th>602_CO_SNTL</th>\n",
       "      <th>937_CO_SNTL</th>\n",
       "      <th>1014_CO_SNTL</th>\n",
       "      <th>939_CO_SNTL</th>\n",
       "      <th>1100_CO_SNTL</th>\n",
       "      <th>802_CO_SNTL</th>\n",
       "      <th>842_CO_SNTL</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2000-10-01</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.00000</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000-10-02</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.00000</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000-10-03</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.00000</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000-10-04</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.00000</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2000-10-05</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.00762</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            335_CO_SNTL  938_CO_SNTL  415_CO_SNTL  936_CO_SNTL  485_CO_SNTL  \\\n",
       "2000-10-01          0.0          0.0          0.0          0.0      0.00000   \n",
       "2000-10-02          0.0          0.0          0.0          0.0      0.00000   \n",
       "2000-10-03          0.0          0.0          0.0          0.0      0.00000   \n",
       "2000-10-04          0.0          0.0          0.0          0.0      0.00000   \n",
       "2000-10-05          0.0          0.0          0.0          0.0      0.00762   \n",
       "\n",
       "            1057_CO_SNTL  505_CO_SNTL  248_CO_SNTL  531_CO_SNTL  935_CO_SNTL  \\\n",
       "2000-10-01           NaN          0.0          NaN          0.0          0.0   \n",
       "2000-10-02           NaN          0.0          NaN          0.0          0.0   \n",
       "2000-10-03           NaN          0.0          NaN          0.0          0.0   \n",
       "2000-10-04           NaN          0.0          NaN          0.0          0.0   \n",
       "2000-10-05           NaN          0.0          NaN          0.0          0.0   \n",
       "\n",
       "            970_CO_SNTL  602_CO_SNTL  937_CO_SNTL  1014_CO_SNTL  939_CO_SNTL  \\\n",
       "2000-10-01          0.0          0.0          0.0           NaN          0.0   \n",
       "2000-10-02          0.0          0.0          0.0           NaN          0.0   \n",
       "2000-10-03          0.0          0.0          0.0           NaN          0.0   \n",
       "2000-10-04          0.0          0.0          0.0           NaN          0.0   \n",
       "2000-10-05          0.0          0.0          0.0           NaN          0.0   \n",
       "\n",
       "            1100_CO_SNTL  802_CO_SNTL  842_CO_SNTL  \n",
       "2000-10-01           NaN          0.0          0.0  \n",
       "2000-10-02           NaN          0.0          0.0  \n",
       "2000-10-03           NaN          0.0          0.0  \n",
       "2000-10-04           NaN          0.0          0.0  \n",
       "2000-10-05           NaN          0.0          0.0  "
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#convert SNOTEL units[in] to SnowModel units [m]\n",
    "for sitecode in CSO_gdf.code:\n",
    "    # overwrite the original values (no use for the original values in inches)\n",
    "    stn_swe[sitecode] = 0.0254 * stn_swe[sitecode]\n",
    "\n",
    "stn_swe.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABzqklEQVR4nO19eZwcR3n2UzOzt3Z1W7csWZZlS5ZP2QYCtoHPGHPZsQmxuUIIMc4XQvjiAOYDEi6TOBDAXPFHMBCIE24HAzK+wNj4lmxJlizJuixpda6OvXeO7q7vj+rq6aO6uqpntD2zquf329/uzkzXvN1d/dZbz3sRSikMDAwMDJofuawFMDAwMDCoD4xCNzAwMJggMArdwMDAYILAKHQDAwODCQKj0A0MDAwmCIxCNzAwMJggMArdoGlACFlICBkmhOTd/x8mhLwva7l04MpcJIQ8kvL4NvcaVAghn6u3fAbNDaPQDcYNhJCXCCFjhJAhQkg/IeRxQshNhBCleUgp3UMpnUQptWuUo5UQ8q+EkF5XOe4ihHzZfe9jhJDVoc9vi3ntevdvSggZccfiPx+RiPABSumlaWSnlJYopZMA3JXmeIOJjULWAhicdHgzpfRBQshkAJcBuB3AJQD+fBxl+BiAVQAuBnAAwKkAuIJ9BMAthJA8pdQmhMwG0ALggtBrp7uf5TiXUrp9/E7BwCAKY6EbZAJK6QCl9B4AfwrgzwghZwMAIeSNhJDnCCGDhJC9hJBP8WMIIYtcazhgiLg0xDFCyErfa6e4u4GZgq+/CMDdlNL9lOElSun33feeAVPg57n/XwrgdwC2hl7bQSndX9tV8GR9iRDyYULIBtfSv5MQMosQcq+7m3mQEDK1Ht9lMLFhFLpBpqCUPg2gF8Cr3JdGALwbwBQAbwTwV4SQaxLGKAH4IYB3+l6+AcCDlNI+wSFPAvg7Qsj/JoSsJIQQ31hlAE+harFfCuBRAH8IvZaKA5fgOgBXADgDwJsB3Avg/wKYAfacfrDO32cwAWEUukEjYD+AaQBAKX2YUvo8pdShlG4A8N9g1EwS/gPA2318/LsA/CDms/8E4DYA7wCwBsA+Qsif+d7/ParK+1VgCv3R0Gu/D435rOsX4D9XKsjsx9copYcopfvc73qKUvqcu1jdDeB8zfEMTkIYhW7QCJgH4BgAEEIuIYT8jhDSRwgZAHATmJUqBaX0KTDr/jJCyJlgHPc9MZ+1KaXfoJT+EdhO4FYA3yGEnOV+5BEAr3RpjpmU0m0AHgfwCve1sxG10C+glE7x/dyndQWAQ76/xwT/T9Icz+AkhFHoBpmCEHIRmEL/g/vSf4Ep4gWU0skA7gBAYg4P4z/AaJd3AfgppbSYdACldIxS+g0AxwEsd19+AsBkADcCeMz93CDYTuJGAPsppbsUZTIwGDcYhW6QCQghPYSQN4Fx3/9JKX3efasbwDFKaZEQcjGAt2sM+wMAfwym1L8f9yFCyIcIIZcTQjoIIQWXbukG8BzAlDwYFfN3YPQHxx/c1+rNnxsY1AVGoRuMN35JCBkCsBfAxwF8CcGQxf8N4DPuZ/4BwI9VB6aU9gJ4FgBFUBGHMQbgXwEcBHAEwF8DuI5SutP3md8DOAXVnQPcMU+BWKGvD8Whf0VVbgODeoGYBhcGEwmEkO+AUSKfyFoWEQgh9wN4OYA1lNJXpzi+DYxfbwHwL5TST9dZRIMmhlHoBhMGhJBFANYBON9w3AYnIwzlYjAhQAj5LICNAL5glLnByQpjoRsYGBhMEBgL3cDAwGCCILPiXDNmzKCLFi3K6usNDAwMmhJr1649QikV1SjKTqEvWrQIa9asyerrDQwMDJoShJDdce8ZysXAwMBggsAodAMDA4MJAqPQDQwMDCYITMciA4OTFJVKBb29vSgWE2uYGWSA9vZ2zJ8/Hy0tLcrHGIVuYHCSore3F93d3Vi0aBF8PT4MGgCUUhw9ehS9vb1YvHix8nGGcjEwOElRLBYxffp0o8wbEIQQTJ8+XXv3ZBS6gcFJDKPMGxdp7o1R6AYGKXFspIx7nz+QtRgGBh6MQjcwSImbfrAWf3XXszg8ZJyKBo0Bo9ANDFJi/8AYAKBYdjKWpLlx++234+yzz8aKFSvwla98BQDwyU9+Eueccw7OO+88vO51r8P+/fsBAC+99BI6Ojpw3nnn4bzzzsNNN90kHXt4eBjvf//7sWTJEqxYsQKXXnopnnrqKQDMKXz11Vdj6dKlWLJkCf72b/8W5XI5dqzR0VG84x3vwMqVK3H22Wfjla98JYaHhwEweuTmm2/2PvvFL34Rn/rUp3Drrbd6subzee/vr371q/jUpz6FL37xi7VcugiMQjcwSImWPHt8Ko5R6GmxceNG/Pu//zuefvpprF+/Hr/61a+wbds2fPjDH8aGDRuwbt06vOlNb8JnPvMZ75glS5Zg3bp1WLduHe644w7p+O973/swbdo0bNu2DZs2bcL3vvc9HDlyBJRSXHvttbjmmmuwbds2vPjiixgeHsbHP/7x2LFuv/12zJo1C88//zw2btyIO++80wspbGtrw89//nMcOXIkcMzHP/5xT9aOjg7v7w9+8IM1XLV4mLBFA4OUKOSY08qym78E9ad/uQkv7B+s65jL5/bgH9+8QvqZzZs342Uvexk6OzsBAJdddhnuvvtufOQjH/E+MzIykspBuGPHDjz11FO46667kMuxxfe0007Daaedhoceegjt7e348z9n3Q/z+Ty+/OUvY/Hixfj0pz/tyePHgQMHcOqpp3r/L1u2zPu7UCjgxhtvxJe//GXceuut2rLWC8ZCNzBIiQK30G1joafF2WefjUceeQRHjx7F6OgoVq9ejb179wJg1u2CBQtw1113BSz0Xbt24fzzz8dll12GRx+Nbx27adMmj+oQvXfhhRcGXuvp6cHChQuxfft24Xjvfe97cdttt+HlL385PvGJT2Dbtm2B9//6r/8ad911FwYGBpTPv94wFrqBQUq05l0L3Wl+Cz3Jkj5ROOuss/DRj34UV1xxBSZNmoRzzz0XhQJTS7feeituvfVW/NM//RO+/vWv49Of/jTmzJmDPXv2YPr06Vi7di2uueYabNq0CT09PVrfSykVWv1xrwPAeeedh507d+L+++/Hgw8+iIsuughPPPEEzjrrLABsQXj3u9+Nr371q+jo6NC8EvWBkoVOCHk9IWQrIWQ7IeQWwfuXE0IGCCHr3J9/qL+oBgaNBWOh1wd/8Rd/gWeffRaPPPIIpk2bhqVLlwbef/vb346f/exnABhXPX36dADAhRdeiCVLluDFF18UjrtixQqsX78ejsDHsWLFikj57sHBQezduxdLliyJlXXSpEm49tpr8c1vfhPvfOc7sXr16sD7H/rQh3DnnXdiZGQk+cRPABIVOiEkD+AbAK4CsBzADYSQ5YKPPkopPc/9+YzgfQODCYW8y6EbhV4bDh8+DADYs2cPfv7zn+OGG24I0Bn33HMPzjzzTABAX18fbNsGAOzcuRPbtm3DaaedJhx3yZIlWLVqFf7xH/8RvNXmtm3b8Itf/AKvfe1rMTo6iu9///sAANu2cfPNN+M973mPkD8HgMceewzHjx8HAJTLZbzwwgsBTh0Apk2bhre97W248847016OmqBioV8MYDuldCeltAzghwCuPrFiGRg0Plo9C735KZcscd1112H58uV485vfjG984xuYOnUqbrnlFpx99tk455xzcP/99+P2228HADzyyCM455xzcO655+Ktb30r7rjjDkybNi127G9/+9s4ePAgTj/9dKxcuRJ/+Zd/iblz54IQgrvvvhs/+clPsHTpUpxxxhlob2/H5z//+dixduzYgcsuuwwrV67E+eefj1WrVuG6666LfO7mm2+ORLvE4XOf+xzmz5/v/dSKxCbRhJC3Ang9pfR97v/vAnAJpfQDvs9cDuBnAHoB7Afw95TSTYKxbgRwIwAsXLjwwt27YxtvGBg0PP7ie8/goS2H8f/edSGuXDE7a3G0sXnzZo//NWhMiO4RIWQtpXSV6PMqFrrIQxBeBZ4FcCql9FwAXwPwP6KBKKXfopSuopSumjlT2BLPwKBpkJtAYYsGEwMqUS69ABb4/p8PZoV7oJQO+v5eTQj5JiFkBqVUbd9hYNCE4LtbGrFvDMYbl1xyCUqlUuC1H/zgB1i5cqX2WPfddx8++tGPBl5bvHgx7r777ppkHA+oKPRnACwlhCwGsA/A9QDe7v8AIWQ2gEOUUkoIuRjM8j9ab2ENDBoFD289jAc3M2deM0ctysL0mgk8nb8euPLKK3HllVfWbby0SKLDRUhU6JRSixDyAQD3AcgD+A6ldBMh5Cb3/TsAvBXAXxFCLABjAK6naaQxMGgSvOe7z3h/N+tUb29vx9GjR01N9AYEb3DR3t6udZxSYhGldDWA1aHX7vD9/XUAX9f6ZgODCYIm1eeYP38+ent70dfXl7UoBgLwFnQ6MJmiBgYpcO78yVjfy1K8nSbV6C0tLVrtzQwaH6aWi4FBCnS3Vxv3Nqk+N5iAMArdwCAF/FZ5s1roBhMPRqEbGKSAX4kbdW7gR8V2sPXgUCbfbRS6gUEK+I3yZo1yMTgxuPXXm3HlVx7B3mOj4/7dRqEbGKSAX4U3cxy6Qf2xdjcr4HVsJL6d3YmCUegGBingt8qNgW7gh1sRIhPfilHoBgYp4H9WjVPUwA+epJXFzs0odAODFDBOUYM45LykW2OhGxg0BRzjFDWIgbHQDQyaDLZjOHQDMbiBnsW8MArdwCAF/G3nDIdu4Aevc7a/f2zcv9sodAODFAgq9AwFMWg4PPMSC1v80I/Wjft3G4VuYJACO/qqXd0Nh27AcWBg/K1yP4xCNzBIge72Ai48dSoAw6EbVHHLz57P9PuNQjcwSAHHoThj1iQApgWdQRV+Ki4LGIVuYJACFYeiNc8en2bm0D/zyxfwbw/vyFqMCQO/g/w9r1g07t9vGlwYGKSAZTtoLTCF3syUy3ce2wUA+KvLl2QsycTAkzuPeX9P62od9+83FrqBgSZsh8Kh8BR6s4YtHhkuZS2CQZ1hFLqBgSY4T9qS5xZ6cyr0J3YczVqECQ2TWGRg0ASwXNK8rZAH0LyUS6FadCSQ+WqQDrf9ZkvWIhiFbmCgi4rFLPQq5ZKlNOlRyFcf/7KVbXTGREDYuZxF9JNR6AYGmqg4rkLPMwu3WcMWC/mqhT5UqmQoSfND5I8wlIuBQRPAstmT2uwWum1XBb/41ocwWDRKPS0ODhQjr2XhWzEK3cBAE9wpWsjlWCGmJiXRLSdIs/SPGIWeFgNj0WtnyucaGDQBKq5lW8gTEDSvhV6xg4ITEvNBg0T8asP+yGumBZ2BQROAW+it+RxyhDRtHHrYQjcKPT1a81FVaix0A4MmgOVZ6EyhN6c6r54HR5OuSw0By6GY0tmCP3z01d5rhkNvQBweKmJH33DWYjQEjo+UsfXgUNZiZI7DQ8wBVsgTgDRvpqgVMiGb9TwaAQ6laMnnMH9qZ+C18YaSQieEvJ4QspUQsp0QcovkcxcRQmxCyFvrJ2K2eOVtv8Nr//X3WYvREHjT1/6AK7/ySNZiZI4HNx8GAMyc1MYaAjepHrRClQHDCt5AHbZDA4laQINSLoSQPIBvALgKwHIANxBClsd87jYA99VbyCxhEi6q2JdBS61GhO046GrN4+x5k0HQvBx62ClqskXTw3IociSs0BvTQr8YwHZK6U5KaRnADwFcLfjc3wD4GYDDdZTPoAHRrLVL6oXth4cxq6cdAJAjzcs9l0LGilHo6eE4NJCoBTRuYtE8AHt9//e6r3kghMwD8McA7pANRAi5kRCyhhCypq+vT1dWgwZBOeMi/lmjkMvh+GgZAEAIadqwxcFixdRzqRMshyJPwgq9MS10UTBTWNKvAPgopdSWDUQp/RaldBWldNXMmTMVRWwMnOxWqR/Oya3PAQDLZncDYKF+zUq5DI5VMLmjxfs/6247zQyHUuRDHHoWs0KlwUUvgAW+/+cDCEfRrwLwQ8JWqBkA3kAIsSil/1MPIRsBtmBLdbKiWRVYvWA5DjoL7NEJ86bNhJGShc62PI66/a5P9vtaCyw7qtCzuJ4qCv0ZAEsJIYsB7ANwPYC3+z9AKV3M/yaEfA/AryaSMgcYzVAQJA+cjLBP8gffdqoPbzNb6BWbBhJiwnHpBuoQWuiNyKFTSi0AHwCLXtkM4MeU0k2EkJsIITedaAEbBSbapYpdfSNZi5ApbFoNUcsR0rRO0bLteE06AMOh1wLLaR7KBZTS1QBWh14TOkAppe+pXazGg1HoVfzZd5/Gun94XdZiZAbLpshxCx3NbKE7aCv4LHSj0FPDFin0BnWKGgAYLFpZi5AaFdvBzT9ejz1HR+synqiy3MkEx2ehkyZO/a8YC71usH1RLt9614UAGpRyMWAYLTevQv/hM3vxs2d78aEfPVeX8ZrUIK0b/NtrFofenBekYtGAQjcWenr4LfTXrZiN2T3tqZ+Tf753C+5+rjfVsUahJ6C7jbFSIyVpRGZD45P/sxEA8Oye/mwFmSCIOEWblI0r2w5aCn4LvUlPBMDf/Wgd/v2RnZl9f5hyqcVZfsfvd+D//Gh9qmONQk9AeytrBPzZX72QsSTpwSMZOt1zSYMDA7Wn/fceH8VVtz+KF/YP1jxWlrADFjpp2hZ06/b2e230AGBTk94XSil+/tw+3Lp6c2Yy2KEol6yqcBqFnoC5k1mK9wsHmnOyA8A58ycDAEbL6XcZT+08VrMc9206hM0HBvGTtXuTP9zA8BdiatYGF8UKmwudrdW4iJYmDcstVrLfWTgCp2gaC71W+q457+A4oseXSdes2Hmk9jBDXrukFnAl0lpo7mlnBSiX5gxb3O8WWls5bzK+855VALKN1rntN1vwmi8+nOrYRnDSW6FqiyRlFc7e47XthBv2yTo0WMSiW36N323JttYXzwS8ePG0TOWoBS11yHCtR1o4V+gdLempn0aAE+LQm9EpymvRLJ01Ca9edgqAbJ3d//bwjtSGRyModDtUbTEt5bK/xoqmDavQN+4bAAB8/4mXMpWD35R6WKhZIVwmNd0YtSv0MZfyadatPQezxtg5NGvHoq0HWdOWyR0tIITUZWHKamHbfTT7RLdwaZC0TlH+rKY1whr2yeIWUNahVHySNmvyCBBMinJSXs96LAp8hFqStCilnqWfFfzWWK5JU/95P9F5UzsAwO2Nmn68VZ97AIs/thov1YHe04Vf7qxi6W0atNAJ0u14uOFEUtYIaniFnnWyA78pzbit5vAr0LSlb+thofN7WbTSK+T/fGoPzvzkb+oSdZMWfmusWcvnDrmJcj3tzEdU68J0ZJhROFsO6gcP+O9lGoPDPzdfedtvtY+vB8Idi9Lu3PjzmUvJkja8Qs/cQndvS7OG6FJKUbZZhx0gvWKuh0Lnk7VUQ1TC6g0HAAA7M6wn47fQmSXWfBp9sFhBayGHdtefUevCNLO7LfBbB1t8fWrTGBz+uXlgoKh9fD1gO9VyEABS95rlTUeIsGp5MhpWoXOOMi1FUC9wqqEZt9VAdYKMuPz15gPpmjzXQ6FXXFlqybrliTBZNtmwHMeX+t+cLUWHihZ62qshi7UuTHysNLelzedTSTPPGqGOu8hCTzMx9h4bdY9PJ0fDKnRuoe85Vp/6I2nB6Ypm3FYDUcXHoxv0x6mfY7WWrFueCFPJqFgapRQORShssfkmx+BYBd3t1ZBcxqGnPw++Y0lDkfqPSOOrqcfcrBWRTFGkMwJ5BFhXm1LdxAgaVqFzn8DhoVKmcvAFpRkfWqC6IL3mTBaaNmOS/pYYiHaITwP+sI7UYqG71lw9nLRpwL+3EKjlkokoNSFsodd6Hlyhp1Fi/mNSWegNUAlVlPqf5no+ufMoAESSlFTRsAq9UXBshFm0zU658FZjaSNM/A/a1M50yVZ8t/D0rvRZp1WFns1DzEMvO10LiqA2yzYrDBVFFnr68bgBlsZC9x+SZn42BOVCgz1F05SEGC5ZuP+FQ97xaZDOrh8HNIJF7J8oTUu5uA/IJFcBpY9yYRdgxqTW1ElWJY9Dt0EpTRWaxRV6Vhw6311wJ3NaSyxrjJRsnNJdza2otfMSv5dpuln5/WRplDOf45csnpaaUqwVth1yikJfZxzxsRFpOxs2rIW+91h2YWkcB/qrHvNmtMIAn0J3t9dpLXR+3NTO1tQRPzv7hr2/09bfaC24HHpWCr3EFLpnoTdpYlHJstHWUn38a/UFcF2WJojB/2ylqTdUshzkCKMTs4xDj4Qtaoqyy5cglTYYpGEV+od+tC5rEfCF+7d6fzepPkfJjfn2LPSUCn1Dbz8Axu2lXdz8h6Xl0bmFXmuKdFpwn86kNtdCR2PsJnVRspxAP9FcjdE6tThF/Yd8+YEXtY8vWTbaW/LI50h2Cj0Utpgm89bf2jFt396GVeiNgF+u3w8AOHteT9Nb6NwBVrbTRZj8bmsfgNq41orteM6esZSVH/lt+MbvdqQTokYcGmS7tq5WbqE352JftpyIhV4b5cJ+1+oUfSlFGn/JYq308jmSWQPzcNhimqhFfh3efO7c1AtTUyj0rGPR2wv5plXoPIOPh0FVrBrLc+bSW6Qly0G3u7CkpUy4Q5aXNR5vcD/AqdO7ADRvLZey5aA1Xy2SVmuUCz82zW31z6c0vpFSxUFbIc+MjYz8o0KnqOYF5X6qjpbcxFboWa26HIxmyFSE1OA1O3jJ2lqvZS3xymMV21PoaTOAOeXyxnPmpDq+Vuw4zPwAnW1Vp2gzLvbFEIcO1GeOWyk0qv97//i8eVrHHh4q4kdr9qK1kEM+l02pEMehoBTI53w7Hug7RbmR096ST/18NIdCz1ibps36agTwSzetqxVA+ms5b0oHrr1gXuoUcduhKFsOuttqC5/kC1JWU4LzpJNafU7RJpsbo2ULFZt6fhWAUya1n0iaxc1fInuhu/NRxcNbGBW4ct5k5HO5TEqF8DkZKCKaYudm2Q4IYUbLhHOK+pGFQuff+bevXYpcrjmtMKDaJ5I7wNKex2jZQldrIXURJ14hkVvov01Z555PhawW+dGyheldrZ5iT5sRmCX290frndRMubjqKw3l8pO11YbIuorsmBum+IU/OQf5jJ5TPhf9Fnqa5uFlm6Ill0MhRya4hZ7BTRrm1eg6WmpOi84S/AHzKJeUE2WkZKOrrZAqHAvwK3RmoX8pRTQDUH3gswpb3LR/EB2+3qxp44WzBN8dLZk5yXutXglStfq7dGV4ft8AWvIEHS155Ek2US5VhV59LU35XMt20JInNUWSNYVCz8Ip2jfMrJju9kLTlkgFqhZ6WyHv/p8u8aNss0iC1BZ6KNomLfh3Wxml/m8/PFxzvHHW4KGswSiX2iz0jftY2Vxd4ytsxeoenycEFZslqeVz6Z2JtcASWuj6maIV20Ehz6J1JraFnsFN8heQSrN9ahTwiVGLhV7xOVbThrfxDlSFGtvh8Qe+klE4w1jZxuIZVZ43a8pl474BvPd7z2j5JPhn2/JhBVQ7dOfXYDGYj6A7PUfLFlbM7QHAFqhamqekBTc4/VObEP2S22WbosVV6JSmM2QbVqH7HTZZUC78IZ3Z3VZznYssYYcVeoprafkKUjELXV8O3smmVqaEi5+FhV62HFgOxQULp3qvZW2h3/zj9fjtlsPY4cvCTcKGXra4hpt114Vy0RxjuBRS6JqTa7hkeTkBnS15lG1n3JW6Z6Hn/VEu+ha6ZTtozRMv/DHNs9qwCt2/Nc8kFMn9yhwhqWmGRgBXfO1coadQhJ5Cz+dSxdcCQN9QCe0tOXzofy3VPtYPPhfShMfVCq58uv20UcZzY7DIGiRzpaYCHuu9dFa391raLvUA8JM1e72/dZ/VoWKwwbPutdzQO4AuN4SUL1DjPTe4zP449DQU1nN7+5HPE+Tz6bNulRQ6IeT1hJCthJDthJBbBO9fTQjZQAhZRwhZQwh5pbYkIfg5pCwUOldaOdK8bcaAqIWeqjyp+4C05Enq3cpI2caktpbUZUE5+MNTrjFBKg248pnkq1JYgx6sC9L0V+ULtH8XXAvl8uGfbvD+1lfobJE8f+EUAHq7P9uhGC3bXm5CLeUHagHXVZFMUU0x2ltyKORyVQv9RCh0QkgewDcAXAVgOYAbCCHLQx97CMC5lNLzALwXwLe1JQnBr9CzoEvDFnqzcuh821bI5dCaz6FSQ8/Glnyuhm7mbDtZK4fOi3plYaFz5dPdHlSEWWr0agMWdSFsn7HCUa8EKd0x+CL5969bxo7XmJ/82MUz3azdHK/JriVCzeAyB2q5pKBcNu4bxIWnTq32Uz5BlMvFALZTSndSSssAfgjgav8HKKXDtKrxulDjFN9zdNSrQw5ky6ETUns3lyxRDaliyjRNM4BHtx0BABzoH0ttoVdsB62FnFey9aJFUxOOEIO3r8uCQxcp9KwzRXm6uM4zYrst9Pzli+vlC9DdAPJryuv161xLToEtmcHCL2up+FgLbIGFzkpkqI/B5/WktkJVoaeY4yoKfR6Avb7/e93XAiCE/DEhZAuAX4NZ6amxYV9/4P9sOHRuxaSnGRoBfodmSz5dJt1H3C31CweGUu9WKrbjbY2Xz+nB5I5W7TEs28Ev1u33xhtvcIuQZ7sC2fcU5Xy4zj2xQt11gPpF6+iOMRhS6Drz80fPMLXEF9i8Z6FnQ7mELXQdOZ7YwToVnTm721sYTpSFLtojR76JUno3pfRMANcA+KxwIEJudDn2NX19fbFfGFbg2XDo7HeOkMytsFpgOyydOOcq9FoaQ1x6xozUu5WyRatcZ8oCX/5Q0ixSvEVO0bRO4npgz9Fqv12d22rbUYWOOi1MujunwTG2SE5xi6594b6tso8HwGunX76MtVespclGLXBo1ELXXegPDbKyzKsWTfMWhhPlFO0FsMD3/3wA++M+TCl9BMASQsgMwXvfopSuopSumjlzZuwXhsOOdLmoesDx8YxZh6bVAstX1rMlT1L1Bl0+h8X5vvOSU5mDOMWaULYdtBSqzqs0i8JopRriVo8ep7roH3Ut9FByVFa7tx1HqqGKuhx6WKHXyxeg+6xuOTgEoNocWQe2QzG5o8XL3OXOxPF+Vvki5m8bpxtI8cIBFko6f2pH1UI/QQr9GQBLCSGLCSGtAK4HcI//A4SQ04m7PBJCLgDQCuCotjQuhsLJBhk6RUmThy3aDvUmWslycHy0knBEFC8cYFmAuVz6a1GxmFMUSB815O9ylEWT6G1upcVwL86sZoY/MUjn4Q/X7gbqSbnofT5Ml+jAcn0BHPzPtDv6YsXG8RH9FnZCCx3QWlkKbpZpe0u+pmidRIVOKbUAfADAfQA2A/gxpXQTIeQmQshN7seuA7CRELIOLCLmT2kN+1AeW/upN7NgmiyUadhCb2aFzifasZEyHnCb0KZF2t2Kn0NPuyjwlHUgmyiXztZgvDOQrjNNvVDwKXQdEYQcep0adeheC9ummN3Tnqq/rB06j1yNHPoff/NxnP/ZB7SPs3yBB54smhuewbEK5k3pAFDNpk6j0JWyESilqwGsDr12h+/v2wDcpv3tMRgqWuhqzWP+1E53/HqNrI5qHDpJTTM0AtbsPo6RlN2B/Hj7JQsBIHXlyYrteH1N0y4KJddC72rNZ5YpyssQc6QpwlQv+ItBaVEuAg49Te0REbQLUjk0dSirZUf7eALpd/Sb3Z1o31AJM7vblI+zBQpdt0TGYNFCj+sY5ueRxk/UkJmiO/qG0d3eAl7rJhML3Z0UzRyHfmykjHV7+2sag1IKQoDpriJLW8uF16kA0lvoPImmq62QSS2XkmUH+nAC9VOEaeB3I2iFLVLqbfH9qIcvQPe+hmkTHdgO9bIqAV/YYo3P6qhmv1uhQofe4jZUrHj0E783ac6jIRX60eEyxiq2tw3LknKpxqGPuwg1g0cQcLzstGk4/ZRJMZ8Wo2Kzbiytvmy8NLdj99ERb4y0iwLn/ye1FTKz0INdftIVYaoX/FtyHYMjTFUA9XP86z4nIvpH51j/wlSvsEVdefrdmuxhC11HjMGi5ZU74TZDmjleWy3TEwRCWGSFt4XKQJkGMkWbtMFF0eWcv/DWcwAAk9paMDCmZ33wMEeuyNJY146bos058LQFvniD5u72AoZKeudRDzy3t18QjZGdU9Sv0HWCfuI59HpQLpoWeig/gYcvqiC8MNXLANRdoHf0scJz0310nG6oM7PQWW2dXA3n0ZAWetltJszvVRZ0h8eh55q3lss9bhLOjEmMD8yniP/mIaR+C133Wuw+xuKlz1swpTpGGn7QnRCzetrHfU5QSrH76Ggkjj9LOs5Ps+hFuThihZ5SjhwBbrpsCTpa9Jup+5VyZ2teq2FIxQ7SNfk6GYC6Dnd+zgund3qv6W46hgIW+gkuzjXeKFssTbxhLPQm5dC/+fAOAMCCacx7niZax1PohfRNkY+NsKSJlfOnAGATNk3yh21X67KP95woudfh2vODSdL1ig5JA/+iqMPji8IW0yZI2Q6FQ1kceZo2dv48iZxm8EEkysX9c6drMaeFriItWSx5rzVUPlfnORkpWehyi6XVEq3TkAq95Cp0vlqbsEV98Ifzf1++BKefUt3K6U7WqkJPz6Fv2s+iB6Z1tnpjpLHQvWYd+dy43w8eYdMRKlObpVPUCnDo6seJOHQWh64vw7q9xwHAM8C0OXSbeuGXhOg5d3cdGQnVT2F/3/Sfa/WEQHBx1I0uKVVstLnNXzh0FnrHobAcGnjGgAmk0KvtzhrBKdqctVy44m73cb65nL4y9tqVFdJz6H1DzELnW9K0Frq/+9J4F2Di/oh2kVM0o7nhr2ejc0+EjsiUCVLX/dsTAFgWcprdm+Wjf1inHvXjWwu5QIOMnA5fE8L63n7v7zQWOm/x6JdFdZRqeWo2t2qhjhpToVsOWvNVhZ5NHDr73ay1XLwazaGwLt3zKAksdN0xKjZFaz7nFWDK50gqD77lK+M73reDh0y2hx5cguxquYz4lJluHfFolIs+rehfVLcfHnYjO3SdoiHKRePwLQeHsOrUad7/+Rq0WS11gooVO7LQ6zQ+4VnPrb6wXmACceg8q7BecaVp0Oy1XLwa5qHGtbqWcdkOKvQ0DuKK282cI5+SwuIPWksGlAsvO9AeinLJstriqC9hTCuxKCb1X/eSjvmaa5QsRzs7EuCJRVXKRfU8jg6zXV9nW/V+EGEdQTX4o2t2HFZv5wfEW+iqF4OXtObPyITj0DmvRmrYetSKZm9BV20bF8yk025cG2oonMaS47XQOfI5fS4fYOeUI+ycxntOeBZ6hHLJbrH3l0LQoaDEYYv6vgC/Qn/vHy1OtXvzLy46Fvre42MAWLnZeuPmn6zX+rzIQtepjTPiJjK1hDn0FPkNDanQK47jtjtj/2eTKco59OZ0inJeLly8SFcZc8qlGoee1kIPKvQ0CpknkmRBgY15Cj1siWUXARXo6lUXykXv+8fcHcK/vPUcrJw/OdXurWw5PoWuvjDxpLklM/US5eJQS4luptDT79yOj7Bz4YqcPyoTwkJ3HJaZmM8RH4eeJeVCXD5s3EUI4NuP7sSXH3hR+fNeF5V8jZSLF4eed8fQn2j+WuiAy6GnMD8s20HB7Ws63lPC26kURJbY+MrC4e9oo0u55HPh89A3WviuhRctS5OcVLKqypAt9IoKnTcb8VW+rIFxqamG+taDQ+KSEIpDHnMzTc+YxRanWuq6N5xC93t8q2Ukx18O/pDyhSXrOPTP/Xozbn9om7Ic/k5FHLkUlnE4bJGkCDkMUy5pqB+gShVkQYFVQr4EjjRURb3gf+B1U/8jHHoKC/3ZPSxkkTuKcynKIDD+Wd8/w0ts93TUJ9m9lqipyZ2t6A+V2dChXHiexhQ3rDdfgyHbcArdr4iyjEPnDwujXLJLHglDtQ64v7Ezh86WluOBFw4CqC0O3b+tBtiWMs0W9+ldx9DdVsiEAhNdTyDbxCI7JeVi+Wrkc6Q5j4MDTBGd68sA1l3cSr76ODr0FadcetrVSwXIUEsHLNtxcHqY+tG4nkWvimi1GikbV1+WxlXogbDF8X9i+HfmvTj07DT61x7a5v2tKoeI89XZ0nLwib5wGoshT2Mdbzow4G3LmRw57e1k/2gZLxwYRE9HSyalGO7dGFzYOFjY4vjKwhHk0HUs9GiFwzRd6vvHyuhuK3ilZtP4V0oV22fh61Eu+RwJzKtaUIuFLioBrBMTzxdmfks4GzYhwharlAvxVU8bfzmqF5nUpEAe33EEZ3zi3lSdUDj+1cedq1oSPKStyxfWldahuWxWt3cv0tA2e4+NhTh0/cm6o4+Fkr3lvLmZ1PhZ75YhXjS9K/B6lg0unJQWerjsLKDfpR5gLfkmdwYbZussLAcGxjBStgPRHarzYkPvADpb8sHsTOVvjqIWDt2yaXTnBvVrQT02oJpg5X9dBw2n0KuUS9Zx6Ow3D1sE0l3gb/5uB8qWgw37Buoil61IufCkk05fqjpJQblUbIqWQjClWdcBBwSjEfJEv7fp0WG2IL7y9BmZ1PjJEYJzF0wRR4eMnxgBBFP/a+TQUzhFH3jhUMApqUvbvHiILdLcGaiTyUwIieyWarkPtVAulh0tdqZDTfr9dfxYYKI4Rd0HvZAnmcah+6st1qJAvCSBOp2E6k3mmW+T2qoKPU1CTzjkUJdD51znMl+8sENZ/Wcd/NfTewCwypFZLPRDJQtnzorGPKdRhPWCnZJysRzqOd440iRItRVy6PbNL93ggSNuSYhz3aJtOnTec7uP4yy3eXk9UCvl0pKPXk/Vc4lQLhMp9d8Lt8vVZhnXSw6/hZ7mweX3+c+/90xd5FIN9xsuMUXq5xjT0CVRha53HQZchc7T/oF01hD/yrlTOjJpfNI3VIrQFEA6qqJe8M+FWqoUAukSpIZKFs5bOMX7X5fSu+P3rBro9En6HHx3eyHC+ddEudTKoYfDQDUWSC9E2hePD6RbZBpOofNJ6neKZmEBccoinyM1KZAwt1YrVB/cX64/ACAYBaBbzQ5wKZdwtmmNCn1WD3uAdSbsaNnCy0+b7skAjJ8idXxGRhTZFW5L3YJO4MRjqf/qY5Qsm/UtaAtRehpj5AjBzO42b27oHE8BzJ3cofxdSahFoQ+5Dlo/dBbIQM4Lauu81HAKnYfltfgSi7Jo8fVfT+8FEE5w0h/Hr8jqAVULvauN1af2O63yqYonBS10XQfxS0dZbeoOf5RLigVytGx7u43xplz4jmJWT3vkPXYqWVEu1blQews6vbMYdikz3gcTSNdH87IzZnr/5zVyHByavnWdCLxBNAC85dy5yscdHykLQ4l1FkiPQychDn1CWOi+sMVM49AdBxctmlpzkbBavOdiudTGs2yK5XODHGOaKJdyyIOvS4Px3Q23ygF/8SENOQLxyuPrWxE1AebIMkfB3xxCNw49yqHr7byGPIVeNRh0d2+sj2b4eLVjHYpAhEut4GG+C6Z1aC0UR93otXCvXh3Khc8vfjr8+UgzrxpOoR8eYn0jbcep6cRqxWjZxvkLpwKoTYHU2szYX4AJUFfoFQGvlyP6q364UqLutaiESgcASLVQ+7n88V7oLUFdHI4snaKOU11sdWRwhKn/utZ11ELXcZjbDsVwyQpkeur07nUcGmnzVouCHxyzMKunTTtwgJc/CO/edBzEXs5LiEOfEFEu/J7M6mnPLGyRUoqy7fi61KeXoxZuDgAGRoMpxcoK3XIinvfqAplOkQL6dIeXYekLfUxDuVR8OwWPAquBitt1ZCTQHEGGRAs9vRg1wR9doWuhRzh0TafokKCWig4HzimboI9Hx0KPUi4Lpqbn1EfKFrpaC9pdvbhCDzcPZ3HoamP4Q6SBdM8HR8MpdN7qq62QzyTeGGATnlJ/l570CqRSYyGaAwPFwP+qq7blOAILXZ+bCydN6DqIy8ISBCkoF9/CUutCTynFq7/4MN7znaeVPu81CxFZ6Clq29QLts9C112kRU2idekSIGqhq16KbYeHAAA9HX7KRY93DmdjnjZzEpbP6cHKeZPVhPBh77FRdLTmtbt67XefT3FZ5XRhi94zNhE4dK4A2lqy49CjfTTTy1Grhf7MS8cC/6tSOCwhKHh7vW7iGhNtX/+YkHJRvRRlK6rQ09zXvqESWuvQAACo1s5Ys/u40uerFrr4ccnKQvcrdNWHv2w5sByKrtaoRakDvrsJOEU1FPJ3H3sJADDfZ1XnNaqBOpRCxLDMmdyeqljacMnCWMXWDsu9+9leAFHKRSfJirrnEs4UnRBx6L1u4frWDGu5RBR6DQqklgw0IFqMS4fqaIlx7qiexnN7uMLzK3T2W5n6CbXXYmPoWSB7j40CqC72tSac8R6nqvCS3YSUS3aci+1Us3hVr8VoOZpBDOgnjImoBq0+mraDmd1teJkbigro7XYcQYGx6hiKQgTkoThzdrc25cKrLC5wax15ckD9Wtg0eC4TqgUdf9imdbVmRrmI2q6llcNKGVrGsetIsB2WygJBKcWm/YORGHjdJrpH3HT7d75sYWQMbQ7d34JO0wI57M6JK1fMdmVgr6dd6G9d/YL3d/9oco0dGYeuS1XUE5ZDvRaDqjJwZ6Y/gxjQPw+u0NsCCl19jKLlYO6UIOetmy4f59NIcz+GihVM62rVDsut2A5etXSGUA51ygWBqKMJ1YLOchxM7mhBV1shM6co5/HDTVtTWeiBJgT6soS/UmXV5oWswttX3fMYGGPKzh9LX6VLlIZAxXbctnHpHas73fOZ3lXNKNSRIYzdR0e9v3lHJhlEDbc5CLKjXPx15lWnJg+zCxsGuhsNkYWuowyLZRvthbDBoRHlEkO5pKkoSinFUNFCd3sL8jk9Q6FsOV7ZWz/YAqk2RtjBe8KdooSQ1xNCthJCthNCbhG8/w5CyAb353FCyLnakrgoVhzPesiqlkvZZpPVs9CRPnwybYlTjnAhJRWFPuzWcbnh4gWB13UUYbFi4/YHWdlev0LXrQRXDkXJAPqOVZ5tWi3Tyl5PG+O/5eCQ97eKDNVyFNHHRdeJVk+MVWx0uc+K6rXg923O5DDnq5d0Vqywhdq/89Lh0P2dijhyOXUOnQqcokC6XItihfkVetpbtGPpy5YTKRLG5VDl8iu2E2oVyYMXlMWoHpv0AUJIHsA3AFwFYDmAGwghy0Mf2wXgMkrpOQA+C+Bb+qIwjFXsQMF7YPw59FKo3RifN7rOlm2HhrDW53hLo9BLtoOF0zrxoxtfBkBNoXsNjQvRPoeqcnztt9s8D/6kUAEmNkay7ABQsaiwPRegngE8MFZBjgCnuAq9liiAMFQ6xfMiUnEWelaUy1i5qtB1LFuguq3n0I1D5300CQkqItUxihUnEhmio4wZ7xx9PQ0FVg3BLIAQAp3UkTiFrtO2MlzOmNsNJ8pCvxjAdkrpTkppGcAPAVzt/wCl9HFKKddcTwKYry2Jiy0HBgMF74HsolzaCuFUc71xrvjyI4H/0zhr+IThykRFoR90FXFbpKGxevjl4Fg1Rjv40LLfWs7ZSLSN+hiUUvy/3+/E5I4WTwlVW3QpiRBA+PqpyNA3zBS6qIwD0XAE1hujZRuT3Hr36gk97Hc4U1THoQkARZGFrcWhi49XD1uMZruyMfR3TIM+hZ7XkAFgYYvhhQlwDQXFYcLFvXSDBvxQUejzAOz1/d/rvhaHvwBwr+gNQsiNhJA1hJA1fX19woPbW/IoutmRmTlFw300PcqlNkHSLExcoevEkPPvCXdz0VHGp05nXvu/f90Zgdd16ZJwpqnuGIcGSyjbDqa6/RaB2iyYcBioyhiyDvM6NEO9MVaxPf5W9eH3Yp5DT76uZfviweFow2wNuqJUcSI7SFULnVKWJyLKDE2TDc0j63o69CgXfs3HylELiSWcKd4TO0ir5mvQeyoKXbQnFX4VIeTVYAr9o6L3KaXfopSuopSumjlzpugj6B+t4PwFU93x2GvjbqFHolzY62nFuPZ8tv6lVuj5nLeCq0S58Djr6V2tgdd1FOn0SezYN54TLFTEG01wx3ESRBy6DuXCSx/8zWtPjx6f4nqGIyNUZOgfjVaM5NClKuqJIOWidgy/ZpFaLtA7j9ZCLpI0RzTkKFpVapVD1cIPZ1YGxkjhFH2+dwAAMHNSm1bYIjc8l80WL/Sq14I3P/cfC5y41P9eAH7v2nwA+8MfIoScA+DbAK6mlB7VlgTs4d3XP4ZpXezBGe8yqRyehe6l/uvL4bfaVriZa2lW3K2HhtBayFWTghQG8fqJxljoKmLs7GNVEsOx119+kLXDu3fjAYVRmD8izDHqUC6lEP0F1OYsD6f7q8jwzO7jaG/JSZxf4w9KKUbLltdiULeZgrDDjsaZrNl9DBeeOjU6hirlUolSLrz8QNIY3qIk0F5pKJcjLqV25uxut56M2nFjZXHaP8B7zareEyfgnyGEaNFPfqgo9GcALCWELCaEtAK4HsA9/g8QQhYC+DmAd1FKXxSMoQSeQMKtjloKvdeCuExRnQnPlar/+DTn0dmax1jF1lLo8U5RdcuWL6an+Kok+qETnhalftTl8OKdC36Okf1OM+GHi/oK3XYcL947jKzi0AeLFhzKEoR0aJ/YmHqi4aQeraBYcXB8JFhnSLXZx1jZZk7RmEzmpLnFr3cc5aJ7Px7cfBhLZnZ5PRhUr+WYIBbfL4eqFGELnR2vl+DkHZf0AUqpBeADAO4DsBnAjymlmwghNxFCbnI/9g8ApgP4JiFkHSFkjbYkqCY9nLtgChMuIw69FObQUzhFeQu4SxZPq6lg/eBYBSvnTdZK2+cNjSPFuTR2GkXLRmshF7CMAeC261YCQMQ6i8OYr6s7h46FLbLQa+m5mMZCPzxYwitOny58LyvK5cmdbBM8tbNVLzokRqGrRPtw/HgNc6mt7+0PvK7KP/M6LpM7g5Sgqo+Hvx0ftqh+Q4oVxgpwH40W5cINJ4FCh8ZOQdTjNU13MQCIRsQLQCldDWB16LU7fH+/D8D79L8+CK8+hBeHzl7PrJZLuLqfhhz7+pmj5W2rFmDUvfG6N8hxKIZKFqZ0tPgs9GQzilsMYQtGxylaElhQALzsPtXHf6xiR7ISdeUAEOBbqxa6ohA+fP232wP/q9wTFk0iblRCfHOjnvW5k8B9GZecNk3LKuWLYFgZ6mzxl7pNnb/4J+cGXldNLNpygCn0sFGguoMMF7MKj6ETv73rCKMWX3PWKWxMDUXKnaEiysW/i0yaF5agnHHajNeGyhTl2+FJ7dXEoiyiCEp2VIkAegr52Ii/hAF7Tfc8hooWKGXe94Kn0JOPGy1ZWBFqbgHoJxaJLA/dXdNY2Y5Y+TqZcNwp6qdcamkJyBdaDpUxRJE6VVnY7/G20nns9Kyedr2ys3EWumbcNFCNhPLGgNoc51TFbEENcSD5WlY59Og90c305GUlLlo0zZVB/XjuFI0NW4TaNd18YDBioet0b/KjoRT6UClaZyJN5let8OLQ88F4eJ3kaJ7deOr0ztTUEW/f1tOuZ6H/buvhCG8NwDOrVSZKnELnV0I53rgS5dCriUEqxwucou7vNEp0VsgnoLK9rtiOMKkI8CkhfVFqwnDJAiFAZ0s+lYVei1O0IiiJzMZQuw782eCRVP7jgeR7wt8Wc+h6lMtxtxQCjwjLEfVsVRnlQhSNuJJlo/f4mJhDnzAWekChZ0i51MChD/hC3dKmqvNCZbMnt3s3XCVssSWf8xzLfugU5xJl8gH6ET9jFTuyJdWhXL7z2C4A0ap+OjL4MbWzFVcsn4UPX7lMeYxwTXg/dBe4emGoaGFSWwG5HOt5q2rNhesUcehY6GXe91cQjqpyHQbGKuhszUeOV/U1ORLKRdcAvPu5fQAQ4NBVncNb3RIS4Xh8v2xJohweZM+4v+ok4FI/zW6hcw69K9BJPDsLvTXc4EJDDl5WkycrAPpRLofcdnyLZ3RVJ3vCGKNlCyXLERb511GkQ6VKDOXCfitHApTtQINoQK+rOX9mF0yrVuarJbGIFWEq4MzZ3cpjVJxoLL0nX0aUy3DJ8nxNeY1npEoThO+tuhPvJZd3ji4Kasrw8FApNusWUI9yqUf1Sz7EVM9CVz+eG1eLZnQJ5FCb4zzHIRoCeuISi8YNIyUL7S25SMuzca+HbrMwQT5h0lhhz7x0DPkccZtMp7MoP//rzQAYD8955yQLne9yauW/H9t+VLyl1Sx9O1aJKnQtLt+y8ZozTwnIUguHPlSseEWYANUwUEdYC90vS5qmCrVguGj5fE3q14JTDOJdk9oYo2789Zwp0aYOSXJUbAe/XL8/0BgjKINKHDr/vhjKRUMTjpRsXLJ4WuB41Wk1MFZBIUe8hdUP1YX+uFu+eUpncIHLaxQq86OhFPpQyYpEE+hyYvXAaNkOZW7pK+Su1gI63YeGW5S6NyhHCM6c3Y2utoJXvCdJAfEwP16Z0A/dqKF5oQeWyaQ+RsV2ULFptN+iIlcKMAtGVMSJyZB4eACOwxoTsyJMUBrj6V3HpO9naaFzapJVfFQT4Iv3szSR9tb0lMtP1/aiS0CZMONLfiy/npcsjoaBqt5Xr8CYYI3Na4b7HR0pBSlejUbVA2MVTO5oES4sqhVa+U5+SkdU752oxKJxw7C7HfYjrVO0WLFxaLCY/EEBth8eDqRGpwmfLNsOTpvJtmJpwh4dh2K4bOF1y2cBgC/KJdkCAuJ4PTU5eGTJ8jnRSBlA3TrmRcLC013Hydx7fCxyzmmjhkbKLBnH72ROGuOgO4de64a1hZFVNvNQycKk9mpGteozUsgRXHjq1EjkkU5mY44AZwrmhopjlddj/9OLFkTeU3eKupSLSJFqUi47+kYCnycazsje42NC6ojLASTv3A65z8iUSEz+CUosGk/4rQ6OtJl4//fu53HJ5x+ClaKo8FjZDnjgdTv9AMz5FE4V17k/w2UWstjdHiyDkES5VGIcVoC6c/epncyKEvXQ1CkfwEPCTj8lWOtCVQ6+sJzSLQ5v053vz+3pB8BKIqiOwecPb64RRlZO0eFixdvqq/K+lFJYDsUrT4/psKP43YPFCi5YOEUwRvLCwinBGZOi1zOnuMhW49DFlIvOrSjkCKb7ZMlrHP/cnuOS6Cf2O2ksniAWXhh0dxre9+ofcuIwXIwq9DS1GQDg588y7/Wj249oH/vS0dFAA4A0D+2eY6O+8rv6YY88a7angz+0aotKXEiZf4yk8+Ac6cuXxG+LVay53uOslMO8qdFWY2wM+fGfuHsjgOo14EibcPb8vgEAwCtPn6FM+3jXU1TzGn5LbHzhN35UAwc4HScOs1PjnksWS9sXOzWT7wmPn58k5NAVKRfXRgvXdGdj6M0LQqp19nWPdxyKZbNFu1h/HHrybuOU7raI8UdIulIhDaXQ2TYymlWYpqWUN2aodocKypaNBVOrSRO6PKnjUOzrH/OiZdKEPfKSrT2uhe6NkTBIWdDDk0PVauB9UEWx7DqVEp/YwayPOZODCl1VIf9qAysAxguFRWTQnBclN2544bROX011xR1PjFM0q0bmQz6nqGrgAC8mJQpHBfRiyGNrwydZ6G78fJdwbrHfiZmiXpRL9D3d+G0rlHav6rMrWTZGyjbOOCVaaRFQX+hfPDSMZW7ElR9pWukBDabQXzoyEvEYp3GKDvnqdYyV9RQ6T7ef77MqdRUI5wkvWzYzcLzOafDU7p6OUOXJhOPkyQ5q55HUFFllDP6ZtkIu4qBVtV4udqMPLl8WLLWclrceGKugp72AfI4oR+tYkh2PH6qLdcmy8cALhzw6KQ1sh7rlCKq7N5X7wTM04xZqlev5q/VskY0LaU1aWHj8vNCRqDk/ZZSLygLnOKyuup9aVA29fMzd+Yt2GnwcIHmO7h8YE77OFqZkOSLH6R9yYkApxVjF9uJkOdLEofOkHqBKH6hisFgBpaHCQZrbap7hGe6BqbMw7etndIW3rYbaGPzc5U5R+XfzxtaiUD0dimGkbEfoFiYHlAZxKEVLnuBtq8K9Uavv62BgrILJnXyBZK8lWXPcQo/jSonOBQFw7/MH8ZffX4N/vneL2gEC7HGrknKlpeoU5c9CR0xTYxUl+M2HWS2cVwh5+GQ5hoqWt+sUHQ8k7/5kceg6i33F/SL/vVXd7fBG469eJnaWe1NcMlbZckApcL5bjNCPpq/lwlO8V8ydHHg9TRy6v16HrkLn3UtETVuV41NdpXqmt5XSpwi43AunMepHPa6VK/QatrQSC12HYhgtWWJrUNE6HilZKORyEWsubT30rYeGPapAddfFH/o4C123tDKnxL772EsefaELfhyvSqrKt8rqd6s6RQeLFq67YD7mTYku1KocethPxqFaJ9+bn0ILXW2MwDgRyiXxUO/YaaHyBWE5ZGLw+zhDEGKsG0/vHad9xAkCP7lwgL1OKi7H7qNVznVMU6HvdxcDP6+lstr6ccxNFqimE/Pj1eV41o3I4AqIKCpSXutFHocu/24eSSN3rMrHANii1NkisAbB5ZAPQgGcJ7Be0jpFd/QNR0qvJnLoVvy1APznoiaD//t+vUGtSUgYPFKER0Cp8q0yykWlfVzJslG2HCye0Sl8nyhkm/I8ABFUF9lqGz2BQtdIfOPzPGC85dTyRcKN5MNQoY9k/oim59DjTi7N1sPvCPU3mlDBEZe7ntYVDVtUlYJb6HxxSpOY9PSuo8iR4KRVsaJePDQMAF4nGz/UOUo2WcUWOpTGABifH84S9cuRNELFdqSOWZ2dW8V2ULYcLJvVHRgjKaq1YjvIEfG1ANRD7TjKPmJ0uJTOQvd3qQd0KBf2XAjvCZLnJ3+u4mKvVZ5VvzM3IoNmYlFcHLr/MzLYdtRCJ4oJPaVQie04OWQj8YY+PTEO5hQR142j0HnR+ykd0T6YujuPTfsHQQgwY1KrNuXywgEW2sapDiYD+626BeI1aTpbg52XdNLDB8csXLBwauA1FSuKL2DibXXywkIpxSd/sQmAnENXttCFCp2PIR+kbIlrqKTZ8XAn8ypeJlVxey8qLuaHroVesapPqai5sArCVUlVczWkLdNIciu9QV99IhFUHKvMQo9fEAD11H9REymd0sxVDr06kGpdnOf2HAch4vIDgJoRxx2iMwUx+XmNjFU/GkahcyfDotB2TicVl4N3Be9sLWhHuTzuhtr5vfi6scajZQsdLXlv5U+TCONQigsXBRW6Smr1ULGCZbO6Y9tzAfIHZp3b7QgITnQOVeoHYApdZA3mFC9oWdCP1H+8Xn16vvMKcegJg4hq0QTg7TbUhPE3Vh6t6IfUAlVL2W+hqzwiO/rY7i1ukU26pwOhUNowVJ7V0bLllcQIgyvjJMojKcoFUJsbvPJke6i9YXLHJIpHtx0R7hA4VGjaZ3f3A6g2DAkeny71X6lj0Xjg4EARUzpbMH9qSKErboH8GC3bOHN2NxxKPWerChyHovfYWKRSoWqYHcdwyQ6Wr9XkfCmlKNtOtJodkq2H46PVSI4wVCa7/yERxwqr00eiWuhsDLhyJEeYiBR6Gg59pBys5Km6SBfL4rrwHKoROxxcoXe3FVDU3D1yhMtMq9KSPGInnHkLqFEu3MIX1tp3R0kOA6WJEUO1RLnozI0Dbtp9W8B4S3ZGcuPg9WfPjv2MCk27vrcfhMQHMOir8wax0I8Ol/CDJ3d7TkQ/0tRyGS0zhVrI5bwkGRUMlSyUbQdXnzc3JIP7h6IcY2ULHa3+VV+PQ7fd+Ngw3UBIsiU4MFrB1BiFrjLZSz5KQGblqzwwo2UxXaHazaUUS7noc+gjERpM7Z5sOTgkVeg6nWkAxqETwuKXdelAjqFiBR0teW8HpRqHzvvExu2akkaIa45RHUMlUSu+FLF2FJbEQqcKjz0PDZ4R8pclXcrdLvd9hVtnSQSVZ22sbOOMU6JJRXyApk39X7v7OAAIa3ir8oN+bNo/gK62AlryxLNKVMCdRtF6MnoP7VglGN2hGyUTV4+FKFAuoxXLU1phqEQRcMV37vzovQDUFVjFdly6QlZaNPnhl8XT60z4Ysi3oBpp091ekIYX6vpHihUbba5S1XXYc4QjRVT9TMWYBZaNUVt0CaC4KDjRDvf+4wGFFnSyKBf3JZVIlce3M3r1bN9cV9nt8IJa/sCJMKrPfPw4luPgfEFNHC5H01Zb5Nvhv7vijMh7KlvBMDpa8iyGOa9noY+U3C1lTFNj9YfWQbvPCtKNkvn9i4cBRNP3Vegn245/YFQmGVcy/3zdOeIxFJUxr7Qo2r6qKuR4Dh3u8eoToxTpQqV2T4qWg7OEVSfhjsNlUZODNdhoQUdLXjuk1j+GP1JEOfVf4uBVSduXRZdwOZLj+uMpF1VHtWynoNM8pWhV6a+qDMmLIw9LPmNWjHUNtWdeRj+lrWHVEAp926F4Z41Or0OOsYqN5XN6UMjpWehcCYW5Y92HdqxiBxwtulEyX32IZeO9amkw5Z0oyGBTGtuMIacw2UuxHW2CYyRNNs4VhystAuqLgsiPwI7Xt9DDccOqMpRcizoO1brXqv4Vlvre0ZLeQt+4fyCohBQt9M0HhiKNzzlUdpHcnxtrMCjwz5btoEUUngJfKGkNTlGd5idly44EEBCS/P0b9w0CiObMBAVhv2TXoyyhn9iz3qQWOi++JCypmSKxaKRso7ONFeCvaARzDhbFyU06kR1ANP5aN0pmqFTBkpldkaI9Kqu27VDJltiVQzJGqSJPmNBxaALy5CTZCBXbgR1jzaWph14OKXTV7X3ZEtM+HKoZvBzDbqZkPkfw+I6jWLtb3kBDhLZCLuDryBG1ZiGT2uJ5e5XzkClSPobseMehcKiMg1d7zjbtZwq1lvIWgHgHmFfYBbe6c1LkzOSIj3+pQtartqkt9NGKjSmdLUJFpMuhHx8po2w56GxhHLqlYaH/fmsfAGDR9GCPQBWqwo9iaGurW9yrb6gkrhGhyHMmtUtTcYrGWuiK1nFFoeqjShZd2YouyGmqLe5y+2Dyh1CVQy9ZjvzB1XR4j5RYUa01rt/o3Xc+rXagDwNjFZzj530VswpHy1ZM0xK1RVZWEoKPIZOjmoEsV+iJNdVdP48o3E+HjivbUYWustsp21SYiR0eJ0kO5iCOUf0pfIdAgyj0vqESLgwl0XDoRrnsdWtwT+1qQUHTQud0Q9jZUZ3saoKEuUqdBcF2WKilrF60DJZDYz+nYr384z0sqSg2pdn9nZgUZAc568AYChQWd2KGF1fAH0MuFSEAbumH68sn3ZN9/WPCc+CoBkCpzY2hkoWutoKnHEdS8Oi89RmHqjX3/L6B2JBDlfsqK1vL5EjgjAWJPOHjgWRqsm+ohDmT22PC/VzaRqW2eyVK6Sn5AWKoQD+Sdjy82Ughln5q4rDFEVl9h5xueBp7QJae0s0sdI3VYGCsgrPn9QiKQbHfqgpkrBxyiiryzoC8QYXKbsWRWOg61kucElNO27dk55G8teb5AyLON00c+uBYBTMmtXnfrTIGVyyyrlf8eVQVZahYiZ3rKhA1mFCNQ28t5BBnE6hQLo5nocfTBCoWevIOMl4GgMWBT48riqX4rFm2gzW7jwsaSySX35Va1lyOhOeEU5Jxz1lT9xQtWXbstpYl06ifGC9d29Gax+rnD2L74WHlNnRhy8cvA6C2YlJKMVSsoL0QtdBVGysD4hoRKpaYJQkL04mZTqwumJQRyMMEY+qGsDHij5fXdU8+Pgx2b/2hfu4YkmO401Lk2PXG0Ug6GypW0Ht8DO0tOW8X+DpJLLMIe4+xdHF/4ppKSYhihS0E4Wqm/jEA+U5DFv/NBpHPLVlZZkA9QuW5PcfRFRuaG5Q1DjzbdnaPuL2hTIS4khR+JBkM/DmPX9z0q4kCDaPQnVjvu+7Wg4eC+RtUlAQ8rAhHR8qxrbUAtYf2wz/dEDkfVasWkK/cKp5vh0oUuuJ5nBMTgw6op8yPujsl0YOnQmFx+kvq+NKYGUeGyxGaApAvTLzyppRy0XB488zEhdO6PAsvrpRsHI4Msz6tfkWUU3j4v/bbbQDEdUMARaeo+2aMgZ5oVXLDKpFySZifna2F2J2GciMYd36dF4oDV5GhYjux5xBG3DCWJGiAydHMFnolPpKAVR1TPzHeSWSaL+tUJcmgWLGx++ioMClHxyL86dpeAPCq+gF6VIeccklO3JBZ6Kq88eVnzIx9TzWEc91e5vQL9wP1yyGjsDjlIu6Mo77T4Hh617HAvVU5j37XMRu24vzQiYDixa2Wz+2J5U6TwC1Lf4kMlYefW/ZXrRSnq1fDL+PHcBKdovLj73pqDwBxcTBAnXIZLVtYGpNhqWq0VOdXiEN3z02mM8oxJSkC43grjngcWatIjhNmoRNCXk8I2UoI2U4IuUXw/pmEkCcIISVCyN/rCEAplVIuMXMnFmWb8Yu5HMFfXb4EQLVMpgyHBpn1dNrMeCeczl7B/+DoPPRlK/5Gs7AwOedLqazUK6RyeB1wJBdddbcxOMYUT7ifKBvDlVdyLg+8cAgAhHyz6raag7q7lhmTBCWRJUPwMrWzJksUuvcdyXL4S0Tz+6vj4wH8FQ/14tCPjZRx/sIpiZUOlZyiEqe77Pg+d3cRVwNF1fA5PlpBp6A8NKAedbShtx8AIguryryoWI4XuhiHJINhTNI9istxQpyihJA8gG8AuArAcgA3EEKWhz52DMAHAXxRV4CKzWJT4+Oe9Tj0NS8d9zjPue6DqGKh8+Sm02fGJ8KoPHtLZnZhxqTWwAKlU+61IokOSbKAVB44IP48kuKMuQxA8uL0wOZDwq42gB5NIQyz0zgeYJEltkOxfG51LBW/xu+2sDBWUY0hDp0sYJ7n0NNewMuXTAcgrlsvQ7UWelUxqzjLj4+WA7vWMFTuiSzlno0hX1gGRll+RVcMzaRCg/FeA3GUn+r8/NsfrgPg7yoWPF52PY+PlhM59KSFYcQrMxJfiuFEUS4XA9hOKd1JKS0D+CGAq/0foJQeppQ+A0C7Yn91pZLUmFCMLukbKmGsYuMU3suTb58UNPFB10JfPCPeQle9vpcsnh74X8cZWZZRLgkOYs9pFRvnC1cO8RhJxZfYGGqOq2Mj5ViFrvLgHh4qYnZPe0yBMPUdDwD8djMrpeDn0FV2GryqnmhOVMdhv1WMju2HmdHQ09GCT7/lbADiyocyhEvnAmoW+ljZjpS08EPlmiY5Rau7FfEY2w4PYYpkUcl7z2vsRzDkNgWJNRY0AhgARHreqtB5Ww4OJY6bZDCEi8WFcSITi+YB2Ov7v9d9TRuEkBsJIWsIIWv6+pj1w2tCx63ajDdWO7MvPbAVAHDVyjkAqh7kJIVu2Q4+8T8bAQCnCPhSnSgVntodOF7RagCiGY1+JFroSVEICREZfOGUW+jJE56fZ5xzVWXHc3iw5DlGI8d73xN/vB+c6vAna6nck97jo5gjoVt0ZfneYy8BYAsLDyHU7sZVstDekgss+DmiVte9PYH3BdR8CvEGg3xuvHhoWHq+Kosjn+O1UEcA484vWDglolCTZOgbYrRR3IISHifeQneDBmIs9BPZJFp099LQO6CUfotSuopSumrmTOZ4G02osaziwefYenAIna15vMlV6KpJBl+4b6v3tzjKRX3VHxa02NKhbLjyEfkUkra0nI+VOa2A+Em288iwe3z8d6g8dDxSJ67WhcqOp6M1j+kxERm6Fjq/plNFbQUlQ+RzRLjA+6FaO4RSipGyjVWnTvWUcV6TTgSAzQcGMakteF1V+k8mNerwdkKSYXgTmrgGFTK6gz+D4S5cwePVdwnxNdXZ76RdfWdrIUDBRWSIOZ6HRb/mTEEmd0AOeSTWaEluyEKhtrwIKgq9F8AC3//zAezX/yoxeHhbrOdbo3PH0ZEyTp3e5VEt/KYnKfSfPbtP+r6qhW077KENW+g6CmjzAVanQuQMJEQtTjg+sUiufN741T8EPieCSsRPUkagyo7HoRRTYlqd6SyQAHMkdrbmA1atigwlywnUypbJknRrueFy+bJqBFGO6PeN3HF4OLJzUYkEC5ejCEOloqjlsJT32LBDSYNm7huSlZxV2f0llR9Qiafn44gijZIs/D3uopaU+p+0c/MsdEk8/Yni0J8BsJQQspgQ0grgegD3aH9TDHgN8lguKadXJ+OsOf5wQbVIgle4DqpPv2WFWAZFDp07OsLKWMdxxvtMiqwHkkC5cEWal8S2AsnnoUK5yDPp5AuLStSPZccXGdPZMQHMGRlum6aaQKIeniYH574Du4QU7RVtSnFpqApnPoFvpZSVk2iroTfqsZEyfrl+P3oUslxF55RUxwVQo0usBFpRdY7H1WVPClvkss1K2LklyVHVezKnqPQrhEi8O5RSixDyAQD3AcgD+A6ldBMh5Cb3/TsIIbMBrAHQA8AhhHwIwHJK6WDS+DyjMC4MSTXKhVKKwbFKwMnEV+Ck43lD6He//FTh+6ocerg1mHe8Iq9Xsmx8+cEX0VrIxVYpVNqO1phY9OKheKdPlcaKP96SOHbZGOy3TAqHxte50PFJAMBze/qj8fAK10JUvCkiC5LHAYDhUjQ6Ja+ZYwGwXUM47T1pYagWXJMlSMkX2f98cjcAYIdbGVUEmRJLmhOAWhBDooXu/k66H5bjxDRBl+9ieQ31uETI6jhyOYYTKJc0ZcMBxZ6ilNLVAFaHXrvD9/dBMCpGG0l9ClW7sYyUbZRtJ8Db8rmTVHGR976M6+CtuurzmxTm0FXDFn+z8SAAcYVBgGeKxh//4iHOgcech2KdC56NKBxDcWsOyHhO96FJeHDbCrXdD4Dd222Hh7E0lL6vkt8gS3jjUKVcBnl0SqiZgq5CHxM03iaESENzeUSFrGpk0iLLKZMbLz1NYYzoKOWELFF2fPJ9TeLQVWvrxJWZzifIUKrwDGZ5uGmSS2JXH6/+GW+0pOHQM28S7YXvtMREuUBta83TtP1bsbyihf7tP+ySvq9q2Q7FWOhJ0SXV72Gfu/Z8cRBRUrLBrzcw10ZcmJ1qBMA7XybeqfhllE02L9tVUkkOkN9Xm9ZuhQHAEztYm7G3nBvsE6vizFSx0FXptGFBuKGKM9MPx6EoWU6EC0+iXLa6Oy5pMkzC9eAZyP/3DWfFDiHjwLlRJZNBh3KJr8mu9qxZMUXskmRQ2e345YhN4gt9TnR8miiXzBX6WCLlora15sr0dF+NZM9CT7PU+aCaCLPDjTOOcuhq38Mn+/teFWMFEflEpZQlU120aFrM4XJlPKWzBW85dy4uF9Vi9yHpnnhFmJIsdMkFdWRFxhTvB8D6ywLAW2Iafydy6PmExB/FRVK0e9OlXEYr4gCCpAYXPBnnosXiecHHABB7UW2HKpeMFXLonl8l2UJXcYrG0nHubymdJ8mo5kXh4nb1pYoNQsTF80RyxJ1LxXZw6vRO8ZtgEXs7JfRWHDKv5ZIctqi2UlWz8PyUCzs9W7flUQiqqfsf+dkGANHJoJqMU+3yI6EapE5RGhsjzI6HVA5bUks9LIfcikqIclHgwGut687BF/qwE0sl2atk2YlcqaosogxP1cYUHPdvYpRcmIdOGueomyAlzRRNuB6WHV+W2RtDEvJXpVzix9CJQ6+l6xGnp0Tn85WHWBGzX24QB/IV3Q5WcZY1R/X9GC6+YgcqsobBk9B0kb1Cdy2XuJNTzRTlD66/RGo1sSj+uKMuX3zTZUtiP+PdGsVnL+zoUI2N5YowbZ9By6GxNAeTQ04PyKxiP5IyE72FqYYiYUwWuRwqynBgrIJTutsiRb6qVn5MNINDUbEVrFL3d9JiL6LjdC10Th9dGaqFwto0xo/jGTsxYaBA8vWwHUdqLDA54sfgc1t2Pb3yuQpOUVlzZUA19DEqC6fGOBUcRrFix3bzCsrBfsdNi2LFSaRt0iBzhb710BBaCzlpiJrKlPeKFrVHS6RaEk3KrZcVgjDB8Diqz164NZYqz1qxkierVKHbTkLaPvsdp3xsqqbQ2SKbTLkkJTglFYKK7+aiyGEB6B+tCBOckpyZ3KJUjWZIurdChZ7Ti0Nfu/s4zp7XE8lSzCcssMWEPrFAsgKK45yDYyRz6CpOUXnynBP4bBhqeRJ8fkbf+8QbmY9g1aliekrFUQ4k73jGFBcGXWSu0NsK+dioDkCfQ/dbIVwxyixjK4HmAJKtF45Tuttww8ULYj3giZRLkoVOkieq0gMTM8scqqYsk5x5ieehQHc4MVEIgH/Hkzwvfv9iH6Z0RKmGpK05d34lWeiqlMvDL/ahvSUXWOR04tBHyxZ2HhkR7mTX9/Z7tYhEKFXsRJogqVGHZScv9rKFWqVcrArl4kjoEtUxbM/giN7b5W4DkDhneNHStdDFcoyWoyVCRFCZ44Hv1fr0CcBgsZLYUEFl0j+75zgICVohKhY6H7vW7EggvqFwrroXlYK3bZMl5Mg7wohjaz05Eiwgx6FKDtwkyiXZKcp+S7NeKUXcs6/cxKBiY6xii7Nu3d9x5yGrqSMeR055rN/b71nKHDqUCy8n/Nqzoh2OeLGouIdfiSZQstAT1IXEV5TU0AFQ47+Td3/J/ioZh87Fi7svoj6kQngLi/jtkZK8WFpYVlVkr9Bj2r5xqGZM8dAtvxXCb1gtThbAV8FNITEorqhWkhwAsPkAezBbpA2a5Ry6zOkki8iwHda09rk9/VIZAR4jK6d+AEkkgsq22KaxvStVlChQbXZyqaBhR5JlXe2YVJsiBIB71okdbLmcPH7cj0e3sXO5/qIFsZ8px/A34ablQlkSdmZqfTTZb1likczgyCcYHIB6YpE8eMClbYRhi67fLXbnpmahV6s+isfZdWQEXZLaOn90Oste181TyFyhD4xF07L9UI3HHClbOG/BlMBr/KbLEouSCloBanHTvHmvuHWcGgf/7J7jAILJJwE5cgn8YkIkguyBOzrCnMNP7DwqFxI8YzX+/UpCmrdKGzvWSk/8nuqO6bBbGe+S06J8aNLWnFvo6nHo8cK8dISFn33m6mBpiXyCM9MPrhCnSmqh1OKAS1okdcogiIbwKBcpj8+DGNJZ14C8ngxHiXcrEsjilfCNC1u01JyZsjBQbhjyEEkRXuWWd9ANRc9coQ8WLWGbMg7VusAbegcipShVms7y96QWuoIC2t/POExxyzT2O4mDp2C9UGOTDRIKlcUVHKrKEX89+KLz3j9aLJWRjZO060nqHSnfa/x2yyEcGCjGNw5XDCPlSllUczxpjLKk0UhgHPd3nChbDw55iWtvv3hh4L28RqZo2XZim1Vf48bYx1mVKpSLipM4WaGz36K5sctd1GQGB8nFH8+haqHLnrWndh0DIO7n6il0iW8lcdcGea4FjwJbOkvcRg+oXsumolxePDSEvqGS1DlAkMyVliwbA2OViMPIs9BVakNItpwqFnr/KIuWkXXY+fEze6Pv+eA4FOdLy4smZGg6jpRykSlSnt25bHZ8h3v/OCrlc2N9Ae7vuDHW7x0AAPyphF5ICMkHUKVN4hRRTjJGSSEyBEgup7DebXX20defGVngdLpxlS0ayz+fPY/5oOIWh6IlL8zFZQFkGZLxbSI5ZI5Vfs9l/VlVHMxJHLpK0to96xkFdrEg0UoWOrnl4CDW7j6O9Xv74wd3ITPixirxDdCrxyfvVoTHaX26ztjqOnPiMhuBZGsQqHLPFy4KKsNqBxSJha7BoasU/58sCJHjPtn1vQPxA0AhNCwhhNNOOF5GM1R7mSqEZCk6RWXROkD8g1txnbtcUYmgsnNLcmzK6DxdCz1uHN71SFT4TcdCr0gs5KT47XKMb8ePpPDLjfsGlTNFRZdizF0g4xpTAOrhrEBtIbGlio2e9oKw3j437ESW8WPbGR2pknlefdai723cx/QAry8vgkpMvgiZKnReQjLpwU1S6LyLSLh4fl5hlVPh0PlWULbF//PvPgNAnI032+1686Zz5sQez+VM4vKTytbGORLZ8fEWUEWhGp6qHNVM0XgrikjGSHTuIjnJCgDWuZZ+fH14WREmtbDFJEV47/MHkCPiTGjmFJUO72GsYsfWQUmiCcYUIjOSmivHnUPwM/Fj6FilNaX+KyQGrN19HKtijEhZxUee7fuj978sfvCqJEwMwcV4x7efAlBtliGCCl0sQqYKfbgkL/IOuFZUQvLFMdehNz/UH1CJcuFhi1ILnUHl2sbVZzilu016nkB8SU+/HDIZbEceiSCL35b1Mg0jl3BPqpmikqxVyOtcyI7lMiTdDt6AV+aTiJOBl4uN7yiDwNiiB9eyHazvHcC8GL9InqhbYE/vOiZRtvHW3OGhItbv7ccf3IifOFTnuPhLhooWls2O53yBaqVDkRIaLVlozccnEAJqIYdeca6Y6ZFkoa/dfRyWQz12IAxZ28rnewfQkic4Z/6UWPnCcoik4BTz12+4IPZ4ItkpSL9X69N1xha3O09cYS6gqoRk+O0W1gQ43A2lmlgkUegJfC+gFskAAEtmdsUqj5Z8zku48aNk2V4EQ7KFnpQpmnw8IHfUqGTB5RPC7Xa7lkdLTPlbLkvc9bTsZAs9qVAZwKzCZRLHE+PhxWPwCoVCn0hQDADihXbjfja/37BSvDPToVy62wqRsswcsvsq29b7ITNsLduB5dBEx6pMIT+3t189MUmBIk3KIo4b4kk3iiupgbnovvR0tHj0axLiFnpKKUbLFj7w6tOF9CyHF8KpWYYqU4XOQ5jkyQZq246Olnyk6xG/KDIL3av2KO23yH4nPXuyfokteSIMn1z2id/gXXc+7cmZtLAkZYomXUtAfD0rCsWTOJLi0LlVK8rQDI4hfs9yHGnGK+CeS8L9GKs4CX0043c8bYUczl0wRWpRsjHiF/utB5lCf/2K2ZH3AHY/VS0wy6GRmu4cXjKMYKyjw4zD/7d3xFuDXBZAfD2KiiVjOUT3tautIEzwEskgT/1PpkiB+J0Gl+Hz164Uvi+jOiyHYt5UeXNojriFfqRsw6HiFpNBOdjvprLQB0YrOG2muHY3hwqHvvvoKM5dEOXhVXgoXpwrrgUelwGQxfmyRUHWL7GQz8VmrPLY7/7RipQDT7JKk2q5yIpzqZQ35UiMQ/ccrHJZ4s6lf7QSW9jLOx7J86JYTuqjGR8GWrEdzJsibzPGxmC/Rbf2Sw+8CABYNF08x/M59Th0WdigjHI57kZfnRvK0QhD5twtKjZ1kGV6jpasSI2jiAwqTlHeZjExU1R8PC8lzLuUhSGjaVWL1/nlCJ/KkEKhNP/xTeUUHUjIEgXc4lwJ57Tl4BBaBZMtKbFo99ERfOqXLwCQ8/gcMgUEAKfGPLgAo3QqEg8YdxBXJNWaZGF2ALfQk7e1sq7satUW5Q9dxXZAiEKylmCIkmXj3o0HPaswDkkZxJRSbNjXL7fQAew6MircXiftdqpjxC+SOUKwfE5PbDJQEnXFYTvUrUcuPhdZNNcAj75KfM7Yb3GWp+sTUa1rI3hve9+w1GhiMhAFh7ucIq3aI+IxBtyG4WkWxyQflx9xi9OQoNGJVI5mstB3HxtRmmjyUq3swT99ZnT1T7LQL/vCw97fMh4/KS2at2wTVfXjaMnnPK5chL3HWMelsyScbVJi0eGhUqKFHbfjSQoHi44R/37ZZsowqRiUSI5vP7oLQDXcL60M9248iGJF7iQeKdt4cPMh3O7WwPaDZd2qhHCy36L7cmCgKIx15mBzIvmBTcpalUW5eFmRiYlF8dRRUtQSRxyl1z9aRv9oRcnaTKKhksKM40KMDw8W0TdUwlDRkipTmVM0yccVkCPGJ8GrwsrCNwG1kGsRMu1YdHS4jEXT5QInrdjH3Qd/8YzoFkolygVgDieZ9SGLDgGqTZVlCVKFPJHKce/GAwDi28cBzPqgMdnC/DrI0omBeEWYtJX1I4lDr9gO2hQ4cJEce48pOvEgd1L/11N7AAAfuyq+ZRrHk4JyByq1S4B4/wqn8mT3oyVPpDsyjiSFLrMqS5aNfI4k3leZc9dSCBwA4uv+80qQr1sRLSwWRmLht4REwLhF5eLPPwSAtXeU5RbkJDpDpYQwR1z9J2789SRZ6M0YtliRpDNzJIWn8cbIIn5PdZVbnMDje5M95n1+82X+gJZcTvrwcgtdZtHJeGO+lVsVSq4SDCKthqcyYWXcM8CUiKxmR3WM6Os/WsOyabd89vXS45MoF9uhWDS9E4skCySHMORQIRYegC/qITjGiBuSK3OUF/K52IJafiRlvMos9LKlVr9b5ieqNv1Wi2UPz699x9ncnixxklfHSE79JyQ+zDgpnr6cEBIr29WnstDDOwU3ZyapdK5KoTIRMlPolLJQuRmCbC0/CBKyvtzJLoqRVUksAtS5wTg5ijG9Hv0oxES5cFRsB7N72hOcmuJFpWw5+M0mZuF3JPCUcck0KjVtqmPI49BHy3ZiEopoYTk0WPRkU6EIZItK2XaUIxLEYZyOmoPYSzqLfj8gb5DRmpcv8hy8TWNcdT6ZY/bxHUeVd11AvCIDkhf7JGpyrpKTOcHhnlCALsmxatly30hVZ4iP1XaKhsd3j0+imr1aLs3iFOXWRJJzgCUWxZ8UrzEtcn6pWuhJW+ukVPUx94GTOeAK+ZxXhVCEYsWWcvBMDvF29Ou/247Pr94CQK0hg9jho6HQc/IIk5FScvF+na5D4uPjd0xHh0tYu/u41AnthzgpKNnBDMRztkntBNl78kWeg0eqxCU5xRkcjkPxwoFBr1mHDLI4dL7opG1woVNWIik2v5hQCjgpIs1KqHcktdCpvoUeHkc1+ECWsSo9TuvTdQQXVMXbK1uxeRy5qJML95rft+kgDoU6uvgf4qSJltSLk1tQsqavLTkidYqqtKRikSFRKXimLAC0SpJ52BhRCuvwUBG/3sAsfPUol+jrxYqNbz68HfdtOqSQJh5vRd36x2cnyiALe/yPJ3YDYNmVKhCdi0osPBBfhEmFwmpRtNB5LHncWHGGyzF3IfjIlcsSv0OW8crHVTV8wtdTJws5iXIpVuzEyCVAnjAmu6+yyLikaqZ+ePMivNArhgfnE1iB2O/V+nQdwa1EWS10IDlU798e3g4AaG+N5xe3HBzCW+94PPC6/zqpTrS4Lf6ze44jnyPSJJQkyuXRbUcSmxAw+in6ur887PO9g/IxBA/MX3xvDe7dyDrKq3Looon22PYj+JffbAWQbFmIFDL/V8V6l9ZhcWm4D752aeI47HvDcrAG0Umx8EC8Ervf7TAkt9BzODJc9hyoceDW/ty47MYYDp2H087sltOagD8OPf77pTkSiOfQ+U5JpdOPigGnYqHHUYJ7j41J76tHdcT4VVQtdEB8LTwLPWFxrO4UFL/ORQNQLiphi+KzemrnUezoY2nmcVmJfEvMnY7h7weSLQ9APtGKFTtRgRViUv8D35FwN+LS5f0P7BtWirMSA2OEhuBp7lzOJMRZ1/56IW+/ZGHk/fAY4SH4mKpxBHGXfKzM6Ku/u+IMpZGidAm3SFUeD7Fl+1U3FFI2xuEhtmv8q/98VvoNYwk+mmqaeGiXoED7cFRDTNM7zOMSi7zCbwm7RzaG3CodLdtSP5GMOlIB39WLaEk7RRx6GElx9OHjmyZsUYdyiXPU/Om3nvT+jw1Firlu/gulosRkzlnLoXjV0hnS4xnlIr85118kV4KEiC0PLtZpM7qkRfO9MULncdqMLq8v5SSFBKs4Lp9TA9tvvUrhmkbH4P+q0OvsM+LrOZjQBSuMqEWZ3F2HI+n5li3SfNu9r38s9jOP7ziC//Oj9QDifTRxlEtS7fCAnBI/kbpTVDyGXiVPuX/mmZeO4TRJ5JJKtum7X7FIKkMhlxOHLaZxioaGUQ0Pbrpqi1zQpAzNuEzRB1446P39noQbJIJ/TOWtYMx7Sa3fADf1P4EvnS4pHQC410LwOt9t/Oj9L5ceD4gfGL+ikCVYVccQ008ly8aZs7uVrXz/GHuPjeJL97NUeZUCSHGUy5M7j+J/1u33Go6oIPzsrt19HIC4RVkYIprBf168k5UInBriv0X490d2en/H+VjiaB9LkfsG5DX/uTM/ObFIPIZqOz+A0Ufhx6R/tIzPr96Miu3Atqm05LBolxD2U/D0+3gZ6hC26P6Oo59kDXX87zeNhc7PUxbWBcTHG/s7sMibU1TH8cNPuShFXEi2ghU72YHWkieRKJewUky0rgXHAICt0ICXQ+TQ9Dswa7GiVPpO+sfwD/GB/3q22gBE6XaIZbje3bUNunH5KghfU14sTZbpWpWDj1F9zR9VIkumuW8T49mPDIsXn4/+dAN+t7XP+z/O0RznQNtxmOVoJHHfgN+6FlMNQLIjL0mJqRlOURlu+81W/PfTe3BKdxuGShb+15mnSI6PWsb3bToY+Mw1582TypAn4kgbnSgXDt56j+NOtx1hUtG3E5pYRAh5PSFkKyFkOyHkFsH7hBDyVff9DYQQeWk3n6DJYXbRk/ruY7vw282Hvf9lrco44vhaoNqcOUmOOBNdpYZKIRe10P1z5ux5PYnOqzgev1ojWsWJFx0jqcaGSA4R9VOykhspVOUI3gMefgqocegyp6gu4sb56dpeBTmiVikPY53Z3aZF/YTBk6wA4NoL5sUutnGUy80/YVSNUvZtjJX/4AuH8N7vrQl8T+wQsU7R2igXfvznfr0ZAKT5BaLdSnjHl1TjPhcTOpnUEcyPXjeZ6mu/3R54fbikZmjweXXv8wcTPhk6LukDhJA8gG8AuArAcgA3EEKWhz52FYCl7s+NAP4taVxVC110gz/9yxfwA7cBwfteuRhnJFi2IvidHuFVVARZlqZlO4mTVRTlwsf7+9edgV/9zauSZYjZJfDXVC30sAWUmASkKEdJ00L3PzP+xUjJMhbw+IO+rfS/XHeOkhxAvBUko0KqcrDf/mvKnZg3KzplhTL5Tu6z15yNL73tvNjP+qNcth8ewju+/SRe8s1plXlRLTIWvBYf+O+qwzYp69b7GgGHniPp6wSF5V91qiSbWsChq85J//cJC7YldASrJ7g64XpOFSqm2cUAtlNKdwIAIeSHAK4G8ILvM1cD+D5ls/pJQsgUQsgcSumBuEGPDJcwEyotvtgNvuJLvxe+f3XC9slPp/jH8N+wy5fNlI7BxgF+srYXD/u2wBy9x8dwoWSSAcw6GS5bARm8qA7FJJscIdh9dDRyLXjiieoD8+vnD3g8MQAcHIjneUXI5wie3XM8Ise2w8N4rWQ77AchwP0vHMQVX+r3jvXeS3F8eIy3Keza/MeJ5td//2VyqzF+675w31Z8y+W7uUWZlFfgR/j7/QrpnQkRQ5xy+fjPn8d+915e/sWHvfdV5OBT58M/2RBY4P07p8RkMXeQD/80OMaR4ZKSXwVg1/O+TQcDjZj99xUApnbJ+pIyGW5/aBv+4/GXAAQXehXkcwS/WLcvUuPn0GARiqcR8DfE6S4Z4vTWtRfMl3+vwtjzAPjb1fcCuEThM/MABBQ6IeRGMAsePXNPwwdfc3rijb5yxSzs7BsOTPBDg0V0t7fg0jNmJNZY/tgbzsQ//GITVsztibSHm9LZgkI+h89dk5zI8tevOd1r7hrG0lmTcN2F8oXlTefMwb7+sYh1vHzuZFyxPLloEQD8yar5sQkTi6Z3KT24N156Gp55KZhws3TWJDy58xg+r5DQAwDvfNmpwqzWpbMm4a0XyiecXw7/AzNjUhue2HkU15w3Fy9fMl37eP79q58/iH+OaV4QxmevXoFP/mITLj1jpteujo/zqqUzhU2Ew5gzuQN/9vJT0ReKJb9g4VS8IuE8fvU3r8SbvvYHnL9wCuZMjqbFDxYtXH7GzMQFf9nsblx/0QIMFis4l7JKk29YORur3e26SlGsFfMm422r5kcogSmdLdi0fxAfff2ZiWMsn9MjHGPprElYMXdy4vEA8P5LT/P6A3CcOr0TD7oU6/86a5Z0YZne1Yr3/tFiHBwMRg6tfv4g3rhyDv7odHk0GgDcdNkSIQ17xqzuRIXKcbHbs/TCU6diVk91Hm07PIxrz5frCgA4c3YPZkxqRT5HAjpuxqSEwAlZPQwAIIT8CYArKaXvc/9/F4CLKaV/4/vMrwH8E6X0D+7/DwH4CKV0bdy4q1atomvWrEk8MQMDAwODKgghaymlq0TvqWwgegH496/zAexP8RkDAwMDgxMIFYX+DIClhJDFhJBWANcDuCf0mXsAvNuNdnkZgAEZf25gYGBgUH8kcuiUUosQ8gEA9wHIA/gOpXQTIeQm9/07AKwG8AYA2wGMAvjzEyeygYGBgYEISgHIlNLVYErb/9odvr8pgL+ur2gGBgYGBjrItGORgYGBgUH9YBS6gYGBwQSBUegGBgYGEwRGoRsYGBhMECQmFp2wLyZkCMBW30uTAYhTMdUxUcaYAeBI4qdOrAyNMkY9ZJgo17MRZAAa43o2yrXI4nouo5SKC1hRSjP5AbAm9P+36jDmhBgjfG2a9TzqdC3qIcOEuJ6NIEOjXM8Guhbjfj1ln28kyuWXZoy6olHOo9YxGuFaAhPnWkyU69ko16JRrieAbCmXNTSmHsHJDnNt6gtzPesLcz3rC93rKft8lhb6tzL87kaHuTb1hbme9YW5nvWF7vWM/XxmFrqBgYGBQX3RSBy6gYGBgUENMArdwMDAYILAKPRxACFkASHkd4SQzYSQTYSQv3Vfn0YIeYAQss39PdV3zMfcpttbCSFX+l5/2H1tnfuj1vNtAkH3ehJCprufHyaEfD001oWEkOfda/1VotoPcAKhztfTzE/963kFIWStOw/XEkJe4xtLb37WGkNpfpTiRucAuMD9uxvAi2ANt/8FwC3u67cAuM39ezmA9QDaACwGsANA3n3vYQCrsj6nJrueXQBeCeAmAF8PjfU0gJeDtTK9F8BVWZ9fk19PMz/1r+f5AOa6f58NYJ9vLK35aSz0cQCl9ACl9Fn37yEAm8F6rl4N4D/cj/0HgGvcv68G8ENKaYlSuguszvzF4yp0A0P3elJKRyhrjxjohk0ImQOgh1L6BGVPz/dRvQcnDep1PQ0YUlzP5yilvMPbJgDthJC2NPPTKPRxBiFkEdiK/BSAWdTt7OT+5tvTuKbbHN91t7OfPBkpAj8Ur2cc5oFdW47wdT7pUOP15DDz00WK63kdgOcopSWkmJ9GoY8jCCGTAPwMwIcopYOyjwpe4/Gl76CUrgTwKvfnXfWVsnmgcT1jhxC8dtLG8dbhegJmfnrQvZ6EkBUAbgPwfv6S4GPS+WkU+jiBENICdnPvopT+3H35kLut4tv/w+7rsU23KaX73N9DAP4LJykVo3k949ALdm05Ttrm5nW6nmZ+utC9noSQ+QDuBvBuSukO92Xt+WkU+jjA3XbeCWAzpfRLvrfuAfBn7t9/BuAXvtevd3m0xQCWAniaEFIghMxwx2wB8CYAG8fjHBoJKa6nEO62d4gQ8jJ3zHcnHTMRUa/raeYng+71JIRMAfBrAB+jlD7GP5xqfmbtET4ZfsAiAiiADQDWuT9vADAdwEMAtrm/p/mO+ThYdMtWuJ5tsOiCte44mwDcDjf65WT6SXk9XwJwDMAwmOWz3H19FZjS2QHg63Czp0+mn3pdTzM/011PAJ8AMOL77DoAp6SZnyb138DAwGCCwFAuBgYGBhMERqEbGBgYTBAYhW5gYGAwQWAUuoGBgcEEgVHoBgYGBhMERqEbGBgYTBAYhW5gYGAwQfD/AeuYEyV1h6zNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "stn_swe.plot(y='935_CO_SNTL', title='Daily SWE [m]');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Export the sites and daily time series to files, for reuse\n",
    "- CSO_gdf to GeoJSON\n",
    "- stn_swe to csv or parquet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note:** Sometimes the statement below produces this warning:\n",
    "```\n",
    "CSO_SNOTEL_sites.geojson: No such file or directory \n",
    "driver GeoJSON does not support creation option ENCODING\n",
    "```\n",
    "\n",
    "Ignore it. The GeoJSON file is still written out without issues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "CSO_gdf.to_file('ClackamasLake_Snotel_SWE', driver='GeoJSON')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7304"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(stn_swe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The file that's written is tiny, only 0.5 MB\n",
    "stn_swe.to_csv('ClackamasLake_Snotel_SWE.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:whw]",
   "language": "python",
   "name": "conda-env-whw-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}