Fetch_SNOTEL_Clackamas.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": 18,
   "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": 19,
   "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' : 45.4562,\n",
    "    'latmin' : 44.8004,\n",
    "    'lonmax': -121.6515,\n",
    "    'lonmin': -122.6089,}\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": 20,
   "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": 21,
   "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": 22,
   "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": 23,
   "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>351_OR_SNTL</td>\n",
       "      <td>-121.856049</td>\n",
       "      <td>45.428699</td>\n",
       "      <td>Blazed Alder</td>\n",
       "      <td>1112.52001953125</td>\n",
       "      <td>POINT (-121.85605 45.42870)</td>\n",
       "      <td>-349122.260934</td>\n",
       "      <td>5.088226e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>398_OR_SNTL</td>\n",
       "      <td>-121.754433</td>\n",
       "      <td>45.096581</td>\n",
       "      <td>Clackamas Lake</td>\n",
       "      <td>1036.3199462890625</td>\n",
       "      <td>POINT (-121.75443 45.09658)</td>\n",
       "      <td>-346153.306893</td>\n",
       "      <td>5.050260e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>401_OR_SNTL</td>\n",
       "      <td>-121.691597</td>\n",
       "      <td>45.188320</td>\n",
       "      <td>Clear Lake</td>\n",
       "      <td>1161.2879638671875</td>\n",
       "      <td>POINT (-121.69160 45.18832)</td>\n",
       "      <td>-339846.412006</td>\n",
       "      <td>5.059786e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1029_OR_SNTL</td>\n",
       "      <td>-122.050003</td>\n",
       "      <td>44.816669</td>\n",
       "      <td>Epa Metolius</td>\n",
       "      <td>917.447998046875</td>\n",
       "      <td>POINT (-122.05000 44.81667)</td>\n",
       "      <td>-373695.245509</td>\n",
       "      <td>5.022342e+06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>651_OR_SNTL</td>\n",
       "      <td>-121.715797</td>\n",
       "      <td>45.320969</td>\n",
       "      <td>Mt Hood Test Site</td>\n",
       "      <td>1636.7760009765625</td>\n",
       "      <td>POINT (-121.71580 45.32097)</td>\n",
       "      <td>-339766.890081</td>\n",
       "      <td>5.074773e+06</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           code   longitude   latitude               name         elevation_m  \\\n",
       "0   351_OR_SNTL -121.856049  45.428699       Blazed Alder    1112.52001953125   \n",
       "1   398_OR_SNTL -121.754433  45.096581     Clackamas Lake  1036.3199462890625   \n",
       "2   401_OR_SNTL -121.691597  45.188320         Clear Lake  1161.2879638671875   \n",
       "3  1029_OR_SNTL -122.050003  44.816669       Epa Metolius    917.447998046875   \n",
       "4   651_OR_SNTL -121.715797  45.320969  Mt Hood Test Site  1636.7760009765625   \n",
       "\n",
       "                      geometry        easting      northing  \n",
       "0  POINT (-121.85605 45.42870) -349122.260934  5.088226e+06  \n",
       "1  POINT (-121.75443 45.09658) -346153.306893  5.050260e+06  \n",
       "2  POINT (-121.69160 45.18832) -339846.412006  5.059786e+06  \n",
       "3  POINT (-122.05000 44.81667) -373695.245509  5.022342e+06  \n",
       "4  POINT (-121.71580 45.32097) -339766.890081  5.074773e+06  "
      ]
     },
     "execution_count": 23,
     "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": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIMAAAD4CAYAAADPeFVkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAANLUlEQVR4nO3dbYxc5XnG8f+VxcFDUrKoBoVd0y6q0KotVnG0stKiRsUldepScIwqoYj0S5vmQ6MASRxlS5WiqpGpTZuoalUJCGqaVKAWHJcCjUtKaF4U0qy9gAGzCkkc8JhgR3Sb0KzA2Hc/nLMP4/Hs7rycM297/aSVd87MM/Os99rzMnPu+ygiMAN4U68nYP3DYbDEYbDEYbDEYbDkrF5PoJF169bFxMREr6cxlPbv3/+jiDi/0X19GYaJiQlmZmZ6PY2hJOkHS93nzYQlDoMlDoMlDoMlDoMlfXk0sZS9s1V275vj6PwCY6MVdmyZZNvG8V5Pa2gMTBj2zlaZ3nOQhRMnAajOLzC95yCAA1GQgdlM7N43l4KwaOHESXbvm+vRjIbPwITh6PxCS8utdQOzmRgbrVBt8IsfG60U+jqreb9kYNYMO7ZMUlkzctqyypoRdmyZLOw1FvdLqvMLBG/sl+ydrRb2Gv1sYNYMi3+dZf7VLrdfstLrDMMaZWDCAFkgyvwPbne/ZFiOdAZmM9ENS+1/rLRfMixHOg5DjXb3S4blSMdhqLFt4zg7t29gfLSCgPHRCju3b1hxVd/uGqXfDNQ+Qze0s1+yY8vkafsMUPyRTjc4DAXoxpFONzgMBSn7SKcbmt5nkDQiaVbSA/ntWyRVJT2ef21tdqz1p1bWDDcAh4Bza5Z9OiJua3Os9Zmm1gyS1gO/A9zZ6gt0Mta6q9nNxGeAjwOn6pZ/SNKTku6SdF6LY63PrBgGSVcBxyJif91dfw/8AnAZ8CLwVy2MbfQ6fyRpRtLM8ePHm5m7FayZNcPlwNWSDgP3AJslfSEiXoqIkxFxCrgD2NTs2EYvEhG3R8RUREydf37Dgh8r2YphiIjpiFgfERPAdcAjEXG9pAtrHvZe4KlmxxYzdStaJ29H75J0UNKTwBXATQCSxiQ9VMjsrKvUj218pqamwrWW5ZC0PyKmGt3nD6oscRgscRgscRgscRgscRgscRgscRgscRgscRgscRgscRgscRgscRgscRgscRgscUXVECiqUYjDMOCKbBTizcSAK7JRiMMw4IpsFOIwDLgiG4U4DAOuyJaI3oEccEU2Cmk6DJJGgBmgGhFXSboF+ACwWBj5JxHxUN2YtcBXgbPz17o3Iv6s5VnasopqFFJ2f4ZXgc0R8YqkNcDXJf17RDzWxlytZKX2Z4jMK/nNNflX/5VwGdCF/gx5C5/HgWPAwxHxrSUe55L8Hiu1PwNAXrZ/GbAe2CTp0iUe55L8Hiu7P0MSEfPAo8B7OpqxlabU/gySzpc0mn9fAa4Eni1i4la8Tt5n2CXpMrIdwsPAByHrzwDcGRFbgQuBz+WHpW8C/jki3P6vT7k/wyrj/gzWFIfBEofBEofBEofBEofBEp/PMGDKvGSiwzBAyr5kojcTA6TsSyY6DAOk7EsmOgwDpOxLJjoMA6Tsi8N7B3KAlH3JRIdhwJR5yURvJixxGCxxGCxxGCxxGCxxGCxxGCxxGCxpOgx5zeSspAfy27dIqkp6PP/a2mDMRZK+IumQpKcl3VDk5K1YZZfkvw58NCIOSPoZYL+khyPimTbmaiUruyT/xYg4kH//E7IwlfNeqnWs9JL8RZImgI2AS/L7VOkl+flzvBW4D7gxIn7c6DEuyS/W3tkql9/6CBd/4kEuv/UR9s5WVxzTzD7DYkn+VmAtcG5ekp+udi/pDqBhQW3evuc+4J8iYk8zP4h1ZrlzJZdTdkm+gM8ChyLir5v6Saxj7Z4r2cn7DLskHZT0JHAFcBNkJfmSFru+XQ68n6zBx5KHoFasds+VbOnkloh4lKz7ChHx/iUecxTYmn//dUCtvIZ1bmy0QrXBL35stMLhZcb5Hcgh1O65kj7tbQi1e66kwzCk2jlX0psJSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGS0rtz5A/7i5JxySdUXFl/aXs/gwA/wD8LfCPrU3NoNyLjdQrtT8DQER8FXi51XH2RgFtdX6B4I0C2mYqqtvRtf4MK3F/hjOVfbGRel3pz9AM92c4U9kXG6nXzJphsT/DYeAesorqL0TESxFxMiJOAXcAm0qZ4SpW9sVG6pXan8E6U/bFRuqV3Z8BSXcD3wQmJR2R9AcdzXgV2bZxnJ3bNzA+WkHA+GiFnds3lHY0oYgo5Yk7MTU1FTMzM72exlCStD8iphrd53cgLXEYLHEYLHEYLHEYLHEYLHEYLHEYLHEYLHHrP7p7Akk/W/VhWK4D+2oLxKrfTHT7BJJ+turD0O0TSPrZqg9Dt08g6WerPgzdPoGkn636Hch2O7APo1UfBmivA/swWvWbCXuDw2CJw2CJw2CJw2BJN0ry3yNpTtJzkj5R1MSteKWW5EsaAf4OeDdwBPi2pPsj4pl2JmvlKrskfxPwXER8LyJeI6vVvKbF57AuKbskfxx4oeb2kXzZGVyS33tll+Q3uvRxw3o+l+T3XjP7DIsl+VuBtcC5eUn+9YsPkHQH8ECDsUeAi2purweOdjBfK1HZJfnfBi6RdLGkN+fj7y9g3laCTj6o2iXpMrLV/mHgg5CV5AN3RsTWiHhd0oeAfcAIcFdEPN3ZlK0sLslfZVySb01xGCxxGCxxGCxxGCxxGCxxGCxxGCxxGCxxGCxxGCxxGCxxGCxxGCxxGCxxGCxxGCxxGCxxGCxxGCxxGCxxGCxxGCxpuz9DzfKPSQpJ65YYd4OkpyQ9LenGDudrJWplzbDYnyGRdBFZ74XnGw2QdCnwAbLS/F8BrpJ0SXtTtbJ12p/h02Sl+kuVZf0i8FhE/DQiXgf+i6wu0/pQ2/0ZJF0NVCPiiWXGPQW8S9LPSjoH2MrpVdmJ+zP0Xlv9GfJf7M3AJ5cbGxGHgL8EHga+BDwBvL7EY92foceaWTMs9mc4TNaGZzPweeBi4Il8+XrggKS31w+OiM9GxDsi4l3Ay8B3Cpq7FWzFkvyImAamAST9BvCxiLi29jF5IKYi4kf14yVdEBHHJP0csB341c6nbWUo/H0GSWOSHqpZdJ+kZ4B/A/44Iv6n6Ne0YrTUrCMiHgUebbB8oub7o2Q7iou3f73t2VlX+R1ISxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSxwGSzq5LFHp9s5W2b1vjqPzC4yNVtixZZJtGxteCdEK0Ldh2DtbZXrPQRZOnASgOr/A9J6DAA5ESfp2M7F731wKwqKFEyfZvW+uRzMafn0bhqPzCy0tt871bRjGRistLbfOdaMk/6a8HP8pSXdLWtvM6+3YMkllzchpyyprRtixZbLZKVuLyi7JHwc+TFZtdSnZhU6va+bFtm0cZ+f2DYyPVhAwPlph5/YN3nksUVNHEzUl+Z8CPlJz12JJ/r+u8BoVSSeAc2jhWtjbNo77l99FpZbkR0QVuI1szfEi8L8R8R+NHuuS/N4rtSRf0nnANWQV22PAWyRd3+ixLsnvvbJL8q8Evh8RxyPiBLAH+LWC5m4FWzEMETEdEevz4trrgEci4tqIuCAiJvLlR4B3RMQP64Y/D7xT0jmSBPwmdTuh1j9KLcmPiG8B9wIHgIP5691e9GtaMRSxVG+u3pF0HPhBm8PXAWc0Delj3Z7vz0dEw52yvgxDJyTNRMRUr+fRrH6ab9++HW3d5zBYMoxhGLQd1L6Z79DtM1j7hnHNYG1yGCwZiDBI+r38nIhTkqZqlr9b0n5JB/N/N+fLz5H0oKRn83G3LvPc05KekzQnaUuv5pzf9ylJL0h6ZZnnXSPpc/n4Q5Kmi5ozEdH3X2Td6SfJelBO1SzfCIzl319K9ikqZB+VX5F//2bga8BvN3jeXyLrZ3022Wct3wVGejHn/PY7gQuBV5Z53vcB99T8nIeBiSLm3LdnR9eKrCE52ccbpy2frbn5NLBW0tkR8VPgK/ljXpN0gOzDtHrXkP3Hvgp8X9JzZNfG+GYP5vxqRDzWaEz9U5N9+nsWUAFeA37c6XxhQDYTTboWmM1/sYmkUeB3gf9sMGYceKHm9pF8Wbc0nPMK7gX+j+z8kOeB2yLi5SIm0zdrBklfBs7oSg/cHBHLnUmFpF8mu5TBb9UtPwu4G/ibiPheo6ENljV9rF3GnJuwCThJdn7IecDXJH15iZ+vJX0Thoi4sp1x+Sl5XwR+PyK+W3f37cB3IuIzSww/wukXQ1lPC6fllTTnlbwP+FJk54cck/QNYAroOAwDvZnINwEPAtMR8Y26+/4CeBtw4zJPcT9wnaSzJV0MXAL8dzmzTfMaZYk5N+l5YLMybyHb6Xy2kMn1+kihyT3z95L9Fb8KvATsy5f/Kdn28/GarwvI/sKD7ESaxeV/mI+5Gvjzmue+mewoYo4GRxzdmnN+3658zKn831vq5wy8FfgXsp3PZ4AdRc3Zb0dbMtCbCSuWw2CJw2CJw2CJw2CJw2CJw2DJ/wNI2/o+tJENAAAAAABJRU5ErkJggg==\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": 25,
   "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": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "351_OR_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "398_OR_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "401_OR_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "1029_OR_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<suds.sax.document.Document object at 0x7f8cf1c4a390>\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Unable to fetch SNOTEL:WTEQ_D\n",
      "1029_OR_SNTL has no data\n",
      "651_OR_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "655_OR_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "687_OR_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "925_OR_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n",
      "1007_OR_SNTL SNOTEL:WTEQ_D 2000-10-01 2020-09-29\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<suds.sax.document.Document object at 0x7f8cf238f6a0>\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Unable to fetch SNOTEL:WTEQ_D\n",
      "1007_OR_SNTL has no data\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": 27,
   "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>351_OR_SNTL</th>\n",
       "      <th>398_OR_SNTL</th>\n",
       "      <th>401_OR_SNTL</th>\n",
       "      <th>1029_OR_SNTL</th>\n",
       "      <th>651_OR_SNTL</th>\n",
       "      <th>655_OR_SNTL</th>\n",
       "      <th>687_OR_SNTL</th>\n",
       "      <th>925_OR_SNTL</th>\n",
       "      <th>1007_OR_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>NaN</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",
       "    </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>NaN</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",
       "    </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>NaN</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",
       "    </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>NaN</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",
       "    </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>NaN</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",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            351_OR_SNTL  398_OR_SNTL  401_OR_SNTL  1029_OR_SNTL  651_OR_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",
       "            655_OR_SNTL  687_OR_SNTL  925_OR_SNTL  1007_OR_SNTL  \n",
       "2000-10-01          0.0          0.0          0.0           NaN  \n",
       "2000-10-02          0.0          0.0          0.0           NaN  \n",
       "2000-10-03          0.0          0.0          0.0           NaN  \n",
       "2000-10-04          0.0          0.0          0.0           NaN  \n",
       "2000-10-05          0.0          0.0          0.0           NaN  "
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stn_swe.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "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>351_OR_SNTL</th>\n",
       "      <th>398_OR_SNTL</th>\n",
       "      <th>401_OR_SNTL</th>\n",
       "      <th>1029_OR_SNTL</th>\n",
       "      <th>651_OR_SNTL</th>\n",
       "      <th>655_OR_SNTL</th>\n",
       "      <th>687_OR_SNTL</th>\n",
       "      <th>925_OR_SNTL</th>\n",
       "      <th>1007_OR_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>NaN</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",
       "    </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>NaN</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",
       "    </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>NaN</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",
       "    </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>NaN</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",
       "    </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>NaN</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",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            351_OR_SNTL  398_OR_SNTL  401_OR_SNTL  1029_OR_SNTL  651_OR_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",
       "            655_OR_SNTL  687_OR_SNTL  925_OR_SNTL  1007_OR_SNTL  \n",
       "2000-10-01          0.0          0.0          0.0           NaN  \n",
       "2000-10-02          0.0          0.0          0.0           NaN  \n",
       "2000-10-03          0.0          0.0          0.0           NaN  \n",
       "2000-10-04          0.0          0.0          0.0           NaN  \n",
       "2000-10-05          0.0          0.0          0.0           NaN  "
      ]
     },
     "execution_count": 28,
     "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": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABdgUlEQVR4nO19ebwcVZn281Z33yU3udkXkkASCCSELWASNtlmYACRTRgFVBREhvlEPhc+xREdFxxhcEVwkBFcEAUZRWGIgIjsawLZISuB3ITs691vd53vj6pTdaq6qvucU327+vY9z+93f7f7dFX16VOn3vOe592IMQYDAwMDg4EPK+0OGBgYGBhUBkagGxgYGNQJjEA3MDAwqBMYgW5gYGBQJzAC3cDAwKBOYAS6gYGBQZ3ACHSDAQMiOoCI2oko475/moiuSrtfKnD73E1Ez2qe3+iOQR8R3VTp/hkMbBiBblA1ENF6Iuoion1EtJuIXiSia4hIah4yxt5ljA1ljBUS9qOBiL5PRG2ucHybiH7ofvYVIpofOn51TNsl7mtGRB3utfjfl0p04VrG2Mk6fWeM9TDGhgK4T+d8g/pGNu0OGAw6nMsYe5KIhgM4BcCPARwL4Ioq9uErAOYAmAfgPQBTAHAB+yyAG4gowxgrENEEADkAx4TaprvHchzFGFtTvZ9gYFAMo6EbpALG2B7G2MMAPgLgE0R0OAAQ0TlE9AYR7SWiDUT0DX4OEU11teGAIuLSEDuJ6AihbZy7Gxgb8fVzATzEGNvEHKxnjP3a/ew1OAJ8tvv+ZAB/B7Ay1LaWMbYp2Sh4fV1PRP+PiJa4mv7dRDSeiP7i7maeJKKRlfgug/qGEegGqYIx9iqANgAnuU0dAC4HMALAOQD+lYguKHONHgD3A/iY0HwpgCcZY9siTnkZwBeI6P8Q0RFERMK1egG8Al9jPxnAcwCeD7VpceAlcBGAMwAcAuBcAH8B8G8AxsB5Tq+r8PcZ1CGMQDeoBWwCMAoAGGNPM8aWMsZsxtgSAL+DQ82Uw68AXCbw8R8HcG/Msd8FcAuAjwJYAGAjEX1C+PwZ+ML7JDgC/blQ2zOha77u2gX435kSfRbxE8bYFsbYRve7XmGMveEuVg8BOFrxegaDEEagG9QCJgHYCQBEdCwR/Z2IthHRHgDXwNFSS4Ix9goc7f4UIpoJh+N+OObYAmPsDsbYiXB2At8BcA8RHeoe8iyA97s0x1jG2GoALwI4wW07HMUa+jGMsRHC3+NKIwBsEV53Rbwfqng9g0EII9ANUgURzYUj0J93m34LRxDvzxgbDuBOABRzehi/gkO7fBzA/zDGusudwBjrYozdAWAXgFlu80sAhgO4GsAL7nF74ewkrgawiTH2tmSfDAyqBiPQDVIBEbUS0QfhcN+/YYwtdT8aBmAnY6ybiOYBuEzhsvcCuBCOUP913EFE9DkiOpWImoko69ItwwC8AThCHg4V8wU49AfH825bpflzA4OKwAh0g2rjESLaB2ADgK8C+AGCLov/B8C33GO+DuD3shdmjLUBeB0AQ1AQh9EF4PsANgPYDuAzAC5ijK0TjnkGwDj4Owe41xyHaIG+OOSH/iPZfhsYVApkClwY1BOI6B44lMiNafclCkT0BIDjASxgjJ2mcX4jHH49B+A/GWPfrHAXDQYwjEA3qBsQ0VQAiwAcbThug8EIQ7kY1AWI6NsAlgG41Qhzg8EKo6EbGBgY1AmMhm5gYGBQJ0gtOdeYMWPY1KlT0/p6AwMDgwGJhQsXbmeMReUoSk+gT506FQsWLEjr6w0MDAwGJIjonbjPDOViYGBgUCeQEuhEdBYRrSSiNUR0Q8wxpxLRIiJaTkThxEUGBgYGBv2MspSLW+7rDjipPdsAvEZEDzPGVgjHjADwUwBnMcbeJaJx/dRfAwMDA4MYyHDo8wCs4WHRRHQ/gPMBrBCOuQzAHxlj7wIAY2xrpTtqYGBQWfT19aGtrQ3d3WVzmBmkgKamJkyePBm5XE76HBmBPglO3g2ONjglw0QcAiBHRE/DSXL0Y6ECjAciuhpOtjoccMAB0p00MDCoPNra2jBs2DBMnToVQo0PgxoAYww7duxAW1sbpk2bJn2eDIcedafD0UhZAO+DU2HmTABfI6JDIjp5F2NsDmNsztixkV43BgYGVUJ3dzdGjx5thHkNgogwevRo5d2TjIbeBmB/4f1kOHmhw8dsZ4x1AOggomcBHAVglVJvDAwMqgojzGsXOvdGRkN/DcDBRDSNiBoAXILiSjB/BnCSm1t6CBxK5k3l3hjUNdZta8eLa7an3Q0Dg7pFWQ2dMZYnomsBPA4gA+AexthyIrrG/fxOxtibRPQYgCUAbAA/Z4wt68+OGww8/MP3HW/W9Tefk3JPDAzqE1KRooyx+QDmh9ruDL2/FcCtleuagYFBvaO7uxsnn3wyenp6kM/ncfHFF+Ob3/wmFi9ejGuuuQbt7e2YOnUq7rvvPrS2tqKvrw9XXXUVXn/9deTzeVx++eX4yle+Env9PXv24LOf/SxeeOEFAMCJJ56In/zkJxg+fDjWr1+PQw89FDNmzEBvby/mzJmDu+++O9arpLOzE5/+9KexZMkSMMYwYsQIPPbYYxg6dCiICF/4whfw/e9/HwDwve99D+3t7cjlcnjwwQcBAEuXLsURRxwBALjyyiuxc+dODB06FNdff33FxtNEihoYGKSGxsZGPPXUU1i8eDEWLVqExx57DC+//DKuuuoq3HzzzVi6dCkuvPBC3Hqroys++OCD6OnpwdKlS7Fw4UL87Gc/w/r162Ov/6lPfQoHHngg1q5di7Vr12LatGm46qqrvM8POuggLFq0CEuXLkVbWxt+//v4Alk//vGPMX78eCxduhTLli0LCP/Gxkb88Y9/xPbtQUrxq1/9KhYtWoRFixahubnZe33dddclGLV4pJbLxWBwYV93X9pdMCiBbz6yHCs27a3oNWdNbMW/n3tYyWOICEOHDgXg+MX39fWBiLBy5UqcfPLJAIAzzjgDZ555Jr797W+DiNDR0YF8Po+uri40NDSgtbU18tpr1qzBwoUL8cADD3htX//61zF9+nSsXbsWmUzGa89kMpg3bx42btwY29f33nsPU6ZM8d7PmDHDe53NZnH11Vfjhz/8Ib7zne+U/M39CaOhG1QFO9p70+6CQY2iUChg9uzZGDduHM444wwce+yxOPzww/Hww47vxYMPPogNG5xQmIsvvhgtLS3Yb7/9cMABB+D666/HqFGjIq+7YsUKzJ49u0hwz549G8uXLw8c293djVdeeQVnnXVWbD+vvPJK3HLLLTj++ONx4403YvXq1YHPP/OZz+C+++7Dnj17tMahEjAaukFVkM0Y97haRjlNuj+RyWSwaNEi7N69GxdeeCGWLVuGe+65B9dddx2+9a1v4bzzzkNDQwMA4NVXX0Umk8GmTZuwa9cunHTSSTj99NNx4IEHFl2XMRbp+ie2r127FrNnz8bq1atx8cUX48gjj4zt5+zZs7Fu3To88cQTePLJJzF37ly89NJLOPTQQwEAra2tuPzyy3Hbbbehubm5EkOjDKOhG1QFYmEs2zZVsgyKMWLECJx66ql47LHHMHPmTDzxxBNYuHAhLr30Uhx00EEAgN/+9rc466yzkMvlMG7cOJx44omxabgPO+wwvPHGG7Bt22uzbRuLFy/2hDDn0NesWYOXX37Z2xXEYejQofjQhz6En/70p/jYxz6G+fMDviL43Oc+h7vvvhsdHR1JhkIbRqAbVAUFQYgXTNlDAxfbtm3D7t27AQBdXV148sknMXPmTGzd6qSDsm0bN910E6655hoATsqQp556CowxdHR04OWXX8bMmTMjrz19+nQcffTRuOmmm7y2m266CccccwymT58eOHa//fbDzTffjO9+97uxfX3hhRewa9cuAEBvby9WrFgR4NQBYNSoUfjwhz+Mu+++W20gKgQj0A2qgrwg0PMFI9ANHLz33ns47bTTcOSRR2Lu3Lk444wz8MEPfhC/+93vcMghh2DmzJmYOHEirrjiCgAOT93e3o7DDz8cc+fOxRVXXFGSJrn77ruxatUqTJ8+HQcddBBWrVoVK2wvuOACdHZ24rnnnov8fO3atTjllFNwxBFH4Oijj8acOXNw0UUXFR33xS9+scjbJQ433XQTJk+e7P0lRWpFoufMmcNMxaLBg5Wb9+HMHz0LAFj6jX/CsCb5DHIG/YM333zTox4MahNR94iIFjLG5kQdbzR0g6qgr+DzmAXDoRsY9AuMl4tBVSAK8bwR6AYVxrHHHouenp5A27333utFZqrg8ccfx5e//OVA27Rp0/DQQw8l6mM1YAS6Qb/jB0+sxNKNvm+u4dBrB3GufQMNr7zySsWudeaZZ+LMM8+s2PV0oUOHG4Fu0K/ozdu47ak1gba84EZmkB6ampqwY8cOkxO9BsELXDQ1NSmdZwS6Qb+iq7dQ1GY49NrA5MmT0dbWhm3btqXdFYMI8BJ0KjAC3aBf0dGbL2rrM5RLTSCXyymVNzOofRgvF4N+xe7O4qRcRkM3MOgfGIFu0K/4xQtvF7UZDt3AoH9gBLpBv2LqmJaiNqOhGxj0D4xAN+hXHDBqSFGb4dANDPoHRqAb9CvsCF9ao6EbGPQPjEA36FdECXTDoRsY9A+MQDfoV3DZffcn5uCOy44BAGw31YsMDPoFRqAb9Cu4hn7I+GGYOMKJervud2+k2SUDg7qFEegG/QrOuFgWmaRcBgb9DCPQDfoVXEO3KDoNgIGBQeUgJdCJ6CwiWklEa4johojPTyWiPUS0yP37euW7Ko+CzTD1hkdx70vr0+yGAfxycxYRchmjPxgY9CfKPmFElAFwB4CzAcwCcCkRzYo49DnG2Gz371sV7qcSevOOJe5rf16eZjcMAHCWxSLCcQeOAgBcdEzyUlsGBgbFkFGZ5gFYwxhbxxjrBXA/gPP7t1sG9QImUC5EhLHDGtGQNalaDQz6AzICfRKADcL7NrctjOOJaDER/YWIDou6EBFdTUQLiGhBf6bsjPJ9NkgHtu1TLgCQtcgEFhlUDO09edz70nqtYhD1CBmBHqVOhUfvdQBTGGNHAfgJgD9FXYgxdhdjbA5jbM7YsWOVOqoCI9BrByLlAgAZ4+1iUEF84+Hl+Nqfl+PFtTvS7kpNQEagtwHYX3g/GcAm8QDG2F7GWLv7ej6AHBGNqVgvFWHkRe2AL67kzrSM0dANKogd7U4d0e4+40EFyAn01wAcTETTiKgBwCUAHhYPIKIJ5NawIqJ57nVTWzLN9qt2wAV6xmjoBv0IU0HPQdmKRYyxPBFdC+BxABkA9zDGlhPRNe7ndwK4GMC/ElEeQBeAS1iKUtXIi9pBmHLJWuTx6gYGSWFmUhBSJehcGmV+qO1O4fXtAG6vbNf0Ybb0tQOPcnE1qIxlGQ3doGLgaiNFmvoGH+oy0sNQLrUDVmQUNQuuQeXgzSQjzwHUqUDvyZv0rLUC323ReZ+xLCPQDQz6CXUp0J9fsz3tLhi4EEP/AeOHbmDQn6hLgX74xOFpd8HAhWcUdVX0DJEpcGFQMRimJYi6FOhGYNQOGGMe3QI4bovm9hhUCmavF0RdCvR93fm0u2DgwmbMo1sAIJsxGrqBQX+hLgX6V/+0NO0uGLiwGQIC3SJ9Dr1gM9z17FqTV93Ag6FcgqhLgb5hZ5f32rgwpgvbZoEovqxFnqFUFX9dsRn/Mf8t3PyXNyvUO4O6gXnMAdSpQG/M+j/LeFSkC5sxZAQSPWMR8gW9e8Jv5Za9PZXomkEdgCsLJiGfg7oU6J8/4xDvta42aFAZhCmXJMm5+FWYUccMQjCKm4O6FOgtDRnvtbG/pQubBSmXRALd08Yq0DGDuoIR6A7qUqCLt9Zo6OmChTT0JBy6m9AT5pYacPC51WcEOoB6FejCvTUrd7oo2GE/dEubQ/cvY+6pgQNuL+vqNa7KQN0KdP+BN6la00WxUVR/kb335Xcq1S2DOgHf/BUMtQpAMn3uQMIPnliJnz27zntvKJd0YTOfKgHc5Fya9+S51U6OHpN8zSAM85w7qDuBfttTawLvjYaeLsKh/5VIzvXuzs6EvTKoF/A86AWjogOoU8pFhFm5gTVb21OruRgO/Xf80JM9fPOmjkraLYM6AXdh1TTL1B3qX6APcg19X3cfTv/BM/jyH5ak8v1RfuhJb8nQprrbWBpoYvOebgBmJ85R9wJ9sPuh87wnL65Np2Z3VOh/0uRcQxuNQDdw0Jh1Yk5MWUMHdS/QBzvlsq/Hcefa09WXyveHvVysCnDog33XZeCDUy4m9N9B/Qv0Qf7w3/uS4+rXm5JnSBHlkiDbIkc9LNKMMZM4rgLoc8lz3diGekPdC/TBvnI3ZNO9xeHQf4v0w4KGN+eca9bBIn35Pa/iU79akHY3Bjy4olIPi3wlUPdk5GDX0GdOGJbq94dD/0GkHbrPBXk9eKhxn3qDZHhr814A9bHIVwJ1r6EPdoGe9lY0HPrPX+vQDVwLKwxwS/d7e7rKH2QgBY9yGeTPOUfdCXQrVMJksFMufKIPETJQVhNhP3QeCKLz/PHFeaBvr03Fpcqgo8fP3zLYn3MOKYFORGcR0UoiWkNEN5Q4bi4RFYjo4sp1UQ2B7T2Mhp628CsuQef819HQbVYflEtjLp3Ftd7wXaFyVdo70VpBWYFORBkAdwA4G8AsAJcS0ayY424B8HilOymLnnyhaOs12FduHhKd1jAwxmAJsyxJTnNPQx/glIvxbqkMfvPyu97rwf6cc8ho6PMArGGMrWOM9QK4H8D5Ecd9FsAfAGytYP+UsL29t6htoGtzScEVl7QmfBHlQuS1q4Ax5i0CA50vNbKn8hjsO3EOGYE+CcAG4X2b2+aBiCYBuBDAnaUuRERXE9ECIlqwbds21b6WRRQ3OdhvNNdm0xqFQijbYpgSk4V4Gwe6R4PRJiuPgb7IVwoyAj3qCQyP3o8AfJkxVtLawxi7izE2hzE2Z+zYsZJdlEdUNORgf3j4RE9rm88YQybotQhA/b6Ixw/0h3ewKxn9gYFOw1UKMn7obQD2F95PBrApdMwcAPe7mtgYAB8gojxj7E+V6KQs9nRFUS6D++GxPYGe0veHKBffKKp2HfE+DvR72t1nhE8lMHPCMLy1eR+A9CKhaw0yAv01AAcT0TQAGwFcAuAy8QDG2DT+moh+CeB/qy3M45C2l0fa8DT0lL7fthHjtqivoQ90gf729o60u1AXmDJ6CABnbpiiJw7KCnTGWJ6IroXjvZIBcA9jbDkRXeN+XpI3ryaiXJcGOt+aFH50JUO+YCObqW7oQTj0n79WvSv1pKF3mPqXFUHBduwzDRkyAt2FVOg/Y2w+gPmhtkhBzhj7ZPJu6YFro986/zDs687j1sdXDviHPylEvrmjp4DhQ6ov0LOC3yI3kDLF50+kSAc6h25qXFcGjDFkLCeFbk/eBGsBdRYpyh/0Ew4ajdNmjANgjKIi5dRTqP6ktxkCfuiWplFUzKE+0O+p2H/jk66PgmufacxaeG93d9rdqQnUlUDnlu6sZXk5uAe9H7pAQ6VhOCo2iroauuJ1OgWX1IEeFRhwwRzYPyVV8CjkPV196EqpxGKtoa4EOn/QMxaBU8WD2Si6dW83fv782977vhQEYTj0X9dtcV+3zzv3DvBVmqF+7AFpwnYTv00ZPST1NNFvvrcXv39tQ/kD+xl1lT6XUy7ZDKGv4HpTDNIHhjGGef/xt0BbGr66jAWzLXocuuJt4RpYY9YKJGUaiAhq6INzflYCvBpWUy6T+sJ49o+fAwB8eO7+ZY7sX9SXhs4FeoByGZwPTJQWm4YxsTjbogNV7rjP/T2jWxoC2vpABKsjF8w0UbAZiAjZCpQ1rBTStonUlUDniaiyFnlCZLBSLjsi8tqkwT1z1zIOXQ6dC/QRQxoGvNufOCUH6/ysBBhzShoyFj3f00DadGBdCXSugWYy5Gnog5VyiQpeSUOL4a5lHLpeLlygD2nI1IFR1O//YJ2flUDBzeT54MK21AUpR3evEegVAxfoOZFyGaQaUNTPrgnKRTN9bm/eOaEW+NKkELs/0H9LmgjPrVoYy7S9bepKoPMbmhEol8GqAY0e2lDUlsaEL/JyATeK6mnojVlrwBsSAxz6AP8tacLxciF86awZAIDuGnBd3NWZLvVTVwK9T+DQB7tRNBcR4p9PwcslNvRfWUN3+t6Uy9SBQPdfmySB+nCUBaDZrQBVC+H/7Sl7YNWVQOcFiS2LkPGMoil3KiVEacCpaOh2dIELVZkc1NAr1r1UYBsNvSIo2I7bYraGlLe0dwl1JdDztp83hIebD1bKhQuKsw6bgB99ZDaAtDh0eLslQHBbVPRz8QR6HXDoYu8H6/ysBJzdH8HiDhA1sDi2p+xSW1cCvWD79SsHu1GUC70Ljp6E6eOGOm2pRIoGKRd+f9Q1dG4UdS4wkAVhPaUCThPcbbGWNPTr7n8j1e+vK4Eubu89P/QauMlpgHOzGcGekIaGzmKMorrJuXiIdy1oY7owfuiVAXdbTPtZzwsuk2mk1xBRVwKdwd/SD3Y/dC4oMhZS1WDsotB/5796YJFzRmPWMYANZEEo2jcG6/ysBLjbYjaTrkDvrAHvGo76EujMN7plUo4U7SvY+Ozv3sCarftS+X4+uS0SNfTqewEU+sEoqnN+LcE2GnpFwHfkfH6llSd/274e7/XMCcNS6QNHfQl0+HytlbKGvqRtNx5ZvAlf+p8lqXw/8zR08gzFafmhE0UYRVUpl4Kj6dcSX6oLw6FXBtzgnknZKPrODj8qe+7UUan0gaO+BDrzBQbg3Oz0NHQ/UVgaCGjomTQ59GDovzblYtvIZWojAvj2p1Zj6g2PalfJMX7olYGTnMtf5NNKCVGooWpadSbQWUAbzBClVuCCC1TO71X9+5kv0NOc8LEFLlQplzxDLmP556coCL/3xCoAQGePrkA3fuiVAGMMGYFySUtDD9a7NblcKgYGFLnIpXWT17nJsUQf7GoiyssljclWHPrP29W9XLIZ8gystSAIdXsQ8EOvgd8xUFGoEaOouEAbDb2CKKJcKL08yV/70zIAQWFWTUR5uaSVnKsSof99BRakXAYw92yyLVYGW/b2YN329tSNoqJykXYm0PoS6AhSLlYNJL7fsjed4rV2hJdLTYX+a0SK5iw/KjDtQgJJYLItJgfPmfLa+l2pG0WDlIsR6BVDlFE07S3t+6ePSeV77QCH7tzmmgr9V+xKvmAjK3DoaVEulVhITGBRcmzc1eW9zqRsFOXPWmPWSsU1WER9CXQgwiiazk0e3VKcvraaEFMJp6mhd/UVAkJQ2w/dZshlhKRraUUFCt+rK9yDgUWJuzQosaerDwBw0wWHe3MiPQ3d+d+QtYyGXkmworwh6WnoPJVnWrxeUENPR4N5ed0OAMD9QjV0rqwrUy55x21Rl4OvFPoq4DZlsi0mx15XoB85eXjqRlFObzoa+gAQ6ER0FhGtJKI1RHRDxOfnE9ESIlpERAuI6P2V72p51JJRlOfvTqs0VkHwcrEsAlH1vVwWrN8JAIGizroVi/I2c71c0tXGeiuQczvoh24Eug72djsCvbUpl3ouFz4XGzLJNfQX1mzH6i360eVlBToRZQDcAeBsALMAXEpEs0KH/Q3AUYyx2QCuBPBz7R4lgBP677/PWOn4ods28wR5Pi2BLni5AI6ni472wBjDD/66Cpt2d5U/OISoREW6FYvyNkNGKC2Ylhzs7E2et+PhxZu812lv0QcquIY+vDmXaiQ04D9ruayVeBf80Z+/gjN++Kz2+TIa+jwAaxhj6xhjvQDuB3C+eABjrJ35T2gL9F10E4GBeQIDSM8PXdTK0zLU8Nvh5bbR9PhZu60Dt/1tNc67/QXlc7vdSMqsaBTVjBR1gkhEDT+dca2EQN8q5P4wlIse9nQ5u75hTVkvJXNq9KZdOQ09KWQE+iQAG4T3bW5bAER0IRG9BeBROFp6EYjoapeSWbBt2zad/pZEkYaeEuWyaMNu73V6lIurobsDkrX0+L3mBie7oU58VE+f89snjWz22nyjqFpfeMSppXl+pdDZ69NHleiBoVz0sLe7Dy0NGWQzVs24LTYMEC+XqEe5aOQYYw8xxmYCuADAt6MuxBi7izE2hzE2Z+zYsUodlQFDsLNWSrlc7n3pHe91evklfC8XwElBoLO48aIYOr+CL2a/+/RxXpuu26JtIyDQa4FyqcTUMhq6HvZ29aG1OQcg3cA5wC9zOVC8XNoA7C+8nwxgU8yxYIw9C+AgIqq6A7aYPhdwtNM0NCBRI01rxfa8XLhAt0irL/wcHbnT1VvA5JHNmDjCHw9Pw1a8Fo84tVKmXLoEga7bh4aM5RXqSFsADFQs37QXrU2OQPcW+ZRD/xsyA8PL5TUABxPRNCJqAHAJgIfFA4hoOrmSlIiOAdAAYEelO1sOYbdFXd44Kbhr27QxLehNOQNcJiGHzs/RoTg6e/MY4lI2HB4HrtgXTqdRyh4NIoWm24fJI5sxY7yTNzvtwLeBiu3tPV60aOpGUYFy0e3D3u4+/NtDSxP3JVvuAMZYnoiuBfA4gAyAexhjy4noGvfzOwFcBOByIuoD0AXgIywFktMJLPLfW5SOH3q3yx2PamlI3cvF8rxc9CzwXOPQGcXO3gKaG4JTzC8SrQYGBot8vjQtOSiOoe7DW2BOkJRzjYp0a9CBCJg3zck9zud42l4uSTT0Gx9aFvB+0kVZgQ4AjLH5AOaH2u4UXt8C4JbEvUkIxoJeLmlp6D19DtXQkLEqEoiiA76eWgk1dC7AdBbG3rztVRji0I0U5Vkb06ZcRNpKtw/5AkM2wysvGQ1dB4z51as8DT2lsbQroKFvrlDOp/qKFEU4fS4hDcajs7eAIQ0ZZDOUWtFYLoh9Lxc9P/QkHDqvyi7Cj/RU93JxOPR0jaKV0NBtQUM38lwPNmOefShtt0Ux9D8tBY6jvgR6UaRoOoaStt2dsBnfgqVzg19yw+6522FG0yi6vb0XANCtUQiXV2UXoUu5cA09bT/0Smjo7+3p9naShkPXQ8FmnndLJmWjqBdYlMAPvW1nZ0X6Ul8CHSEvl5Qol2Ub92LN1nZHQ8+nM8mGNTpsWlNOEOgauwX+0AxtlGLnAghXKwIgpL9VuxZjTk3RtD0axB2XjkbIQ9afX7MdQHo7jYGOvJCWOW2jKHfGyGX0dsEAMK61qSJ9qS+BHk7ORenVFAWcFbsvJQ09bzNMG9Pivdf1Q+dbSAoJZhmEqxUB+hWL+OKQduh/PqGXS3inYzR0Pdg28+ZC6kZR2ymFl7X0NfRK1cGpM4EeUSQ6RRUol6JRNG/bgTzkGc1IUX/81M91ilsE23RD/23bWVTSp1z879VZq/k5PL2yMYrqocB8yiVto2jB5fOdXbDe816pHWd9CfRQxaL0BboezVEJ5AsskEMlq+vl4vmhq/fBZgzhmqrJQv+RerZFUaDrCBC+wP/LKQcCMJSLLgp2sVE0zfS5GTdNtW4fRlaofkJ9CfQIDT3NLW02VQ2deXmiAX2jKJ+gOuNos2KqRjf0nwcWWZpuj5WCqIE9u0o9HxGfDw2u26LOuF796wW46lcLlM+rJ3CaA0DqRU8KtltIJgGH/tZ7+ilzRdSfQK+B5Fwcjh96eppkVnAxSayha/QhmnLRqynqBBal74cu3s+7nl2nfP5Db2wE4Hi6AHoa+hMrtuDJN7eon1gnYIw59hnu5ZJyLhe+ezQaeoURTp+bNuWStSg9Db1gByiXjK4futt/HY4vinLxKhZpBhalHfq/Zmu791qnDzwacO22DgCGQ9cBH3Y+v/mc2CakJa5ufziH7tipdO5pvmDjA0dMwKQRzZg4XN/jpb4EekhDT8MoKd7MSiS818WLa3dglVD5JKmGruPlUmAsgnLR81LhgUVph/4/uvQ977XqLgNwSukBQEM2+e9YsWmv/skDGJw6DCsL1a7I5X+vz6EDeruu3oJTYvHYA0d5Ow8dqDsX1zDCnK2uq14S8EAcAMhZhD7bdt0pK+SXpIC9Qum3jGYulyTjx6LcFjUjRfm10qZcRFwwu6gsQFn4laT0OXSOfa5P+0DEog278erbOzBxRDMyRDj7iP2kz+VyW5xbY4c1Bnbn1YSvoXPqx0bGypQ5KwheMxdItsjXlUAHWJFRtNoc9n8/5/Cqo1oakMtYYMyNasukM9k4kmroOrDdKkMitN0Wi7xctLtVMYh+/rL42LFT8P2/rsKn3j8NjyzelOh3hDXUgYQL7ghWwFp/8znS5/JFUaQUhzRk0JNPXk1KB2ENXec56y0wp0BGwU60yNeVQC+iXKzqh97/yTV6DW/OeQmYHI+TqnYDk0c2Y97UUd57xwKvNhaMMXS5FXp0eMGCXRwpqlpTtDdvoyFreYFFafuhi9BZ7HhRhvGtjQCS/Y5KFKweiOBFV0RqojFreVlOqw2buV4uCYyzfQUbDRkLts0SzYm64dB78gVs2NUZEOhpUC68XuR3LjzcS8CURhk6W/DTBfQ09O89sRLfe2KVdh8YQxEf6FMu5c9/d0cnDrnxL3hwwQYvsCjt0H8A+IeZ4wDoVaPiD3tDBbItplXeMG14tJUwtRqz6Wnotu3YdzwNXWNe9BVs5DKO0T/J1K4bgf7Vh5Zh1Zb2QImwbAqUC8dx00Z7nFgahtECY4FMhzpeLr968R3vtV7of7HbokrFIm7UfWzZ5uJcLinJ84aMhUPGDwORXjUqbrjLen7o+n15dtV2/ZMHMMLlFfnr9ErQMdcP3bmnOuk+el0O3aJki3zdCPRn3CCPjp6wITAdLcayyBPoabguFmyHZuHQ0dDbhbHU2QZGUi4KlAk/4tW3d2LTnm43l4t+fyqBvO24g+a0Uyk4/5MEFnHc88Lb2uemjSMnD9c+1xfovvhKGkT4+PLNeHmdXpG1pBy6bTPkbeYK9GQael1w6AWbeT6ookaesSpTmV0X3BCajkC3Qxq6mgDaui+YcF+HuhKDPzhUIkX5A7rPXVjEEnRpCHQe0JKxCNmMXt6OQsjlrhaMu2lgneuHr4PdXY4nWa9AsSQNIvyXexcCUDPOcnAvF69YtcKOvK9g498fXg7AyaduUbK5XRca+sZdXd7rnR2+22BaJeg4GjwNPQXKxQ4G9WQVEwd965EVgfc62mgU5UIKlMvabe2B95ZFqYb+c4GR5YmYEmjoWYsSb6/FPg0kbNzdFdj9AVCam/e9/C4A4IEFbV6bZeklS6sEPA09o66hP7F8C377ivN7GjKWw6EnuKd1IdB5bukwkhoYdDBtTAvOO2oiAF9DT4P2ydssUaRoWMvQ09DjKRcZQXbb31YH3je4HCO/drXBxy+Tceg0Pb9+Zy7wxSnp7xiIni6/jKCKuhQKqJw6YywA4EtnzfDaMlZ6qbJ78zYYmEcBqTxnoiE3l3HmRJKfURcC/ZcvRnOJldCAVNFXsD1Bzjn0NLwRuLsfR07R4yf8gBU0QprtSLdFBzKXmjF+WOB9VtDQ09BMizV0DcpFSIeQlC8F9AyzaWNYU66oTWUXy8dsTEuj12almLfp7yu3YdWWdi0OXaRjc4ZycRA3gJV4YHT6wm9sLqPOqVUCfQUbeZuhOec7v1uKGrpIXXGoUgxRBS58L5fy1zrX3elwZDOWdsWjSsDT0C0LOc0KUDwzH+DsVpJq6ANQnmPkkCiBLv9Der30CZUxim7a3VX+IAmIkaKy2NvlU085d34PerfFeIFe/a153va3Xp7bYpWfOq5d83qigCtIFYZixoRhRW2qAsxJzhVs87xcJIYk7CrZkEk39D+goWumShWN1Trb6xdD9GKaFbl0MXxIcWZBJYFecOZ3TkwPnUBDP+Hmp7TO45g4vAnHHThKS0Nvbfb9UgjJF/m6EOjnHBmdB4LcB6aatIuoofP0tdXOAtfdWyzQCWoTZfLI5qI21YUpkkPnkaIS54eTLWVdty7n2kpdqQjEpFBJ3BazHuWiHiB12c9fCV1v4An0qOdRhXKJ0tCtFDOrDmnMYlRLg1akaItQq/d/FrYZDh0AJo8cEtmehkeEyKH7XibVzbnBOfuc4KdrESm5cBbsYH1WQENDtyMKXCgYRcNcazaTbui/GNCiW26sYNt+pR1FSjBK+NdCCgRVtO0qpjhUxrLXnYcByiVFjzZuK9IpVi0atYc0ZKrDoRPRWUS0kojWENENEZ9/lIiWuH8vEtFRuh1a0rYbU294tMhlrRRKUS5AdSe9qKHz7VS1Jxr38V3UtttrU93K9RWYRxl5bQoa+r7uPvQWbPx1xeZAu0rof/iYBkFDTyOPOF/QuEDX0Qh/9dI72NPlZEl0WDA9Q3WSRFBpQ7TtcKg4DvS449AkXCdJtaCk4PVNPQ1dQfHpEQT6N847LLHnU1mBTkQZAHcAOBvALACXEtGs0GFvAziFMXYkgG8DuEu3Qz960nFVe3GtfNQWH4B7PjkHj173fq+da0LVvM8ih55Waay/r9wKAHjELaYA+PSTLHiyIBEqE5XTTGtDASQqFYvC/LBF6Xq58HnGvW1UH7x3dhSPhe6zy4XXQBToUXy5ytzq5gJdyHiXSei/nQT5ghtYpOGH3iMs0qNaGqqSy2UegDWMsXWMsV4A9wM4XzyAMfYiY2yX+/ZlAJN1OzTU5ZSGNcoHsfIBPHr/kThs4nCvPY3tecFmnrHGX1CqO9FE4x2H7y4o1xeeLOiVf/tHfO70g702WXBt5dsXHB5oV6lYFO7r31du9RbJdDj0oIau2odTbn068F51ex2lgQ5EyoXTDJ84forXpjK3evI2iEJG0RT90G03b5KOl4uooYv5/nV3oDICfRKADcL7NrctDp8C8Bet3sDnxVQCJvhkyGWDP8d/+KtzoxljgQjNtDR00b2OQ9We0FdgyGYsjG9twoFjh7pt8veE96G1Kbgwq1QsCo/bhp2dIPcnrdxcmaK6KvAXSnV/4SiOWFXLFzXQT54wNdCngQQ+j46ZMtJrU6FceNCcaJ+xiNKNFLX0crmIAr0xm9zoLyPQoyx6kV9HRKfBEehfjvn8aiJaQEQLtm2Lrpje6Arlx5Zvjvw8CtxCngtVU6i2R0Q+pBlnUuI5efrOgIauuFsRKZcGngY4L/87ojLiif2QolxC4za8Oectkg8s2BB1Sr9C5NBJ0U0uSmCpbq+5Bvrv587C0QeMADAwc8H0FhgaMlaAA1ehXKIC1jJWZZ4zHc04XLFIxWOnJ19AxiKsv/kc14vLv6YOZAR6G4D9hfeTAWwKH0RERwL4OYDzGWORBDhj7C7G2BzG2JyxY8dGftl9bl6Dp97aKtE1B30RXh1On5z/1dLQw1ng0qJc+DAME7RjbysneQ1OuQB6/vRRtA+QzCj6nQuPSLVKTzhSVOW2Ru04VSOZ+fc3ZC0hudfAk+h9BSeKWRToqru/8LxKQrm0CO69Hb3qOdX9bIvqXi75AgspXsnuq4xAfw3AwUQ0jYgaAFwC4GHxACI6AMAfAXycMaZfEUETfW6F+3BmP49mqNJWrKuXW9/DRtHqfD/HSQc7i+W3zvf5a9WJ0lfwaxx6KQwUaLCoFKeAWsWi8AN6+KThRQ9yNeH5oWcIC9/ZhefXbFeqvMQxqsUJrFGlCbwxJUqNzqsEvGIOgTa1hS3qWdc1iorzTLdMI8/A6byXv6k2CybRS+pqXdbyyBjLE9G1AB4HkAFwD2NsORFd435+J4CvAxgN4Keu4MgzxuaodkZ8OI5SyJcc5WIHVN9tkWsZ3A7AZVm1jTVf/P1iAKEoNAXNGPA5dEAvJ43P4wfbVWqKivft+S+fBsAvDHHGrPHSfakUonYd3X12IIArDiJX+nnXyKzKw/NjLUF5GYgCnRdzEBc5FQ09LASBZBq6bTtaekdvQSu2wE7AoYfTY/CfpXtfpVxJGGPzAcwPtd0pvL4KwFVaPRAgTvrRQxtLHBmEM0GKNbdqUx55QYMS/1fbZ5r7K4sPjK8Zy13D4dCdc3QM1X7e76BEV9FANu9xcrL/5lPHBoLHxrc2YnRLcfh4fyN8fwHHN78Z5QU6Xww/ecJUfOw4x7tDlUPnY2YJGvoAZFywrzuPYU3ZQOF0VTovQ8Uauq4QtBlDQ9YV6DrRv7xikUakqM2CAXxWFSiXqqFL4K9UBkWkB0T4NEPyvskgbAhM02caAKaO9qvS+xy6XF+eXrkNG92kRY1aAt35X8Shu/9lJuyvX3oHAHD8QaND10gnKjDK0Cu7zW/vdpIwHX3ACG9eWpYeh26Rv/vTHYfHlm3G0rY9WucmxZ6uPrQ25wLj2KdgcOdGSBEZKxnl0uj6tOsUo7FtaEeKspCG7tv9lLsBoMYE+noh8CKcx6MU4gR6Up9OVcQJ9GrL83OPmogDx7RgpKDFqkwUnqN5e7uTcTGnUagjH6rME+6Hyi0J0+aOMVH+/KTY3t6DG/6wxKtXK2qWsg/v3c+/DQD43avvem3Kboss6GUD6NN51/xmIc69/Xmtc5Nib3cfhjfnILLoKnRelIauS7kwxsAY0JjTr/9bcJPQ8UVWVUO3IjR0XZlVUyXoRIGuMrD5AkMuG0G5VFtDZyGB7t7galMuBdsuEqQqE4Vrkxye+6WGNlks0HmkqDyK88FUNy3ynJueBOBTgiKNJCvQP3DEfnh48SZcceI0r005l4t7LAUol4HHuSxp24Nzj5roCVFALZeLmIKYQ9cPnd8/vgtVTUDnxZ4IUcwq2ka47q5VTxr6zo4+77XKtuWPb2zEhp3FCX+qbRStFcolXyg2GnHIdCWsYfjGHh0vlzgNvXxHjjtwVGQ7Kbr7VQpDXOOnSCPJLnL83DFDg7smlbnJf7NF4tySPr2m0Je3MWfKSPyfUw9y3it5uRQrLNkEBUcA306kk/MfcBZ5nefdZsG6u0ndUWtKoK/f3gEi4KSDx1Qk0U61CwqLbmVAeoWACzYL0AKAwNNJ9CXMI+okHYr1Q3f/y9ySqaNbMG5YsXFcNXNkpTDU9esXhYnsw+vPwaCLmsrU5F8lhojrzO0/vbFR+ZxK4dEl7wEAdnX2gohw3T+6aSUUhPGabe2B0m0AvMIQ6lW1nP+cVlRVvnzFBVppKViIckkqs2qKctm2rweMOYNbCa222ulzw1ppWqleHU5PP8gqLLh1kg7xBTlJxaIo9zTnGukE1HDDnbhIycoh3tsgX6rptkjJPLieWRUdpV0NvLTOKdCxequTTdWjRRXm1rKNe4vaRGEa4fAWC09D1xTooispaRiqwzUDksqsmtLQMxZh0ohm5YLGh09qxTFuKLSItCmXjMZkrVQ/wpPaF6Tlwbeut158JABouWN5GnqoIyrG2YJdvCA416h+aUHAr5STsQi3XHQEAHnKxadLgg+vytzkY0oCX6vDGx87LZrKqgY8O4D7vlK7WO4Toathc8pFW0MX7onqriuaQ68DyqW3YKO1OYesRUp8bUPGClT+4EjdKJrQE0EXTjhxAg3dHTA+pjruWPGUi/ykd9zTituryaGLizF328xaFpobnLGRnaf8MHF9SuSHnsBtkYfcN2Sr//gfMs5J9HbjBw8FUDmlS3fHwu+vroYuPvM6vyXeD12pG/75eqf1D3gwi6qGbrNiTwggzVwu3MslPQ491rtExigaSu6VREOPo36kKZeI+5q0TJcKRG6XC3RL4EtljZI+5RLUxtS0OZ9yScK18nPSiI9ocP29TzhoDADhdyj2hXPvHLqpELhA1uXQbVsU6Doceowfen9GilYL3J/c0dDVrP9RTh3Vrm4T5dlhUXUpl33dfXh1/c6idt8YWb4vXrIzd5LreLlERVUCan7oYZcu7xqoZjoHVvQ6a1nKW/yo/joLk/zvEO9LkohC3uU0BHohIj7BsSXInR83XjqutYD/bOY0KRcxR76OAunMcf99XXHoPMeDZZFaOk0WneO36pRLhBBzCiFU78Hh4fJhqORQ8dIAh2qjKu2auFE0nMtFITkXYyiKCASqq6GL/tE8+MUJ83Z+mOy9FekSDlW3Rb5DaMxagn1G+nQPaVX2AaIXepVnZNGG3QCA2/62OtCuY1wFfNnAU4foLghiNS2VPoSjXpNGANeWQC8w5LIaGjqiNbkkPKMOojR0oupWUon7JhWNTsz7DfgaulKwF9d8MjGUi7SGXtyuKgiTQNTQfQ7d57BlFzm+gCXh0J92vVMacxnvOjpzK82Uu94zkgk+I7LjsLuzL7Jdt/YAH79GTT/0IIeenHKpq1wufXnOoVuK2mAch653k3XBb24gUKCK2iQQLyhV/L+5lwsXxnocenTov8qkL7AYykXDy+WXL7yNDTs71U5CMHKQC/RMRr22aRyHrvLg3vXsOgCuhm7p04niIrBpd3FAXn8iyliuMg4ThjcBAG6/7OhAu5WUcnHnuqqGH/Ryca+p7Lbov0+af6q2BLqb+D6rSFOEB4Wj2n7odtxkreIWN07AKLkthoyiRLzKvUp+neKxAISFRaInrIQfuoog29XRi288sgJX/PI16XM4xF3JW5sd/2exwrvsPLUjNHQV6ki8r0lLlYnnpFceMeS+KdkPz80wtPPTpaAKIYGuHCnqfp/l5tdR9cAq9kN3/vdnTdGqgRtFMxYp5XYIb1s4qu2HHhVMY1WZcolNcqRgPb//NSeBlEiXZBU9j8qH/stdI/q+qkWKtvc4uWnWuMEsKhBTOu9yt/vi9lr21vpRnn6bimbKfwPAjaLOax2BLAoLneyCSSDWZeXIKOy44uIbPCO14rMWDv1XVb74Do4rLjr5ecK1UXm7DmpMoDMtL5ewLydHtY2idpT2YVWXcunpiy6hFSUY4/D48i0Agr8ja5FXq1QG8Ry6/E6hEGsUlReEPfkCtuyNNhTLnh9G1rKUPRp8Dj3MHcud/+lfLfBeN2atRJGi4rP1X0+vVT4/CXylx29TsYl4tGZMBLKyUTQ0T2969E2180M0q3L0b5GXS/C6qqgpgd7LNfSMmjbIEBdR6PyvmlGURQh0qu62lhe3+K+PHhNo5z3SHQvV2IA4Dh2QDwzauKsTO9p7ij9Q0ILO+8kLuPjOl+QOjsA9z68vahOTY6l6uYijoeKux11Rj5o8HONamyritggADy5sUz4/CQq2Uy4yoJUqUKxicJcIv/xbMqPoRkWbAt/gcMpH1b4TTm9RV7lceGBRpTV0HT7q9Xd3YURzDgeOHSp9TiGCcqm222K3K9CnjmkJtPupfEufL46VeA+yivl14vzQAUeoyQzJ2m0dke0qHPrKLfukjovDnq5irwox9F721nKbQREdp0j6XnXSgQD0OWPnnOry5iJ4/U0RKikQdnY4OfrDybx0M5v6HLr8DjbqfL4RVbXvxFEudeGH3pfnHLrj5SI7MPEcuj7l8qGfvoh/+P4zSufYERq6yra6EuAaulhRHfD9v8v1RUzcJD4c6ho6CySREuFw4PpjUk0/9DlTR8b0wfkvbRSNCP1X5VsBn+tNsvtM1W0xVOUeUBsHbgwd0xLMwqmbdjacbVEVfvSuz6GrKqOVpFxqRkP/+XPr0NFbQE/eRmuzL4hlFs54Dd39vFpui6HtF3+to0XpoqvX+bLmsEDnxsgy54vGt0PGD/Neq+bXidLExL7I3JJZ+7Vi4ojmonbdbIthzwgZ9IXK7k0d7dQ2VXUv44clCSwC/LmVhEMX+3zqjLHK5ydB3i4uH6ei1UZ5CwHJQ//Di4ws1m13dpFckVJdpOP90LW6UzsaOjdGPL58s+D3LCdA4jT0atcUjYqOtKi6ybn+7aGlAKIEuhz9JNYNFavZ62joYZ7T6wvkNOw4d1TdXU9vwcbvF2xQOifsBeILZr+PMogL/Vedm3znpZN7O9yXUS0NGD+sSf0CCeDMi2INXT6Fgn9O4BqagUWvv7MLQHzAUjn84oW3AfgRrKqLdLEfut+ug5oQ6Jf998ve6xs/OEvIHSL/sJTS0KuVy8ULmRcEmYrBp5IYPiQXeO+PRenzuvuiF1FVu0ZfwY7Veojkk3PFuaPKdOWRxZuK2r70P0vKnyjgtqfWAACu+4fpXp8A+QXSAxdEIYO57PlzXernxOmjvXMBPYP7rY+vBODYMnZ29iqfnwRRufpVjMO+V0mw3V/g5Mfjyl++hm/97woAwPod6kFnAPDGu7sBAJ88YarTD0stP48juyI0dE0ttCYE+otrdwBwHvTzjpqoHJkYy6Fb+lqMDgoRk00laKJSOGr/EUVtPode+tzuGLdHLQ09hi9zvFzKX8NmxQ8uwDX88hdQ1cZLYVyro8nyUoeyCySHmCmRQ8UYOKQhi6MmD/cefj63kygrOzp68dcVW7TP10Ekh66g9Hh2qtDzrhP6/9RbW73Xl8zbX/q8KPBU08p+6HbxnADqgHIBgPdPd1Jqehq6pN9zXM6PagcWeX7oRV4uVfl6AMCE1ibMGF/smeMJoDKacbfrd/309acG2rOWpeyHHtbEOGQFcqyGbskJ0sYK5Ps+67AJAIAd7UFNVvXB48eJle5Jwb4STolc7RiLSiHOy0U1QCuc6kPHpnDyIb79YOzQ4lKH5dDZ69ubuLFa2Q89JlJ0wBtFAX/7klEMw928tzsyWX9qNUUDXi7V5dALMeHyHjdXRoBwDv2AUUMC7coaeoQmxmFJauilvJdk7mklCjg05SxMGT3EK/LMocp1RmvocgFWQJRA99sHEgq2HVHvViEnTsQ4Aur56Z1z/NcjhzRg3rRR6OqN3qFGgZckBHyDu6ofet5mGCIY65PKLKkZT0RnEdFKIlpDRDdEfD6TiF4ioh4iul61E6NbnEro3AClyqEDfli2iKrncglFjfE+VLNCvR0TLu9HaJbui2d0Cj0x2YxiLpeIyuxiX2RubWmjaPnzKzHsf1q0Ce/s6CxaHFRjHJLOjXCisiReLhwHjXViFao5P5P6oUfFegA+NaeWpdVHLktobcoqnS+m2eACXdUPPVwuUpXKC6OsQCeiDIA7AJwNYBaAS4loVuiwnQCuA/A9rU64v4I/NKpeLhmLMG10S1F7tSkXnpAqTbfFuMLKvEUmsChKiFaUQ0dyo6jMQyOee+WJ08oeXwofnuNwrP99+ZzAtaUplwg6zrJU/NiD45k0zSrgB25VU8uP9HJRSI/Buxqe4zpG0UDgnGUpLSxA0APKD/1Xcx4ohKjJpAu1jIY+D8Aaxtg6xlgvgPsBnC8ewBjbyhh7DYCW78+oIY6G/uQXTgGg7lNaK7lcbnnsLQDJKZdP3PMqjvzG41p9KJXQCigv0OOEqKqXSz7iweWQNorGpEWWNTzNmtjqvZ42ZkiJI8ujuSGD9TefgzNmjXf74PZRWsOGe16IQ5cc0nzoviZxWwyjmqxNlG1FhXeO9UPX2NWLX5m1SDlgjQt0kYtXNYqGFzh/XslfQ4SMQJ8EQHQXaHPblEFEVxPRAiJasG2bH5HYZ9s458j9MN71JFDNy8BY/IMPVD8yThToqm5MgBOtubc7j+Wb9qBDCPSRgR3DO8tyvk4ocnG7iobe0ZPHso174v3QJakGZ7cRcT7U72mlZ4BqjENUjILK7yjO+eG3q2LM0AaMbmnwdi3VfD4cjTTYpqLV+hx6Mi+Xd3d0BpK2WRYp7ZgA3970z++b7LWRIuWSt+2iYh9A/2roUWqW1rcxxu5ijM1hjM0ZO9Zf1fIFhlxICAJyN8fLYhfxmZ+/pMoCvShfh973n3Pb8zjl1qeVzol6YHg/gPI3Lm5xzFoW3pX01f2XexfinR2dsbQKSRoD43YLJKlJ8XE/6eAxgeOXtO2W+PbSvsCqMQ7RidvU8qFnQ+cCQNsu9QIV41ubcPQBIzC+1fHsqKZAj9LQVXYqcYFFqgL95Fv/jtWhdMqq1cV6Q7V3eT9UaZ+o+9qf+dDbAIhOmpMBFEdsJADPg86hYhRlMTdYbJO1fJ/8n3/HCd/9m9zBJRCsEZjMbXF7VLbBEogTgpDU6OI49K37ugORo6XAMwNS5DIrL8jCiYv88+WEEN9R/OKTcwMPSNgFMQ6in3JxH9Q0qShjngrVEKbSuAD77SvvSp0fvpaYYKyaHLodQcVlLIWF0RvHYLv3WxIsTqrVxbjNTHSPVaVcwkZiP12J/DVEyAj01wAcTETTiKgBwCUAHtb7umIwxvDenm4s3bjHa+MruMxEi3IH41DlOd/d2YlNMUWWZTC6pQEfO+6Aoj5UUwMKF531+5GMQz943LDIMY4CPyyqH/xzmTGJW1xkFoS2XZ1eIeFsxsIHj5rofRbnfRNGKaO8H9EndanoXPmK3h1RD74OGHNLpllqv6ESyEd4P6kUgfEol6JFgf8W/WdN9VnlBVPEnahq6L/j5VI5Y3dZgc4YywO4FsDjAN4E8HvG2HIiuoaIrgEAIppARG0AvgDgRiJqI6LW+Kv64BVh3trspzn1ihJLCXTnf6RSqjk4P3pyldLxHL2hnQbgPDjV9iKIS1kLlNeESnHwsr+Cnx4nc2SvFWvgleA6v/9E8B6OEQJHhjbJhV/wufm3L55S9Jkqh+1RLmGjqGxgUWih5nM7Kiq4HJxF3/fDrjaHXlSWsJKUS4Lfourlcv2DiwEAS9v2Bq6h0oWwN1hSRw6pmc0Ymw9gfqjtTuH1ZjhUjDbOi9CgZPye+epYyiiqeo9/9ORqtRNcOPncQxZ8Re+QMCYOl0+exBhzBXLxZ7IceimPIdlx5FRL1HWcdnnKJTpIqvyD15SL11U+9NMX8Z8XHYkPzy0d7t3pBpmEg4qcPjj/pd3tPKNoUMtWoRrCgnC/4U04ZJx8vn6vL8ylXCogBFWRtxmacmENXcOfP45ySaKhW3ouxqOHNgj9UOtD0c4rod0v9dB/3m/RxczT0CVCzfn50ULI+V9NP/QoDT3J9+cUoh3jyr4BYqRoOQ49zh4hP47l6ACC3ISNW1wI5QXp/qNKuyk+IJHn5R3XCDykoVjvUebQI+IDVPjWqJ1XVrGyFwe/x0kTQekgamFSeUa4sAwrcLr50EXo0qMNRRy6ipdLMCvpgM/lErXiqnm5wD0/XkOvxnx9a/NeFOwIgZ5QQ+/NK0Rncqt7ZBoE57+Mhh6n4cv+jjjuXOyLzJyPX1zKPzStTU62yQevOd5r+8Un53qvZSrU3PmMU28zSkOX3fEATs6PPyzcWCSQVQOLwuOasyytIs/8HqeRDyZfiPJDV3FbdP7HBRbJDEfcAqYqjDnEnqiG/sfZRgZs+lwvUEAYFhU/9FJG0SS+uqo460fPAQBaGoMPv2URFHJaFSEu+2EUeG6JqEIOshplfEIsecqlnNFRtmJR/OJSXpByAXGQUEJQTCnckJXz2AGidzwqD97X/7wcm/d2B0LFAbWHv8CKNdtshqR2sWFwOwn/WVXNNRTJoaskOYuhXNzfIrPbCN8Hvx+klAtm/1FO8ZXjDhzt90OBPvrrii1o78kHxmPA1xTl3RZliI6XS5SLXFKfTh2EXfuylnz63Dff24vb3fzbHF0KAr23hIYu7+UST1/JTjKuLS3ftDfyc4LcA1wq6rVcX7jmKi4uOUEzbJDQ0D9wxASs3tIe+ZlKYNHGGF9xVQ69SEPPWNLpMcLXIlFDr6KKHuflIvs7KmEU7YnZ9aq4TwLAkZNHIJexAvSeih/6p3+9AEDYrqJn9+NIX6C7YyveIBUvl6gFgUN1S9mUs9DdZ2PEkJx2BZOmkObnTFa5Dnz+gUUBbx/AKTjBGIs0+obBhViUsJLVKOO+S2U72l4mulXWKBoX5CTjHeLbE/zzo6qrl4KTeiD6M5XAond2xBW7VnRbLOLQLfRqaOiMOddKI5I6TDEAzn2RTXJYiRJ0PW6K6DFDGwNxHqqUy6NL3itqU6VcAMSE/g9wyiWaQ5fwcolYEDhUB4dvX/cf6a+4qtr9OzuD0ZQZS14DihMycRpFGH0RkWvha5frCovxklGZqOX667gtyu2+oqNey9+XqFTGonAPl+iLAkNMkBbUtNu42AZlo2gRh07Ia3HonHJJ7hmiiijqiBQ8QzxvoSJ7hNzi1N1XwPZ9TmDZ9f90SPAaKu6TsTy8fIwFRyXz3NeOQBd+FL/hMmWh4lZsp01+cAo287Q60dCkOtnHCC5MgGsUTegdIpujuZRAl13ckmY4lIGsUTS+L+UfPL44izRLULiXn/qlvkP2wVvatif2M1KkXMICXZ9Dd/zQK5GCVxWFQlyBC0mBzo2i4d2KpGfcBXe8gA/c5ti7MhZh7LBGL0aBSF752tsdvYOXde8VlZ6/velHJCe1+9WAQHf+i7eH3/Cb//JW2fP5zy6loctMFtGbRBTosnTJsdNGAQAuP35qqA/yHLr4G/7jwiPwzfMOK+pPKfTmOc1QikMvJ9BLeZZIdaOkDzigUrEojnIpP+H57s6KEeIyO4Q4LxsAIG6EK9GPvoKNc29/3nv/oaODOe2UcrlEuD3mMhb6NDh0Pq7+Iq98CW3k7eK0ypaGUTR8W2QXJ5HSzGUsPHLt+/HwtScCUHOf5JTsTRccHuyH5G5DNMyu3ur3qRq5XPoVUYFBcVn6olBKQ1fZFosCneeJBuQ1dCJgnivURShp6MIDe9T+wz0Da5xVPgyPQ8/GG4jLPf/xgUXyWsMHj3SCxH74kaMiP5fxUomrTAO4HHyZ86PS9wYEohSHzyJrmjr9Km+8Cs+dg0JBQCpj6tSeLNZKdTR0nlJBNU11JRDHocuOw4J33DxBYaOoxm/JWIQJw5swcYTjraKSd4k/a8Obg8XYZYLeAKBHKMY+tNG/xoCnXKL8yMV0knu6Shsn/cCiCCFmyQ8ON5SEIes2GOWOBaiF/ounD2nIepyv7EMblf2NQ3YrF6eVksJOI2sRxg1rxIVHRwcPy/DxcZVpnDYJDj1Co40ruBEHm7ESCcb8Y0qdHzynmGpI4raYy+j5ofPqR16kqIL0uOf5t/GLF95W/k6OcCANwN0F5frwwpodke06of/R7pNy5/PEXOF4BkuSThSVNHEu16dRVJj4979aOptcKU1OZXDiDHnlPDY4okprATycWJ1yaWnIeIJZ9qHty5cX6OV6UiqwSF57KQ6wCvQFcgIZiPZpl+nLz55ZV3RPcwo7P8AZqzi7howmFRZS4bVJVptjjMW6LWoFFrnuoDouct/63xX45iMrsGbrvvIHRyBKQ5cVgoDv+x2GjgtmVHCStAOFzctlFgdJyWnovqI4rtVP76GaZz+M1N0WPQ5dGFtxnLfsLZ0+1ufg42kGVSOFTHsYJTV0DaNoc0PGmyx9FdDQk3LoGYWoRicdcrw2TBKUC9+VRI2pboh2RllDR7zfondMfD/CAn2/UF4eWUEWawhMGPqvElj0/lueCtCSf160CV/8pxnK35237eISdAqCdO6UUXAqXgbhe+zI9yX8nKgqLUDxrk+WPvrj6xu91/dddazQB+e/bmxA6hq6V6BCmKxNQnBO+CEoOh+lNHR5Kz6frIfuF0wSKUu5RIU0A26kqOQkE8dgSEPW48KlNfSI/MxePyRX/lLJuVS0l5IauoRXA98ZdUZ4+KjQPyJEDV3m7Lj0vYDcAikK9F9eMTeQgA6Q3+Lz64SFR9ay9L1cSO35aNvVha37egLvdRC107AUOPRwsWzvGuR/LouonYK8i3O08iTr3nv7350AwkvnHYBRLWJyL/l7EoUaEOjOf/EmtTbl8J8XHwkAOGTCsJLnR2n4HD5vXL4fnENvDnloJNbQFTTbV9/2NY+MRZ6GLhtFVwm3xbjAIpVUr715hmwJgW5ReeMsX2AnRCzoMm6PLQ0Z/OPMcYE28QGW8SJgLD4FsO9BFX8+F8Q3nnMoTp0xrmhcZRfJOHtCLkNalEtHbwHtPflEkaITR8hnARURZaxW8v9mxTsVwPfaUfktUWl8GZObG6vcXOjFv0XNQ+XcI/cLnT/AjaJxuRkOdj0Cyt0g/nnSmqJcgIQz64nW6FII1wbkUDGKHhDKEMgF8zcfWSG1U+B53KPojqSBRUrai22XDK2XyeXCqYSoACAZQdiUyxQtBuLDJ6WhywQWlRLobh9bGqOZTVltzrcnBNud0H+1J//N95x0DPe/tkE6sCjq8+k6aXtt5hTXSKAZx+0gATWPMqB4DvgZG8uf+7U/LQOAIsVFpVgHAJwwfUzgPXeHHbBui2VzM2gYFP029zukwoEdwd0UEiBJvVxUjKIjhwRdoLhgXtK2R6rU2Co378jolsaiz2R98uPyp6g8LH0Fu6SGLpMEifuRRxtFSwtkxhh2dPRi0YbdwfMCGnrp7wdKh/7LeA1xOiQuWZnKPXGOL+bQ+xSycfbkCzj7x05QzfmzJ0prg1HPQHu3WvFywF+YknDodsz89K5T4sd09gb7/OzqbYH3Oh4mUV4uSSpAJaVcasAoyjXsYLuXsF4ishGIKRKtsH3xNfSgQJelXNbv6MT6HZ348SVHB9pVjKK7XRdNXrxXpE7K0S6bhfDy4aGFAZCbKL15G0+s2BL5mSVsR8vlQenLs5JGUZkkSH0ljaKlBQBPaBaXHAyQ09BfWhftIgc4ixKV2V7zPkZRBEBwfpay13ocesLAoiVC1OqE4U3SAixq/u6T9P4KXMdLxxBhjJQt5RfhjspRLlV1eLd9ydxguUidTIfRHHr58983ZWTM7hNuH6S7EDxf77TKIcooCsjXCPQ4+IhfIut7DfgeImGBrpK+Ngoq+dA5zfTkF04BoOY3/fya7SU/lwksKuXzr7Jj6itjFJXZlkblYuEoR1Xwzz5x/JTYYyqRxqAc9xtnzPTPd/6Xj3qNHgvVwCLx/MkjmqXvadQzqKOh82cpineWFaK7OvuwL+a7yylPYeVsRCgoyJc5Ul0BoO+HHuW+6Zw/wI2ia7Y6UZlxUX0yXhlAtNuijBbFwUOCw5SLrIYOAJcde0BRm4oFv7O3gPdNGYlhbnEGMa95uUtEFWEQobK4RUElcKMvoraqCKlcLO4BcQbekt4lrpCbMrqlZB+SopwgKhUcBchrhFH5jgCHv83bTHpxGjXE96Y48/AJ0uHyUfeqQ1FDZ4zhIz97GQCwvSPoiqxCubz69k5s3B2TjrgMvRkOHgynuk4a1ONcQ97QHbX79OJFNLuQOuVyzW8WAoiOogNkKBfnf9zzKWtB5xXiw+5YvLJ3yT64XzB+WLHlX8Uo2tlbwDChgLHIQ5cbh5HuwxoXeCEzUUpNRNnUAQCwbONeTB0d/8BLRXqW4NAdd7/4c72gjxI7HJl7kssQ/umwCbGfl9splAqOAuTjJGIpF/d93i5NcYX7Azj8vqw2GDVWqpTL/a9twMotjtIU1rCdyk1Kl4uEVWZeiMrZIeOHFilvOtpx+FBZBS42EDGB5xFQAxo6xxMrNgfey1IunPuO8r0G1INQwqv4PRJhzvkSW2ueH0JGi+rszaNF8LIRH9JCma01d1/74YdnR37ul0yLv04pFzjZIBSeGbJUpkyZkmOlA4tKPzTlNGNAzrd/eHOuaFse7Efp+yprFJUVqEVui+6cl6Fd+go2vvPom977US0N0mXbovqnSkW+K6SVzhVRLvKOA5NGNBf583OUM9yLfX7flOK8SzpRmuFjZROuFWw7Wl4MdLdFjo8eG+Q7vUCBMr9sZ4eT23jkkIbIz1UTzn/1nEOL2spNtlJ8b0bhBnX2FgLUSdAoKifQ46gOmYlSSjB4O6YSF3hl3Q5cfs8rAIoDtALXktBi8iXGtFyBjDhvisAxZcbz/z24GNvbe0suCuUWJhWjaCmU4tCB8snb7n1pPQ7+6l/w1FtOmtY7LjsGTbmMZ3f69K8XYE+Jgi5R9yoq4EsW4aupFvqIy+ZZbmEQNfSoPPJcvpZ73h96o817feDYIK0nq0A6Gnqy0oZRqBmBPmlEkCqQpVzae5yJ2BqjSak6+k8b04I7P3YMfvbx93ltcUYYjlICxNNsS0wSxhjWbWtH264urBYoHlE4/9ilhOJQXqA7/0u62ZXgU2R2TP963+t4bf0uAMAlc/ePPa6ca5dtM+/3RLk/lnL3W7+9A8d/96lAn6NQboF8cKHz0Jai3LJW6dB7TxDH0CFU4ndw7Gjvwanfe9q5ToSXCxAtnER87c/LA+/5kIqL1fxlxdV3OKLmbhKBHkY5qkRE3o53iVXxcpkyekjR57I2hc8/sBgA8I1zZ0WmD5CxM8Vz6JyG0xPoqXLo4uCHg0BkKZfuvmj/cQ7VslINGQtnHb5f+QMFcDok2n+7fH3UGTc+5mlZSzf6rmUqXi7czS8qdS4gt53k+dSj8MjiTQCA9Ts6MLIlZjckvC4lTDNW6Uo7B/7bfO91HOUCRLv7Pfmm73ZZWqDLGbtLuS425jIljeb8syQa+hJhPvSGvkummHqUYOBzQabIBxC9+Cbx/gp3SdbVD3DGIKoIOuCM5zs7OmNda0U6ddLIYluTKuXCnReC15CzM+UjCn1wqCxwRefqnVYZ8Al6w9kziz6T9XLhPsdxJcVUwoqB6IjTcr6+Uka4Uv7fgnA74SC/gng4O2Apz4JNruW/nIZeauUvlVny9Xd3AwDufGZt5DU6e/PY4dJfTj9KUxWyD3A4h7hzvvO/3DVKCfQX1uxInAe8IWPFpl0GgJseXQEA2NXZG/m5bGARx+PLgzECOS95W/z8/E1EQNrTKx3qRbQ7lVIdwnO3IWsFfNplsCnGMwVwE1pJ3ouo1AEcG3d34dX1O3FfTBAeV/5uPOdQXDB7UtHnPj0q15c4o6bM/YzT0Pl1+5VyIaKziGglEa0hohsiPicius39fAkRHSNz3d8v2AAgWnuVpVy4ES5OoJOiUTQKUcVgRfzhdWd7HnWDf/bsWgDAc6u2FX0WxidPmIpfXDHXex8WimJBWxGvrd+J77rVneKEGDegvRiTTxrw/dB/K2R/4zj9UCcvyuPLt3h0hIiv/HFp4H0UP8hBRJB1n26JcMeUdfcrJ7DvfGatXCdi0JizSmroyzY6gU17Yyg7lXz9AHDU5OGB9zl3N1bK9vH6O7uK2rg7Z0OMI0EY4XHmihhPIyCDUik0VDTSvoLtzeU4vP5u8W8GfA399EPHx6QKcf5LKxtjo5UNmd8S5+UC6BWa9r6/3AFElAFwB4CzAcwCcCkRzQoddjaAg92/qwH8l8yX//vDDrd36+Mriz6TpVy4ht4YYyjJWOWtznxFver900r2Mw7/Md8VphGThJeqivOdBYB/ft9kjG9txDfOOwyNWV+AhW/4zo5oTe8jP3vJex2XaneYm0/kseWbIz8HgO/+xfGCiNqOXnmiPzZrtxXzym+9F8yPXTJStIRdoziHuL6vrki/RGGdUJlKhKy2uG5bR9nFvhRUIxM/dlzQcUAmedsbEcLtUjdCUp5yie5f3M5DFbI7NsaYk2u/xM6rFPjiGycrVD1MDhpXHOcg+1t684WSu+n+5NDnAVjDGFsHAER0P4DzAawQjjkfwK+Z04uXiWgEEe3HGIud7au27AN3HLrw6Pjtz0+eWoPfvPxObOd2dPSCqJTbIuGhNzbihRKRlJxqKKWxnPa9pyO3SOKwR90gJyMew02PvokHXtsQee3Ne7oxamgxLx0WZtf+9o3IACJxAu4fIYzD1zrjB89EXIN5Aq41ghsUH4KfPbMOTwmFbQEEjLlAabrDIsJbm/dF9kNmGvOyX++/5e9F+W/EfoiLEMe1p033Upf+4fU2LGnbrdUHEVG/Q0Tc1po3f/jOl2LHSzTIh20XfNG88pcLYud/2H309EPHe6khxPl+wx+X4u7n3468Rtwu5Iu/X4yhMYnHwogz9gPO3OzJ22XHkd+XUnmCACfXeFRxbp5aQ1Sagv1w/n/0v18uudg15SzMmTKqKJEf4Oy6dnT0lv0te7vzsWNnEeGBBRs8ryQRH40IXhQhczcmARAlURuA8J486phJAAICnYiuhqPBo3XigV77LRcdWfSlrc1ZXH3ygWjbFe/PDDhbghnjW2Pzi/zrKQfhjQ3RWzARx2asIg3oOxcejq8+5GRVO3S/+DS+GSK0NGZwYihzGgA88flTcNr3nsbZh0+IDX46ePxQnHBQ8bkAcP0/HYJxrU14491dsaH53fkCdnX04evnzio52S+ddwAeW/YeDh4fnSlve3sv5kwZiRERuWAO3a8VH5mzPx5YsAFnHja+SACNGdroGRA/dPSkyPqqYj/itCTA9yz5w7+eEPn56YeOx4L1uyL56ymjW/Dkm1swblgj3jdlZNHn1585A7P3H4Grfr0Apx86LnYR533gaRii8C+nHIhfvrA+djwnDG/Cc6u340PHFCssAHDS9LE4f/bEsj7x85duxudPP6Sofc7UUfjQMZNKGih3dvRiR0cvxgxtwPz/exLGCcFvLQ0ZfPTYA3DfK+/ijFnjS+6qxg1rxMbdXTj5YKfPl/38FRx9wIiS/RYxbUwLnlixBRe/bzI+f8bBgc/OOWI/bNjVKaWVzpwwDGfGBHvx5/XE6aOLan1yTBzeXKQEcJw4fQwumD2xrBvoweOH4ryjou/pBbMnYU9nX9lsojMmDMMHjoh2vvjMadOxfFO0jWLU0OLEeyKo3CAS0T8DOJMxdpX7/uMA5jHGPisc8yiA7zLGnnff/w3AlxhjC+OuO2fOHLZgwYKS321gYGBgEAQRLWSMzYn6TIZEawMgOhVPBrBJ4xgDAwMDg36EjEB/DcDBRDSNiBoAXALg4dAxDwO43PV2OQ7AnlL8uYGBgYFB5VGWQ2eM5YnoWgCPA8gAuIcxtpyIrnE/vxPAfAAfALAGQCeAK/qvywYGBgYGUZAyUTPG5sMR2mLbncJrBuAzle2agYGBgYEKaiaXi4GBgYFBMhiBbmBgYFAnMALdwMDAoE5gBLqBgYFBnaBsYFG/fTHRPgBiEpfhANRSuBWjXq4xBkDpqs/934dauUYl+lAv41kLfQBqYzxrZSzSGM8ZjLHo0HXGWCp/ABaE3t9VgWvWxTXCYzNQf0eFxqISfaiL8ayFPtTKeNbQWFR9PEsdX0uUyyPmGhVFrfyOpNeohbEE6mcs6mU8a2UsamU8AaRLuSxgMfkIBjvM2FQWZjwrCzOelYXqeJY6Pk0N/a4Uv7vWYcamsjDjWVmY8awsVMcz9vjUNHQDAwMDg8qiljh0AwMDA4MEMALdwMDAoE5gBHoVQET7E9HfiehNIlpORP/XbR9FRH8lotXu/5HCOV9xi26vJKIzhfan3bZF7t+4NH5TmlAdTyIa7R7fTkS3h671PiJa6o71bRRX+qqOUeHxNPNTfTzPIKKF7jxcSET/IFxLbX4m9aE0f1J+o/sBOMZ9PQzAKjgFt/8TwA1u+w0AbnFfzwKwGEAjgGkA1gLIuJ89DWBO2r9pgI1nC4D3A7gGwO2ha70K4HgABOAvAM5O+/cN8PE081N9PI8GMNF9fTiAjcK1lOan0dCrAMbYe4yx193X+wC8Cafm6vkAfuUe9isAF7ivzwdwP2OshzH2Npw88/Oq2ukahup4MsY6mFMesVu8DhHtB6CVMfYSc56eX8O/B4MGlRpPAwca4/kGY4xXeFsOoImIGnXmpxHoVQYRTYWzIr8CYDxzKzu5//n2NK7oNscv3O3s1wYjRSBCcjzjMAnO2HKEx3nQIeF4cpj56UJjPC8C8AZjrAca89MI9CqCiIYC+AOAzzHG9pY6NKKN+5d+lDF2BICT3L+PV7aXAwcK4xl7iYi2QevHW4HxBMz89KA6nkR0GIBbAPwLb4o4rOT8NAK9SiCiHJybex9j7I9u8xZ3W8W3/1vd9tii24yxje7/fQB+i0FKxSiOZxza4Iwtx6Atbl6h8TTz04XqeBLRZAAPAbicMbbWbVaen0agVwHutvNuAG8yxn4gfPQwgE+4rz8B4M9C+yUujzYNwMEAXiWiLBGNca+ZA/BBAMuq8RtqCRrjGQl327uPiI5zr3l5uXPqEZUaTzM/HaiOJxGNAPAogK8wxl7gB2vNz7QtwoPhD45HAAOwBMAi9+8DAEYD+BuA1e7/UcI5X4Xj3bISrmUbjnfBQvc6ywH8GK73y2D60xzP9QB2AmiHo/nMctvnwBE6awHcDjd6ejD9VWo8zfzUG08ANwLoEI5dBGCczvw0of8GBgYGdQJDuRgYGBjUCYxANzAwMKgTGIFuYGBgUCcwAt3AwMCgTmAEuoGBgUGdwAh0AwMDgzqBEegGBgYGdYL/DyY75jYUhCgfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "stn_swe.plot(y='398_OR_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": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "CSO_gdf.to_file('ClackamasLake_Snotel_SWE', driver='GeoJSON')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7304"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(stn_swe)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "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
}