NWM_A&A_retrospective.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Purpose:\n",
    " - Download retrospective hourly analysis & assimilation NWM files from the Google Cloud archive\n",
    " - Obtain and plot the discharge time series of a stream of interest\n",
    " - Run the cells in order to import libraries, select parameters, define functions, download files, and plot a discharge series\n",
    "    \n",
    "### Inputs:\n",
    " - (1) begin date for data retrieval\n",
    " - (2) end date for data retrieval (inclusive)\n",
    " - (3) desired product (channel routing or forcing)\n",
    " - (4) directory for downloaded files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import python libraries\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import pandas as pd\n",
    "import requests\n",
    "import sys\n",
    "import xarray as xr\n",
    "\n",
    "from matplotlib.dates import DateFormatter\n",
    "from time import sleep\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# select input parameters\n",
    "\n",
    "begindate = \"20190915\"  # data is avail. beginning 20180917\n",
    "enddate = \"20190930\"  # inclusive\n",
    "product = \"channel_rt\"  # 'channel_rt' or 'forcing'\n",
    "destfolder_name = \"test\"  # directory will be created in current working directory if it does not already exist\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define functions\n",
    "\n",
    "\n",
    "def get_netcdf(filetype, begin_date, end_date, output_folder_name):\n",
    "\n",
    "    output_folder = os.path.join(os.getcwd(), output_folder_name)\n",
    "\n",
    "    if not os.path.exists(output_folder):\n",
    "        os.mkdir(output_folder)\n",
    "\n",
    "    if filetype == \"channel_rt\":\n",
    "        prodstr = \"\"\n",
    "    elif filetype == \"forcing\":\n",
    "        prodstr = \"forcing_\"\n",
    "    else:\n",
    "        print(\"Product error. Choose 'channel_rt' or 'forcing'.\")\n",
    "        sys.exit()\n",
    "\n",
    "    t1 = pd.to_datetime(begin_date)\n",
    "    t2 = pd.to_datetime(end_date)\n",
    "    dates = pd.date_range(t1, t2)\n",
    "\n",
    "    for i in range(len(dates)):\n",
    "        date = dates[i]\n",
    "        save_dir = date.strftime(\"%Y%m%d\")\n",
    "\n",
    "        if not os.path.exists(output_folder + \"/\" + save_dir):\n",
    "            os.mkdir(output_folder + \"/\" + save_dir)\n",
    "\n",
    "        for hr in range(24):\n",
    "\n",
    "            url = (\n",
    "                f\"https://storage.googleapis.com/national-water-model/\"\n",
    "                f'nwm.{date.strftime(\"%Y%m%d\")}/{prodstr}analysis_assim/'\n",
    "                f\"nwm.t{str(hr).zfill(2)}z.analysis_assim.{filetype}.tm00.conus.nc\"\n",
    "            )\n",
    "\n",
    "            filename = os.path.basename(url)\n",
    "            write_file = os.path.join(output_folder, save_dir, filename)\n",
    "\n",
    "            for attempt in range(30):\n",
    "                try:\n",
    "                    r = requests.get(url)\n",
    "                    with open(write_file, \"wb\") as f:\n",
    "                        f.write(r.content)\n",
    "                    break\n",
    "                except Exception as ex:\n",
    "                    if attempt != max_attempts - 1:\n",
    "                        sleep(0.5)  # Give NOAA time to wake up\n",
    "                    else:\n",
    "                        m = \"Could not download file.\\n\" + str(ex)\n",
    "                        raise Exception(m)\n",
    "\n",
    "        print(f\"{save_dir} done\")\n",
    "\n",
    "\n",
    "def get_series(comid, begin_date, end_date, datafolder_name, filetype):\n",
    "\n",
    "    data_folder = os.path.join(os.getcwd(), datafolder_name)\n",
    "\n",
    "    t1 = pd.to_datetime(begin_date)\n",
    "    t2 = pd.to_datetime(end_date)\n",
    "    dates = pd.date_range(t1, t2)\n",
    "\n",
    "    df = pd.DataFrame()\n",
    "\n",
    "    for i in range(len(dates)):\n",
    "        date = dates[i]\n",
    "        date_dir = date.strftime(\"%Y%m%d\")\n",
    "\n",
    "        for hr in range(24):\n",
    "            filename = f\"nwm.t{str(hr).zfill(2)}z.analysis_assim.{filetype}.tm00.conus.nc\"\n",
    "            nc_file = os.path.join(data_folder, date_dir, filename)\n",
    "            data = xr.open_dataset(nc_file)\n",
    "            Q = float(data.sel(feature_id=comid).streamflow.values)\n",
    "            # time given in \"minutes since 1970-01-01 00:00:00 UTC\"\n",
    "            # below line returns a numpy.datetime64 object in UTC time\n",
    "            timestamp = data.sel(feature_id=comid).time.values[0]\n",
    "            df.loc[timestamp, \"Q\"] = Q\n",
    "\n",
    "    return df\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20190915 done\n",
      "20190916 done\n",
      "20190917 done\n",
      "20190918 done\n",
      "20190919 done\n",
      "20190920 done\n",
      "20190921 done\n",
      "20190922 done\n",
      "20190923 done\n",
      "20190924 done\n",
      "20190925 done\n",
      "20190926 done\n",
      "20190927 done\n",
      "20190928 done\n",
      "20190929 done\n",
      "20190930 done\n"
     ]
    }
   ],
   "source": [
    "# download files\n",
    "\n",
    "get_netcdf(product, begindate, enddate, destfolder_name)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract a discharge time series for an individual stream from the downloaded files\n",
    "# this example shows a reach of the Sabine River in Southeast Texas during TS Imelda\n",
    "# compare: https://nwis.waterdata.usgs.gov/nwis/uv?cb_00060=on&cb_00065=on&format=gif_default&site_no=08030500&period=&begin_date=2019-09-15&end_date=2019-09-30\n",
    "\n",
    "comid = 8331880  # Sabine river near USGS gage 08030500\n",
    "sabine_timeseries_df = get_series(\n",
    "    comid, begindate, enddate, destfolder_name, product\n",
    ")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnAAAAGPCAYAAADcP+3yAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABqBUlEQVR4nO3dd3gc1dXH8e9RlyzLTZZ774Vq0wwEmWJKAgRIoZcQSkhIIBDeEJJgSCNACIFAAqETwJDQOxgwhI5NMS5gY9x7ldXrff+4I7NeVFbSSrur/X2eZx9JM3dnztXujo5uG3POISIiIiKJIyXWAYiIiIhIyyiBExEREUkwSuBEREREEowSOBEREZEEowROREREJMEogRMRERFJMErgRCRiZtbFzP5pZmvMzJnZPbGOSUQkGSmBS1JmVmBm15jZfDMrNbMyM5sbbOvXQPm9zew/ZrbOzCrNbJWZ3WtmYxsoWxj8cXdmdlYj538k2F8Rtn2WmX0Wtm1ZyPHqzGybmX1qZreb2T4trPe04BybzKzIzD40s/PMLCWs3FVm9l5QrsLMFpvZ38ysT1i5bmZ2t5ktCI5XEvweLzazjLCyucFxnzOzjUF9fhmFWFPN7AIz+9jMis1sg5nNNLNDGznuKUGMFWa2PIgpPcJf4aXAecDdwGnAbRE+L2bMbJfg/bYsqPMaM3vdzKbHKJ6I3wdmNj3kvR/+SIvgXAeb2R1m9lnwGV9uZg+b2ehGyo80syeD99v24PsRDZT7vpndFxy3LvwzG1Z2uJnNMLPNZlZuZh+Y2bebiz147sDgd7B7JOVbwswyzOynweeqvr6fm7+u7Rvt80UYU0TXnZDyrfosm1lPM7skuL6sD+r+oZn9yMxSGyhvZvYzM1tk/vq/KPjdWVi5fub/hrwS/E6dmZ3YRBznm9m84JjrzOw2M+sZye9KAOecHkn2ACYDG4AK/B/iH+H/KN8CbAYWhZU/H6gFlgFXAmcDfwTWB8c4Jqx8IeCAcmBmA+fPC/aVAxVh+2YBn4VtWwZ8ApwaPH4E3ASsDc7z1wjrfXxQ/j3gp8AFwMvBtmvCyj4X/D4uAn4IXA8UAV8CuSHlBgFvAdcExzsfuA+oAx4JO+bQ4FwrgZeC738ZhVj/FmyfEbyOvwA+D2I4MqzsGUHZZ4Bzgt9jLXBXhL/DWcDHsX4Pt+C9vj9QGbyHpgev5W+D17cmRjG15H0wPdj/45D3f/3DIjjX7OA9ey3+c/sbYB1QCuwRVrZv8JlaAfw8eKwA1gB9GngflACvB8/5rJHzD8RfJzYGv/efAG8GdfpeBPHvG5Q9sx1eh2eC9/7DwIVBbDcHr8v0GL03IrruBGVb/VkGvgVUA08CF+OvW48Gx3uogfJXBvvuC+K6P/j5t2HlCoPtS4L3hgNObCSGPwb7n8Rf0/8MlAFzgMxY/P4T7RHzAPTo4BccugOrgovqhAb2dwP+FPLzfkANPkkJv4D0ARYBxcCwkO31H+L/BheU/mHPOwv/R/UpIk/gXmgg1mzg8eBcF0RQ91fxf4yyQralBBeM9RE8/7jgXKdGULY+qRoYsi2z/nfBV3/EG/vDHVGsQAY+iX4s7Pm9gSrg4bDzbyAsqearJGGPCOq1AJgVQblUICMO3u/PAluA/Ab29e/oeFrxPqh/bfq28lzfAFLCtg0P/lCGv2f+FnwuR4RsGxW8j24MKzsYSA2+/9pnNqTcLfjrx8Sw98bs4P2d3kz87ZLAAXsFx726gX0prf19t9P75WvXnbZ+loFhwJAGtv8jeP6uIdv6BdeYO8PK3hNs7xuyrSvQK/i+kEYSuOCY1aHXp2D7scFzfhzr33siPNSFmnzOAwYAlzjn5ofvdM4VOecuD9n0W8CA051zJWFl1wfHy8W3+oR7Bv9f+klh208Bnsf/YW0151w5vhtvC/Cr8Ob8BuQBW5xzO7ptnXN1+BaJsghOuSr42q01ZZ1zlc65NRE8tyWxZuAv5uHH3Yy/uIaWnYpP7P4eVvYWghaRxoKxoFscGAccFNKNV2hmQ+u7Ac135S7GJwJTgudONLOnzHd9l5nZO2Z2ZEPHN7OTzOxXZrbSfHf0E2bWy8zSzOxPZrbWfJf/DDPLbeZ3CDACmO+c2xS+o6HXwswOM7PXzHdFlwbfTwkrU9+tOd7MHgjqtdV8V2VecwG18H0QclrLs7Du8wjO9Ubwvgnd9iXwETA+rPj3gOedc0tCyi4GXgS+H3aMFc652ghCOBD41Dk3L+S5tfjW4n74P/INMrNC4J3gx7tD3nPTQ8ocYGavBu+VYjN72SIbVlHfLfx6+A7nXJ1zbl1YLHlmdr35bvhKM1sRdBVmhpVz5seIft/8sIqKoIvzmxHE1JiGrjut/iwDOOeWOueWN7DrseBr6HvjWPw15uawsjcH248NOW6xc25zU+cO7AekAQ+GxfUk/m/GyREcI+kpgUs+x+D/sD/SXEEzywEOAf4XelEP5Zx7DVgOHN3A7gp8s/wpIcfsh7/4PNDiyBs+fwm+FW4AX/+DFO41YEKQCIwys2Fm9nPgcHxz/k7MLMXM8oNxHQfhL461+Nax8LKZQdnBZvZdfEK7DN9C2RoRxRrUfzZwlpmdHpx/AnAXPvH+W8gx9wy+vh96IufcRnwXzZ40biE+WV4LfBZ8f1qwvd6pwGX4bvmLgbXmx1q9hW9J+Su+C68r8IyZHdfAeS7D/0G4Dj++7ujgeDcDewO/B/6NTyh+30S89ZYBu5vZbs0VNLOTgBfwr/Gvg0c+8GojScEMoCdwBb61+ezga3tYhO9KKzE/hm1gG4/XD9iR1JpZf3wX6vsNlH0P6GsNjI2NQAYN/3NUv21SE89diG9RAridr95zjwGY2TeAV4Ah+PfCn4AxwOsRJHHLgq8nWzNjxswsG9/KeDZfdbc+A1wS/Bxuf3xL1iP490Ym8EQQb7MivO605bPclP7B19B/ePbE/0M2N6zsR/jW2dacq358cEPvjXJgj5b+s5KUYt0EqEfHPvCtVR9HWHY3/H9zf2um3JNBua7Bz4XBzyfiE0AHjA/2XYL/Q5RF0AQfdqxZRNiFGrL/ouAcxzQTZy7wH/zF0AWPSuCsRsqPDCnngjhOaKTsD8PKvkcDXdQh5YfSdNdZxLECo4EPw86/BpgcVu7vwb6UBo7xBr6lpLn3xGeEdaGG1KUE6Be277/4rpKxIdvy8En/8vpYQt4zCwnpesUnoi54X6SEbH8O2E4z48CAg4PfYW3wmlwPfJOQrumgXBf8Z+OesO05wFLglZBt04OYng89P3B1sH1aCz6Pzb0PfoZPXk8BTsAnwVXB765XpOcJO+bJwTkvDNk2OdjW0PvrB8G+SY0cbxaNd6E+GbxO3cO2Pxwc8+ZmYm20CxX/j8tmoHfItgH4IR1vNnNcwydEDt8V+XDwux7ZQNlf4RONcWHbLwief2DItvrP35SQbb2C99bbEb4+zV53iMJnuYHnZQDzgNWEjEHDJ6vLG3nOCuDpRvYV0ngX6h400IWNb+Gvr3er3t/J9FCGm3zy8Be4SHQNvm5vplz98RrqPnoNn0zUt8Kdgh97U9FA2daq79rt2mQpnwB9gR8XdRq+FedZ4HYz+04D5VcDh+FbLa8KztO9kWM/F5T9HnBHsK3Z7rQoxVqMn7RwB/6P/A/wg8afsZ1nCWcD1S6sSy1QEexviyecc2vrfzA/m+0I/AV+xyxF59x24J/4cVS7hB3jPudcVcjP7wVf7wmL+z386927qYCcc6/iu/GexrfQXoL/g7Tedp4hfRjQA3ggaP3IN7N8fAI3EziwgZaav7vgr07gpuDrUU3F1BLOub855y50zj3gnHvUOXcxvqVzML6Vs0XMbAxwK77l5h8hu+pf+8oGnlYRVqYlbsG/Tv81s73MbISZ/Qb4dhuOiZn1xbfe3et8qxMAzrnV+G65KWbWo7HnB6/bt/Atwpvxn9sbgcVm9kJYa+P38K3IG8PeGzOD/QeHHX62c+7tkHNtDmLar6mYQkRy3WmPz/JNwATgR8650PdBNg2/L1p9LufcR8DbwKXmZ9YPDVobH8b/w1d/XmlCs9PQpdPZTvOJTr2mErNQ9cf7WmLonKszsxn4rooH8P95XRbh+SNVPxaqucT0Fnz3xm7OuZpg2yNm9hLwDzN7xu085qycry7ST5vZLOA1M1vvnHsm9MDOj2mqH9f0n2CczktmNsqFjaeJUESxml9KYia+xeG8+ieb2WPAYnxXZH33djmQbmYWlniAbxEtb0WcocK72XvjW7YaWmJiQfB1KH6Gcb0VYeWKgq8rG9neA9+C0qjgj+m3g4RyV/wf7kuAu8xseZDk1S+r8VITh+qOT4zrLQ47zyYz24rv0ms3zrlHzOxa4FB8N29Egm7SF4BtwHEh7yv46rXPDH8e/r0RWqYlsb5kZufj34f13X1r8a1d/yDyfybDDQ2+NvbeMnySu7WJ2MrwXa+/N7Pe+Akf5+GHKcwADgqKjsb3Rmxs6DhAQdjPixsoUz+UosmYgrgiue5E9bNsZr/C1/3XzrmnwnaX0/D7olXnCvEd4CH8P3PgW93ux//jehytf28kDbXAJZ+FwBgLW6OsEYvxM8iaGz+0G7AqaFlpyAP4C+7f8Rfv1yILNWITg69fNFYgqO9ZwDNhf7gAnsCPdframnahnHOz8BfxU5oqF3gEn1geG0HZtsT6DXzL0hNhsRbhu7YODNlc3zrW0Fim/nx9IkRLteRC3tiEk8YGxje2vbmJKzs452qdcx85536HX6YFfGsWfHUtPBPf+tHQo4idhf/hbFE8bbQCP/4uImbWC78MTQ6+izf8tW7uvQGtfH84527Dz1jfN3gMxXdLQ+vHiDalxa+Bc25j0MI5Dd8F+Q0zGxzsTsF3tzb2vggf3B/V90Uj152ofZbN7EfAH/Azjf/QQJG1QJ/wMWnBz31acq5Qzrm1zrlC/ISSg/CzYs/A12lDcA2TJqgFLvk8hZ8d+F2amUjgnCszs1eBQ8xsuPOz13YSNHsPYefumPDjfGh+oc+p+DXbIpm9FhHzMxGPw7fQLGyiaC/8+72h93xa2NemZBLZLNT65v9IyoZrSax9wn4OLxu6fU7wdW9CEr6gO2g4fsxdNG3ErzfWUGJcv21ZlM8ZqfrWoPrkpL71cKNzbmYD5RsympDWluD32B0/Pq29DePrrZUNMj8z9kX8mmyFzrmvJU3OudVmth7/3gi3D7AutHu8pYJWpfrucMxsWvBtc7/rhpIh+Op909h7yxHh76cBH+D/MeofHGMJfnxvS94X4UYFX1sbU/h1JyqfZTM7Dd/afz9+3b+GzMGP8d0V+Dhk+x74cXMfRnKuxgR/V74M4umJ7xpvdpKdqAUuGd2GH2PxFzMbF77T/HT50BmZV+P/e7zXzLqEle0dHK8U30XSlJ/jx3M0mui1VDA77H58S8QfGuhKCLUB33X07eB59cdIwS9zUkHQrRf8Dr7WZWBm38d3J88O2dbYGKwfBl9nN7K/KRHHih/7BmHT7oMxPAfz1YUefIvcJvyisKF+gn+No5rABYn6C8C3LGTlfzPriu+uWQF8Gs1zhjN/J4KGrnP149Tqu+DquxZ/3chr39Dr/BOznZau+Wnw9flWhvs11sCq9Gb2A3wy9kIEz8/Gj/mbgJ/k81ETxf8DHGkhd14ws1H4LsWovTeCcXjnAk86v0xJU0qDrzuNHQuGJcwBTg+Slvpj98e3VL3tnGu0q9L8zO6vdXUHrd+H4Ft862ObAexlZsc0UD7bvr6czWQz2y+kTC/85/PdZmKK+LpDFD7L5u+GcTd+fOgPmrh+PomfOPOTsO0XBtvDu1zb4s/4dQL/GsVjdlpqgUsyzrltwQf3OeBDM3sQ/x9nHX5A+Un4Qb2/Csq/ZWY/ww9wnWdmd+PXJRqBn1bfHfi+c24pTXDOPU/b/rD1M7P67q5cfLfhd/FLH/wl6KZp6vy1ZvZn/FIDHwT1qMbPlN0XPxuqfkr7nsAMM3sY38Vj+P90T8L/5x+6NMePzS+H8WywLw84Ep88PR2Mr9rBzH6C/511DzZNta9uiXSz8+vwRRxr0Lr5InBi0NLybHDsH+HHJv4u5HdQYWb/B9xpZk/jL8y74v8I3Oeca9N/0o24At/N9D8zuwU/IPss/Fig7zYyCDuabgJyzexxfAttCv71PQ3/Pr8R/PpVZnYufkzOJ8F4zXX4RGkqPpE4MuzYA4DnzOwZ/DCCc/CzVV9sLqhI3gfB90uC8YwL8DMhp+ATlIXs/D5szAP4bvQHgUEhnyGAEufcEyE//xH/mXrVzG4Mtl2Mn0H5p7D4v4FvpQLfAp9rZvXj8d5wzr0RlBuCn4n8FP4fx9H4Vf834N+jzVmMH7f7IzMrwY+Lmuf8unI/x7fgvWtm/8J/Tn8EpONv+daU3YCHgs/OG/jW4j74z/iu+GtK/Xpm9TOXHzOzf+NbbzPxS5Z8D//PwLshx56Hn0B0cxDvufjrQuj6mg2J+LrT1s+yme2FT0yL8QnciTv/L8Lb9T0uzrm1wfXoN+Yn8ryB7/I8DX8tWht27Pr3wbDg63FmNjI41u9Dyt2E76n4ONh0Av6z9n/OudB/PKUxjU1P1aNzP/ADb//MV38YyvHr/PyBsNvmBOX3xa/pth7/X9dq/G1VxjVQtpBGpo+HlbuHyJcRqZ9aXocfizQPvzbU3i2s93fxs5/qF7r9BDg/rMxA/IzOz/F/uCvxF9QbCFvRHzgAvy7ViqBcCT4hvghIa+D8oXUJfwxtaaxBuSzgl8HrV/9H7lXg4EZ+B6fhW74q8V3PvyPCuybQ9DIijS2FMRH/R6IoeJ+9AxwVyXsGn7Q6fNdf6PaLgu1jm4n3COBfwft8e1DnpcHrO6yB8lPwLVZbgt/5UvwfukNDykwPzj0BnxgV4Vvv7gK6Rfh7jOh9ELzH5wWxV+HHef6FsGU5WnmeZQ2UHx28VtuDx9PAqAbKTW/iuNNDyvXAr9O4Joh/BT6pjniJCPwknE+D54cf/0D8mNpS/Ht/JrBvBMcswE+mehV/LasKXsP/4f/BsLDyOUGdPwveQ5vxn/Pfhr4WQXz/xM8aXxiUnQccHUFMEV932vpZxo/1bOz1c4Qt24JPJi8O3n/1M+QvDv89hfwOGnyElTsDv5Zc/TXrDZpZCkqPnR8W/CJFRCQC5mcYX4lf8641M4ylkzJ/t5LbnHPnxzoW6fw0Bk5EREQkwXRYAmdmqWb2OzNbav7+cEvN7Pch4z4wb7qZrTGzcjObZf62QKHHyTSzm81sk/l7FT5lbb+tjIiIiEjC6MgWuP/DD7D8KX6a98+Cn0MHdl6GX2DzQmAv/EDXl4NZa/VuxA92PAk//iEPP2A0tZ3jFxEREYkLHTYGLpiptdn5hfrqt92LH8z6rWA6/hr87Wn+EOzPxidxlzrnbjOzbvjZQmc55x4IygzCr7t0pItg9peIiIhIouvIFrg38VPlxwKY2Xj8UgvPBfuH4ZeE2HErG+cXf3wDPzMM/AJ/6WFl6hdwrS8jIiIi0ql15Dpwf8avS7XAzGqDc//BOXdrsL9v8HV92PPW49dbqi9Ti1/AMLxMXxoQrO10LkB2dvakQYMGtaUOzaqrqyMlJXnnhiRz/ZO57pDc9U/muoPqn8z1T+a6Q/vXf9GiRZuccw0uGN+RCdz3gdPxK1LPB3YH/mZmS51zd4aUC+/TtQa2hWu0jHPudvxaSkyePNnNnt2ahfEjN2vWLAoLC9v1HPEsmeufzHWH5K5/MtcdVP9krn8y1x3av/5m1uit+Toybb4OuN45N8M596lz7n78AoX1kxjq11MKb0kr4KtWuXX422zkN1FGREREpFPryAQuB9/9Gao2JIal+ATtsPqdZpaFn2n6drBpDv6WQqFlBgLjQsqIiIiIdGod2YX6NPBLM1uK70LdA38vu/vA32MjuP/eFWb2Gf4WIr/G32bjwaBMkZndCVxnZhvwtzO5AX8LoZkdWBcRERGRmOnIBO5C/H3absV3ea7F36Pw6pAy1+JvbnsL/h567wHTnHPFIWUuBmqAh4OyrwCnO+fCW/dEREREOqUOS+CCJOyi4NFYGYe/YfD0JspU4JPBC6MZn4iIiEii6MgWOBEREZEG1dXVsWnTJrZt20ZtbWJ0qnXr1o2FCxe26rmpqal0796d/Pz8Vi1FogROREREYm7VqlWYGUOHDiU9PR1/g6b4VlxcTNeuXZsvGMY5R3V1NevXr2fVqlUMHjy4xcdI3tX3REREJG6UlpYyYMAAMjIyEiJ5awszIyMjgwEDBlBaWtqqYyiBExERkbiQbHd1aEt9k+s3JSIiItIJKIETERERSTCaxCAiEXHOsWprOSu3lJGXnc6I3rlkZ6TGOiwRkaSkBE5EmlRRXcsNLy/i8Y9Ws7G4csf2nIxUjpzYj58dMorBvXJiGKGISOytXr2aq666iueee44NGzbQu3dvjjrqKK688koGDhwY9fMpgRORRm0preLUO95jwdrtHDmxL/uPzGdYfhe2l1fzxuJNPPHRap7+ZA2XHTGGEc7FOlwRkZhYunQpU6ZMYdiwYdx7772MGjWKJUuWcMUVV7DXXnvxzjvvMHTo0KieUwmciDSoorqWc+6bzRcbS7jj9MkcOr7PTvuP3MW3vv32yXn8/tmF7FmQyp77VNMtJz1GEYuIxMaPf/xjUlJSmDlzJjk5vkdi8ODBzJw5k1GjRvHjH/+YZ599NqrnVAInIg267sXPmbN8K7eesufXkrd6fbtlcdtpk7jrrWX88dkFHP33N7nvB3szNL9LB0crIp3RVU/PZ8Ga7R16zvH987jy6AkRl9+yZQsvvPACv//973ckb/VycnK44IIL+M1vfsPWrVvp0aNH1OLULFQR+ZrP1m3nnreXcfI+gzlql35NljUzzj5gGJfvnUVJZQ3f+ec7zFtd1EGRiojE1uLFi3HOMW7cuAb3jx8/Huccixcvjup51QInIl/zh2cXkpeVxi+mjYn4OSN7pPLIeftw+p3vcdLt73LHGZPZZ3ivdoxSRDq7lrSExVpjd49wwfjgjIyMqJ5PLXAispP5a4r43+JNnHfQCHp0adkFZ2RBLo9eMIWCvEzOvPsDPli2pZ2iFBGJD6NGjcLMmD9/foP7Fy5cSFpaGsOGDYvqeZXAichO7nxzKTkZqZy0d8tvrgzQr1s2D527L/26ZXHW3R/w0YqtUY5QRCR+9OzZkyOOOIJbb72VsrKynfaVlZVxyy23cNxxx9GtW7eonlcJnIjssKmkkqc/WcP3Jg+iW3brZ5MWdM3iwXP2pVduBqff9T7z12hMnIh0Xrfccgu1tbUceuihvPrqq6xcuZJZs2Zx2GGHkZ6ezk033RT1cyqBE5Ednv5kDdW1jpP3aV3rW6i+3XwS1zUzjTPv/oCVW8qaf5KISAIaNmwYs2fPZsKECZx22mkMHTqUqVOnkpKSwscff0zfvn2jfk4lcCKyw+MfrWZC/zxG9+kaleMN6J7NvT/Ym8rqWs64+322llZF5bgiIvFm4MCB/Otf/2L16tXU1tZy00038d577/HGG2+0y/mUwIkIAF9sKGHuqiKO22NAVI87qk9X7jhjL1ZtLeeH982msqY2qscXEYlHF154If/+979ZsGAB5eXlUT++EjgRAeD5T9cCcMxu/aN+7L2H9eTG7+/OnOVbueLxeTum1YuIdGbf+973uPzyy8nOzo76sZXAiQgAsxZtZNeB3SjIy2qX4x+1Sz8uOnQU/52zijvfXNou5xARSRZK4ESEbWVVfLRiK4Wje7freX568CiOnNiXPz63kFmfb2jXc4mIdGZK4ESE/y3eRJ2Dg8YUtOt5UlKMv3xvN8b0zePChz5iycaSdj2fiCSWZBte0Zb6KoETEV5ftJFu2ensPqh7u58rJyONf50+iYzUFM65bzYllTXtfk4RiX/p6entMtg/npWXl5Oe3ro1N5XAiQjvLd3MfsN7kZrS8L38om1gjxxuOWVPlm0q5VePfZp0/3WLyNcVFBSwevVqysrKOv01wTlHWVkZq1evpqCgdT0fupm9SJLbUFzByi3lnLHf0A49777De3HJtDFc9+Ln7DO8J6fsM6RDzy8i8SUvLw+ANWvWUF1dHeNoIlNRUUFWVusmfqWnp9OnT58d9W4pJXAiSe7D5dsA2GNwjw4/948OGsH7S7dw1dML2G1gdyYOiO69AkUkseTl5bU6oYmFWbNmsccee8Tk3OpCFUlyH67YSkZqChMHdPxFMyXF+Ov3d6dnTgY/m/ERFdVa5FdEJBJK4ESS3IfLt7LLwG5kpqXG5Pw9u2Rw3Xd3ZcnGUq594fOYxCAikmiUwIkksaqaOuauLmLSkI7vPg114KjenLHfEO56aylvf7EpprGIiCSCDkvgzGyZmbkGHs8G+83MppvZGjMrN7NZZjYh7BiZZnazmW0ys1Ize8rMBnZUHUQ6myUbS6iqqYuLsWe/PHIcw3t34dL/fML2isQYwCwiEisd2QK3F9Av5LEn4IBHgv2XAZcAFwZlNwAvm1nXkGPcCJwAnAQcCOQBz5hZbPp+RBLcZ+u2AzCub9dmSra/7IxUbvje7qzbXsE1z38W63BEROJahyVwzrmNzrl19Q/gKGA78B8zM+Ai4Brn3KPOuXnAGUBX4GQAM+sGnA38wjn3snPuQ+A0YFfg0I6qh0hn8tnaYjJSUxiW3yXWoQCw+6DunH3AMB58bwXvfrk51uGIiMStmIyBCxK2s4F/O+fKgGFAX+Cl+jLOuXLgDWBKsGkSkB5WZiWwMKSMiLTAwnXFjOqTS1pq/AyH/flhYxjcM4dfPjpXs1JFRBphsVjt2MymAS8CezjnPjazKcBbwBDn3IqQcncBA5xzh5vZycB9QLoLCdrMXgUWO+fOa+Rc5wLnAvTp02fSjBkz2q1eACUlJeTm5rbrOeJZMtc/Eet+0WtlTOiVyjm7Zrb5WNGs/4LNtVz7QQVHDUvne2MyonLM9pSIr300qf7JW/9krju0f/2nTp06xzk3uaF9sVrI9xzgA+fcx2Hbw7NJa2BbuCbLOOduB24HmDx5sissLGxRoC01a9Ys2vsc8SyZ659odd9SWsW2F16mcPdRFH5jeJuPF836FwIrmMsjs1fyo2/uwy4DYz/JoimJ9tpHm+qfvPVP5rpDbOvf4f0mZlYAHAv8K2TzuuBr37DiBcD6kDKpQH4TZUQkQvUTGMb2i/0EhoZcftQ48nMzuezRuVTX1sU6HBGRuBKLgS9nApVAaF/mUnyCdlj9BjPLws80fTvYNAeoDiszEBgXUkZEIrR4fQkAo/vEZwLXLTudq4+dyMK127n9jS9jHY6ISFzp0AQumLzwQ2CGc664fnswpu1G4JdmdryZTQTuAUqAB4MyRcCdwHVmdqiZ7QHcD8wFZnZkPUQ6g2WbS8nJSKWga9vHv7WXIyb25ahd+vK3VxbzxYaSWIcjIhI3OroFrhAYxc7dp/WuBW4AbgFm49eKmxaa6AEXA48BD+MnPZQARzvnNFVNpIWWbSplSK8u+P+r4tf0YyaQk5HKz2Z8RGWNPuoiItDBCZxz7jXnnDnn3m9gn3POTXfO9XPOZTnnDgrWgwstU+Gcu9A518s5l+OcOzpYSkREWmj55jKG5efEOoxmFXTN4rrv7Mb8Ndv503Na4FdEBHQvVJGkVFNbx8qtZQzpFR8L+DbnsPF9OHPKUO59Zxkfr9wW63BERGJOCZxIElqzrYLqWsfQXvHfAlfvkmmj6Z2byW+emEdtXcevXykiEk+UwIkkoWWbSwEYmiAtcABds9L59bfG8+nqIh58f0XzTxAR6cSUwIkkoeX1CVyc3AM1Ukfv2o8pI3px3QufsamkMtbhiIjEjBI4kSS0bHMZ2enxvYRIQ8yMq4+dSHl1Lde98HmswxERiRklcCJJaPnmUob0yon7JUQaMrIglzOnDOWROSuZt7oo1uGIiMSEEjiRJLR6WwUDumfHOoxWu/CQUfTMyeDqpxfg1wEXEUkuSuBEktDaonL6dc+KdRitlpeVzqWHj+H9ZVt49tO1sQ5HRKTDKYETSTLlVbVsK6umX7fEbYED+N7kQYzrl8efnvuMimrdoUFEkosSOJEks6aoHIB+3RK3BQ4gNcW48ujxrN5Wrpvdi0jSUQInkmTWFVUAJHwLHMC+w3tx1C59+cesJawNElMRkWSgBE4kyazZ5hOd/gk8Bi7U5UeOo9Y5/vy87pMqIslDCZxIklkbtMD1TfAu1HqDeuZw7oHDeeLjNcxZvjXW4YiIdAglcCJJZm1ROfm5GWSmpcY6lKj5UeEICrpmMv2p+bpPqogkBSVwIklmzbaKTtP6Vq9LZhq/Ce6Tet87y2IdjohIu1MCJ5Jk1hVVdIoJDOG+tWs/Dhrdm+tf/HzHOD8Rkc5KCZxIkllTVE7/TtYCB/4+qb//9kRqnWP6U/NjHY6ISLtSAieSREorayiuqKFvJ2yBAz+h4aJDR/PSgvW8OH9drMMREWk3SuBEksjG4koACrpmxjiS9nP2AcMY27cr05+aT0llTazDERFpF0rgRJLIphKfwPXuxAlcemoKfzx+F9Ztr+AvL30e63BERNqFEjiRJFLfApef23kTOIA9B/fg1H2GcO/by5i7aluswxERiTolcCJJZGMStMDV+8URY8jPzeTyxz6lprYu1uGIiESVEjiRJLKxuJIUg55dMmIdSrvLy0pn+jETmL9mO/e8vSzW4YiIRJUSOJEksrG4kl65maSmWKxD6RBHTuzLwWMLuOHlRazW2nAi0okogRNJIhuLK+ndyce/hTIzrjpmAs7Bb5+Yh3O6zZaIdA5K4ESSyMaSSvKTYPxbqEE9c/j5YaN55bMNvDBPa8OJSOegBE4kiWxKsha4emftP5Rx/fK48qn5bK+ojnU4IiJtpgROJEk459hYUpkUM1DDpaWm8Kfjd2FjSSV/eVFrw4lI4lMCJ5Ikisqrqa51SZnAAew+qDtn7DeU+95dzkcrtsY6HBGRNlECJ5Ik6hfxTdYEDuCSaaPp0zWLyx/7lGqtDSciCaxDEzgz62dm95rZRjOrMLMFZnZQyH4zs+lmtsbMys1slplNCDtGppndbGabzKzUzJ4ys4EdWQ+RRLQjgUvCMXD1ugZrw322rpi73lwa63BERFqtwxI4M+sOvAUY8E1gHHAhsCGk2GXAJcH2vYJ9L5tZ15AyNwInACcBBwJ5wDNmltq+NRBJbF/dhaHzL+LblMMn9OHQcX3468xFrNxSFutwRERapSNb4C4D1jrnTnfOve+cW+qce8U5txB86xtwEXCNc+5R59w84AygK3ByUKYbcDbwC+fcy865D4HTgF2BQzuwLiIJZ0tpFQA9uyRvCxwEa8MdO4EUM377pNaGE5HE1JEJ3LeB98zsYTPbYGYfm9lPgsQNYBjQF3ip/gnOuXLgDWBKsGkSkB5WZiWwMKSMiDRga1k1ZtAtOz3WocTcgO7ZXDJtDK99vpFnP10b63BERFrMOuq/TzOrCL79K/AIsDtwM/BL59zfzWwKvot1iHNuRcjz7gIGOOcON7OTgfuAdBcSuJm9Cix2zp3XwHnPBc4F6NOnz6QZM2a0S/3qlZSUkJub267niGfJXP94r/t9Cyp5f20Nfz+kS7scP97rH662znH1uxVsq3T88YBsuqS3/vZiiVb3aFP9k7f+yVx3aP/6T506dY5zbnJD+9La7axflwLMds5dHvz8kZmNAn4M/D2kXHhGaQ1sC9doGefc7cDtAJMnT3aFhYUtDLtlZs2aRXufI54lc/3jve7/WfMhBWXb2y3GeK9/QwpGF3HsLW/ydkk+fzhul1YfJxHrHk2qf/LWP5nrDrGtf0d2oa4FFoRtWwgMDr6vv8dN37AyBcD6kDKpQH4TZUSkAVtLq+jRJbknMITbZWA3zpwyjAfeW8Gc5VtiHY6ISMQ6MoF7CxgTtm00sDz4fik+QTusfqeZZeFnmr4dbJoDVIeVGYif0VpfRkQasKW0ih45SuDC/XzaaPp3y+JXj83T2nAikjA6MoH7K7CvmV1hZiPN7LvAT4FbAIIxbTcCvzSz481sInAPUAI8GJQpAu4ErjOzQ81sD+B+YC4wswPrIpJwtpVV07OLJjCEy81M46pjJ/L5+mJuf+PLWIcjIhKRDkvgnHMf4Geifg+YB/wB+A1wa0ixa4Eb8EndbKAfMM05VxxS5mLgMeBhfKteCXC0c662nasgkrCcc2wpUxdqYw4b34cjJvTlplcWs3xzaazDERFpVofeicE596xzbjfnXJZzbrRz7qbQ2aTOm+6c6xeUOShYDy70GBXOuQudc72ccznOuaODpUREpBFlVbVU1dTRU12ojZp+zATSU1O44nGtDSci8U/3QhVJAvWL+GoMXOP6dsvisiPG8OYXm3ji49WxDkdEpElK4ESSwNayIIFTF2qTTtlnCLsP6s7vnlnI1iDpFRGJR0rgRJLA1rJqAE1iaEZqivGn43dhe3k1f3xuYazDERFplBI4kSSwVV2oERvXL48fHjic/8xZxTtLNsc6HBGRBimBE0kCX93IXglcJH52yCgG98zhisc/paJaE9xFJP4ogRNJAlvLqkgxyMtSF2oksjNS+cNxE/lyUym3zloS63BERL5GCZxIEthSWkX3nAxSUlp/w/Zkc+Co3nx79/78Y9YXfLmxJNbhiIjsRAmcSBLYWlZF9xy1vrXUr745jqy0VK58ar7WhhORuKIETiQJFJVX0z1bCVxLFXTN4ufTRvO/xZt4cf66WIcjIrKDEjiRJFBUXk03JXCtctq+Qxjbtyu/e2YhZVU1sQ5HRARQAieSFJTAtV5aagpXHzuR1dvKueW1L2IdjogIoAROJCkUlSmBa4u9h/Xk+D0GcPsbX2pCg4jEBSVwIp1cXZ2juLJGCVwb/fKosZrQICJxQwmcSCdXXFGDc5CnBK5NCrpmcfFhmtAgIvFBCZxIJ1dU7u+Dqha4tjt9P01oEJH4oAROpJNTAhc9oRMa/qE7NIhIDCmBE+nklMBF197DenLs7v257Y0vWbG5LNbhiEiSUgIn0sntSOB0J4aoufzIcaSlGFc/syDWoYhIklICJ9LJqQUu+vp2y+InB49k5sL1zPp8Q6zDEZEklNbYDjP7eSuOd4dzbnsb4hGRKFMC1z7OPmAY/5m9iqufXsCv9tSyIiLSsRpN4IDrgVVAbYTHGgQ8ASiBE4kjReXVpKca2empsQ6lU8lMS+W33xrPWfd8wMvLMzg01gGJSFJpKoEDmOyci6h/wMyKoxCPiERZ/W20zCzWoXQ6U8cWcPDYAp5cvIFLtldQkJcV65BEJEk0NQbuKqAl94z5I7ClbeGISLRtL6/WIr7t6LffGk9NHVzz/GexDkVEkkijCZxz7irnXMRz5J1zf3LObYtKVCISNbqRffsamt+Fw4em89hHq5mzXP/DikjHiGgWqpmlmFlKyM99zeyHZrZ/+4UmItGgBK79HT0inb55WVz51Hxq6zShQUTaX6TLiDwLXAhgZrnAbOA6YJaZnd5OsYlIFGyvUALX3rLSjMuPGsu81dt5ZPbKWIcjIkkg0gRuEvBq8P3x+JmmBcA5wKXtEJeIRIla4DrGMbv1Z++hPbn2hc8oKquOdTgi0slFmsB1BbYF308DHnfOVeOTuhHtEJeIREFdnWO7ErgOYWZcecx4isqr+evMRbEOR0Q6uUgTuBXA/mbWBTgceDnY3hPQzQBF4lRZdS11DrpmNbdikETDhP7dOHHvwdz/7nIWrdfKSiLSfiJN4G4A7scv7LsaeCPY/g3g03aIS0SioLjCd+V1zVILXEe5dNoYumSkcvXTC3BOExpEpH1ElMA5524D9gV+ABzgnKsLdi0BfhPJMcxsupm5sMe6kP0WlFljZuVmNsvMJoQdI9PMbjazTWZWamZPmdnAiGoqkoSKK2oAtcB1pJ5dMrj4sNG8+cUmXl6wPtbhiEgnFfHN7J1zc5xzjzvnSkK2Peuce6sF5/sc6Bfy2CVk32XAJfjZrnsBG4CXzaxrSJkbgROAk4ADgTzgGTPTPYJEGlCfwOVmKoHrSKfuO4RRBbn8/tmFVFRHejdCEZHIRXxVN7M9gKn42ac7JX7OucsiPEyNc25d+Ebz9/i5CLjGOfdosO0MfBJ3MnCbmXUDzgbOcs69HJQ5DVgOHAq8GGldRJKFulBjIz01hSuPnsCpd77HnW8u5cdTR8Y6JBHpZCJdyPcyYA6+dawQ3/pV/zigBecbbmarzWypmc0ws+HB9mFAX+Cl+oLOuXL8WLspwaZJQHpYmZXAwpAyIhKivgUuT12oHe6AUflMG9+HW177gnVFFbEOR0Q6GYtkkK2ZrQWmB2PhWncisyPxy5F8hm/F+zUwFpgAjAHeAoY451aEPOcuYIBz7nAzOxm4D0h3IUGb2avAYufceY2c91zgXIA+ffpMmjFjRmurEJGSkhJyc3Pb9RzxLJnrH491n7WymnvmV3FDYTY9syIeMdEq8Vj/jtJY3TeU1fGrN8vZq28q5+3aeW90n8yvPSR3/ZO57tD+9Z86deoc59zkhvZF+m95CvBKW4Jwzj0f+rOZvQt8CZwBvFtfLOxp1sC2cE2Wcc7dDtwOMHnyZFdYWBh50K0wa9Ys2vsc8SyZ6x+PdV/0xhKY/xnTpn6j3cfBxWP9O0pTdV+e+hm3vLaEX3x7N/Yc3KNjA+sgyfzaQ3LXP5nrDrGtf6T/kv8DOCuaJw4mQ8wHRgH14+L6hhUrAOqnca0DUoH8JsqISIjiihpSDLpkaJ5PrFxQOJL83Ez++OxCLSsiIlETaQJ3FbCnmX1sZveb2V2hj9ac2Myy8F2oa4Gl+ATtsLD9BwJvB5vmANVhZQYC40LKiEiI4ooacjPT8POEJBa6ZKZxybTRzF6+lRfnf20Ol4hIq0SawP0BfwutaqAH0Dvs0Swzu97MDjKzYWa2D/BfoAtwbzCm7Ubgl2Z2vJlNBO4BSoAHAZxzRcCdwHVmdmgwK/Z+YC4wM8J6iCSV4ooazUCNA9+dNJDRfXK55vnPqKqpa/4JIiLNiHRQzAXAyc65h9twroHAQ/gu0I34cW/7OueWB/uvBbKBW/BJ4nvANOdc6P1oLgZqgIeDsq8ApzvntNCSSAOKK6q1iG8cSEtN4fKjxnHW3R/wwHvLOWv/YbEOSUQSXKRX9nLgo7acyDl3YjP7HTA9eDRWpgK/lMmFbYlFJFn4FjglcPGgcHRvDhiZz99eWczxew6kW7ZaRkWk9SLtQv0rcJFpII1IQimurNZdGOKEmXH5UWMpKq/m1te+iHU4IpLgIr2yH4i/cf03zWwBfizcDs65Y6IdmIi0XUlFDcPzk3eNpngzoX83TthzIHe/tYxT9x3CoJ45sQ5JRBJUpC1wm4DHgFfxs0U3hz1EJA6pCzX+XDJtNCkpcN2Ln8c6FBFJYBFd2Z1zUV0DTkQ6hmahxp9+3bI558Dh3PzqF/zggGHsPqh7rEMSkQQU6b1QJ5jZrg1s39XMxkc/LBFpq8qaWqpq69QCF4fOO2gE+bkZWtxXRFot0i7U24GJDWwfH+wTkThTfyN7JXDxJzczjYsPG837y7bw0gLdSEZEWi7SBG5X4P0Gtn8A7BK9cEQkWpTAxbfvTx7EqIJc/vz8Z9TUanFfEWmZSBO4WqBbA9t74G8mLyJxpiRI4HIzNQYuHqWlpvCLw8fw5aZSHv1wVazDEZEEE2kC9zpwhZntuCO2maUBVwBvtEdgItI2xRV+tR+1wMWvw8b3YfdB3fnbzMVUVOuGMiISuUgTuMvwa8F9EdzM/n5gMXAA8Iv2Ck5EWm+7ulDjnplx2eFjWFNUwYPvrYh1OCKSQCJK4Jxzn+PHwT0I9AR6AQ8AuznnFrZfeCLSWvUtcHlaRiSuTRmZz/4je3HLa19QUlkT63BEJEFE2gKHc26tc+4K59w3nXNHOed+7Zxb057BiUjr1ScDupVW/Lt02hg2l1Zx95tLYx2KiCSIRhM4MytoyYHMLF/3ShWJH/WzUHPVhRr39hjcg8PG9+H2N75kW1lVrMMRkQTQVAvc2hYmcV8Cw9oYj4hESXFFNdnpqaSnRtzQLjF06bQxlFTV8I/Xl8Q6FBFJAE39a27A+WZWEuGxNNBGJI6UVNao9S2BjOnblW/vPoB7317GD/YfRp+8rFiHJCJxrKmr+wqgJfdAXQdUty0cEYmW7bqRfcK56NBRPP3JGv7+6hf87tsN3fxGRMRr9OrunBvagXGISJTpRvaJZ0ivLpy49yAeen8F5xw4nMG9cmIdkojEKQ2OEemkiiuqyVMLXMK58OBRpKYYN85cFOtQRCSOKYET6aRKKmq0hEgC6pOXxZlThvL4x6tZtL441uGISJxSAifSSRVrDFzCOv+gEeRmpHH9i5/HOhQRiVNK4EQ6qeKKao2BS1A9umRwzjeG89KC9Xy8cluswxGROKQETqQTqq1zlFbVqgs1gf3ggGH06pKhVjgRaVDECZyZ9TGzS83sH2aWH2zb38y0eK9InKm/jZa6UBNXbmYa5x80gje/2MTsZVtiHY6IxJmIEjgzmwR8DpwCnA3kBbsOA/7QPqGJSGvpRvadwyn7DqZXlwz+9sriWIciInEm0ha464G/Oef2ACpDtr8I7B/1qESkTXbcyF4tcAktJyON8w4azv8Wb2LOcrXCichXIk3gJgH3NrB9LdAneuGISDTU38heXaiJ79R9h9CrSwY3zlQrnIh8JdIErhzo0cD2scCG6IUjItFQ34WqWaiJLycjjXO/oVY4EdlZpAnck8CVZpYZ/OzMbCjwZ+DR9ghMRFpPLXCdy2n7DaGnWuFEJESkCdylQE9gI5ADvAl8AWwDft0ukYlIq+1I4LSMSKewcyvc1liHIyJxIKIEzjm33Tl3APBt4P+AvwFHOOcOcs6VtmN8ItIKX7XAqQu1szhtX98Kd5NmpIoILVzI1zn3qnPueufctc65mW05sZn9ysycmf09ZJuZ2XQzW2Nm5WY2y8wmhD0v08xuNrNNZlZqZk+Z2cC2xCLS2RRXVJOWYmSla63uzqJLZhpnHzCM1xdt5NNVRbEOR0RiLKL+FTP7bSO7HFCB7059wTlXHuHx9gXOAeaG7boMuAQ4E7/u3G+Bl81sjHOu/q7ONwLHAicBm4EbgGfMbJJzrjaS84t0diWVNeRmpWFmsQ5Foui0/Ybwz9eXcOusL/jHqZNiHY6IxFCkA2S+CwwGugBrgm39gVL8uLhBwAYzO8g592VTBzKzbsAD+AWBfxuy3YCLgGucc48G287Az3I9GbgteO7ZwFnOuZeDMqcBy4FD8evSiSQ93ci+c8rLSueM/YZyy6wv+GJDMSMLusY6JBGJkUj7V/4CfAAMdc4Nds4NBoYC7wFX45O5RcBfIzjW7cB/nXOvhm0fBvQFXqrfELTovQFMCTZNAtLDyqwEFoaUEUl6xRXVdM3U+LfO6Kz9h5KZlsI/ZjX5v7KIdHLmnGu+kNlS4Fjn3Nyw7bsDTzjnhgbdok865xpd2NfMzgHOB/ZzzlWZ2SxgnnPuJ2Y2BXgLGOKcWxHynLuAAc65w83sZOA+IN2FBG5mrwKLnXPnNXDOc4FzAfr06TNpxowZzda3LUpKSsjNzW3Xc8SzZK5/PNX9T+/50QyX75PdYeeMp/p3tI6u+wMLK3llRQ1/PjCb3jmxH+eYzK89JHf9k7nu0P71nzp16hzn3OSG9kXax9IHyGpgeyZQEHy/Hr/ESIPMbAzwR+BA51xVE+cKzyitgW1fO3xjZZxzt+Nb/Zg8ebIrLCxs5lBtM2vWLNr7HPEsmesfT3W/9pP/0b97FoWFe3XYOeOp/h2to+s+Zo9yZl37GnOrCvjdURM77LyNSebXHpK7/slcd4ht/SP9120mfgzaXmaWEjz2Av4BvByU2QVY2sQx9gPygXlmVmNmNcBBwAXB95uDcn3DnleATw4B1gGpwXEaKyOS9Iorq7WESCfWr1s2J+w5kIdnr2RDcUWswxGRGIg0gfshPkF6D38z+0rg3WDbOUGZYvyCv415Ap/k7R7ymA3MCL5fhE/QDqt/gpllAQcCbweb5gDVYWUGAuNCyogkvZKKGnK1iG+ndv5BI6iprePON5v6v1lEOqtmr/BmloJv8foOfrLCWHyX5ULn3KL6cs6515o6jnNuG/7ODaHHLgW2OOfmBT/fCFxhZp/hE7pfAyXAg8ExiszsTuA6M9vAV8uIzMW3EookPeecZqEmgaH5XfjWrv359zvL+dFBI+iekxHrkESkA0XSAueAj4G+zrlFzrmnnHNPhiZvUXQtPiG7Bd861w+YFrIGHMDFwGPAw/hJDyXA0VoDTsSrrKmjps6RqwSu0/tR4QhKq2q59+3lsQ5FRDpYswlcMNvzc6B3tE/unCt0zv0k9FzOuenOuX7OuazgVl3zwp5T4Zy70DnXyzmX45w7OlhKRESA7RXVgG6jlQzG9cvj0HEF3P32Ukora2Idjoh0oEjHwF2G77bc3bS0u0hcq78Pap5a4JLCBVNHsq2smofeX9F8YRHpNCJN4B4B9sFPIqgws+2hj/YLT0RaqiRI4DSJITnsObgHU0b04vY3vqSiWiNJRJJFpFf4nzRfRETiQX0LnLpQk8ePp47klDve49EPV3HKPkNiHY6IdICIEjjn3L3tHYiIREdJpR8Dpxa45DFlRC92G9Sdf76+hO9PHkRaauzvziAi7avFn3Iz62tmg0Mf7RGYiLTO9h0tcErgkoWZ8ZOpI1m5pZxn5q6NdTgi0gEiSuDMrJuZ3Wtm5cBq/B0XQh8iEieKlcAlpUPGFjCmT1duee0L6uqav8e1iCS2SFvgrgd2A74NVAAnA78AVgHfb5fIRKRVNIkhOaWkGBdMHcHiDSW8vFB3FhTp7CJN4I4ELnTOvQjUAnOcczcAvwTOa6/gRKTliiuqyclI1TioJPTNXfoxpFcOt7z2BX4JTxHprCK9wncH6pf6LgJ6Bd+/A0yJckwi0gYllboParJKS03hgsIRzF1VxMyFG2Idjoi0o0gTuCXA8OD7hcCJwYK+xwNb2iMwEWkd3Qc1uZ2w50CG5Xfh+hc/11g4kU4s0gTuHmDX4Ptr8N2mVcB1wJ+jH5aItNb2impytQZc0kpLTeHnh43m8/XFPPXJmliHIyLtJKIEzjn3V+fcTcH3rwJj8ZMXdnfO/b0d4xORFiqprNFttJLcN3fpx/h+edzw8iKqa+tiHY6ItINWjXJ2zq1wzj3mnPs02gGJSNuoC1VSUoxfHD6GFVvKePiDlbEOR0TaQcRXeTPbBzgEKCAs8XPO/TTKcYlIKxVXVGsSg1A4pjeTh/TgplcW851JA8lKT411SCISRZEu5HspfsbpmcDuwC4hj4ntFJuItEJJRY3ugyqYGZcdMZYNxZXc+/ayWIcjIlEW6b/pPwN+qvFuIvGtts5RWlWrFjgBYO9hPTlodG/+8foSTtxrMN1ylNiLdBaRjoHLA55rz0BEpO1KKnUbLdnZL48cS1F5NX9/bXGsQxGRKIo0gXsIOKI9AxGRtiuuqAYgT12oEhjXL4/vThrIvW8vZ/nm0liHIyJR0ui/6Wb285AfVwJXmdn+wFygOrRscFstEYmx+hvZ56oFTkJcMm0MT3+ylj+/8Bm3njIp1uGISBQ0dZW/MOznEvxts8JvneUAJXAicUBdqNKQPnlZnHfQcG6cuZjZy7YweWjPWIckIm3UaBeqc25YhI/hjR1DRDpWfReqJjFIuHO/MZw+eZn87tmFusWWSCcQ6TIiGWaW1cD2LDPLiH5YItIa9V2oWkZEwuVkpHHptDF8snIbT8/VLbZEEl2kkxj+A1zQwPbzgUeiF46ItMVXCZxa4OTrTthzIOP75XHtC59TUV0b63BEpA0iTeD2B15qYPvLfH1MnIjEiBI4aUpKivHrb45j9bZy7n5rWazDEZE2iDSBywFqGtheB3SNXjgi0hYlldWkphjZum2SNGLKyHwOGVvAra99waaSyliHIyKtFGkCNxc4qYHtJwPzoheOiLRFcUUNuZlpmFmsQ5E4dvlR4yivruUvLy2KdSgi0kqR9rP8DnjCzEYCrwbbDgG+CxzXHoGJSMsVV9So+1SaNbIgl9P3G8rdby/llH0GM3FAt1iHJCItFFELnHPuWeBoYAhwU/AYDBzjnHum/cITkZaob4ETac7PDh1Fj5wMrn56Ac5pWRGRRBNpFyrOuReccwc457oEjwOcc8+3Z3Ai0jLFFdW6jZZEpFt2OpdOG8P7y7bw7KdrYx2OiLRQpOvA9Taz3iE/72JmvzezhsbFiUiMlFSqC1Ui9/29BjG+Xx5/fHYh5VVaVkQkkUTaAvcIvgsVM8sH3sCPffunmV0SyQHM7MdmNtfMtgePd8zsmyH7zcymm9kaMys3s1lmNiHsGJlmdrOZbTKzUjN7yswGRlgHkU6vuKJG90GViKWmGFcePZ41RRXc9saSWIcjIi0QaQK3K/Bu8P13gC+ccxOA04HzIjzGKuD/gD2ByfjJEE+Y2a7B/suAS/D3YN0L2AC8bGahy5TcCJyAnxF7IJAHPGNmWjNBBN+FqhY4aYl9hvfim7v245+vL2H1tvJYhyMiEYo0gcvG38we4FDgqeD7D4FBkRzAOfekc+5559wXzrlFzrkrgGJgP/NrHlwEXOOce9Q5Nw84A7/G3MkAZtYNOBv4hXPuZefch8Bp+OTy0AjrIdJpOecoqawhN1Nj4KRlLj9yLM7Bn55bGOtQRCRCkSZwi4HjzWwQMI2v7srQB9jW0pOaWaqZnQjkAm8Dw4C+IcfFOVeO76qtv9PDJCA9rMxKYCG6G4QIlTV1VNc6tcBJiw3skcP5B43gmblreX/plliHIyIRsEimj5vZ8cBD+HXjXnHOTQu2XwHs75w7KqKTme0CvANk4Vv0TnHOPWtmU4C3gCHOuRUh5e8CBjjnDjezk4H7gHQXErSZvQosds412JVrZucC5wL06dNn0owZMyIJtdVKSkrIzc1t13PEs2Suf6zrXlTp+NlrZZw2PoNDBnd8K1ys6x9LnaHulbWOy/9XTm66MX1KFiktWAy6M9S/LZK5/slcd2j/+k+dOnWOc25yQ/si+lfdOfeYmQ0G+gOfhOyaCTzaglg+B3YHuuPHst1rZoWhpworbw1sC9dkGefc7cDtAJMnT3aFhYWNFY2KWbNm0d7niGfJXP9Y1/3LjSXw2uvsucs4Cvfo+Lk9sa5/LHWWulfmr+GnD33E+i4jOGnvwRE/r7PUv7WSuf7JXHeIbf1bsg7ceufcR865upBt7znnPmvBMaqCMXCznXOXAx8DFwPrgiJ9w55SAKwPvl8HpAL5TZQRSVo7bmSvMXDSSkfv2o+9h/bkuhc/Z1tZVazDEZEmNJrAmdlNZtYl5PtGH208fyawFJ+gHRZy/iz8TNO3g01zgOqwMgOBcSFlRJJWSaVP4LSMiLSWmXHVsRMoKq/mzy98HutwRKQJTV3pd8FPGqj/vjER3YPFzK4BngVW8tXs0kLgm845Z2Y3AleY2WfAIuDX+HFyDwI454rM7E7gOjPbAGwGbgDm4rtyRZLa9vJqAE1ikDYZ1y+Ps6YM5Y43l/LdyQPZc3CPWIckIg1o9ErvnJva0Pdt0Bf4d/C1CJ94HemcezHYfy1+uZJbgB7Ae8A051xxyDEuBmqAh4OyrwCnO+e0hLgkve0VPoHrlq0uVGmbiw4bzTNz1/Lrx+fx1E/2Jy014tE2ItJBOuxT6Zw70zk3xDmX6ZwrcM4dGpK84bzpzrl+zrks59xBwXpwoceocM5d6Jzr5ZzLcc4dHSwlIpL0isqVwEl05Gam8dujx7Ng7Xbue2d5rMMRkQY0m8CZWbaZXRncBqvEzIrN7BMz+7WZZXdEkCLSvKLyalJTjNxMdaFK2x05sS8Hje7NDS8vYv32iliHIyJhmkzgzCwNf8urX+EnGtyM7+JcDvwWmBmUEZEYKyqvJi8rDWvB+l0ijTEzrj52AtW1dVz99IJYhyMiYZpLvs4FRgJ7Oufmh+4ws4nAa0GZW9snPBGJVFF5DXnqPpUoGtKrCxcePJLrX1rE0fPWccTE8JWeRCRWmutC/Q7wh/DkDSAYn/anoIyIxNj28mqNf5OoO++gEUzon8evn5jH1lKtDScSL5pL4Cbgu1AbMxOYGL1wRKS1ipTASTtIT03huu/sxrayKq5+Rl2pIvGiuQSuB7Cxif0b8bfFEpEY215erS5UaRfj++fxk4NH8vhHq3l5gW58IxIPmkvgUvHrrjWmLigjIjGmFjhpTxcUjmRs365c8finFJVVxzockaTX3CQGA/5tZpWN7M+Mcjwi0grOOSVw0q4y0lK4/ru7cewtb3HV0/O54fu7xzokkaTWXAvcvcAa/G2rGnqsAe5rzwBFpHnl1bXU1DnyspTASfuZOKAbP5k6ksc+Ws2TH6+OdTgiSa3JFjjn3FkdFYiItJ7uwiAd5cKDR/LWF5u44vF57DGoB4N75cQ6JJGkpBvciXQCSuCko6SlpnDjibuTYnDhjI+orq2LdUgiSUkJnEgnUD+oXAmcdISBPXK45oRd+WTlNv7y0qJYhyOSlJTAiXQCaoGTjnbULv04ae/B/PP1JczbVBvrcESSjhI4kU5ge4Vf7UcJnHSk335rPKMKcrl9bgVri8pjHY5IUlECJ9IJqAVOYiE7I5V/nLonVbVw/v1zqKhWS5xIR1ECJ9IJ1CdwuVnNLe0oEl0jC7pyzq6ZfLKqiN88MQ/nXKxDEkkKSuBEOoHt5dV0zUojNcViHYokoUl90rjw4JH8Z84q/v3eiliHI5IUlMCJdAK6C4PE2sWHjmbqmN5c9dR8Pli2JdbhiHR6SuBEOoEtpVX06pIR6zAkiaWkGDeeuAcDe2Rz/v1zWLmlLNYhiXRqSuBEOoGtZVX0UAInMdYtO507ztiLmjrHmXe/r5vei7QjJXAincCW0ip65iiBk9gbWZDLbadNYsWWMs7/9xyqanSnBpH2oAROpBPYWqoWOIkf+w7vxbXf2ZV3vtzM5Y99qpmpIu1Aaw6IJLjKmlpKq2rpkaNJDBI/jttjIMs3l3HjzMUM6pnNRYeOjnVIIp2KEjiRBLctGGekFjiJNz87ZBQrt5Rz48zF5Odmcuq+Q2IdkkinoQROJMFtKa0C0Bg4iTtmxjUn7MK2sip+8+Q8uuek861d+8c6LJFOQWPgRBLc1iCBUwucxKP01BRuOWVP9hrSk4sf/pg3Fm2MdUginYISOJEEt6UsaIFTAidxKis9lX+dMZmRBV05/99z+GjF1liHJJLwlMCJJLgdLXDqQpU41i07nXt/sBe9u2Zyxl3vM291UaxDEkloSuBEEtyWUj+JobtmoUqcK+iaxQM/3IeuWemccsd7zF+jJE6ktZTAiSS4rWVV5GWlkZ6qj7PEv4E9cphx7r50yUjl1DveY+Ha7bEOSSQhddgV38wuN7MPzGy7mW00s6fNbGJYGTOz6Wa2xszKzWyWmU0IK5NpZjeb2SYzKzWzp8xsYEfVQyTebCmt0vg3SSiDeubw0Ln7kpWeyil3vMfn64pjHZJIwunIf9kLgVuBKcDBQA0w08x6hpS5DLgEuBDYC9gAvGxmXUPK3AicAJwEHAjkAc+YWWo7xy8Sl3QfVElEQ3p14cFz9iU91Tj5X++yeL2SOJGW6LAEzjl3uHPubufcPOfcp8BpQG9gf/Ctb8BFwDXOuUedc/OAM4CuwMlBmW7A2cAvnHMvO+c+DI6zK3BoR9VFJJ7oPqiSqIbl+yQuJcU46V/v8dk6daeKRCqWg2a6Buevn08+DOgLvFRfwDlXDryBb7UDmASkh5VZCSwMKSOSVHQfVElkI3rn8tA5+5KaAife/i6frNwW65BEEoLF6ibDZvYIMAqY7JyrNbMpwFvAEOfcipBydwEDnHOHm9nJwH1AugsJ3MxeBRY7585r4DznAucC9OnTZ9KMGTPatV4lJSXk5ua26zniWTLXPxZ1d85x7stlHDI4nRPHxjaJ02ufnHWH6NR/Q1kd131QQXGV4+JJWYzpmTijYpL59U/mukP713/q1KlznHOTG9oXk1tpmdkNwAHAAc652rDd4RmlNbDta4dsrIxz7nbgdoDJkye7wsLCFsfbErNmzaK9zxHPkrn+saj79opqql98iUnjR1L4jeEdeu5weu0LYx1GzESr/t/Yv4JT7niXv35Uzj9PnUThmIK2B9cBkvn1T+a6Q2zr3+FdqGb2V/wEhIOdc1+G7FoXfO0b9pQCYH1ImVQgv4kyIkljw/ZKAAryMmMciUjb9e2WxSPn7cfw/FzOuW82z3+6NtYhicStDk3gzOxv+AkJBzvnPgvbvRSfoB0WUj4LP9P07WDTHKA6rMxAYFxIGZGksbHYJ3C9c5XASefQKzeTh87dl10GdOPHD37IIx+sjHVIInGpI9eBuwU4C9/6ttXM+gaPXIBgTNuNwC/N7Phgjbh7gBLgwaBMEXAncJ2ZHWpmewD3A3OBmR1VF5F4saG4AlALnHQu3bLTuf/sfdh/ZD6XPTqXm19ZTKzGa4vEq45sgbsAP/P0FWBtyOPSkDLXAjcAtwCzgX7ANOdc6AJBFwOPAQ/jJz2UAEc3MJZOpNP7qgUuK8aRiERXl8w07jxjL47bYwB/eXkRv35iHrV1SuJE6nXYJAbnnEVQxgHTg0djZSrwC/1eGK3YRBLVxuJKMtJSyMuOyXwkkXaVkZbCX767GwV5mdz2+pdsKqnkbyfuQVZ64sxQFWkvunmiSALbUFxJ79xM/DrYIp1PSopx+ZHjuPLo8by0YD2n3vEe28qqYh2WSMwpgRNJYBuLKzX+TZLCWfsP4+aT9mDuqiK+8893WL2tPNYhicSUEjiRBLahuEIzUCVpfGvX/tz7g71ZX1TB8be+xbzVRbEOSSRmlMCJJDC1wEmy2W9EL/7zo/1INeO7/3yHF+eva/5JIp2QEjiRBFVVU8fWsmrNQJWkM7ZvHk/8ZH9G9+3K+f+ew22vL9EyI5J0lMCJJKhNJboLgySvgq5ZPHzuvhy1Sz/+9Pxn/N+jc6mqqYt1WCIdRmsPiCSoDboLgyS5rPRUbj5xD0bkd+GmV79gxZYy/nnqJLrnZMQ6NJF2pxY4kQS1rsjfhaFPnrpQJXmlpBg/nzaGG7+/Ox8u38Zxt77Nko0lsQ5LpN0pgRNJUPXLKAzokR3jSERi79t7DODBc/Zhe3k13/77W7yycH2sQxJpV0rgRBLU6q3lZKen0iMnPdahiMSFyUN78uRP9mdIfg5n3zubv81cTJ1uvyWdlBI4kQS1elsZA3pk6y4MIiEG9sjhv+dP4fg9BvDXmYs49/45bK+ojnVYIlGnBE4kQa3ZVsGA7uo+FQmXlZ7KX763G1cePZ7XPt/At295iy82FMc6LJGoUgInkqBWbyvX+DeRRpgZZ+0/jAd+uA9FZdUc8/e3ePyjVbEOSyRqlMCJJKCyqhq2lFapBU6kGfsO78UzPz2Aif27cfHDn3Dpfz6hrKom1mGJtJkSOJEEtCaYgTpQLXAizerXLZsHz9mHCw8eyaMfruLom99k4drtsQ5LpE2UwIkkoFVbgyVE1AInEpG01BQumTaGf5+9D0XlNXz7lrd44L3lugWXJCwlcCIJSGvAibTO/iPzef5nB7L3sJ5c8fg8zr53Nuu3V8Q6LJEWUwInkoBWby0nLcUo6Kq7MIi0VO+umdx71t789lvjeXvJJg674XUe/2iVWuMkoSiBE0lAK7b4NeBSU7QGnEhrpKQYPzhgGM/99EBG9enKxQ9/wrn3z2FjcI9hkXinBE4kAS3ZWMrw/C6xDkMk4Q3vncsj5+3Hr44ay+uLNjLtr6/z6By1xkn8UwInkmDq6hxLN5UwondurEMR6RRSU4xzvzGC5356AEPzu3DJfz7h+7e/y+frtPivxC8lcCIJZk1RORXVdQxXAicSVSMLuvLo+VP40/G7sGh9Md+86X/88bmFlFZq3TiJP0rgRBLMko2lAIzorS5UkWhLSTFO2nswr15SyAl7DuT2N77k0Bte5+lP1qhbVeKKEjiRBPPlxhIAtcCJtKOeXTL483d25dEf7Uf3nAwufOgjjrv1bT5YtiXWoYkASuBEEs6SjSXkZaWRn5sR61BEOr1JQ3ryzIUHcO13dmVtUTnf/ec7nHPfbJYE/0iJxIoSOJEEs2RDKcN752KmJUREOkJqivG9yYOYdelUfnH4GN5Zsplpf32DKx7/dMdt7UQ6WlqsAxCRllmysYQDRuXHOgyRpJOdkcqPp47kxL0GcdMri3nw/RU8Mnsl+/dPZcSuZQzqmRPrECWJqAVOJIFsLa1iQ3ElY/t2jXUoIkmrV24mVx07kVm/mMqJew3mzVU1TL1+Fpf99xOWby6NdXiSJNQCJ5JAFq7dDsC4fnkxjkREBnTP5nffnsieWRv4pKoPD72/gkc/XM0RE/vyg/2HMWlIj1iHKJ2YEjiRBLJACZxI3OmRlcL0IyZwQeEI7nxzKQ++v4Jn565l90Hd+cEBwzhyYl/SU9XhJdHVoe8oM/uGmT1lZqvNzJnZmWH7zcymm9kaMys3s1lmNiGsTKaZ3Wxmm8ysNDjewI6sh0isLFxbTH5uJvm5mbEORUTCFORlcflR43j38kO4+tgJFJVX89OHPuLAP7/Gza8sZm2RJjxI9HT0vwS5wDzgZ0BD7+TLgEuAC4G9gA3Ay2YWOuDnRuAE4CTgQCAPeMbMUtsvbJH4sHDtdsb10/g3kXjWJTON0/cbyis/P4i7zpzMyIJc/vLyIva/5lXOvucDXpq/jurauliHKQmuQ7tQnXPPAc8BmNk9ofvMr4lwEXCNc+7RYNsZ+CTuZOA2M+sGnA2c5Zx7OShzGrAcOBR4sUMqIhID1bV1fLGhhANHDY11KCISgZQU4+CxfTh4bB+Wby7lkdkr+c/sVbzy2QZ6d83k+D0HcNweAxjbV0MipOXiaQzcMKAv8FL9BudcuZm9AUwBbgMmAelhZVaa2cKgjBI46bRmfb6Rqto6jX8TSUBDenXhF4eP5eJDRzPr843M+GAFd/xvKbe9/iWj++Ry7O4DOGa3/lqKRCJmsbq3m5mVAD9xzt0T/DwFeAsY4pxbEVLuLmCAc+5wMzsZuA9IdyGBm9mrwGLn3HkNnOdc4FyAPn36TJoxY0Y71gpKSkrIzU3eWxwlc/3bs+7vrqnhtrmVFOQYv943m64Z8beIr1775Kw7qP6trf/2KscH62p4d00Ni7f5LtUR3VKY3DeNSX1SKciJ/4kPeu3bt/5Tp06d45yb3NC+eGqBqxeeUVoD28I1WsY5dztwO8DkyZNdYWFhW+Nr0qxZs2jvc8SzZK5/e9b9vns+YHCvEp776YF0yYzHj61e+2StO6j+ban/McHXlVvKeHruGp7+ZC0Pf76dhz+HsX27Mm1CXw6f0Ifx/fLi8u4reu1jV/94+kuwLvjaF1gZsr0AWB9SJhXIBzaGlXmjvQMUiZVPVxdx4Kj8uE3eRKRtBvXM4YLCkVxQOJIVm8t4acE6Xpq/nptfXcxNryxmYI9sDh3Xh8Ixvdl3eC+y0jVvL9nF01+DpfgE7TDgAwAzy8LPNP1FUGYOUB2UeTAoMxAYB7zdwfGKdIgN2yvYWFzJxP7dYh2KiHSAwb1y+OGBw/nhgcPZVFLJKwvX8+L89Tz0/grueXsZGWkp7DOsJweN7s1Bo3szskD3Rk5GHZrAmVkuMDL4MQUYbGa7A1uccyvM7EbgCjP7DFgE/BooIUjWnHNFZnYncJ2ZbQA2AzcAc4GZHVkXkY4yb00RABMHKIETSTb5uZl8f6/BfH+vwVRU1/Le0i28sWgjry/ayO+fXcjvn11I/25ZTBmZz37De7HviF4M6J4d67ClA3R0C9xk4LWQn68KHvcCZwLXAtnALUAP4D1gmnOuOOQ5FwM1wMNB2VeA051zte0dvEgszFu9HTMY31+zT0WSWVZ66o5Wt98Aq7eV+2Tu843MXLie/85ZBcCgntnsO6wX+yqh69Q6eh24WfgJB43td8D04NFYmQr8Qr8XRjc6kfj06eoihuV3IVfj30QkxIDu2Zy092BO2nswdXWOz9cX8+6Xm3n3y828vHA9/wkSur55Wew5pDt7Du7BHoO7M6F/N42h6wT0F0Ekjjnn+HRVEXsN6xnrUEQkjqWkGOP65TGuXx5n7T+MujrHZ+uKeX/pZj5csY0PV2zluU/9XMGM1BTG98/bkdBNHNCNIT1zSEnROLpEogROJI4t21zGuu0V7K0ETkRaICXFGN8/j/H98zhzf79tQ3EFHy7fxkcrtvLRim088N5y7nprKQC5mWmM7+fLT+ifx8QB3RhZkEt6avyvRZeslMCJxLG3vtgEwP4jesU4EhFJdAVdszhiYl+OmNgX8Lfn+3xdMfPXFDF/zXbmr9nOwx+spLzaDynPSEthTJ+uTOifx+g+XYNHLr27ZmrWaxxQAicSx95esol+3bIYlt8l1qGISCeTnprCxAHddprhXlvnWLqplPlriliwZjvz1hTxwvx1zPjgq+VZu2WnM7pPLqP6dMW2V5PxxSZG9elKfm6GErsOpAROJE7V1TneWbKZQ8b10UVRRDpEaooxsiCXkQX+/qzgx+JuLKlk8foSFq8vZtEG//XZuWspKq/mgYXvAZCXlcbQ/C4M7dWFofldGJafw9BeXRiW34XuORmxrFanpAROJE7NX7OdrWXV7D9S3aciEjtmRkHXLAq6ZrH/yPwd251zPPnia+SP2JVF64tZuqmUZZtL+XDFVp6Zu4a6kBtcdstOZ1DPbPp3y2ZAj2wGdM+m/45HFvldMjWJooWUwInEqWc/XUtailE4uiDWoYiIfI2Z0T0rhQNG5XPAqPyd9lXW1LJySznLgqRu6aZSVm8rZ9nmUt76YhOlVTsv3ZqRmkKv3Ax6dgl75GTQvUsGXTJSyclIJScjjZyMVLIzUklLSSHF/ISNFDP/fdBbUVPnqK6to6qmbsf3/uGoCb6vqnXU1tVRU+uorXPU1Pmv1bV1O/1c5xypZqSkGKkptuP7o3fr12G/64YogROJQ845npm7hgNG5dOji7oeRCSxZKal7uiKDeecY3tFDWu2le94rN5WwaaSSraWVrG5tIrlm8vYWlpFcWVNDKKPzOQhPWJ6fiVwInHo45XbWLW1nIsOHR3rUEREosrM6JadTrfsdMb1a/oOM5U1tRSVVVNWVRs8anZ8X986VucczrHjZ4dv0UtPTSEt1cgIvqanppAefE1L+er71BQjLdW3rqWlBD+n2I6vZrbTuWrrHHV10CUzlTeXdMzvrCFK4ETi0OMfrSYjNYVpE/rEOhQRkZjJTEulIC/2d41IjcPxeVqhTyTOlFTW8NiHq/nWrv3Iy0qPdTgiIhKHlMCJxJnHP1xFSWUNp08ZGutQREQkTimBE4kj5VW13PHmUnYb2I3dB3WPdTgiIhKnlMCJxJE/PLeA5ZvLuOyIsbEORURE4pgmMUTRn1/4jPveKiXvnVfITk8lK92vVbPz9yl0yUyjS0aa/5qZGnyfSpfMNHIy0sjN9Ovc5GamkZOZSmZa7AdwSvtas62cm15ZzIwPVnLuN4bvtFimiIhIOCVwUbTbwO7s2y+NXgX5lFfXUV5VS0V1LeXVtWwpraKiupbSqhrKKv3X0FWqm5Keajsldjsnfjt/n5Phk0WHX2unrs5R59gxzbrOOWqdo7rGUVPnFzWsrq3zCxvWOaqDRQ+rgm01tfXf+/JVIYsg1tQ6quvqdjpWXU01OW/PJD01Zcc07vS0+unbKWSnp9IlM5XsdB93doaPPXSBxvo61u/LTk8lMz2FrDT/NTMtJWq3lnLOUV3rqKippaLKv1ZlVbW8OH8dT32yBsPPgspIS6FrVhq7DexOt+z0YIFH/7vyvxvHsuWVvF48P1gIcuf9OxaQDBaNrN6xkGQdG4srWb+9khSD8w8awSXTtHSIiIg0TQlcFB0xsS9Zmz6jsHC3Zss656iorqOksoayqprgay2llTWUBgleabCtpLKGssoaSiprdyq7uaQsKOefV1lT16J4zfzNjNNTjLT6ZKt+jZxUIz3FJ1+h6+VkZ6SQkeq3ha6vkxYcZ+Xq1RT0KaAqSFiqa75KVKpr69hWXs2abeU71vMpraqlqoVxA2SkpZCVlkJmeipZ6SlkpqUGiV3973fn8rV1jqqaOip3PGp3/NyYA0flk5eVvqP8ltIqbp31xdcS7xSDtNQUzNWRsW5VsH6Q/52lmIX8Tv22tOD3nZ2RSl5qOmP75jG2b1cOHlvA8N5fX/RSREQknBK4GDEz36WakQpkRuWYNbV1lFbVUlldi4XcViTFDEthp1uN1C9eGG2zZm2isHDXFsddVl1LWeXOizTWt1aWV9dSWVNLRXXdTl8rQ75WBF9DhTbSmRmZaSk7WtMyQx4Zwfb67u7sjFRG9M5lTN+uX4u1orqWmjq3Y5HH+luq+LrPorCwsMW/MxERkZZSAteJpKWm0C07BbITa+2wtNQU8lJTEmLNs6x0jUcUEZHY0yxUERERkQSjBE5EREQkwSiBExEREUkwSuBEREREEowSOBEREZEEowROREREJMEogRMRERFJMErgRERERBKMEjgRERGRBKMETkRERCTBKIETERERSTBK4EREREQSjBI4ERERkQRjzrlYx9BhzGwjsLydT5MPbGrnc8SzZK5/Mtcdkrv+yVx3UP2Tuf7JXHdo//oPcc71bmhHUiVwHcHMZjvnJsc6jlhJ5vonc90hueufzHUH1T+Z65/MdYfY1l9dqCIiIiIJRgmciIiISIJRAhd9t8c6gBhL5vonc90hueufzHUH1T+Z65/MdYcY1l9j4EREREQSjFrgRERERBKMEjgRERGRBJOUCZyZ9TazW81smZlVmtl6M3vFzA7rgHMfb2YvmtlGM3NmVthAmVnBvtDHjCidP27rbmZDG6h3/eMXUYohbusflBlhZo8HZbab2SNm1ieKMcSk/maWbmZ/NrO5ZlZqZmvN7EEzGxxW7lwze83MtgW/o6FRjCHe6/4vM1tiZuXB6/+kmY2LYhzxXv9Od92LpO6d+brXgte+Xa97IecZYGa3m9kqM6sys9XB525gC45R/3o1u3SImfUws/vNrCh43G9m3cPK/M3MZptZhZkta0l9kjKBAx4F9gbOBkYD3wKeB3p1wLm7AG8DP2+m3N1Av5DHeVE6fzzXfSU717kfcAHggP9GKYa4rb+ZdQFeAgw4BNgfyACeNrNofVZjVf8cYE/gD8HXY4FBwAtmlhZW7iVgejvEEO91nw2cCYwDDse/D2aaWXqU4oj3+kPnu+5FUvfOfN1rtv4ddN3DzIbhP2MTgTOAkcCpwATgg2j+sxjiQXy9jwSOCL6/P6xMCnAvcF+Lj+6cS6oH0B3/wTi0mXIZwJ+BVUAp8AFweMj+wuA43wI+BiqAOcCkCOPID55f2MC+WcDfk7HuDZR9GXgpGeoPTAPqgB4h27oF25qMOZHqH3Kc8cFxdmlg3+Rg39DO9NpHUveQMrsGZcYkQ/3p5Ne9Fr72ne6611j9aefrXsgxnwNWAzlh23OC7c+GbDPgEmAxUBn8Tv4U7HNhj1mNnG9csH//kG0HNPaZBi4FlrWkTsnYAlcSPI4xs6wmyt0NHAScDOyCz5CfNrPdwspdD/wf/g/Ol8CzZpYThThPNLNNZjbfzK43s65ROGai1B3Y8R/TIURvmna81z8T/+GuCNlWgb+QHdCG49aLt/rnBV+3tuA5rZVQdQ9aJc4CVgDLWnDcxiRK/ZPhutfca9/Zr3vh9W/v6x5m1hPfAnaLc64sdF/w863AkWbWI9j8R+A3wJ/wLXTfxbeUgm/JJDheP+D4Rk67H/73/nbItrfwyfGUttQnNPikewAnAFvwb5J38G/IfUL2j8C/eQaHPe8J4Fa3838jp4TszwW2AT+MIIamWuDOxXeh7AKcCCwFXk6GuoeV+yOwAUhPhtce6A0UAX/Hd7d2AW4Oyt7WWeoflM/AX8yeamR/VFvgEqXu+K6zkuAcnwEjk6X+dPLrXiTv+6BMp7zuNVZ/Oua6t09wvOMa2X9csH/voE4VwPmNlB0alJ3czDl/BXzZwPYvgcsb2K4WuEg45x4F+gNH48cBTAHeNbNfBUX2xDehLjCzkvoH8E38Gz3UOyHHLQE+xTcRtyW+251zLzrnPnXOzQC+DxxqZnu25bjBseO67vWC8RFnAvc456qjcUyI7/o75zbi/9M7EijGX9S6Ax8Cta09btg5Yl7/4LX9N75uZ7WpQi2QIHV/ANgD3xKyCPhPtFq1473+nf26F8n7vjNf9xqrf0dc90K4xsIL2T8e3yr4Sjudz5qIo0XCB5AmDedcBX6cwcvA1WZ2BzDdzK7HDyp0wF5A+IeovEMD9Wbj38ij8G/qNkmQuh+Nb56+I9oHjuf6O+deAkaYWT5Q45zbZmbr8K0R0TpHzOofXMQfwreyFDrnNrf1mC0R73V3zhXh/4AtNrN38d1MJ/D1gc+tEu/1D9NprnstqHunvO41V/8OuO4txtdvAr5VMVz9eLUl+PdbNKwDCszMXNDEZmaGb3FcH40TJG0C14AF+N9HFvARPkvu65x7rZnn7YtvEq0ftzKR1swmadouQCqwNsrHrRePdT8HeN05tyhKx2tK3NXfObcpOO7BQAHwVDSO24gOqX8wm3JGUK7QObeu7aG3WTzX3YJHZgRlWyue698prnstrHunu+61pP7tdd1zzm0xsxeBC8zsry5kHFzQwv1j4Pmg3AL8xIVD8IlfuKrga2ozp30H3x27H1+Ng9uPr1YjaLto9C8n0gM/bfpV/PThXYFh+ObbdYSMt8A39S4HvgMMx4/JuRQ4PthfiM/YFwCH4TP7h/GZdZcmzt8T2D3k+T8Mfu4b7B8B/DY431DgKGAh/j/Q1M5c95Byg/H/eZ/SlvomYv3xXQv7Be+DU4HNwF8Svf74PxRP4Gd77Qn0DXlkh5TrG/xOTg7OcVTwc8/OXHf8kgb/B0zCv/+n4P94bSXs89FJ698pr3uR1D2kbKe77kVaf9rxuhdyjhFBrG8DB+OXMynEj8lbCwwLKftn/GfvrOB5ewM/CqlTWfB+7QN0a+Kcz+O7mPcN6vcp8HRYmZH4a9wNwJrg+92BjGbrFM1fUCI88P/N/hE/RXpr8EIsDn55PUPKpePXovoSn3Gvw19QJ4W9mY8B5uIz9g+BvZo5/5l8fRqyA6YH+wcBrwdv4ErgC+BvtPEPWCLUPaTcVfgBt1nJ9NoHZa4JzleFHwP1c4J7Fidy/flq4G9DjzNDyk1vrkxnrDv+c/88fvB6FX7G2wPA2GR47emk171I6h5SttNd9yKtP+143QuLZxDwL3xCWY1PmO4ABoaVSwF+GfK7WAn8IWT/D/EzxGtpZBmRoFxPfGK8PXj8G+geVmZWI7+foc3VRzezbyXzq+i/BvR2QbNvskjmuoPqn8z1T+a6Q3LXP5nrDqp/PErKWagiIiIiiUwJnIiIiEiCUReqiIiISIJRC5yIiIhIglECJyIiIpJglMCJiIiIJBglcCIiIiIJRgmciCQ9M7vHzFzwqDazDWb2mpn9OLgVUKTHKQyOkd+e8YqIKIETEfFm4m8kPhSYBjyNXx3/f8E9H0VE4oYSOBERr9I5t845t9o597Fz7gb87YP2BC4DMLNTzewDMysOWun+Y2YDgn1D8SvVA2wMWuLuCfaZmV1mZkvMrNzMPjWzUzu6giLSeSiBExFphHNuHvACcEKwKQO4EtgN+BaQDzwU7FsZUm4CvjXvZ8HPvwfOBn4MjAf+BNxmZt9s5yqISCeVFusARETi3ALgUADn3F0h2780sx8BC81soHNulZltCfZtqL9fZND9+nNgmnPuf8H+pWa2Nz6he7ZDaiEinYoSOBGRphngAMxsT3wL3O5Az2AfwGBgVSPPHw9kAS+YWeitb9KBZdEPV0SSgRI4EZGmjce3tnUBXsRPdjgN2IDvQv0fvmu1MfVDVY4GVoTtq45uqCKSLJTAiYg0wswmAkfgx7CNxSdsv3LOLQ32Hx/2lKrga2rItgVAJTDEOfdq+0YsIslCCZyIiJdpZn3xLWa9gUOAXwFzgOuBHHwi9hMzuwUYB/wu7BjL8d2t3zSzp4Fy51yxmV0PXG9mBrwB5AL7AnXOudvbv2oi0tloFqqIiHcosBbfzfkKcAx+HbhvOOdKnXMbgTOAb+Nb1a7ET07YwTm3Otj+B2A98Pdg12+A6cClwHzgZfyM1aXtWSER6bzMOdd8KRERERGJG2qBExEREUkwSuBEREREEowSOBEREZEEowROREREJMEogRMRERFJMErgRERERBKMEjgRERGRBKMETkRERCTBKIETERERSTD/D0ggwd11D+SyAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the discharge time series\n",
    "\n",
    "plt.rc(\"font\", size=14)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 6))\n",
    "ax.plot(sabine_timeseries_df, color=\"tab:blue\", label=\"Q\")\n",
    "ax.set(\n",
    "    xlabel=\"Date\",\n",
    "    ylabel=\"Discharge [cms]\",\n",
    "    title=f\"COMID {comid} from \"\n",
    "    f'{pd.to_datetime(begindate).strftime(\"%b %d %Y\")} to '\n",
    "    f'{pd.to_datetime(enddate).strftime(\"%b %d %Y\")}',\n",
    ")\n",
    "date_form = DateFormatter(\"%b %d\")\n",
    "ax.xaxis.set_major_formatter(date_form)\n",
    "ax.legend()\n",
    "ax.grid(True)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.5"
  },
  "vscode": {
   "interpreter": {
    "hash": "f495720fd4e62ac3edfedb951329873889a5dcf6e6c0cb1e5a7b0d2770d51a8a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}