{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/eolson/anaconda3/lib/python3.5/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.\n",
      "  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')\n",
      "/home/eolson/anaconda3/lib/python3.5/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.\n",
      "  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')\n"
     ]
    }
   ],
   "source": [
    "# imports\n",
    "from sqlalchemy import create_engine, Column, String, Integer, Float, Date, MetaData, Table, type_coerce\n",
    "from sqlalchemy.orm import mapper, create_session\n",
    "import csv\n",
    "from sqlalchemy import case\n",
    "import numpy as np\n",
    "from sqlalchemy.ext.automap import automap_base\n",
    "import matplotlib.pyplot as plt\n",
    "import sqlalchemy.types as types\n",
    "import numbers\n",
    "from sqlalchemy.sql import and_, or_, not_, func\n",
    "from sqlalchemy.sql import select\n",
    "import datetime as dt\n",
    "import os\n",
    "import re\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "import netCDF4 as nc\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# First, calculate river concentrations from EC data at Hope and Gravesend\n",
    "# NOTE: am assuming all forms of N are reported in grams of N, so that mole conversion is \n",
    "#  by division by 14.006720 (mol. weight N), rather than weights of full molecules. This \n",
    "#  is consistent with NO23 values in Susan's notebook.\n",
    "#  This is consistent with report in 'precdetectlimits.pdf', although it is rather old.\n",
    "#  However, assume Si in grams of SiO2 (mo. weight = 60.08). This is consistent with \n",
    "#  Si data in Susan's notebook (http://nbviewer.jupyter.org/url/bitbucket.org/salishsea/tools/\n",
    "#     raw/tip/I_ForcingFiles/Rivers/FraserRiverNutrients.ipynb?at=default&fileviewer=file-\n",
    "#     view-default)\n",
    "# Also appears to result in rough consistency with Debby's river Si and NO3 values\n",
    "mwN=14.006720\n",
    "#mwSiO3=76.083820\n",
    "mwSiO2=60.08\n",
    "#\n",
    "#  nh4 muM  (N) = nh4_cst from mean of data at Hope/Gravesend\n",
    "#  don muM  (N) = don_cst from mean of data at Hope/Gravesend\n",
    "#     - only 1 data point available. Should use it or set to zero? Sutton (2013) suggests DON may be unreactive on timescales\n",
    "#       of exchange of water in SoG, but Clair (2013) reports that organic nitrogen dominated total nitrogen in \n",
    "#       Canadian rivers.\n",
    "#       -> for now, using available data point\n",
    "#  no3 muM  (N) -> daily from smoothed data\n",
    "#  sil muM (Si) -> daily from smoothed data\n",
    "#  dia muM  (N) = 0.001\n",
    "dia = 0.0 #0.001\n",
    "#  phy muM  (N) = 0.001\n",
    "phy = 0.0 #0.001\n",
    "#  mes muM  (N) = 0.001\n",
    "mes = 0.0 #0.001\n",
    "#  zoo muM  (N) = 0.001\n",
    "zoo = 0.0 #0.001\n",
    "#  pon muM  (N) = 0.0\n",
    "pon = 0.0\n",
    "#  bsi muM (Si) = 0.0\n",
    "don=0.0\n",
    "bsi = 0.0\n",
    "#  oxy ? empty\n",
    "oxy = 160\n",
    "# other rivers const vals:\n",
    "const_si=55.69\n",
    "const_no=6.26\n",
    "const_nh=4.09\n",
    "const_tur=1.0 # default river turbidity (NTU)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# load EC Rivers database\n",
    "dbpath='/ocean/eolson/MEOPAR/obs/ECRivers/ECRiversDB'\n",
    "engine = create_engine('sqlite:///'+dbpath+'.sqlite')\n",
    "Base = automap_base()\n",
    "Base.prepare(engine, reflect=True)\n",
    "Profs=Base.classes.profiles\n",
    "session = create_session(bind = engine, autocommit = False, autoflush = True)\n",
    "# relevant data types: \n",
    "#('NITROGEN DISSOLVED NITRATE',)\n",
    "#('NITROGEN DISSOLVED NO3 & NO2',)\n",
    "#('NITROGEN DISSOLVED ORGANIC (CALCD.)',)\n",
    "#('NITROGEN NITRITE',)\n",
    "#('NITROGEN TOTAL',)\n",
    "#('NITROGEN TOTAL DISSOLVED',)\n",
    "#('NITROGEN TOTAL ORGANIC (CALCD.)',)\n",
    "#('SALINITY',)\n",
    "#('SILICA DISSOLVED',)\n",
    "#('SILICA EXTRACTABLE',)\n",
    "#('SILICA REACTIVE',)\n",
    "#('SILICON DISSOLVED',)\n",
    "#('SILICON EXTRACTABLE',)\n",
    "#('SILICON TOTAL',)\n",
    "#('TURBIDITY',)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define function to plot results of query as repeating fxn of year-day if output rows contain Year, Month, Day, value:\n",
    "def plotYMDV(qout,pcolor,ptitle):\n",
    "    Y=[]\n",
    "    M=[]\n",
    "    D=[]\n",
    "    Amm=[]\n",
    "    date=[]\n",
    "    for Y0, M0, D0, Amm0 in qout:\n",
    "        Y.append(Y0)\n",
    "        M.append(M0)\n",
    "        D.append(D0)\n",
    "        Amm.append(Amm0)\n",
    "        date.append(dt.date(Y0,M0,D0))\n",
    "    Y=np.array(Y)\n",
    "    M=np.array(M)\n",
    "    D=np.array(D)\n",
    "    Amm=np.array(Amm)\n",
    "    date=np.array(date)\n",
    "    YD=0.0*Y\n",
    "    for i in range(0,len(Y)):\n",
    "        YD[i]=date[i].timetuple().tm_yday\n",
    "    plt.plot(YD,Amm,'.',color=pcolor)\n",
    "    plt.plot(YD+365.0,Amm,'.',color=pcolor)\n",
    "    plt.title(ptitle)\n",
    "    plt.plot((366,366),(0,np.max(Amm)),'k--')\n",
    "    return(YD,Amm)\n",
    "    \n",
    "# definte function to smooth data with gaussian filter, plot result, and return dict of days, values\n",
    "def gsmooth(YD,Amm,L):\n",
    "    allt=np.arange(1,367)\n",
    "    fil=np.empty(np.size(allt))\n",
    "    s=L/2.355\n",
    "    sdict={}\n",
    "    for t in allt:\n",
    "        diff=[min(abs(x-t),abs(x-t+365), abs(x-t-365)) for x in YD]\n",
    "        weight=[np.exp(-.5*x**2/s**2) if x <= 3*L else 0.0 for x in diff]\n",
    "        #weight=[np.exp(-.5*x**2/s**2) for x in diff]\n",
    "        fil[t-1]=np.sum(weight*Amm)/np.sum(weight)\n",
    "        sdict[t]=fil[t-1]\n",
    "    plt.plot(allt,fil,'k-')\n",
    "    plt.plot(allt+365,fil,'k-')\n",
    "    diff_ex=[min(abs(x-t),abs(x-t+365), abs(x-t-365)) for x in allt]\n",
    "    weight_ex=[np.exp(-.5*x**2/s**2) if x <= 3*L else 0.0 for x in diff_ex]\n",
    "    plt.plot(allt,weight_ex/np.sum(weight_ex)*np.max(2.0*Amm),'k:')\n",
    "    plt.plot(allt+365,weight_ex/np.sum(weight_ex)*np.max(2.0*Amm),'k:')\n",
    "    return(sdict)\n",
    "    #plt.plot(YD,weight,'*')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "units: ('MG/L',)\n",
      "don_cst = 17.1346325192 uM N ....DON\"T USE\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEKCAYAAAD0Luk/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGcNJREFUeJzt3XmUHGW5x/HvLwREVkERJDFhVcTjggdjNGpGOUJEMB69\nKqDgcsWooKhXBfV4A+65inAVEYOIgktUFAmKGlTmXhYhQcKeQFgSsxFAQlhy8IbkuX+875Ci6Znp\nme5U9aR+n3PmTFfVW1VPv91TT9Vbb72jiMDMzOppVNUBmJlZdZwEzMxqzEnAzKzGnATMzGrMScDM\nrMacBMzMasxJwDYJSeMlbZDk71gLJN0t6fXtlpU0WdLSzkZnmzP/gXYRSYslrZW0RtIDkq6QNE2S\nGsq9StJfJD0kabWkiyS9oLB8cj4An9Gw3uWSjinr/QDDeghF0hsk/TW/v/skXSfp05K26nSAm6lS\nH/5plpQkvUfS5WXGYcPjJNBdAnhTROwIjAe+DpwInNNXQNIrgT8BFwLPAfYEbgSulLRHYVuPAkdL\nGteJwBoT0aYi6e3Ar4CfAOMiYhfgncBY4Ln9rLNFGbHZkPlJ1BHASaD7CCAiHo6I35EOgO+RtH9e\nPgP4UUScERGPRsSDEfEF4Grg5MJ2HgR+1DCv9SCkyyR9OV+NPArsKWkHSedIWiFpqaQv9SUHSaMk\nfTOfud8BvGk4+wVOBU6OiB9GxIMAEbEoIk6IiDvzvqZL+pWk8yU9SKqfl0u6Kl8ZLZf0HUmjc/kz\nJX2j4f39VtLH8+vnSLpA0r2S7pT00UK5l0ual6/OVkr6ZmHZRElX5n3OlzS5of6+mOvvIUl/lLRz\nYfnR+crvPkmfG0Y9TZB0i6R/5s+k6VVSviLcqzB9rqQvFqYPy7GvzrG+aBixDErSfrlOVku6SdLh\nDTF9T9KcXFeXFU9e8rpz8ntdkE8UrFMiwj9d8gPcDby+yfwlwDTg6cDjwOQmZd4LLM+vJwP/AJ4N\nrAH2zfMvB45pMZbLgMXAfqSThdGkq48zga2BZ5ESz7G5/IeAW4HdgWcAfwXWA6Py8u8Cq4EHCr/7\nXl+fy+yX1xk3SGzTgX8Bh+fppwEHABNISXQccAvwsbz8NcCSwvrPANYCu+by1wKfB7YA9gDuAN6Q\ny14FvCu/3gaYkF/vDtwPHJKnD8rTzyzU3yJg7xzfZcBX87L9gYeBScCWpMT3f80++wG+JzcW6voK\n4IvFz75Qdj2wV2H63ELZA4BVwIG5Ho7O294yL7+4yWfW93v2QN9b0vfxf/Pr0bkuTsyvXwc8xMbv\n5bmk72lffZwOXF6o838Ax+QYXwLcC+xX9d/r5vLjK4GRYQWwc/4ZBaxsUmYl6cD8hIi4FzgL+GKT\n8q34UUQsjIgNed9vBD4REY9FxP2kP9Yjctm3A6dHxIpIZ/Bfa4jluIjYKSJ2Lvzue/3SXOyZ+fc9\nfetJ+nk+e3xU0rsKm/xbRFyct/2viJgfEXMj+Qcwk3RAJCIuB0LSq/O6/wZcFRGrSInjWRHxlYhY\nHxGLgR8U3tc6YB9Jz4yItRExN89/N/D7iPhT3sdfSMnk0EKM50bEnRHxL+CXQN/7fBtwcURcGRHr\ngC8w9KaT7xTq+ivAkf2UG6gZ71jgrIi4Ntfb+aTkOjG/p8ObfGZ9v9/csK3f5vtYD0h6gJT0+7wS\n2DYiZkTE4xFxGfC7hph/X6iPzwMTJY0BDgPujojzcow3AL8hfd+sA5wERoYxbDwL20C6F9DoOaQz\n0UYzgEMkvXgY+y32MhlPOktbmf/QV5MSzC55+e4N5ZcMY3//zL+feH8RcWRE7ARcRzpTbxYbkvaV\ndHFusuk7MBaT4i/YeNA5Cvhpfj0OGFM4gK0GPku6igJ4P/B8YKGkayT1NXONB97RsN4kYLfCPu8p\nvF4LbJdfP6muImJt4b23alnh9ZK8zaEaD/xHw3sYO8xtTS0k9p2BjxSWPYeGzyvHPKYwXayPR0nf\n9d1zjBMbYjyKJ9eztWF01QHYwCS9nPTHcHlErJX0N9JZ0P80FH0H8OfG9SPiAUmnA19i6GebxfJL\ngcdIzR3NtrOSJ9+4Hd/wPr5HOntuXFfA4oh4EXAbsBx4K3DaEGID+B4pUbwz19MJpDPuPj8H/iRp\nBvAK4C2F93VXRDy/6U7SfYij8nt4G3BBbttfCpwXEdMGibOZlaSmL/J2t2HjVVCrGut6RT/l1pKa\nVPrsxsYD7lLgKxHxtaesleK6hNSU1uzzvjwiivd9BrriWMFTb+qPI33efZ5YLmk7YKe83lKgNyIO\nGWD71gZfCXQpSdtLOox08Do/Im7Ni04i3Qg9XtJ2knaS9GXSJfwp/WzuNOBVQLEbaV8//pZ6D0XE\nPcAc4LQcmyTtJem1ucgvgY9JGiNpJ1L7b3H9D0fE9hGxQ8PP9jkBkJPLp4Dpkv5d0jNyrPuS2u8H\nsj3wUE4A+wEfbtj/9aSz7R8Af4yIh/KiucDDkj4jaWtJW0h6oaQD877fJanvimIN6YC4gdR76XBJ\nByvdFN9aqWtuK2fRFwCHKXX13ZLUXPfEQTRvZ8Mg2zgu1/XOwOeAWf2Umw8clWOcQm4iy84GPiRp\nQt7vtpIOlbQtQEQc2s9ntkNDAhjMNcDaXMejJfWQmnl+XihzaK6PrUgnLFdHxHJSs9HzJL07r7ul\npAPzZ2wd4CTQfS6WtIZ0M+yzwDdJTRIARMSVwCGks9yVpJtyLwEmRcRdzTYYEQ8D/0Vq1+8zjnTj\nd3k/cTQ7+zsG2Ip0A/gBUlfOvsvys0ldV28gtY3/euC32c9OI35Juqo5GviHpPtIB7jv5/3151PA\nuyQ9lMs2Oyj+jHQDt68piHy/4zBSe/3dpJuOZwM75CJTgFvydk8jXWn8KyKWAVNJB+D7SM0bn2Lj\n31S/V105oR9HOgiuICWnYvPOc4ErB3ivkd/LHNJN7EWk5q9mPg68mdS8ciTp5n5fHH8n3Rc4I7fj\n3w68Z4D9DhRP/wtTO//hpPsl9wNnAEdHxKJCsZ+RerL9k3TD+t153UeAg0n3aFbkn6+TvofWAWp+\nZd9QKJ1BnE76gp8TETMalh/FxjO/h4GPRMSNedli0hnUBmBdREzoWPQ2bJI+D9wbEWdXHYs9maSZ\nwK8i4tKqYymDpHOBpRHxn1XHUkeD3hNQeuz/DNIZ1ApgnqSLImJhodhdwGsjYk1OGDPJPQxIB/+e\niFjd2dCtHRHR35mjVSwiPlh1DFYfrTQHTQAWRcSSfFk3i3QZ/ISIuDoi1uTJq3nyXX+1uB8zqyc/\nWVyhVnoHjeHJ3buWkRJDfz4A/KEwHcClktYDM938YGZFEfH+wUvZptLRLqKSXge8D3h1YfakiFgp\naRdSMlgQEVd0cr9mZjY8rSSB5aSeJH3G0qRHSX4YaSYwpdj+HxEr8+/7JF1Iuop4ShKQ5EtCM7Mh\nioi2Bndspa1+Humx+fG5D+8RwOxigdzX/Nekbl93FuZvkx/8IPc9Phi4ub8dRReMo1H8mT59euUx\nOKbNJ6ZujcsxjdyYOmHQK4GIWC/peFKf5L4uogskTUuLYyZp7JOdgTMliY1dQXcFLsxn+aOBn0bE\nnI5EbmZmbWvpnkBE/JE0fkpx3vcLr48lPXTSuN7dbBw0y8zMuoy7bg6gp6en6hCewjG1phtjgu6M\nyzG1phtj6oSWnhgug6TolljMzEYCSUQJN4bNzGwzVd8ksGIFzJyZfpuZDcVmdPyo5/8TWLEC9t4b\nHnsMtt4a7rwTdh/O/9Ews9rZzI4f9bwS+N3v0gcI6fcll1Qbj5mNHJvZ8aOeN4Y3s0xuZiXqouNH\nJ24M1zMJQPogL7kEDj3UCcDMhqZLjh9OAmZmNeYuojYinXzyyVWHYGaZrwSsdPnspeowzEY8XwmY\nmVlbnATMzGrMScDMrMacBMzMasxJwEo3ffr0qkMws8y9g8zMRij3DjIzs7Y4CZiZ1ZiTgJlZjTkJ\nmJnVmJOAlc5jB5l1D/cOstJ57CCzznDvIDMza4uTgJlZjTkJmJnVmJOAmVmNOQlY6Tx2kFn3cO8g\nM7MRyr2DzMysLU4CZmY15iRgZlZjTgJmZjXmJGCl89hBZt3DvYOsdB47yKwz3DvIzMza0lISkDRF\n0kJJt0s6scnyoyTdkH+ukPTiVtc1M7PqDNocJGkUcDtwELACmAccERELC2UmAgsiYo2kKcDJETGx\nlXUL23BzUE24OcisM8pqDpoALIqIJRGxDpgFTC0WiIirI2JNnrwaGNPqumZmVp1WksAYYGlhehkb\nD/LNfAD4wzDXtRrw2EFm3WN0Jzcm6XXA+4BXD2f9YtfBnp4eenp6OhKXdRd3ETUbnt7eXnp7ezu6\nzVbuCUwktfFPydMnARERMxrKvRj4NTAlIu4cyrp5me8JmJkNQVn3BOYB+0gaL2kr4AhgdkMg40gJ\n4Oi+BNDqumZmVp1Bm4MiYr2k44E5pKRxTkQskDQtLY6ZwBeAnYEzJQlYFxET+lt3k70bMzMbEj8x\nbGY2QvmJYRuRfGPYrHv4SsBK54fFzDrDVwJmZtYWJwEzsxpzEjAzqzEnATOzGnMSsNJ57CCz7uHe\nQWZmI5R7B5mZWVucBMzMasxJwMysxpwEzMxqzEnASuexg8y6h3sHWek8dpBZZ7h3kJmZtcVJwMys\nxpwEzMxqzEnAzKzGnASsdB47yKx7uHeQmdkI5d5BZmbWFicBM7MacxIwM6sxJwEzsxpzErDSeewg\ns+7h3kFWOo8dZNYZ7h1kZmZtcRIwM6sxJwEzsxpzEjAzqzEnASudxw4y6x7uHWRmNkK5d5CZmbXF\nScDMrMacBMzMaqylJCBpiqSFkm6XdGKT5c+XdJWkxyR9smHZYkk3SJovaW6nAjczs/YNmgQkjQLO\nAA4BXggcKWm/hmL/BD4KfKPJJjYAPRFxQERMaDNe2wx47CCz7tHKlcAEYFFELImIdcAsYGqxQETc\nHxF/Bx5vsr5a3I/VxCmnnFJ1CGaWtXJwHgMsLUwvy/NaFcClkuZJOnYowZmZ2aY1uoR9TIqIlZJ2\nISWDBRFxRQn7NTOzQbSSBJYD4wrTY/O8lkTEyvz7PkkXkpqXmiaBYltxT08PPT09re7GzGyz19vb\nS29vb0e3OegTw5K2AG4DDgJWAnOBIyNiQZOy04FHIuLUPL0NMCoiHpG0LTAHOCUi5jRZ108M14T/\nn4BZZ3TiieFBrwQiYr2k40kH8FHAORGxQNK0tDhmStoVuBbYHtgg6QRgf2AX4EJJkff102YJwOrF\nYweZdQ+PHWRmNkJ57CAzM2uLk4CZWY05CZiZ1ZiTgJlZjTkJWOk8dpBZ93DvICudnxMw6wz3DjIz\ns7Y4CZiZ1ZiTgJlZjTkJmJnVmJOAlc5jB5l1D/cOMjMbodw7yMzM2uIkYGZWY04CZmY15iRgZlZj\nTgJWOo8dZNY93DvISuexg8w6w72DzMysLU4CZmY15iRgZlZjTgJmZjXmJGCl89hBZt3DvYPMzEYo\n9w4yM7O2OAmYmdWYk4CZWY05CZiZ1ZiTgJXOYweZdQ/3DrLSeewgs85w7yAzM2uLk4CZWY05CZiZ\n1ZiTgJlZjTkJWOk8dpBZ93DvIDOzEaq03kGSpkhaKOl2SSc2Wf58SVdJekzSJ4eyrpmZVWfQKwFJ\no4DbgYOAFcA84IiIWFgo8yxgPPAWYHVEfKvVdQvb8JWAmdkQlHUlMAFYFBFLImIdMAuYWiwQEfdH\nxN+Bx4e6rpmZVaeVJDAGWFqYXpbntaKddc3MbBNz7yArnccOMuseo1sosxwYV5gem+e1YkjrFg8O\nPT099PT0tLgbG0lOOeUUJwKzYejt7aW3t7ej22zlxvAWwG2km7srgbnAkRGxoEnZ6cAjEXHqMNb1\njeGa8AByZp3RiRvDg14JRMR6SccDc0jNR+dExAJJ09LimClpV+BaYHtgg6QTgP0j4pFm67YTsJmZ\ndY4fFrPS+UrArDM8lLSZmbXFScBK57GDzLqHm4PMzEYoNweZmVlbnATMzGrMScDMrMacBMzMasxJ\nwErnISPMuod7B1np/LCYWWe4d5CZmbXFScDMrMacBMzMasxJwMysxpwErHQeO8ise7h3kJnZCOXe\nQWZm1hYnATOzGnMSMDOrMScBM7MacxKw0nnsILPu4d5BVjqPHWTWGe4dZGZmbXESMDOrMScBM7Ma\ncxIwM6sxJwErnccOMuse7h1kZjZCuXeQmZm1xUnAzKzGnATMzGrMScDMrMacBKx0HjvIrHu4d5CV\nzmMHmXWGeweZmVlbnATMzGrMScDMrMZaSgKSpkhaKOl2SSf2U+bbkhZJul7SAYX5iyXdIGm+pLmd\nCtzMzNo3erACkkYBZwAHASuAeZIuioiFhTJvBPaOiH0lvQL4HjAxL94A9ETE6o5HbyOSxw4y6x6t\nXAlMABZFxJKIWAfMAqY2lJkKnAcQEdcAO0raNS9Ti/uxmnAXUbPu0crBeQywtDC9LM8bqMzyQpkA\nLpU0T9Kxww3UzMw6b9DmoA6YFBErJe1CSgYLIuKKZgWLZ4g9PT309PSUEJ6Z2cjQ29tLb29vR7c5\n6MNikiYCJ0fElDx9EhARMaNQ5izgsoj4RZ5eCEyOiFUN25oOPBwR32qyHz8sZmY2BGU9LDYP2EfS\neElbAUcAsxvKzAaOyUFNBB6MiFWStpG0XZ6/LXAwcHM7AZuZWecMmgQiYj1wPDAHuAWYFRELJE2T\n9MFc5hLgbkl3AN8HPpJX3xW4QtJ84Grg4oiYswneh40gvjFs1j08dpCVzmMHmXWGxw4yM7O2OAmY\nmdWYk4CZWY05CZiZ1ZiTgJXOYweZdQ/3DjIzG6HcO8jMzNriJGBmVmNOAmZmNeYkYGZWY04CVjqP\nHWTWPdw7yErnsYPMOsO9g8zMrC1OAmZmNeYkYGZWY04CZmY15iRgpfPYQWbdw72DzMxGKPcOMjOz\ntjgJmJnVmJOAmVmNOQmYmdWYk4CVzmMHmXUP9w6y0nnsILPOcO8gMzNri5OAmVmNOQmYmdWYk4CZ\nWY05CVjpPHaQWfdw7yAzsxHKvYPMzKwtTgJmZjXmJGBmVmNOAmZmNeYkYKXz2EFm3aOl3kGSpgCn\nk5LGORExo0mZbwNvBB4F3hsR17e6bi7n3kE14bGDzDqjlN5BkkYBZwCHAC8EjpS0X0OZNwJ7R8S+\nwDTgrFbX7Wa9vb1Vh/AUjqk13RgTdGdcjqk13RhTJ7TSHDQBWBQRSyJiHTALmNpQZipwHkBEXAPs\nKGnXFtftWt34oTum1nRjTNCdcTmm1nRjTJ3QShIYAywtTC/L81op08q6ZmZWkU11Y7itNiozMyvH\noDeGJU0ETo6IKXn6JCCKN3glnQVcFhG/yNMLgcnAnoOtW9iG7xSamQ1RuzeGR7dQZh6wj6TxwErg\nCODIhjKzgeOAX+Sk8WBErJJ0fwvrAu2/ETMzG7pBk0BErJd0PDCHjd08F0ialhbHzIi4RNKhku4g\ndRF930DrbrJ3Y2ZmQ9I1o4iamVn5SnliWNJYSX+VdIukmyR9LM/fSdIcSbdJ+pOkHQvrfFbSIkkL\nJB1cQkwfzfOnS1om6br8M6XEmJ4m6RpJ83NcX83zq6yn/mKqrJ4K+xmV9z07T1dWTw0xzS/E1A31\ntFjSDTmuuXlepXXVT0yV1pWkHSX9Ku/jFkmvqLqeBoirc3UVEZv8B9gNeGl+vR1wG7AfMAP4TJ5/\nIvD1/Hp/YD6puWoP4A7yVUsJMU0HPtmk/As2dUx5P9vk31sAVwOTqqynAWKqtJ7yvj4B/ASYnacr\nrad+YuqGeroL2KlhXtXfqWYxVf239yPgffn1aGDHqutpgLg6VlelXAlExD2Rh5GIiEeABcBY0oNj\nP87Ffgy8Jb9+MzArIh6PiMXAItKDZ5s6pr5nGJrdpJ66qWPKsazNL59GulJbTYX1NEBMUGE9SRoL\nHAr8oGHfldVTPzFBhfVU2H/j33qlddVPTH3zG23yupK0A/CaiDgXIO9rDdV/p/qLCzpUV6UPICdp\nD+ClpDPKXSNiFaSDMvDsXKzxIbPlbMKHzAoxXZNnHS/pekk/KFz+lRJTX3MCcA/QGxG3UnE99RMT\nVFhPwGnAp4HiTa2qv0/NYoJq64kcz6WS5kn6QJ5XdV0VYzq2ML+qutoTuF/Subl5Zaakbai+nvqL\nCzpUV6UmAUnbARcAJ+Sz78Y/ltLvUjeJ6Uxgr4h4Kemgd2qZ8UTEhog4gHSl9BpJPVRcTw0xvVbS\nZCqsJ0lvAlblK7mBuhaXVk8DxFTp9ymbFBEvI12lHCfpNVT/t9cY06uptq5GAy8DvpvjehQ4ierr\nqTGutTmujtVVaUlA0mjSwfb8iLgoz16lNMYQknYD7s3zlwPPLaw+Ns/b5DFFxH2RG9eAs9l4KVVK\nTH0i4iHgEuBAKq6nhph+DxxYcT1NAt4s6S7g58DrJZ0P3FNhPTWL6bxu+D5FxMr8+z7gtzmGSr9T\nDTFdCEyouK6WAUsj4to8/WvSwbfqv73GuC4ADuhkXZV5JfBD4NaI+O/CvNnAe/Pr9wAXFeYfIWkr\nSXsC+wBzy4gpf9B93grcXFZMkp7Vd1kn6enAG0g3eSqrp35iur7KeoqIz0XEuIjYi/QA4l8j4mjg\nYiqqp35iOqbKegKQtE2+2kXStsDBwE1U+51qFtPNFX+nVgFLJT0vzzoIuIWKj1H9xHVrR+uq3TvX\nrfyQzpLWA9eTDmrXAVOAnYE/k3rmzAGeUVjns6Q72wuAg0uM6Tzgxjz/t6Q2wbJielGOYz5wA/Cp\nPL/KeuovpsrqqSG+yWzsiVNZPQ0QU6X1RGpT7vuO3wScVHVdDRBT1XX1EtIICdcDvyH1wqn8O9VP\nXB2rKz8sZmZWY/73kmZmNeYkYGZWY04CZmY15iRgZlZjTgJmZjXmJGBmVmNOAmZmNeYkYGZWY/8P\nqliT4q62AfwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe058e40940>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# calculate mean DON from data at Hope and Gravesend:\n",
    "\n",
    "qdataG=session.query(Profs.Year, Profs.Month, Profs.Day, Profs.value).\\\n",
    "            filter(and_(\n",
    "                Profs.variable_name=='NITROGEN DISSOLVED ORGANIC (CALCD.)',\n",
    "                Profs.station_name=='Fraser River (Main Arm) at Gravesend Reach - Buoy'))\n",
    "if qdataG.count()>0:\n",
    "    YDG, valG = plotYMDV(qdataG.all(),'r','DON, red=Gravesend, blue=Hope')\n",
    "    \n",
    "qdataH=session.query(Profs.Year, Profs.Month, Profs.Day, Profs.value).\\\n",
    "            filter(and_(\n",
    "                Profs.variable_name=='NITROGEN DISSOLVED ORGANIC (CALCD.)',\n",
    "                Profs.station_name.like('%Hope%')))\n",
    "if qdataH.count()>0:\n",
    "    YDH, valH = plotYMDV(qdataH.all(),'b',)\n",
    "else:\n",
    "    YDH=[]\n",
    "    valH=[]\n",
    "\n",
    "units=session.query(Profs.unit_code).\\\n",
    "            filter(and_(\n",
    "                Profs.variable_name=='NITROGEN DISSOLVED ORGANIC (CALCD.)',\n",
    "                or_(Profs.station_name=='Fraser River (Main Arm) at Gravesend Reach - Buoy',\n",
    "                    Profs.station_name.like('%Hope%')))).group_by(Profs.unit_code).all()\n",
    "for row in units:\n",
    "    print('units:',row)\n",
    "val=np.concatenate((valG,valH))\n",
    "don_cst=np.mean(val)\n",
    "print('don_cst =', don_cst/mwN*1000, 'uM N ....DON\"T USE')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "units: ('MG/L',)\n",
      "nh4_cst = 4.42713321781 uM\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEKCAYAAAD0Luk/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xu8XHV57/HPF5JIA1sEithESRUUXiRea2NOMbArR07E\nC57aU+P92GJzerTa01pBT9sEe3oqvXjp0baGQz2tRfFWMCpKaGUXKEGjBjQ3iCAYSERAJFs4tAl5\n+sdvDVmZzOyZ2bNmzZpZ3/frtV9zW5dnZj9rnrV+67d+o4jAzMzq6bBhB2BmZsPjImBmVmMuAmZm\nNeYiYGZWYy4CZmY15iJgZlZjLgJWOEmLJO2X5PzqgqTvSXpRv9NKOlPSzmKjs3HnjbRiJN0h6R5J\nP5V77tckXZN7vF/S05rmWy3p4y2Wd2Y2/XsHG/khZnUBiqQXS/qqpD2S7pX0LUm/K2le0QGOqVIv\n/GlVlCS9SdJ1ZcZhs+ciUD1B+r/8VovnW91vNw2S5gAfBG7sJyBJ6mf+HtbzX4DPAH8PnBgRxwOv\nBp4MPKXNPIeXEZv1zFehjggXgWr6U+B3JD2+zevdfin/DnAVsL2XlUu6RtL/knS9pIeAp0p6vKRL\nJO2StFPSHzaKg6TDJP1Ztuf+XeClvawv58+BNRHxNxHxY4CI2BER74iI27J1rZb0GUkfl/Rj4E2S\nfl7SDZIekHS3pP+TFUAk/aWkP216f1dI+q3s/s9I+qykH0q6TdJv5qb7eUkbJT0oabekP8u9tkzS\nv2Tr3CTpzKbP773Z57dH0lckHZt7/Q3ZEd+9kt4zi89pqaQtku7P/ictj5KajxglfSx/RCjpZVns\nD2SxPnMWsXQk6dTsM3lA0nckvbwppr+StD77rK6RdGLTvOuz97ot21GwArkIVNM3gCngd2e7AEmL\ngDcD76X7opH3euA8YAL4PvC3wL8CTwOeC7w4ex3g14FzgGcDzwd+uSmWj2RfAD/K3Tbu35RNcyqw\nEPiHLmJ7BfDpiHgCcCmwj3TkdCzwH4AXAf89m/aTwK/kYnkCcDbwyayIfQHYBPwMcBbwDkkvzib/\nEPDBiDgaOAn4dLaMBcAXgfdGxDHAO4HPSTouF+NrgDcBxwOPy6ZB0mnAXwKvAxYAx2XvuxevJX3+\nJwGnAL/XZrq2e+OSngtcAryF9Ll9FFgnaW72+hda/M8at+s6xPdYvmXF+AvAV0ifxduBSyU9ven9\nXEj6LG4m/U+RNB9YTzoy/GlgJfCRLFesKBHhvwr9Ad8jfYktBh4gbRi/Bnw1N81+4MfAj7K/B4CH\ngb/LTXMF8MvZ/Y+RvrC6jeEa0h554/ETgUeAx+WeWwn8U3b/n4Bfz732YuBR4LAe1nl6Ns+83HOf\nzN7bQ8DrsudWA1MdlvUO4HO5x3cAL8zunwf8Y3b/BcAdTfNeAFyS3f/nbH3HNU3zLuBvm577CvCG\n3Of3ntxrvwFcmd3/feATudfmk4rri3rIj7fkHr8E2JHdPxP4flOePC33+LE8IBWiC5uWvR1YPot8\n3ZPLxR9l/69rs9eXA7ua5vkE8Ae5mPKfx5HAXlJh/BXgn5vm/Wvg98vYFuvy5yOBioqILaS9zXe3\nmeS5EXFs9ncMcFHjhexweyIiPttHCPleJouAucDuxh4haWM8Pnt9QdP0d85iffdntz/TeCIiXpO9\nt28B+bb/g3rASHp6tue6O2si+iPSnmPDp0h75pD2Oi/N7p8ILMwfmZA+7ydmr/8qaU97u6SvSWo0\ncy0CfqVpvtOBJ+XW+YPc/YeBo7L7B31WEfFw7r13667c/TuzZfZqEanJMf8enjzLZZ2by8VjOXAU\nBun/2dxj6U4OPvrJfx4PkQr/gizGZU0xvpaDP2fr05xhB2AzWkP6AvzzFq/N1MTzIuDnJO3OHh8N\n7JP0zIj4z12uO9+UsJN0JHBcZLtjTXZz8InbRQcFKv0VqXmpeV6R9sSfCdwC3A38EvCBHmID+CvS\n5/TqiHhY0juAV+Ve/yRwlaSLSHv/r8y9r9sj4pSWK0nnIV6bvYdXAZ/N2vZ3ko66VnWIs5XdwGPN\nGVmTx3HtJ2+p+bPe1Wa6h0lHGg1P4sAX7k7gjyLij1vNKOlK0l58q//3dRGRP+8zUy7u4tCT+ieS\n/t8Nj70u6SjgmGy+naSjvv80w/KtTz4SqLDsS+hTpHbUXvwe8AxSG/2zgXXAxaRzBPl+/Ce2X8RB\ncfyA1Db7AUkTSp4m6Yxskk8Db5e0UNIxwPlN8/9GRExExOOb/iayAkBWXN4JrFbqEvuELNanAyd0\nCHEC2JMVgFNJzS/59d9E2tv+v8BXImJP9tLXgWlJ75J0hKTDJS2W9Pxs3a+T1DiieJD0hbif1Eb9\ncklnK50UP0KpK243e9GfBV4m6Rey9veDztlky9nfYRlvzT7rY4H3AJe1mW4T8NosxhWk5qKGi4H/\nJmlptt4jJZ0j6UiAiDinzf/s8U0FoJOvAQ9nn/EcSZPAy0iFueGc7POYB/whcGNE3E06En6GpNdn\n886V9HyfEyiWi0D1NO95vZe0N9dNF9H0YsRDEfHDxh/w/4GHIutxQ9rzuoO0591NDABvBOYBW0nt\nvp/hwGH5xaReSDeTTmp/bqb4Zoj706R24DcA35d0L+kL7qPZ+tp5J/A6SXuyaVt9KX6CdOK30RRE\nROwnfSE9h9S2/cPsvTR6Za0AtmTL/QDpSONfI+Iu4FzSF/C9pOaNd3Jge2r7/4mIrcBbSV+Cu0jF\nKd+88xTgX2Z4r5G9l/XAd4EdpOavVn6LdBL9AVJz2OW5OL5JOin8YUk/Am4lncjuVadc3Au8nNRx\n4D7gw6RzJztyk32CdNR7P6nTweuzeX9COom/kvRZ7QLeR8pDK4haH903TZT2Ij5ISvJLIuKiptdf\nQarg+0kn994VEV/NXruDtBe1H9gbEUuLfAPWO0n/E/hhRFw87FjsYJLWAp+JiKuHHUsZJH0M2BkR\nfzDsWOqqYxFQuvT/VtJe1C5gI7AyIrbnppmfneBCqa/x5RFxcvb4duDnIuKBwbwFMxtVLgLD101z\n0FJSF7Q7s0O7y0iHwo9pFIDMUaTDvgZ1uR4zqx9fWTxk3fQOWsjBXbzuIhWGg0h6JfDHpHbi/Nn8\nAK6W9Ciw1k0QZtYQEb867BjqrrAuohFxBXCFpBcCHyf1rwY4PSJ2SzqeVAy2RcT1Ra3XzMxmr5si\ncDepX2/Dk2nfq4SIuD7rznVcRNwfEbuz5++VdDnpKOKQIiDJh4VmZj2KiL4GeOymrX4jcHLWt3we\nqbvWQWOHSDopd/95WWD3S5qfXfxB1v/4bGBzuxUN+/LpTn+rV68eegyO03E6TsfZ+CtCxyOBiHhU\n0ttI/ZIbXUS3SVqVXo61wKskvRH4N9K4Ia/OZj8BuDzby58DXBoR6wuJ3MzM+tbVOYGI+AoH2vgb\nz300d/9PgD9pMd/3SBfimJlZBbnrZg8mJyeHHUJXHGexHGexHGe1dHXFcBkkRVViMTMbBZKIEk4M\nm5nZmHIRMDOrMRcBM7MacxEwM6sxF4FeTU/Dhg3p1mwQnGNWIheBXkxPw/LlcMYZ6dYbqRXNOWYl\ncxHoxebNsGUL7NsHW7em+2ZFco5ZyVwEerFkCSxeDHPnwmmnpfs2MtasWTPsEDpzjlnJfLFYr6an\n097Z4sUwMTHsaKwH2YU1ww6jM+eYdamIi8VcBKw2RqYImHXJVwybmVlfXATMzGrMRcDMrMZcBKw2\nVq9ePewQzCrHJ4bNzEaUTwybmVlfXATMzGrMRcDMrMZcBMzMasxFwGpjJMYOMiuZewdZbXjYCBs3\npfUOkrRC0nZJt0o6v8Xrr5B0s6RNkr4h6UXdzmtmZsPT8UhA0mHArcBZwC5gI7AyIrbnppkfEQ9n\n958JXB4RJ3czb24ZPhKwgfKRgI2bso4ElgI7IuLOiNgLXAacm5+gUQAyRwH3dTuvmZkNTzdFYCGw\nM/f4ruy5g0h6paRtwJXA23uZ18zMhmNOUQuKiCuAKyQtBz4OnNLrMvK9NyYnJ5mcnCwqPDOPHWQj\nb2pqiqmpqUKX2c05gWXAmohYkT2+AIiIuGiGeW4jNQU9vdt5fU7AzKw3ZZ0T2AicLGmRpHnASmBd\nUyAn5e4/DyAi7u9mXjMzG56OzUER8aiktwHrSUXjkojYJmlVejnWAq+S9Ebg34CHSF/2becd0Hsx\nM7Me+WIxM7MR5aGkzcysLy4CVhseO8jsUG4OstrwFcM2btwcZGZmfXERMDOrMRcBM7MacxEwM6sx\nFwGrDY8dZHYo9w4yMxtR7h1kZmZ9cREwM6sxFwEzsxpzETAzqzEXAasNjx1kdij3DrLa8NhBNm7c\nO8jMzPriImBmVmMuAmZmNeYiYGZWYy4CVhseO8jsUO4dZGY2otw7yMzM+uIiYGZWY10VAUkrJG2X\ndKuk81u8/lpJN2d/10t6Vu61O7LnN0n6epHBm5lZf+Z0mkDSYcCHgbOAXcBGSZ+PiO25yW4HzoiI\nByWtANYCy7LX9gOTEfFAsaGbmVm/ujkSWArsiIg7I2IvcBlwbn6CiLgxIh7MHt4ILMy9rC7XYzZQ\nHjvI7FDdfDkvBHbmHt/FwV/yzc4Dvpx7HMDVkjZKekvvIZoV48ILLxx2CGaV07E5qBeSfhF4M/DC\n3NOnR8RuSceTisG2iLi+1fz5PbXJyUkmJyeLDM/MbKRNTU0xNTVV6DI7XicgaRmwJiJWZI8vACIi\nLmqa7lnA54AVEXFbm2WtBqYj4v0tXvN1AjZQHkXUxk1Z1wlsBE6WtEjSPGAlsK4pkBNJBeAN+QIg\nab6ko7L7RwJnA5v7CdjMzIrTsTkoIh6V9DZgPaloXBIR2yStSi/HWuD3gWOBv5QkYG9ELAVOAC6X\nFNm6Lo2I9YN6M2Zm1hsPG2G1sWbNGvcQsrFSRHOQi4CZ2Yjy2EFmZtYXFwEzsxpzETAzqzEXATOz\nGnMRsNpwzyCzQ7l3kI2G6WnYvBmWLIGJiVktwlcM24wKyLGyuXeQ1cP0NCxfDmeckW6np4cdkY2b\nGueYi4BV3+bNsGUL7NsHW7em+2ZFqnGOuQhY9S1ZAosXw9y5cNpp6b5ZkWqcYz4nYKNhejrtnS1e\n7HMCNhgF5FjZijgnUOjvCZgNzMQELFvWeboZrF69uqBgbCwVkGOjyEcCZmYjyr2DzMysLy4CZmY1\n5iJgZlZjLgJmZjXmImC14bGDzA7l3kFWG75OwMaNewfZ+Jmehg0bajV2i5XI+XUIFwGrjhoP4mUl\ncH615CJg1VHjQbysBM6vllwErDpqPIiXlcD51VJXRUDSCknbJd0q6fwWr79W0s3Z3/WSntXtvGaP\nmZiA666Da69NtwUP4uWxg2puwPk1qjr2DpJ0GHArcBawC9gIrIyI7blplgHbIuJBSSuANRGxrJt5\nc8tw7yAzsx6U1TtoKbAjIu6MiL3AZcC5+Qki4saIeDB7eCOwsNt5zcxseLopAguBnbnHd3HgS76V\n84Avz3JeMzMrUaG/JyDpF4E3Ay+czfz5KzonJyeZnJwsJC4zs3EwNTXF1NRUocvs5pzAMlIb/4rs\n8QVARMRFTdM9C/gcsCIibutl3uw1nxMwM+tBWecENgInS1okaR6wEljXFMiJpALwhkYB6HZes7J4\n7CCzQ3U1dlDW4+dDpKJxSUS8T9Iq0l79WkkXA78E3AkI2BsRS9vN22YdgzsSmJ5OF4osWeJuYTU2\n0LGDnGM2BEUcCYz/AHKNS8UbPyDt/sG1NbAi4ByzIfEAct3wpeI2aM4xG2HjXwR8qbgNmnPMRtj4\nNwdBOlxvHKr7ML22Bn5OwDlmJSuiOajQ6wQqa2ICli0bdhQ2ZAMdO8g5ZiOqHkcCZmZjyCeGzcys\nLy4CZmY15iJgZlZjLgJmZjXmImC14bGDzA7l3kFWGwO9TsBsCNw7yMzM+uIiYLM3PQ0bNqRbs6I5\nv0rhImCz0xg584wz0q03VCuS86s0LgI2Ox450wbJ+VUaFwGbnREcOXOgYwdZsUYwv0aVewfZ7Hnk\nTBsk51dH/mUxM7MacxdRMzPri4uAmVmNuQiYmdWYi4DVhscOMjuUTwxbbXjsIBs3pZ0YlrRC0nZJ\nt0o6v8Xrp0i6QdIjkn676bU7JN0saZOkr/cTrJmZFavjD81LOgz4MHAWsAvYKOnzEbE9N9n9wG8C\nr2yxiP3AZEQ8UEC8ZmZWoG6OBJYCOyLizojYC1wGnJufICLui4hvAvtazK8u12NmZiXr5st5IbAz\n9/iu7LluBXC1pI2S3tJLcGZmNlgdm4MKcHpE7JZ0PKkYbIuI61tNmO+9MTk5yeTkZAnhWV147CAb\ndVNTU0xNTRW6zI69gyQtA9ZExIrs8QVARMRFLaZdDUxHxPvbLKvt6+4dZGbWm7J6B20ETpa0SNI8\nYCWwbqa4cgHOl3RUdv9I4Gxgcx/xmplZgTo2B0XEo5LeBqwnFY1LImKbpFXp5Vgr6QTgG8AEsF/S\nO4DTgOOByyVFtq5LI2L9oN6MmZn1xheLmZmNKI8iamZmfXERsNrw2EFmh3JzkNWGxw6ycePmIDMz\n64uLgJlZjbkImJnVmIvAbExPw4YN6dasaM4vK5GLQK+mp2H5cjjjjHTrDXVkjMTYQc4vK5l7B/Vq\nw4a0ge7bB3PnwrXXwrJlw47KxoXzy3rg3kHDsGQJLF6cNtDTTkv3zYri/LKS+UhgNqanYcuWtIFO\nTAw7Ghs3zi/rUhFHAi4CZmYjys1BZmbWFxcBqw2PHWR2KDcHWW147CAbN24OMjOzvrgImJnVmIuA\nmVmNuQiYmdWYi4DVxkiMHWRWMvcOMjMbUe4dZGZmfXERMDOrsa6KgKQVkrZLulXS+S1eP0XSDZIe\nkfTbvcxrZmbD0/GcgKTDgFuBs4BdwEZgZURsz03z08Ai4JXAAxHx/m7nzS3D5wTMzHpQ1jmBpcCO\niLgzIvYClwHn5ieIiPsi4pvAvl7nHRu7dsHatenWKmnkxw5yjtkAdFMEFgI7c4/vyp7rRj/zjo5d\nu+Ckk2DVqnTrjbSSLrzwwmGHMHvOMRuQOcMOIC+/pzY5Ocnk5OTQYunJF78IjzyS7j/yCFx5JZx3\n3nBjsvHiHDNgamqKqampQpfZzTmBZcCaiFiRPb4AiIi4qMW0q4Hp3DmBXuYd3XMCjb20Rx6BI46A\n226DBQuGHZU1GelRRJ1j1kJZ5wQ2AidLWiRpHrASWDdTXH3MO5oWLEgb5cUXe+O0wXCO2YB0dcWw\npBXAh0hF45KIeJ+kVaS9+rWSTgC+AUwA+4GfAKdFxE9azdtmHaN7JGAjYaSPBMxaKOJIoKtzAhHx\nFeCUpuc+mrt/D/CUbue1NqanYfNmWLLEPzA+ALUfO8j5ZS147KCqmJ6G5cthyxZYvBiuu84bqhXH\n+TWWPHbQONm8OW2g+/bB1q3pvllRnF/WhotAVSxZkvbQ5s6F005L982K4vyyNtwcVCXT0wcO10f1\nUN3tztXl/Bo7RTQHuQhYcTq1O3sDtn50c16jZjnmcwKjZHoaNmxIt+NqpnbnxgZ8xhnpdgifw8iP\nHdTJuOdYp/MaFcixUeQiUIZhJmeZXwwztTtX4MTkSI8d1Mmwcqwq+QWVyLFR5CJQhmElZ9lfDBMT\n6RD92msPPVT3icnBGkaOVSm/wDk2Sz4nUIbGxrJ1a0rOsvpob9iQNtB9+9KGce21sGzZ4NfbzpBP\nTI71FcPDyLGq5RcMPcfK5hPDo2QYyTms4lNRY10EoPwcc34NnYuAdTYOe0YF9fgY+yIwDOOQXzCy\nvYrcO8g6m5hIh+gjlNgHKbDdufZjBw3CqOcX1L5XkY8ErNqq2O5s42WEc8xHAjb+3OPDBq3mOeYj\ngbKMaJtjJYxLu/OgOcdmb0RzzCeGR4WH8bVBc47VkpuDRoWvZLRBc47ZLLkIlGHc2hxHdIyasR47\naJxybETza1S5Oagsw7pYrOg24hFudhj76wTGIcdGOL+Gwc1BVdRuL6bs/tSD6vvsZofhG+ccc36V\nzkWgSFW66GRQG9M4NTuMonHPMedX6VwEilSlvZhBbUydRnK0wRr3HHN+lc7nBIpUtQG1+mkjHmSf\n86KW3eNyxuKcgHOs3GVX/NqLIs4JEBEd/4AVwHbgVuD8NtP8BbADuAl4bu75O4CbgU3A12dYR4yF\nPXsiNmxIt6Noz56Iq6+OWLIkYs6ciGc/+9D3smdPxA03zO497tmTltlu2QNczurVq2e3rqpxjnVe\n/pByrGzZ92ZX3+Pt/ropAIcB3wUWAXOzL/lTm6Z5CfCl7P4LgBtzr90OHNPFegb4UVlX8kkP6W/u\n3PSF02qa2WwYN9xwYPnNyx7GcqxczrFCFVEEujknsBTYERF3RsRe4DLg3KZpzgX+Lvsm/xpwtKQT\nsteEzz2Mhnx7M8Dhhxf/M5FFtSP7BOJoco5VTjdfzguBnbnHd2XPzTTN3blpArha0kZJb5ltoNZB\nERfY5JN+yRK46qrifyayqBN/PoFYPufY7JZTcXNKWMfpEbFb0vGkYrAtIq5vNWH+is7JyUkmJydL\nCG8MFHWBTSPpZzrR18003ayniKF6i1qOdeYcq4SpqSmmpqYKXWbH3kGSlgFrImJF9vgCUjvURblp\n/hq4JiI+lT3eDpwZEfc0LWs1MB0R72+xnugUy9BVtadAN+OhVzV2O1hV/0+dcqyqcY+5sq4Y3gic\nLGmRpHnASmBd0zTrgDdmQS0DfhwR90iaL+mo7PkjgbOBzf0EPDTdXqQzyHFP2i270+FzlS4wGqLK\njx00qjnm/Bpt3Zw9JnURvYXUBfSC7LlVwK/npvkwqRfRzcDzsueeSupNtAn4TmPeNuso8qR58brp\nKTDILmWdlj1Tt8ER6OVQBudYB7PNMefX0FBGF9Gy/iq/gTY2kLlz2298g9wY+ll2N7HPNG83/bX7\n6dddEudYB7Nddj/51Zh/THKsbC4CZet0kU6/G0Ondfe7ofV6gVG3e50jcFFNxAgUgYjRzbHZXsA2\nZjlWtiKKgIeNKNogh/Ntt+xBnZTr9ge4R+SHusdi2AhwjlU4x8rmn5e0wY6/3u04NVUbz6aNsSkC\nZXOOVZZ/T8Bmd3Vlp94ljdehu4tlRuSimtWrVw87hNHkHBtrPhIYdb3uIXXaq/MvO1kz51hl+Uhg\nHLTbY+q2L3ive0id9uqqNF69FaNVLvVyrYFzbKy5CAxTu4tser34ppefFex0YVlNBs2qjVa5NJuL\nu5xjY8vNQWVq7mHRrsfDoHtCdOpdsmsXfOlL8NKXwoIFrWO36mn1P2qVSxGD72nTa445v2altB+V\nKeOPUejD3Y9W/Zzb9cseVF/wbi62mSlO99Gurnb/o1a5NOhrDXrNsbvvdn7NEr5YbIS0uxqz3UU2\n/Vx802oj7PaLvFWcYzIswNj8slgrM/2PWuVSP79OVnSOrV07Fvk1DC4Co2SQe1/N62i1EXb7RV72\nnmOJxjrHyvofDSLHGkcCI55fw1BEEfA5gTLN9krPbttLZzqX0Es3v1ZxDvIq1ZKM/cViZfzo+6By\nbAzyaxh8xfA46LTx9dKnutNGWPMNbeyLQDvOsbHlIjDqZtr4GhvuQw/BS15yYM/rIx+BlSt724s3\noKZFoF2O5QvD5s0H7907x0aGi8Coa3dond9wTz01TXvLLelHuffuTRuur7LsWS2LQKscW7z44MJw\n5ZVwzjlp7945NlJ8xfCoa3fRTP6KyltugQ9+MO2d7d0Ljz46+Kssd+2CtWvT7Rip5dhBrXKs+Yrd\n738/feE7x2rJRwLD1u4kbHO7K5QziuKuXXDSSfDII3DEEXDbbQcuGLPR1OokbKtcKmukTudYYdwc\nNM6G1UNn7VpYterA44svhvPOG8y6bHhm+t0A59jIcBGoizIvqfdeWj05x0aSzwnUwWwG++rHggVp\no7z4Ym+cdeEcqzUfCVSdf1bPBs05NrJ8JFAHS5akbqKHHw6nnOJhd/uwZs2aYYdQTc6xWuuqCEha\nIWm7pFslnd9mmr+QtEPSTZKe08u8M+nlty/KmH5o1N9osQYXXnhhy+cHmTMjk1/gHKurToMLkQrF\nd4FFwFzgJuDUpmleAnwpu/8C4MZu580t45DBkXodwXjQ019zzTUzTzAIsxjBcyhxzkLZcc42x/Jx\n9pIzZY/APevPs+RRYp2fxaGAAeS6ORJYCuyIiDsjYi9wGXBu0zTnAn+XfZN/DTha0gldztvWTL9C\n12oPq9dfret1+qmpqW5DL84sfoVpKHHOQhXi7CbHrrpqqqvpe1n2IMz68yz5l76q8H/vxqjE2a9u\nisBCYGfu8V3Zc91M0828bbXLzXadGXrN5ZH4lbtef9/VetJNjn3sY7PLsZHIL3CO1dycAS23kMbF\nRm42X7vSag+r8fOnrabvdfmV0/h9VytcNzl2772zy7GRyS9wjtVYxy6ikpYBayJiRfb4AlI71EW5\naf4auCYiPpU93g6cCTy107y5Zbh/qJlZj6LPLqLdHAlsBE6WtAjYDawEXtM0zTrgrcCnsqLx44i4\nR9J9XcwL9P9GzMysdx2LQEQ8KultwHrSOYRLImKbpFXp5VgbEVdKOkfSd4GHgDfPNO/A3o2ZmfWk\nMlcMm5lZ+YZ+xXC/F5MVHMslku6R9O3cc8dIWi/pFklXSTo699q7swvktkk6u6QYnyzpq5K2SPqO\npLdXNM7HSfqapE1ZrP+7inHm1n2YpG9JWlfVOCXdIenm7DP9eoXjPFrSZ7L1bpH0gqrFKekZ2ef4\nrez2QUlvr2Cc784+w29LulTSvMJj7PdCg37+6OFispLieSHwHODbuecuAt6V3T8feF92/zRgE6lJ\n7Wez96ESYnwS8Jzs/lHALcCpVYszW/f87PZw4Ebg9CrGma3/fwB/D6yr4v89W/ftwDFNz1Uxzv8H\nvDm7Pwc4uopx5uI9DNgFPKVKcZK+F28H5mWPPwW8qegYS/ug27zJZcCXc48vAM4fckyLOLgIbAdO\nyO4/CdjeKlbgy8ALhhDvFcB/rHKcwHzg61mSVi5O4MnA1cAkB4pAFeP8HnBc03OVihN4PHBbi+cr\nFWdTbGe+wVHcAAACvklEQVQD11UtTuCYLJ5jsi/2dYPY1ofdHNTXxWQleWJE3AMQET8Anpg93xz7\n3ZQcu6SfJR253EhKikrFmTWxbAJ+AExFxNYqxgl8APhdIH+CrIpxBnC1pI2SGr/CUrU4nwrcJ+lj\nWVPLWknzKxhn3quBT2T3KxNnRDwA/Dnw/Wx9D0bEPxYd47CLwCiqxJl0SUcBnwXeERE/4dC4hh5n\nROyPiOeS9rSXS5qkYnFKeilwT0TcxMwXOQ798wROj4jnAecAb5W0nIp9nqQ91ucBH8lifYi0h1q1\nOAGQNBd4BfCZ7KnKxCnpaaRmykXAAuBISa9rEVNfMQ67CNwNnJh7/OTsuSq5R2kcJCQ9Cfhh9vzd\npDbEhtJilzSHVAA+HhGfr2qcDRGxB7gSeH4F4zwdeIWk24FPAi+S9HHgBxWLk4jYnd3eS2oGXEr1\nPs+7gJ0R8Y3s8edIRaFqcTa8BPhmRNyXPa5SnM8H/iUifhQRjwKXA79QdIzDLgKPXYgmaR7pYrJ1\nQ45JHLxHuA74r9n9NwGfzz2/Mjtb/1TgZFK7dxn+BtgaER+qapySfrrRa0HSTwEvJp20qlScEfGe\niDgxIp5Gyr+vRsQbgC9UKU5J87OjPyQdSWrH/g7V+zzvAXZKekb21FnAlqrFmfMaUvFvqFKctwDL\nJB0hSaTPcmvhMZZ5AqbNyY8V2ZvdAVww5Fg+Qeol8K+kdrg3k07K/GMW43rgCbnp3006A78NOLuk\nGE8HHiX1pNoEfCv7DI+tWJzPzGLbBNwMvDN7vlJxNsV8JgdODFcqTlJbe+N//p3GtlK1OLP1Ppu0\ng3cT8A+k3kFVjHM+cC8wkXuuUnGSzlVtAb4N/C2pF2WhMfpiMTOzGht2c5CZmQ2Ri4CZWY25CJiZ\n1ZiLgJlZjbkImJnVmIuAmVmNuQiYmdWYi4CZWY39O6RxQe2yEhdDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe080294978>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# calculate mean NH4 from data at Hope and Gravesend\n",
    "\n",
    "dataG=session.query(Profs.Year, Profs.Month, Profs.Day, Profs.value).\\\n",
    "            filter(and_(\n",
    "                Profs.variable_name=='AMMONIA DISSOLVED',\n",
    "                Profs.station_name=='Fraser River (Main Arm) at Gravesend Reach - Buoy')).all()\n",
    "YDG, valG = plotYMDV(dataG,'r','temp')\n",
    "    \n",
    "dataH=session.query(Profs.Year, Profs.Month, Profs.Day, Profs.value).\\\n",
    "            filter(and_(\n",
    "                Profs.variable_name=='AMMONIA DISSOLVED',\n",
    "                Profs.station_name.like('%Hope%'))).all()\n",
    "YDH, valH = plotYMDV(dataH,'b','NH4, red=Gravesend, blue=Hope')\n",
    "\n",
    "units=session.query(Profs.unit_code).\\\n",
    "            filter(and_(\n",
    "                Profs.variable_name=='AMMONIA DISSOLVED',\n",
    "                or_(Profs.station_name=='Fraser River (Main Arm) at Gravesend Reach - Buoy',\n",
    "                    Profs.station_name.like('%Hope%')))).group_by(Profs.unit_code).all()\n",
    "for row in units:\n",
    "    print('units:',row)\n",
    "val=np.concatenate((valG,valH))\n",
    "nh4_cst=np.mean(val)\n",
    "print('nh4_cst =', nh4_cst/mwN*1000, 'uM')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "units: ('MG/L',)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEKCAYAAAD0Luk/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXt8XFW5//95Zk+SlhJ6oW3SpE3TENomHRsFbcu1kaKU\ni6LAEcQLB+WiIAjqsUX9noIcVDj+jrZegFYUFaUgCo0K2ERbbgm3Ugq9pS3QC9Nyr90pvSb5/P5Y\ne8+s2bNnMpNMJpOZ9X698srsmX1Ze+/PWs9az1rrWUISBoPBYChMAgOdAIPBYDAMHMYIGAwGQwFj\njIDBYDAUMMYIGAwGQwFjjIDBYDAUMMYIGAwGQwFjjIChXxCRFSLypYFOx2BARBaIyO8zsa+IvCYi\np2UudYZ8xxiBHENEtorImyIyVPvuyyKywrPff4nIJhF53znmByJSrP1+nYi8IiK2iLwhIr8WkSOz\neS+9QUTKRWSxiISdtG9x0j55oNPWz6QzYSdnJvckMkoi0i0iNQORJkN6GCOQexDqvVzn8z0AQER+\nBuAyAJ8HUArgTABzANyv7b8MwIdJHgVgKoCJAL7bmwSJiNWb43pxnVEAWgEMBXCSk/bjADwG4GMD\nmTZDUvyMUs4YKkNyjBHITf4XwDdF5CjvDyJyLICvAriY5LMku0luAHA+gLki0ggAJF8juds5zALQ\nDWBXKhcXkUtE5EkR+T8ReQfAAuf7L4nIehF5V0QeEZEq7ZiPicgGEdntGCnpxX1/A8Aekl8gudW5\nD5vkb0n+wrnORKeW+SUR2Qbgn87394vILuf6K0Wk3vl+hvN9JD0i8mkRWeN8FhGZ77Q43haRpSIy\nwvmtRER+LyLvOOd9RkTGOL8dJSK/EpGdIrJDRG52r+E8vydE5H9F5D2nRTZXu361k8Y9IvIPAKPT\nfE5DnXTaIvK8iEz320lEfiMi39e2Z4vIDm17nIg8ICJvOWm8Js10JEN/3sUi8lOndfe6iPxERIr0\nNInIDc7zf1VELvYc+2MR2ea8x1+KSEkG01nwGCOQmzwPYCWA//L57TQAO0iu0r8k+TqAp6HVmEXk\nsyKyB8BbAN4iuSiNNMwEsAXAWAC3iMi5AOYD+BSAMQCeAHCvc53RAP4M4DtQBdorAE7S0nGSU4i+\n5/zXP78nIic6u84B8GCK6TsVqoVzhrP9MIBjnPS+AOAPznN5FsBeqOfm8lkA9zifrwXwSQCnAKgA\nsBvAL53fLgFwFIBKAKMAfAXAfue33wI4BKAGwIegnvtl2jVmANgA4Ggoo36X9tsfATwH9az+x7lO\nOnwSwH0ARkK9g4fSaBERUMYPwF8BrAYwDurZf11EPub8/tke3tn4NNL7PajnMR1Ag/P5e9rv5VDP\ntwLAfwJY7FR2AOBWALXOsbVQ7+K/07i2oSdImr8c+gPwGlSBNQ2qQDoawJcB/Mv5/bsAWhMcey+A\nO32+PwYqs1+XYhouAbDV893DAC7VtgMA3gcwAcAXvGkCsAPAl9K8980ArtC2P+E8AxvAo853EwF0\nAZiY5DwjoFo+pc72zQDucj6XQhmF8c72egAf1Y4dB1W4BwBcCuBJAB/wnH8sgAMASrTvLtLe0SUA\nNmm/DXXSM9Z5XocADNV+/wOA36X4jBbozxqqxr0Tyn0W0Y/z+TcAvq/tOxvAdufzTJ93PN99Tmm8\nswUADgJ4T/vb7byjGmefLQDO0I75OIBXtTQdAjBE+/0+AN91Pu8FMEn77QT3WPOXmT/TEshRSK4D\n8DcAN3h+egeqoPJjnPO791yvAPgRgC+mkYQdnu2JABY6tcD3ALwLVaushKrBeff3bqfCu9DujeRf\nSY4EcD2AYs++r7sfRCQgIj9yXDr/hioIiaib5Y8APu24IM4DsIqq5eTe14Pafa0HcBhAGYDfA/gH\ngKWOG+NHTo17IoAiALvc2jGAOxDr1nlDuw+39XAknNaG9h0AbEvjGQHas6UqGV93zpsOVQAq3ft2\n7uEGKEOVLveRHKX9jUSsO7ACwHZte5snvbtJHvD+7rjejgCwSns/j0BVjAwZwhiB3OZGAJdDFbQu\n/wIwQUQ+rO8oIhMAzALQkuBcRQD2pXFtb8fedgBX6hmd5JEkn4bqa6jy7D9BS9vJItLh+LD1P/c7\n13X0Tyh3U7rpuxiq1XAayREAqqEKIQEAqj6TbQDOgnIF/dFzX2d67msYyV0kO0neTHIagBOda3wR\nqhA+AOBo7ZgRJH198x52ARgp2ugvxD+7ntCfrQAYDyDss9/7UIWoi1552AFVo9bvezjJTzjnvbiH\nd5aOO2gnlOF0meh85+L3PHZCVWj2AZimpXMEyeFpXNvQA8YI5DBODf4+KL+1+91mAHcC+IOIzHRq\nwdMAPABgOckVQGRYqduJWQ/V1P+zex5R4/jT8a3eCeA7Eu1wHS4iFzi//R1AvYh8SkQsEfk6VE3a\nTfOTJEtJHuX5c797ytn1/6AKhN+LM7xQREoBfNCTFm+ncymUS2K3iAwD8EPEG7E/Avg6lO//T577\n+oE4ndwiMkZEPul8bhSRkIgEoNwShwF0kXwDwHIAPxGRUlHUiMipPT1Ektuh+nxuEpEiETkZyrhE\nb06N9U/WajvefdZQraQDAJ7x2e9FAGeJyEgRKXfu3+VZAB0i8m0RGeK8t2lu5YLkH3t4Z6/7XC8R\n9wL4noiMdvqP/h9UKytyy9rzOAXA2QDud1o5SwD8VNNypYh8PI1rG3rAGIHcw1t4fR+qNhf5nuTV\nAH4F1bnZAeWv/xeAC7TjTgLwsojYAP4C4Lckf6L9PgHK351aosiHoFxKSx2Xy0sA5jq/vQvgP6A6\n8d6B6oN4KsGpkl3jXajWzAEATzppfwHKjfJVfVfPob+DqtGHAayFGmbqZSlUZ/I/Sb6nfb8Qajjt\ncqcTvRWq4xJQHZYPANgDYB2AFYh2KH8RykW1HsoP/idn/4S3p33+nHOf70IViL91fxA112MUVCd/\nIpYBuBDK9/45AOeR7PK5zu+h3tNWAI9CPQO1E9kN4BwoA/sa1OCBJVAd4ZlAT8f/QBm+lwCscT7f\nov2+y7mXnU6ar3QqOwAwD6pP4WlHd8sB5Puckawiytj2sJMa3vZTKKNxF8lbE+z3EahMdCHJv6Rz\nrCF7iEgllB/35IFOiyEWxzV2FcnPDXRasoGIzAbwe5LpusQMGaJHI+A0hTdBDSHbCTW07SKSG332\na4YaQvdrkn9J9ViDwVCYGCMw8KTiDpoBYDPJbSQPQzUpz/XZ7xqopvNbvTjWYDAYDANAKkagErHD\n/V5H7GgViEgFgE+RvB2xnXY9HmswGAoXko+ZVsDAkqmO4Z9CdeAYDAaDYRARTGGfMGLHMfuNSf4w\n1KgRgZowc6aIdKZ4LABAREzAKYPBYEgTkr2J0xVzgp6mhVtQQ7QmQg2JexFAXZL9fwM1ZC2tY1VS\nNGybbGggi4rUf9vmQLNgwYKBTkJKmHSmeL0UNWaeZ2Yx6cwcTrnZv2EjqMYffw1qfO46AEtJbhCR\nK0XkCr9Dejo2JetUWgo88QTw+OPqf2lpSocZDIm46aabYr8wGjMYUnIHgeSjAKZ4vrszwb5f8mzH\nHZsypaXArFm9OtRgSAmjMUOBY2YMp0FjY+NAJyElTDozi0lnZjHpzC1SmjGcDUSEuZIWQ34iIjAa\nM+QTjqb71DFsWgIGg8FQwBgjYCgYFixYMNBJMBhyDuMOMhgMhkGKcQcZDAaDoU8YI2AwGAwFjDEC\nBoPBUMAYI2AwGAwFjDEChoLhxhtvHOgkGAw5hxkdZCgYzGQxQ75hRgcZDAaDoU8YI2AwGAwFjDEC\nBoPBUMAYI2AwGAwFjDEChoLBxA4yGOIxo4MMBoNhkGJGBxkMBoOhTxgjYDAYDAVMSkZAROaKyEYR\n2SQi83x+/6SIrBGR1SLyvIicpv22Vfvt2Uwm3mAwGAx9o8c+AREJANgEYA6AnQCeA3ARyY3aPkeQ\n3Od8/gCAB0nWOtuvAjie5O4ermP6BAwGgyENstUnMAPAZpLbSB4GsBTAufoOrgFwOBLAO3o6U7yO\nwdCvmNhBBkM8qRTOlQB2aNuvO9/FICKfEpENAB4GcK32EwE0i8hzInJ5XxJrMPSFm266aaCTYDDk\nHMFMnYjkQwAeEpGTAfwewBTnp5NI7hKRMVDGYAPJJzN1XYPBYDD0nlSMQBhAlbY93vnOF5JPikhQ\nRI4m+S7JXc73b4vIg1DuJV8joDfXGxsb0djYmELyDAaDoTBYuXIlVq5cmdFzptIxbAFoh+oY3gXg\nWQCfJblB2+cYkq84n48D8CeSx4jIEQACJPeKyDAAywHcRHK5z3VMx7ChXzGhpA35RiY6hntsCZDs\nEpGvQRXgAQB3kdwgIleqn7kYwPki8kUAhwC8D+BC5/AyAA+KCJ1r/cHPABgMBoNhYDBhIwwFw403\n3mhGCBnyiky0BIwRMBgMhkGKiR1kMBgMhj5hjIDBYDAUMMYIGAwGQwFjjIDBYDAUMMYIGAoGMzLI\nYIjHjA4yFAxmspgh3zCjgwwGg8HQJ4wRMBgMhgLGGAGDwWAoYIwRMBgMhgLGGAFDwbBgwYKBToLB\nkHOY0UEGg8EwSDGjgwwGg8HQJ4wRMBgMhgLGGAGDwWAoYIwRMBgMhgLGGAFDwWBiBxkM8ZjRQYaC\nwcQOMuQbWRsdJCJzRWSjiGwSkXk+v39SRNaIyGoReV5ETkv1WIPBYDAMHD22BEQkAGATgDkAdgJ4\nDsBFJDdq+xxBcp/z+QMAHiRZm8qx2jlMS8DQr5iWgCHfyFZLYAaAzSS3kTwMYCmAc/UdXAPgcCSA\nd1I91mAwGAwDRypGoBLADm37dee7GETkUyKyAcDDAK5N51iDwWAwDAzBTJ2I5EMAHhKRUwD8HsCU\ndM+hj95obGxEY2NjppJnMJjYQYZBz8qVK7Fy5cqMnjOVPoFZAG4kOdfZng+AJG9NcswrUK6gY1M9\n1vQJGAwGQ3pkq0/gOQC1IjJRRIoBXASgyZOQY7TPxwEAyXdTOdZgMBgMA0eP7iCSXSLyNQDLoYzG\nXSQ3iMiV6mcuBnC+iHwRwCEA70MV9gmP7ad7MRgMBkOamMliBoPBMEgxoaQNBoPB0CeMETAUDCZ2\nkMEQj3EHGQoGM2PYkG8Yd5DBYDAY+oQxAgaDwVDAGCNgMBgMBYwxAgaDwVDAGCNgKBhM7CCDIR4z\nOqi/6egA1q4FQiGgtHSgU2PIR4zGChYzOijX6egATjkFOPVU9b+jY6BTZMg3jMYMfcQYgf5k7Vpg\n3TqgsxNYv159NhgyidGYoY8YI9CfhELAtGlAURFQX68+GwyZxGjM0EdMn0B/09GhamfTphl/raF/\nMBorWEyfwGCgtBSYNctkzhwgb2MHGY0Z+oBpCRgKBhM7yJBvmJaAwWAwGPqEMQIGg8FQwBgjYDAY\nDAWMMQIGg8FQwKRkBERkrohsFJFNIjLP5/eLRWSN8/ekiEzXftvqfL9aRJ7NZOINhnQwsYMMhnh6\nHB0kIgEAmwDMAbATwHMALiK5UdtnFoANJPeIyFwAN5Kc5fz2KoDjSe7u4TpmdJDBYDCkQbZGB80A\nsJnkNpKHASwFcK6+A8mnSe5xNp8GUKmnM8XrGAwGgyHLpFI4VwLYoW2/jthC3stlAB7RtgmgWUSe\nE5HL00+iwWAwGPqLYCZPJiIfBXApgJO1r08iuUtExkAZgw0kn/Q7Xp/R2djYiMbGxkwmz2AwGAY1\nK1euxMqVKzN6zlT6BGZB+fjnOtvzAZDkrZ79pgP4M4C5JF9JcK4FADpI/p/Pb6ZPwGAwGNIgW30C\nzwGoFZGJIlIM4CIATZ6EVEEZgC/oBkBEjhCRI53PwwB8HMDaviTYYOgteRs7yGDoAynFDnJG/CyE\nMhp3kfyRiFwJ1SJYLCJLAJwHYBtUR/BhkjNEZBKAB6H6BYIA/kDyRwmuYVoChn7FxA4y5BuZaAmY\nAHKGgsEYAUO+YQLIGQwGg6FPGCNgMBgMBYwxAgaDwVDAGCNgKBhM7CCDIR7TMWzIfTo6gLVr1aLq\nZglFQ38wSDVmOoYN+U9HB3DKKcCpp6r/HR0DnSJDvlHgGjNGwJDbrF0LrFsHdHaq/8+aaOSGDFPg\nGjNGwJDbhELA1Knqc2cncN11BVdTM/QzBa4xYwQMuU1pKfCTnwCWpbbb21VtzWDIFAWuMWMEDLnP\nzJmqtlZUBNTXA9Om9eo0JnaQISEZ0thgxIwOMgwOOjpU7WzatF6P3jBhIwxJyYDGso2JHWQwpIEx\nAoZ8wwwRNRgMBkOfMEbAYDAYCpjCMAIdHUBbW0EN+zJkGaMxwyAl/41Agc8GNETpt9hBRmOGQUz+\ndwy3tanM2dmphn89/jgwa1bmr2MoXIzGDAOE6RhOhVBIDfkqwPG/hixhNGYYxOR/SwAYlON/DYMM\nozHDAJC1loCIzBWRjSKySUTm+fx+sYiscf6eFJHpqR6bFUpLVfPcZM7CYCA6aY3GCos8GgjQoxEQ\nkQCAnwM4A8A0AJ8Vkame3V4FcCrJBgD/A2BxGscaDJnDdNIa+ps801gqLYEZADaT3EbyMIClAM7V\ndyD5NMk9zubTACpTPdZgyCh6WOD162MCgZnYQYaMkERjg5FUjEAlgB3a9uuIFvJ+XAbgkV4eazD0\njSSdtDfddNMAJsyQN+TZQIBgJk8mIh8FcCmAk3tzvF5Ta2xsRGNjY0bSZSggSkuBJ54wnbSG/mMA\nNbZy5UqsXLkyo+fscXSQiMwCcCPJuc72fAAkeatnv+kA/gxgLslX0jnW+W3AA8gN0mVGDSmSCwHk\njMYMmSRbo4OeA1ArIhNFpBjARQCaPAmpgjIAX3ANQKrH5gp51tdjyEGMxgy5SI9GgGQXgK8BWA5g\nHYClJDeIyJUicoWz2/8DMArAL0VktYg8m+zYfriPPpNnfT2GHMRozJCLFMZksRRwa2nr16u+niee\nMM31fKOxsTHj/tR0MBozZBqzqEyGMZM+Df2N0ZghkxgjYDAYDAWMCSBnKDzyaLq+IUcpMI0ZI5AG\nBaaN3KMAhtcYjQ0wBaAxL3lpBPojI+naOO44YOfOzJ3bkCI5NLzGaCxPySGNZYu8MwL9Zch1bWzZ\nAsyeXRCVhNzCO12/qgpoaVF/KbyMTM1ANxrLY/qosUEJyZz4U0npO62tZDBIAmRREdnWlpHT0rbJ\n2lp1XkBdI1PnNqSBbasHHw6ToVD0hYRC6rckGI0ZUqIPGss2jqb7VPbmXUugv2I7lZYCjz0G1NYC\nwaA67yCPGzU4ceP2b9sGbNwY/X7Dhqw13Y3G8pwc0Fg2ycshov05FtuM884ROjqAE09UPhRAlcyt\nrUlfSiZjBxmNFQC90Fi2MfME0ojGZQJ35SEdHcCzz6rPM2b0+GJ7ZQSMxgqbNDWWbQrbCLi9c26V\nKckc/DR2NeQxaRsBozFDjlPYk8XSGMpVgKO+DD7Mnj07vQOMxgwFwOA1Amn0zvWmI89M2sk/0g4e\nZzRmKAAGrzsISKsHLZ3OtnSb9sYXnMcYjRlymMLuE+hH2trURKDOTlWze/xxNWLMD+MLNvQGozFD\nJijsPoF+JJ2mvfEFG3qD0ZghVzBGwAd3HenHH09Q69KcudHMTNSXvYNpo3ZlNC3Gb5yfJNWY56X3\np8aMvgwDHi7C/UOGpvT3O7ZNNjSoOf0NDaRt027fybbiU2njSHLIENrtO9na2vcZ5j6XMvSB2bNn\nD3QSeibBS+8PjRl9DX5gwkYMAD5t89KVf8WsQ4+jFHvRcSCIU+YUZyS4mHEDZJbHHntsoJPQMwle\nen9ozOjLAKToDhKRuSKyUUQ2icg8n9+niEiriBwQkW94ftsqImv0BegHNX7O3HPOAYYMAQCsLT4e\n694YFZ+xetHu7q8YNTmL8U0kfun9oLGC0xdgNOZHT00FKEOxBcBEAEUAXgQw1bPPaADHA7gZwDc8\nv70KYGQK18mICyUruFEG9cSGw7R/djebl77NUEhFl4w0sfvQ7va7VF6SBd/EoNFYopfeDxorGH2R\neen/QgbcQakYgVkAHtG25wOYl2DfBT5G4DUAR6dwndx/N7bNRKWIrq9QiGxp0Xbrr9jD+UQ/PyPb\nVhnGaKyAycNnlAkjkIo7qBLADm37dee7lBsbAJpF5DkRuTzZjm7zdunSHGyteVYS6djZEdOq1P2r\n7e3AsGHaiI+CbHenST8/IzcQZGenek/P5qJj0misfzHPyJdgFq5xEsldIjIGyhhsIPmk346BwI3o\n7gauvBK49dZGrF7dmNakmIzOqvSeTMuBHeu245TZAazbCkydCvzkJ1FNrV/voy93PKCJD5wY5xl1\nPLsBazkNIQxDJp9SKAQcccQ3sW+fKkSvu653UYGNxgYxpaXoePgJrP37NoTOnojSQfiMVq5cmX74\nk57oqakA5Q56VNtOyx2U6u8AuHgxaVmptda8reaMuvv8TuZ+V1TE1soLGAx2RxYcsiz1Uzgc71/1\na90nafEXNP3tsm1uTl1fbnqMxvKHPOwSyFqfgIVox3AxVMdwXYJ9FwD4prZ9BIAjnc/DADwF4OMJ\njtXzQNKX5PcyM+ruS3SycJh2TQOb5XSGSjbFZFK/aybL5/kkxD7jvMDW5r396rJNVV/6vtnWmN2+\nk62VFzAslWwYstFoLFPYNlvvfCnyPPOkSyA7RkBdB3MBtAPYDGC+892VAK5wPpdB9Rv8G8B7ALYD\nOBLAJMdorAbwsntsgmuQTG20gl/+SSWDp1w7SnAyu/lpNuBFBnGQIbzEptvWx4/S6CGdedg31Te0\nEssOncCGUGdKhXRfLpfKaJjeaCyt2rfPyWybbKjdyyAOsgGrGbYmsGXRWqOxvuI8a9sawYYhG1lU\n1J03xjFrRiAbf64RSIVEmTFZBveOrGhu7kEEtk275Rm2Nu+N7NfavJdBHFKZCwfZ1rI3pWvq6Uyn\nNloQeEosu+UZtrTEv59suzfS1Vja+iJph222Ln6ZdljtqB6FU1PFAbbVfj5iHIzG+oCmMTs4km1L\nXmY4nB8utII1AmT645v1sgboObMmbGaHOlkU7GJDqDOla9s24wq1ghqb3ROeEssO2ym7N/o706bz\nnvz05frw/dKYvEugmw21eyPGIZV0Go0lIYc11lcK2giki/uS3Y5Bb4dbKs1s9zxtLXtpNz+dkjKM\nfzYFtBIrVfdGb55rf8YO8tNXURFZW+ufxkzpS7+20VgSsqSxbFPQRqA3FtqtMYVC8bW2RB1ucU3q\nFJXhpq+5ucD8s32sOqXq3uiN37u/NabryzUAgYC/xjKlr94+i0FLBqrm/amxbFOwRqA3FlrXjm2T\nTU1kSYl6AkOGqGa73zFxTerWVoYD43knLmPYmkC2tcXq0rZVB3KoM+JyStaxl1dkqOrk99y93/XG\n792fGvPqq62NbG9X2kqksVT1FXf+sK06kIPdEZdTQfQBZLBq3l8ayzZ5aQS8ht62VW1a93ema6H9\ntKOfIxhUNbhUaF9lsxgHCXRziOxn+yrNvxjqpB06ga3WydEO5GAXW5r2KnGFc9zB2FeyXHXSM20q\nFUR3GLJXX62t8b77dG4lUdnUm8fh1Ve43Y49f6iTzTWXM4iDznm72bb4JdphO2030qBjAKrm6Wos\n2+SdEQiHoz7UUIhctow89lhG3DZ1ddGX4bp0amr8a/E6ra1RX20gEI25EgrFnzsZtk1WVpKAO3a7\nm/Pna7oMdrHNOok2jmQDXmQRDrIBL9IOnRCtruWyg7GvDFDVKdUKIoAYH317e1Rvrj5090CqGtPL\nJssiFy1Sx4fD0dZmSUnPOvXT15IlnrIv2MWWwOlswGoW4QAbitfTtkZEm5xGX/166Vx7vHlnBGpr\no4WyW2Dr20C0AJ88OfpdKBTfrNMttp4Z9QJ/2bL4cyejuTk2TcXF5KpV0YJkfEUn2yefo4Y61jSw\nLXCiWgSkqIhcvFiNhceRbLVOpt3yTIqveZDh6+PoX1Lpd7FtlWH0Ttvx4+P1JaKOD4fJiRNT05i3\nQuHun66+Wlvj9RUOq7/qamVg6qeq1qYdHMm28RfQDgyPWh+jr3675J135mb/QN4ZAe/IHb+/piaV\n6b3GQpvUGzciw1t4W5ba33ueZJnUm9ErKpQBaGhQBYf7/ZCSboabnovU/O3gSDbXXM7me99muG6O\nM9nsEBtCnQmHD+Y1GW5Te9+Lt7B2aW0lgf8X2W/8eP9KhmVFWwje75NpbNmyeI0tWpSeEQiHo30I\nwaDSl7fCU1xMhtvtqKVqaGDYmsA7KxawvfZMoy8yoxpzWwCWpd5NrvUP5J0RcDNWTY2qresjePSM\n5C28a2qizW/VnGbEYrujNby1NLcGV1enXnBP7iA97ozbh6C7mfS/JUvUMXbYZl313kjzvmZiJ4NW\nVyRtccMHc9HpmEn6oU3tdcUkKmhdY2FZSi9+Bb17DqfR5quZdDQWDsdeM5k7yK1temMbebWu64tU\nBqGkqJNAN4uC3YWtLzLjGvP2HS5ZkluPL++MgGUp4YbD5LvvHuLtt7/AsWN/ReBaAp/mkCEzeeyx\nU1hdPYnB4HgC1Swuns4ZM07m3Lnn8KijriDwfQJ3EXiU48evY1NTR0wh8bOfxXYKupk0UQ1S38/P\niHgLErcJT8ZnYBG1v5tBYzJ8y97cdDpmkhQ79tIpq1J1E/u9602bOlhZ+TiBn1LkcgJnc+jQD3HS\npGNZXFxNoJLBYC1raj7Ek0+ezU9+8jwOH34NgR8RuIfACk6c+Ar//vcDCTXmthoSzUfR7yEYVLVN\nbytW11BRUawxWbgw9veKigLWF5lxjeX6CKFMGIFshJJOma6u9/HKK3/B2Wf/FVu3NqOiogKnnXY8\nhg9vQFnZKTjllEqMHz8SRUVFOHiwCOvXH8a4cR3o6rLR1rYby5fvhFru4HEAYbz11g5ceOEOkCUA\nqlBSMgHPP1+F99+vQlVVFWy7Chs2VKGraxza24NYtw6YNSs+XWvXAhs3qs/BIPDTn0Yj9T72GDB7\nNrB1K1BeDvzzn0BFhf/9TZqk9t++HaiqAs46SwsLTJ8FX/0SM5hx47n7xkJWuCH13YjITzyRICpy\nRwc6nln5QyV0AAAgAElEQVSPtQjh4YeHYfv25BGU3XfY1bUT69bdj1NO+RteeeUZTJ06DTNnHofj\nj/8gioo+gZkzK1BWdiQOHizCq68GUVV1AN3dNjo6OtDa+g7+9rcwgDCANQBeR2fnDnz60zsBjAJQ\nhWHDqtDeXoVf/7oKEyZMwL//XYXXXqtCV9cYrF8vvq9VXydABLj9duDCC9W9zJypHtvGjcDYscCK\nFbH6qq6OPdf//R8wcWKB6gvIuMbwzHr8+PshyLBhmDEjPyN055QRACaAPBEvvngepkxZiKefHhf3\n0PUQ7PX10e8OHVLbGzcqnQOqbvTPfxKlpe/hrLO2Ixzejubm7RgxYgdWr16N117bjkBgO7q63gQ5\nCpdfXo5x48pQXl6OsrLo/9LSMtTUlOPVV8tQXz8aM2ZE1+KpqABeeME/jLuegauqlNgqKqKZOCb8\nO5IFi88TUoh577f4eVxZ1dGBjhPPwClrb8c6FKO6pguPPWHFnE7XyZFHEn/+880oKnoanZ1Pg/w0\n1q27BuvXP4Rjjz0yLg3usXPmRJPY0QF0dalC99VX1XeWBTzwADB1ahdmzXoDW7bswPDh21Fevh2v\nvPIKVqxYgW3bdkDFU7QBjMVVV5WjoqIsRl/Dh5dj4sQybNtWjilTyvCZzwxHaalEHllra+JH9tGP\nKrm0twNTpqiC392n4PQFDJjGgsH9eOSRR3Deeef17/31A6JaFAOPiFBftCwYVO9SfzleC/7ww+ol\nXn+9KminTgV+8APgO9+JZopbblG19G98Q2XiYBB49FGVwd1ztrV1Yteud7B//xsoK3sTtv0G3nzz\nTWzf/gY2b34TBw++gbfeehM7d76Bzs4DOP744zFz5kycccYZOPXUU1FUVJRwsZGOjlg9uvtNnAhs\n2+bZ37tzAeK+Y7escmtpMc93bRvaTvk2Tu36JzpRDICorRW88EJ031NOAdauJaqq/o5hw/4bGzfu\nRGfn1QCuhwpwCyxZAlx2mf/1E2ns2GOBgwcRaXk8/LAyBK6+ioqAX/wCuOiiWAPy5JMHsHPnW9i/\n/w2MHfsm9ux5A9u3v4kNG95AV9ebePPNN/D662/i3/9+AyUlJZg5cyZOOOEEnHPOOZg+fTpEJC2N\nPfOM+m3mTKMvL5nV2EGUly9Bd/cPcfzxx+H+++/H0KFDs3YvIgKS0qeT9NWflKk/ADG+zfr6+Ak9\n3qGAbseX7i91R4+1tKjOXvc3fYio7v9P5u/3iwq5e/duLl++nDfeeCM/8pGPcNSoUbz88qs5ZcqG\n6ISxBBN2vL7fZH7iQibRzM2IOzts0w6dwFpsinS666N3VIf9CwROJBDiV77yZ1pWF4EFMX037e3R\n8ycKweCnsZaW6OAcb4dwSUn8e/XTmHfaiD6SZ8eOHXzggQd43XXXcdKkSaypqeGCBT/gtGlvp6Sx\nVEZLFTp91dhTT3UzEHiAQBVFzuZvfvP8gNwH8q1jePJksrxcdaAmKozdEAzeIX7ezl3vyB3Liu7v\nZmSvYdFfdKKokN6x4n/5yzZ+7nPfIzCWwKcZxMtss06K2dnPiOmFyuLFPcxIbG8nv/3taKk12Elz\nlIpvX59ts33pKhYXqwzqhmXYs2cPL730a7SsMgYCSzh9elekwIUzAm30aKUxPeSCXhjrk8S8Gkqm\nsUAgdY3pI5CCQf9Ac7atCpuVK5/j2WdfSmAEgW8yiDdS1pg7YspvVnQMRmMpa2zLli2cPfvjLC6e\nRstaOaAVubwzAjU1sYWutzC2LHLp0ujIB3ckRUmJytT6MDxvjai+PmpAQiG1HQyq/3qLwWt8vIWA\nOwRRH21SV0dOnbqPAfkhLYzgNSjiW9YIti5+OaaQ0Y2Ym/YhJd20Ap0MVdsM1XXGG5v29uhEBJHB\nn0l7MYQv0QgNb8a9444nWF09iaNGfZmW9W5kpJl7Dj+NeQvjRYtihw277yvTGtPj/eitDVdjuv5q\na9WcgerqnRS5jEGM5C9Qwj0JNFZfH71WXR2VrqxODinupGX5LKhiNJaSxoLBbt5wwxIeffRojhv3\nv7SsQzEaGwjyzgjoNZjaWv/ZmHrtLBgkr78+tkVQUaEm7rjHtrREM5XbBPTO5Gxqit3Pxc+t5BoJ\n79A9NaOzm1OL23iODGWJHEPLej5mqF4wqK7luhIWL9xHy4kxZOEwAzgcP7Lt29+OvdD8+X1XzkDi\nWeCjdfHLKQ3Ta272fz8qrx/m2LHfZVlZOW+9dVmMYXBrwa4R8GpML0BdF53eSlu4UP31h8bcfb1u\nJdf9qKeluFilobamk8snzeAHYPGoQCMt6w1fjbnXal4WXQjJdWsUFXXHjpw0GktBY2/zqKPO5fTp\nH+Qf/rAuocayTd4ZAT0Er25dvX0B1dVRvU6ZojKlrmE/X6jbOgyH48dWpxIuwttU9BqB6Jjsbt45\nbzMDgXsIjKFl3RNT+9TTFb73MdZgE4EuAl0swT4Gra6CaAnYwZGRNXR7Gt/v1zfj8tprb/P440/j\nnDlncMuWN9jcHJ1o6NbM3WMTacy24yeIBQLRgre+Pvsaa2nxn8zmjvl/7JcvUOTbBKoYDK7y1Zht\nk823Ps86vOxorJtAF0M1e2Oft9FYRGN+Nfsnn1zN8vKJ/MpXvsXHHjsY4zLUNTYQbqG8NAJeS0xG\nJ9y4zWzdCLgZ1utr1ztxvB2ygUC0A6+uLn5FJjdqqetD9QvV67ZQdCG4v0drly/Tsmp49tk3MBCI\nXeA6HCZrK953WgLOLE8c5JJF+/z7BObPH/yZ08VWyyqmsuh3bN9MNy0ruqrbmjVrOGnSJM6bN4+7\nd3dG3nFJiXrHNTVeXfhrzH3n7jusq4sv9DOlMd1Q6NFx/dwRuu69k8ii7qL7GQiM5tVX/zlmclhk\n3YxAF6vxSqTFGcQhtjTtjX/QRmMRjdWO3xdZ1W3p0qUcPXo07777vpjKiOsy1DU2EDGFsmYEoBaa\n3whgE4B5Pr9PAdAK4ACAb6RzrLZfTCepi+4brakhb73VP5yEHr/HzaDt7ep4fTq+++f6f3Ur7m2W\n6yM9wuFo4e2NG+/9TKprFxWRwDsUmcERI65gMNgZOZeq5bnRIrsYxCHWVh8aUP9iNknkg020n8rM\nUWP5kx8+xNGjR/OPf/wjyfiOfLfgdgOvqe8WREJ+xJ9fZeyFC5UrJ1WNrVqVuNPfqzE90KffaDXd\nPeTVmKs//TmtWuVecxWBcRw37u7I84x1R3WxBptZhINGYwn20/NjEAfYesznePP3vseJEydy9erV\nceEj9KjE1dUDN6M4K0YAQADAFgATARQBeBHAVM8+owEcD+Bm3Qikcqy2b6T5nahjxn3obq2ookIF\n10rUQqipiWa6SEesVqPyNsH9Ysa4mV3vrEul6RfrDrBZU9PIM874HHfv7vQEtOvm2JEHWF7WGRl9\nksoi5fmA13Am26/ptnVOQXaAVfgRx4wYyaeeeipmHzfD6wVsSYkaTODtLE2kL/cduJ261dVq/0Qa\nKy5O3Onv1ZgT6DNOX4FA7AixVDRm294oqBs4cuQE/td/3U7b9hqBbn73km2sKO+MCc1SCKSqsXCY\nrKlQfXTTsJpflRJOr63lzp07I+dx9RUKxQb1KylR/TB52ycAYBaAR7Tt+Ylq9AAWeIxAOsfG1KD8\nmtl6Rh0/Ppph3QI6HI7NGHqmKypSwZ/0GpW3puBtCRQXRz+7hXOqTT9vx+Cf/rSPp512Gv/zPy/n\ntGndke9VayHe6Jj5A1HCYdUhGsAhjsGNLA+W8CXPw3ddOk1N5FVXxT7PJUsSh3X205c+PNQdyeNq\nxqsx/T3q8wf8NKYbCq+hSlQxSc2FobTa1vYqJ0yYwDvv/F1CHbt/7sALgyIcJmuqOynYxxE4hycc\nUcrd27fH7bN4sWrl+2lsIMiWETgfwGJt+/MAFiXY12sE0jk2oUDdERh6p16ijOJdmMZtCfjVfvxG\nBOijPXS3gN6092v66b5e9399vRMDvt4tEGxOmzaTIt+i8m/7hzN2r+d1jQ1W3GfTm3uJNtVJ4GcE\nJvGB370cc07dZejWwnUXTnt78rDh6eiLjA0K54YXTtSC82pMd++0tChN6S2FpqbELgxvf4I+Z8Z1\n5a9fv55HH13OQOAvkXP6acx1i+VLUNHMaKyTwEUUnMEVf3srocZ0A+4a8YFqWeWdEQgEFnDkyAX8\nxjcWcMWKFZEb1cWvhVFP6IfTm4B+sd/dffya3d6OYW90R7fl4fXb+s0E9vPjbt36LocM+QADgdsi\nHUyumFzXg9dtNZgzaTL3hl/G9X4XrfEuIVDFiRNfi5vgtWyZf0hvvVWpd+T7zaLVC+tU1uz1FuZ+\nI0QS3Xs4rPqpwuFouvQWZ6L+J6/GEg2kePzxVQwGx9CyHo/TWFFR9Nh8WYgsWV72Mwx+GlMzyi8h\ncDpravbHzb1YuNBfY4GAMtzZYsWKFVywYEHkL5vuoEe17XTdQake6+u7S/aCU/H1JVpL2G+GoDdD\n1tVFm4De4Xt6mhLNBE5U09qwYQdHj67kggV/iRQibuHT3BxbOxzsLYJEkX39CmW/d23bZFXV7wlU\nsqpqU6Rg1M+pj+QpKfEfTeNeM1WNpaov7z32pDF94Rg3dIWfC8ibpnQ19tBDyzlyZBnvvntzjMba\n25We9FZusnUYBgPJ8rLX6Pu96z17ujlq1BUUOYU1NXvjNOY+I9f4JgpBMxBkywhYWudusdO5W5dg\n3wUAvtnLY31vsq9rS/vVtLwFjusq8lskRu8QdmuGfmO63d/1FkOimpZamOR5AqN5zDHPx6VHP3aw\ntwi8z869Bz/3jP783Rr8vfc+xFGjynnXXeviMnJRUay/250olWg0zezZs33TmE2N3Xln7H1XVPi3\nPLxp6o3GRo++ncAU1tW956sxv0mQgxE/jSVy//lp7Nprv8WpU2eyqcmO0Zg3LpTer+g3uXQgyIoR\nUNfBXADtADYDmO98dyWAK5zPZQB2APg3gPegYucemejYBNfwvclEhUg66LMv9UzudRW5a7nqL95t\n6uk1Q7806e4BdzJZohpwdO3avxCo5L337vDN9N4WSK6sa5oufrVqv0yq15KHDCEfeOA5WtZoWtZz\nMc/ZdQ0uXhzr73Z944l8w7mgsfb2+I5atzDxC2jWd41dT+CjfOSRg3Ea01ucg1lfZPzzS2QEvBr7\n7//+BUtKptCy3o3TmHfJ0Jqa2PlDA20AyCwagWz8eTOot1OmL2tLJ6uNejNuqouD95SmVGvAwG08\n5pgPcteuvb6Zvq+FU67i1rT0oGyxbpWtHDGigoHAQzG1fL2Jr/tti4tVJ3CyllO2NWbb/ktG3ntv\nNjXWSeATPOusy7lnT3fcaLhC0hcZm+ct6+886qhyWtaWmMqX32xgd8lbffBBLjyzvDUCifoB+oJf\nTStREz5ZB2JfrknGG4GJE7t58cWX8IILLuC//90V3d8poeyw3afCKRdI1kGn31u0lvZvikzjDTf8\nJOYdlZfHPjt3bPbixf6TAb2FazY1pp/fb7ZvNjVWVWWzrm4af/azn0X3DxeevmLz/Gpa1mjee+9T\nkZZBSQn5/e/HPrv77otvdenuoYFuPeWtEeirjzYVEnXmkanXCr01yVTXLNVrF+EwuX//fs6aNYs3\n3XRTZCc7dAJbrZNph04Y1Dk0ncK2tZUMBA4R+BhFrmZra3dcIDX9r6YmvsBN1mmXTY15z+9doDzb\nGnvllVdYVlbGFq2Ty7ZGsLX2C5EQCYORdPWl3snrBMbzu99d6hsOvCeNuTPH3dDSA0neGoFsuEF6\nukZPmU/PaMmCSCU61lsA7Nq1ixMmTOADDzxAu/lpNuBFBnGQDXiR4abncsYHmS7pFLZ79nRz5MjL\nKHIWp007HHn2foHUVHM+dtKX687z882TsUagvzWWyvl70liiOS+91diKFSs4duxYbrrvPtrWCDZg\ntdJY7d6c8nOnQzr6Unm2g8AHOWrUDyJDdL2TBZNpTO+H0ie1DhR5awTIvvtoUyHRNXRhJMp8iaOI\n+s92TqWgWbVqlQpWdXtrJPxvEAdZW+OsM5BkRalcJZ3C9vvfv5VDhjQwELBjo61qhWFdXWznvV8k\nz0TX844O6m+NJTt/TxrzGr9gMDMau/322zl18mQun3QhgzjonLs72nk9yDSWjr46Ozs5Z87ZLC39\nEgOB7phnrU8WTKaxbFRQ0yGvjcBAkihYlJ75vEagoiI+I6btcrBt3nfzzawaX8X6qTtZFOyKGACA\ntHCILYHTc0N9aZBKYfunP/2JY8aMZyCww7eg02fd9hRaIRsViL7Sk8a8boqKimgMo75q7KrzzuPH\nTj2N04/Zw6Ki7pgZ0oNRY6m+72uuuYYf+tAcx+XYe43lkr6MEegn/Jrh3sznuoPckMH6WHC3ea6P\nvugxMJzmp7105Af4keNm8bHHDmjxjJxY8FhDOzhy4NuhaZDI7eF+bmpq4/Dhozlp0guRQs8N8uZX\n08212li6uIWMvgqZV2PuPVpWbJA6PfxEbzX2bmA4jxtaxq9++eueRW3yS2N6CJfrrlvIyZPrOXny\nbt/BIINVY8YI9AO6IPRCPZHbqKfZxMnCCsTQ2hrx01rYz+HDPsEvfOFL7O5WnaNBy4mLjgNsq/18\n7qrSQyK3h/vZsl4lMI4if4tpWVVXR8N3JBoLnyu1sXTwPg9vTCGva8tvroi30OqdxnayODiJd9zx\nW5LMS425M3yLipYRGMcJE16L0djSpfFhwAebxowR6AfS7WjSa3R+s4ndpn2P57RtttZ+QfPT2qwd\nP5k/+eEPtVpJNxtq90ZGc/h1COYa3jWidRdIIPAegalUgeEY9+cuDOPnBhms9FVf7jvPiMasNRw5\nYjTb2tryUmPq7zkCowk84zvyx12nZLBqzBiBfiDVZmCyGp13kop+zkSx3G1brQcbqtmrMuKQjXw5\nUMryYJD/ePDBuFpJuh2C/U6C0kK/d93tUV9/gEccMZuBwDdihnUWF8cPCa2pyY0p+pmgr/pyf+uT\nxqptFjkjz+6rmsqKceO4Y8eOvNNYcfFrBCpYVPSXuJg/1dXxs4EHo8aMEegn3NECyfyr6QRGI2PD\nD/sN8XMzW00N2XTbetrWCBLg45bFsSNHst2z7F/KNcoeqnIZqen5+Se0k+qFix22+dQdL/KCT13A\nc889n0891cV7741do3nR9a+wekJnJINmahx/othB2cZ9HsmGZSZ7v73VWGRIc/X7bAmcThtHkkVF\n/NFVV/H444/n+++/n3IaYk7cg4CyqrGwTbv5af7jwR089tg6XnfdTxkOqxAZkZhBVhfvu3svKyuZ\ncY1lG2ME+olUakCJanSJMk6yDOWNWFhT3cnmmstpB0fSDp3Aedf/jDU1k/mPf+yOmenZY42yhxtJ\nVmikhffm/GJ3O+mxQyfwizKeM4YeyX1vvknb1lf96maopJ3hwHjWFW+mu9yfvhJYXxhMGkvWYuiN\nxmJHs3WzqfprbLVOZrhuDp9a3sHPfOZzPP/8i/jUU92payyFjJJVjTn6Whk4kScfcRSvu+qqSBpq\nahjp+K7DWoZK2mlZ3TGhtQdbK4A0RqDfSLWWnagzzy/jJMvUth0/ISoQ6GZd9V7WVHcyECAt6xqK\nnMHpJWtVK6Ghoecp/0luxHvNPk18cW7ODo5ka+UFtAPD/a/Z/DQr8V0Cx3AqHmPzonVxyyA2ySfY\nilm0HL+121R3XR59qVEONo0lG5CQrsbihjSXdzIQ6GJJsVrcaMqUfRw69CMMBG5mw5CNqWmsh5vI\ntsbs5qc5HS9QcCGPwhw++tOXtPW8nTTgEBfi6ki/iNv35MaiyvX+Dy/GCPQTqfptkx2fKPMmylDR\n2kqiv8METmcA17INM1Nrvya5EW/rI9Fyg34Fr+93YZsNtXvVxKPirQxbE+Ku+d/zfkWgksAWAl20\nrG7PPXezpeYy2sGRrC3eGmkJ+I266k1GLWSN6f0IRUX+GhMJE6ikhftT01gPN5GKxhIZ9kQaC9Xs\nZdDqZqhkkxrGql33qeUdFFxO4FQCe2lZag6Evp63q82GIRtpWbHLvCZqwOYyxgj0I8kK7P4iWZwc\nVSC+zWJU8ScoibRf4zKL94sEN6LnX315Qr99/MZPJ5sYB3SztmIv7WX/iuzw4IMPsqysnFNqX2LQ\n6oop4PX1fO2wSm+43WZtbfJRV+lS6Bqzbf9ge5FastVNC620MIIrMZQMhWiH7fhCWtdYkpvoSWOJ\ntJToe69Lq+W2VZGOu+7ubl5zzTc5dOiHaQV2x+jLLdxra8lwu0qvHbbjlhRNFG4klzFGYJCSrPbj\nZpqaGuUrjx3u1sU7MJuVAJd84zuR4YPBoBrWF161SwUEc5ryPZUuPflr0xmKGNf0x0G2WSeRDQ18\n5M9/5pgxY/jYY89z2bL4KfmrVpHz5vkXEsnWcEgXo7F4jVVVxRas47Gdv8YIlgH85//eHdWHE07C\nbt+ZMY2l27cRZwRqLiODQXZPn87/njePoVCIa9e+y9tui67ZMGSI0lVTk1oi0m+dcbeTfjBMDvNi\njMAgJJUOQb3gW7gw1gjchm/xXzUzGQyWOYuJdzuF7gHWBl9TsYbQrtwxPVRnevJL+xW8tk021B9i\nkdXJhvpDMbW3ZffuZXX5PlqBTobwEsMo5w8CUzlmxEg2N7dG7lv3xf7mN7GLfCSLypioDyZVP26u\njA7qb9LRWDisaumRd4JDXCTX8BcV02hZFQTaI0a9JXA6G4rXZ0xjSfs2fDQWDqtBEwHpYk35XoYD\n47kHw3ixVLK++hhu2fJGnMYsSxmAVDSWzJWWq30FxggMQtKN9WLb+jKA3QS6WVPdSctaRWAMgWW0\nnEwZ7exSvs+eQgSnUrt2V/CKZJxwmHbJaLZhJu2S0WrbVjVFC4dYgn20cJh1xZtYiVsIlPPYY55z\n3Dnd2n2oe9LXBwZUyOVU6amwK1TS1Vhsf1RUY4HArwhUEWhnLTaxGadlXGNx+nK+TKaxYuxnEIc4\nrWgDR+E/CHyY06a+5auxkhLy1lvzV2PGCAxCeuPW8K5EZVluE/sZBq0y3jB6NsMoZy3aoy2DYHdK\ncweS1X6i7iYtrd5FcpcscYKddcUUIiL/S9UJvE517DbtZcOQjSzCARbjYCSdgUBs0z2d+OzZWHdi\nMJJZjS2mZVVw+cQP08aRGdVYQldRShrbR+BCqk7gPTEaUxF4oxWNpUtTb216yXWNZXuN4Y0ANgGY\nl2CfRc46wi8C+JD2/VYAawCsBvBskmv046PKLdJpdobDsTMdXT96ezu5eOE+/utn93L80Ufz5yJs\nxzEcj+0MWl1xwexiOo57qNrou+ijJ9raqHzCxbPVRCMnR4XD5JCSbifj7SfwZYp8gMDWSHrt5qcZ\nDoznPPyAAW14XnGx6hNwF/BO9zkORj9uNsikxm7/r1s4dvhwPhUIsB3HsAI7aAW6Eg+rTFFjiYaP\n9qyxMIEZDAQudIxBrMauwsIYI9DUpNKZjxrLihEAEACwBcBEAEVOIT/Vs8+ZAP7ufJ4J4Gntt1cB\njEzhOv34qHKfRPnGWym6+mq1cEqovpNBHGIIL/F3485nNYIcjf+kYD+Lgl2Rsc/usMBIKIEUqjbe\nkT6W1cWGUKe2pm83G8a/Q7t9p2f/MIHZBM6hZdmRzN3SokZl1BZvZQAHWYwD1Edv9KV2NRCjuAYr\nvdGYhUOsxSbeM24aRyHAsfgpATW6a/Lk2JnLkRn2vdBYbU2nZ/RZIo21EaiiyILIEE+vxsTTEki0\nhnM6zy1XNZYtIzALwCPa9nxvawDAHQAu1LY3AChzPr8G4OgUrtNvD2owkCjfRNfdVbU1d7HrqMi7\naOEQx+IlAqcROIPAu5EmfXQ4YDfHjznA8Kpdyas2tpp23xDqZFFRN0PFG9kip6vlLpv3JhwZNGrU\n3wiUEbiJQCfHj48OvWtvd2t9rguoO/J7Ltau8pXea6yb5Xid5fgbgYkE5hNQkxhjKgsBtf5FuL2H\n6rOrsfpDLMJB1mITw5MbSdtOmMYdO7o4bNgPCYwl8CABJtWYqrwM3pnAqZItI3A+gMXa9ucBLPLs\n81cAJ2rbLQCOY7Ql8AKA5wBcnuQ6/fiocp9kzc72dnL+fOW31TOdWyOLbh8i8HUCkxgIPM5QSJ8u\n79S4ireqoaSLno0Zx+8mwl3bODy5kW0V56kmuWNR7EW/YUOoMybM9nvvvcdLL72CwWAVgccjBcmq\nVVF/r75oiZt529t7jp+TaQpldFAieqexbs/ntwg0EvgoJ03a5rQ0XS1Ga/XhdrtnjVXNVJ2/rsYW\nLYqEqtA19vzzGzhs2GwCpxDYTte3n0xjFRXKDeS2LnJ1dE9fGSxGYJzzf4zjSjo5wXX68VENDhIN\ngdSjSUZHCtGTQfXvl/Goo8bx2mu/xS1b9nH8mKj7xcIh1pS/zyAOsQEvxixk713bODI133XeB4MM\n181hbU0nLauLEyb8jmVl5TzvvKtoWbvpdvQ2NcWH93VXXnOjhLoxwLI58sJorC8a0z8fJvADDh8+\nmrff/hs2L3yZNdgS+T1oOSvipaIx1wC4QmloiEwUtKwOjh37PR511NEMBBYS6KRlKWPlVh5cjQUC\nav6Jih46cBrLNtl0Bz2qbafiDtrouoM8+y0A8I0E1+GCBQsifytWrOiv5zao8C5DeNtt3glk8X+1\nteSffvcaz5s9h1UVFfz5BZdxUvA1reWgLR5inRRpc7c2742sbVyEg2yruVjlqsrKyEWXy+kMBB4g\nEKLIDP7qV09rNcxoLHrbdleriq2d6bMy/RZM6U+MEfAnVY3pM43rJnfyjm8/wIYp9Tyhupr3V85k\nMfYT6GZxsXINpaSx6gtjTxwM8m+3PEmR/49AGUUu4tKlr0dbMNoayLHDp6PDQQdSY/3NihUrYsrJ\nbBkBS+sYLnZq83Wefc7SOoZnuR3DAI4AcKTzeRiApwB8PMF1+vfpDVK8hWldXXSquz6iY/x4cuxY\nUqehUqwAABMdSURBVIQsKemmhUNswGo+iqH8MMBqHEmRJYSTUS2nlhaum8PW5r3RTrlQJ4uCnWyo\nsVWnXFsb2d7Od0pK+AMUsxjVBD5E4K+cNq07OsEnbLOl5jI2Bz4WOefSpfGFiB4KItuzNI3G/Emm\nMZGocZg0SW2Pr+xkXfFmBnGQ0/E8b0cJx6CYwMcI/IMiXays6GQQBxnCSyoirrZITZzGWlrI+nq2\nA7wGRbQwisCnCKyJWS+hrWWv0pYT/bS1eS9vuilWY7fdFqupwToTOFWyYgTUdTAXQDvUEND5zndX\nArhC2+fnjrFYo7mCJjlGYzWAl91jE1yjf5/WIEaPKeQGU9NjwFgWOW6cnhmis4gX48vcg2FchiEc\nVjyHwAgCn+fYMX/kip//jdPqDsfEpbfDNttqLuYrUsklFWfwhzf+kI0f+jCPgMU5KKOFFQS6GQx0\nsum7rWxt3stwmGxeuJYhvMQgDnII3qdldcWEhwBUcz0cjp0RrS/u3d8YjSWmJ43FrsrV7YzAiWrs\nVZRxwoiFBKYTGE+R61k+9q+cMmGHCipYGx2e6WpsvVRwYflp/PZ13+b0SbUcgSJ+BlW0sE6dO9DJ\nhVdtiLQs/TQWG/ZCzQnQXV7Z1li2yZoRyMafyaCJ8evQ07/Tg1/pncRBHHT8sqsZDlaxZuJhAjsJ\n/Jwin2JpaZljFKYRmM3p02dzVl0dKxCgoITAdI4cdgmrxvyRFt5jCC+xvngzgzjEeqxlCGto4RCH\nlHTFdQ76uanq6z2hAbLsqzUaS0xPGoud2d3NssCbDOAwS7DP8e+vZvsvm52FWtYSuJFqtFop1aix\n4zh0aCNPPnk2jzv2WI5GgIIjCJzE0aXXclL5o7SwT2kM6xiMzD4/yFBxuxoSnYLG9NFAuT7bNxMY\nI1BAJOrQc0fY6M35Y46hk2Fc3/9BLr55V4yhqK0lH3qom8A7BNYQ+Bd//ON/cfmPFnEi/kngAN2O\nZMup9QWkk2WlexjAIdZgU0wIAfe/ZXVxSEm3b6TKVBdA6S8KfXRQTyTTWHt71P1YUkIWF7v60jS2\ncF/Me6+pISdN6qSaP/IsA4F/8uc//xdX/PiXrMYTVB3MHo2hk2V4nYLDkQLfwuFIH0OuayzbGCNg\niGDbqsnb0uINOtfN8RWdEUPhjqJYtixqPPTx1HrHHdDNGmxmHV6mt1M5gAORjkBBJ4NWl1oa05md\n6Ybp1d0I+pT9XJ+JaYjHnXUbG4tHdQK7kwm9GnPH7+u18VQ1prbVX0mRGppcX+9MDDMaI2mMQEHh\nHevs+jr91kFub4926AEqUqde4JeUxC7Vqtf+wmGytvqwEylyE8PFE7kscC69QwWPwu5I+IcADnJ8\nRWdcDJhwONaNYFnxUSR7momZz2O8c41UNeYNM1FVFV+pyITGRuCtSAsh0eSvQteYMQIFgte36XX/\neDOGCrYV/d2deONtPsc0kbWZwsEgOb6ik6t+s5qty97ists2xGVQ9dfJAA6yptx/JrE3/ntNxb4e\no04mu+/BmEkHC+lqbNmyNDXmlLT6ZLBUNGZhf0wYiJ7WTq6p7kxLJ4NdY8YIFAhe36Y+9jlR7aeh\nITZDuouI+HWi2WGbrbVf4DL5JC3HTwtEJ92EQowb6aNn1CnHdkaGFCZeCaqLTfKJtHJaIfh0c4V+\n1Vg4Okt42cSvxSzrmKrGiou7fQtqpbGo+6ip+mtpleSDXWPGCBQIXt9mT7U0txm/dGm02T5kiAoJ\noA8rdYfNNdTupYWDLMY+j082mjmamrwZPLqPey6/TsVQSM0gDeElNTs0jZxWCD7dXCEdjbnuk3BY\n6aJHjTmzhAM4yCIteKCfxmINQazGlizxDUOk1h3GIYawRq07nEZJPtg1ZoxAnqP7Kr2+Tb0jOFEI\nAL157kZadH+rqVFNev/FXmIzaUlJNCx1Swv5y1+SRx+d2Ah576GtZa8KHdCLnJaKTzdVzOigeNLV\nmNd9os8vcDXm9gvU1Tl9Csv2xvj242v7sRpraiK/+12yrCxFjYVtttV+Pm7h+XSeQa5GCe0JYwTy\nDG+G7I2v0hsCQJ+ha9veVaTIqVPdmlx8C8D9CwRiI4a66aqoiE7O6bFzzZPTBqIzzmis7xrzuk/c\nioW3BaF3DtfXk0XBxPrKiMZyQF8DgTECeYRfDas3vkrXBeP6Wb0jM7ydxpblnQ0aiRUX2a6riz1e\n/62yUkVzrKmh7+iNVO41WxnVaKzvGvPqy9uC8OrDrYzoo9UApZeiov7R2GDv7E2HTBiBAAw5wdq1\nwLp1QGcnsH49IAJMmwYUFQH19epzOoio/6WlwKxZ6j8AhEJAXV10v2OPBaqqAMsCpkwBbrsNOOYY\nlf0qKoClS4Fnnok9vro6enw4DMycCbz6KtDVpe7j2WfTu9d169K7N0PvyKTGXH0BsRoLhdR5gkFg\nyBB17smTgYkTgUBAaWfpUuCII4Du7v7RmNFXehgjkCO4mcfNkDNmAE88ATz+uPrvZpCeWLsW2LhR\nZYD29tgM0NGhMtsttwBNTcB99wEHDgBbt6rMFQgAH/oQsHmz2n77bZV59WuXlgKPPQZUVka/6+xM\nnJ6ODqCtTf1PdK/pGjhD78iExnrS19q1wP33A7/8JfDUU8CPf6y0tHWrKvSHDgWOPlqdo780ZvSV\nJn1tSmTqDwXeVCfT66BK5PP0G+XhjuTwRor0DudzR3OkMlrCXSTcG81Ub9Yna5YPRGec0Vjqzz1V\nfbnuoObm6FKTQ4YoLQ0ZEu8GypbGBnNnbzrA9AnkD+l0ZPXk83QzgL6ghncij2XFZ1A3c9ntO9k2\n78HI+q7J0uFex2+kUq6NwS700UGpaixVfemdy35xfPz+Iv0IeaqxbGOMQJ6QbkdWKsK3bbWAuD58\nT6/5e1sCIuSiRSpzRhac1QOx9OG+BusY7HwiHY2lqq/W1tjOZVdnQ4ZEV5HTR/9cf70jJ31RY6Ox\nPmGMQJ6Qbm2mJ+HrGd7NkK5rSK9Nuc3tmP3Gvx275N+SJX26t0Jpluc66WgsHX2FQtEFaEKhaHA3\nN/JonL4aSHvhr2ObB0ZjvcYYgTyhN7WZZML3Zni/mZb6eWKX4OtmW/GpGamlGXKHdDWWjr78Zovr\n54lb4rHprYy1BAqdvDMCXV1dkZs7dOhQpp7ToCCTtZmEGT6BUzhu//adynLkQeb87W9/y3379nHL\nli288847I98fPnw4Rm+FQKY0ltSg+GjMd383LvUg19ju3bu5dOnSAbt+3hkBF9u2WVlZWXCZNJPE\nZfgenML50qS+8847+cYbb0S2582bx3feeYdbt27lvffeG/n+4Ycf5n/8x38MRBLzAl+9JNFYvuhr\n3759/PGPf8zu7m6S5Ntvv8358+cPWHry1giQ5Pvvvx/5/MILL/ArX/lKnx5WwZOnwyj27t3Ld955\nJ7K9ePFivvbaa777ekcH6Rq74447eNddd/VHEguHPNXYrl27Ip6Jrq4u3nTTTTx48OAAp0qRNSMA\ntdD8RgCbAMxLsM8iqIXoXwTwwXSOpY8R0Nm7dy/XrFkT2V69ejVXr17dp4eXC2Q1vkmeDqO48cYb\nY9w8yUimsZ07d8YYj0cffZS7du3qa/IGlKzHz8lTjZ111ll88cUXBzoZvmTFCEDNKt4CYCKAIqeQ\nn+rZ50wAf3c+zwTwdKrHaudI+cYffPBBPvDAA5Ht/fv3p/3wesOKFSsydq7+jG+SMJ051ibvzfN8\n6qmn+PWvfz2y7TbLUyEdjd18881sb28nqdKZLY31Bf15Dlj8nBQ0lsl81B/8/Oc/59133x1JZzoa\nyzaZMAKphI2YAWAzyW0kDwNYCuBczz7nAvidU5I/A2C4iJSleGzafOpTn8L5558f2T7vvPPw+OOP\nR7Y7k80x7wMrV67M2LnSjW/iF34hEQnT6Q0kNMCk8jzfffdd3HLLLZHt+vp6fOUrX4lsix7EJoN8\n73vfw+TJkwEALS0tqK2txf79+wGoilN/aawv6M+zN/Fz0tFYQlLQWCbzUSZ44YUXcM8990S2zzzz\nTJx55pmRdPaXxnKFVIxAJYAd2vbrznep7JPKsX3mr3/9K0466aTIdigUwvbt2yPbf//733HgwIFM\nX7ZPpBPfpKMDOOUU4NRT1f8+ZdIcgyQ2b94c2T506BA++9nPuq1DlJaWorS0NLI9YsQITJ06Natp\nDAaDePXVVzF06FAAwOuvv45p2gvr6OjA8uXLs5qmnkg3fk4+a+zgwYPYtm1bZHvt2rVYsGBBZPuI\nI47A8OHDI9s1NTUYO3ZsVtM4kPRXALmsmk7LsmBZVmR7zZo1mDBhQmT7D3/4Q6QQAYCxY8eiQ1P5\niSeeiPfffz+yffbZZ0dqfQDwyU9+Mmb7vPPOi9m+4IILYrbPP//8pPt/+tOfRjC4PxK8q7LyXASD\nia935pmfwNq1+yO1urlzz4n53Zvee+65x3P8mTHbc+fOxb59+yLbH//4x2O2P/axj8Vsn3766THb\nc+bMidk+7bTTYrYbGxtjtmfPnh2zfeyxx8Zs68+nuLgYn/nMZ9DV1RXZvvbaawe8NlZcXBz5PGHC\nBKxZsyay/c477+CRRx6JbL/wwgs48cQTI9vr16/HxRdfHNlub2/HV7/61cj2pk2bcPXVV0e2N2/e\njK9//esx29ddd13S36+99tqY7e9855qIvn7960244YavJbzepk2bcMklV0VaDuvWbcIll0TT197e\nHtP6am9vx5VXXhmzfcUVV8RsX3755ZHtjRs34rLLLot5Xl/+8pdjfv/Sl74Us33ppZdGtjds2BC3\nfckll8Rsf/GLX4xsv/TSSzjzzDNjzvetb30rsl1eXo7GxsbI9tSpU/GJT3wChYrohaPvDiKzANxI\ncq6zPR/KD3Wrts8dAFaQvM/Z3ghgNoBJPR2rnSN5QgwGg8EQB8k+1ZCCKezzHIBaEZkIYBeAiwB8\n1rNPE4CrAdznGI1/k3xTRN5J4VgAfb8Rg8FgMKRPj0aAZJeIfA3Acij30V0kN4jIlepnLib5sIic\nJSJbALwP4NJkx/bb3RgMBoMhLXp0BxkMBoMhfxnwlcVEZK6IbBSRTSIyb4DTcpeIvCkiL2nfjRSR\n5SLSLiL/EJHh2m83iMhmEdkgIh/PUhrHi8i/RGSdiLwsItfmaDpLROQZEVntpPUHuZhO7doBEXlB\nRJpyNZ0islVE1jjP9NkcTudwEfmTc911IjIz19IpIpOd5/iC83+PiFybg+m8wXmGL4nIH0SkOONp\n7OtEg778IY3JZFlKz8kAPgjgJe27WwF82/k8D8CPnM/1AFZDudSqnfuQLKSxHM6MbABHAmgHMDXX\n0ulc+wjnvwXgaQAn5WI6netfD+AeAE25+N6da78KYKTnu1xM590ALnU+BwEMz8V0aukNANgJYEIu\npROqXHwVQLGzfR+ASzKdxqw96AQ3OQvAI9r2fCQJLZGlNE1ErBHYCKDM+VwOYKNfWgE8AmDmAKT3\nIQCn53I6ARwB4FlHpDmXTgDjATQDaMT/3875vFhZRnH88wUVHX/QqISETCnhLkSREAdJkqQSZlsR\nkq5btDFI/4cQF21aJCLoIk2coIVJBCIEqTM1OBJiQjo0k0Ni2MKFfVs8z+tcp8mN432P3PPZ3Pse\nLtzPfe773vO85zzPnUkCET1vAKtmxUJ5AiuA63PEQ3nOctsFnI/mCfRXn/76wz78NK71tstBXdlM\n9oQ8b3sKwPYk0Owime0+QZfdJb1EuXP5gXJShPKsJZYRYBL43vZ4RE/gEPAx0Nkgi+hp4FtJP0pq\nFt5H81wHTEs6Ukstn0vqC+jZyTvA8fo8jKftO8CnwG/1/e7aPjffjm0ngWeREJ10ScuAk8BHtu/x\nX6/WPW3/Y3sTZaa9XdIOgnlK2g1M2R7l8ZscWx9PYND2ZuBt4ENJ2wk2npQZ62bgs+r6N2WGGs0T\nAEkLgSHgyxoK4ylpPaVM+SLwArBU0vtzOD2RY9tJYAIY6DheW2ORmFL5HyQkrQH+qPEJSg2xoWvu\nkhZQEsAx22eiejbY/gv4BtgS0HMQGJL0K3ACeF3SMWAymCe2f6+PtyllwFeJN563gJu2L9bjU5Sk\nEM2z4S3gku3pehzJcwtwwfafth8Ap4Ft8+3YdhJ4uBFN0iLKZrLhlp3EozPCYWBvff4BcKYj/m7t\n1q8DXqbUvbvBF8C47cNRPSWtblYtSFoCvEFpWoXytH3Q9oDt9ZTz7zvbe4CvI3lK6qt3f0haSqlj\njxFvPKeAm5I21NBO4Eo0zw7eoyT/hkievwBbJS2WJMpYjs+7YzcbMP/T/HizfthrwCctuxynrBK4\nT6nD7aM0Zc5Vx7PAcx2vP0DpwF8FdnXJcRB4QFlJNQJcrmO4MpjnK9VtBPgJ2F/joTxnOb/GTGM4\nlCel1t5852PNtRLNs77vRsoEbxT4irI6KKJnH3AbWN4RC+VJ6VVdAX4GjlJWUc6rY24WS5Ik6WHa\nLgclSZIkLZJJIEmSpIfJJJAkSdLDZBJIkiTpYTIJJEmS9DCZBJIkSXqYTAJJkiQ9TCaBJEmSHuZf\nTxHrNFYnpJwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe056634ac8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# NO2+NO3 at Hope and Gravesend\n",
    "dataG=session.query(Profs.Year, Profs.Month, Profs.Day, Profs.value).\\\n",
    "            filter(and_(\n",
    "                Profs.value<2.5, # remove single high outlier\n",
    "                Profs.variable_name=='NITROGEN DISSOLVED NO3 & NO2',\n",
    "                Profs.station_name.like('%Gravesend%'))).all()\n",
    "YDG, valG = plotYMDV(dataG,'r','temp')\n",
    "\n",
    "dataH=session.query(Profs.Year, Profs.Month, Profs.Day, Profs.value).\\\n",
    "            filter(and_(\n",
    "                Profs.variable_name=='NITROGEN DISSOLVED NO3 & NO2',\n",
    "                Profs.station_name.like('%Hope%'))).all()\n",
    "YDH, valH = plotYMDV(dataH,'b','NO3, red=Gravesend, blue=Hope')\n",
    "\n",
    "units=session.query(Profs.unit_code).\\\n",
    "            filter(and_(\n",
    "                Profs.variable_name=='NITROGEN DISSOLVED NO3 & NO2',\n",
    "                or_(Profs.station_name.like('%Gravesend%'),\n",
    "                    Profs.station_name.like('%Hope%')))).group_by(Profs.unit_code).all()\n",
    "for row in units:\n",
    "    print('units:',row)\n",
    "    \n",
    "YD=np.concatenate((YDG,YDH))\n",
    "val=np.concatenate((valG, valH))\n",
    "no3dict=gsmooth(YD,val,60)\n",
    "#print(i,np.max([v for v in no3dict.values()]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "units: ('MG/L',)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAEKCAYAAADkYmWmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXmYFOW1/79nFlDCqMCAyC5hnRlAMCiGba5bECNqorji\nrphE45W4IDEB4o0xN0FRExNBrjEmLtH8ENwFZRQFV3BhkSgoyKAo28yALDPT398fVTVTXVPdXdVd\n3V09cz7PM890dVdXnar61unznve8bwlJKIqiKLlJXrYNUBRFUZJHnbiiKEoOo05cURQlh1EnriiK\nksOoE1cURclh1IkriqLkMOrEQ4SInC8iL2Rwf0tE5LJM7S+XEZHpIvJwEOuKyGcicnxw1iktGXXi\nGUZERonIGyKyS0S2ichSETkaAEg+QnJctm1MhIh0FpE5IlIpItUi8qmI/J+I9Mu2bWnGz6CK0AzA\niPWjIiIREemdDZuU4FAnnkFEpAjA0wDuBtAOQFcAMwHsD2Db+aluw+N+2gNYBuBgACNJHgJgGIBX\nAZyUTduUuLj9qITmh0ZJHnXimaUfAJL8Fw32k1xMchUAiMjFIrLUy4bMdV8XkTtFZBuA6eb7l4nI\nGhHZLiLPi0gP23dOEpG1IrJTRO4FIEkcwxQAVSQnkfwcxgFVk3yI5J/N/fQ0o7zLRGQjgJfN9/8l\nIl+a+68QkRLz/WPM9xvsEZEzReQD87WIyFQz4v9GRB4TkcPMz1qLyMNmq2aniLwlIh3Nzw4RkQdE\nZIuIfCEit1n7sM61iPxBRHaIyHoRGWfbfy/TxioReRFAsc/zdLBpZ7WIvCsig91WEpEHReQ3tuWx\nIvKFbfkIEXlSRL42bbzWpx3xsJ/vViIy22xdbRaRu0Sk0G6TiNxinv8NInK+47t/FJGN5nW8T0Ra\nB2inEgd14pnlPwDqReRvIjLOckQO/ERHxwL4FEAnAL8VkdMBTAVwBoCOAJYCeBQARKQYwL8BTIPh\nkNYDGGltSERGmk5wh/nf/nqHiHzfXPUEAPM92jcGwAAAPzCXnwPwXdPeFQD+CQAk3wawG4A9T3we\ngH+Yr38OYAKA0QC6ANgJ4D7zs4sBHAKjVdMewNUA9pqfPQTgAIDeAIbCaClcYdvHMQDWAugA4A8A\n5tk+ewTAOzDO1f+Y+/HDBACPw2hxPQrgKR8tEgLGjxeMlttKAEfAOPfXichJ5ufnJbhm3XzYeyuM\n8zEYwBDz9a22zzvDOL9dAFwCYI6I9DU/+z2APuZ3+8C4Fr/2sW8lFUjqXwb/APQH8H8ANsFwMAsA\ndDQ/uxjAax63czGAzx3vPQfgUttyHoA9ALoDmARgmWP9LwBc5tP+TwBcZVs+DYZTrQbwgvleTwD1\nAHrG2c5hACIAiszl2wDMM18XwXDq3czlNQD+y/bdI8xzlwfgUgCvAxjk2H4nAPsAtLa9dy6AV2zn\n7z+2zw427elknq8DAA62ff5PAH/3eI6m2881jIh3C4z0EwB8BuB48/WDAH5jW3csgE3m62NdrvFU\n6zz5uGbTYaTsdtj+dprXqLe5zqcAfmD7zskANthsOgDgINvnjwP4pfl6N4AjbZ8dZ31X/9L/VwAl\no5BcB+AyABCjI/CfAGYDuCCJzX3hWO4J4G4RmWUuC4yoriuMCMq5vnPZC9thOFEAAMmnAbQTkcvR\n9Bg2Wy9EJA/A7QDOghHd0vwrBlADI/J9Q0SuBvAjAO+RtL7fE8B8EYnYjqsWwOEAHgbQDcBjInIo\njOj9l+Z3CgF8aWVQzL9NNvu+sh3HXnO9tjBaMTtJ7rWtu9Hcj1cazi1JishmGNfADz0AdBWRHeay\nwPjhes3ndgDgcZIX2d+wnU+YttnPzUZE27uT5D7n52bqqg2A92zZsDwkl6pTkkDTKVmE5H8A/A1A\nWbKbcCxvAjCZZHvzrx3JtiTfBPAlDKdgp7v1QoyqmRozh2v/s96zUi8vw0jX+LXvfBhR+/EkDwPQ\nC42OFSTXwnAM42GkUh5xHNcpjuP6DskvSdaRvI1kKYDvm/u4CIYT3Qegg+07h5F0zU07+BLGD9PB\ntvec5y4R9nMrMH4AKl3W2wPDCVocYXv9BYyI1n7ch5I8zdzu+QmumZ8fnS0wfvgseprvWbidjy0A\ntgH4FkCpzc7DSB7qY99KCqgTzyAi0l9EpohIV3O5OwyHtTzG+ktExE9u8X4A06Sxw/BQETnL/OxZ\nACUicoaI5IvIdTAiWQAAyddJFpE8xPFnvfeGueqdMG7oh8UsTxOj6uYop/mO5SIYTfqdIvIdAL9D\n0x+hRwBcByP3/YTjuG4Xs5NWRDqKyATzdbmIlJmR/m4YEXo9ya8AvATgLhEpEoPeIjIm0UkkuQnA\nuwBmikihiIyC8ePQeHBGrfdFrhswONo61wCuh/GD8pbLeu8DGC8i7USks3n8Fm8DqBGRm0TkIPO6\nlYrI90w7H0lwzTa77C8WjwK4VUSKzf6TX8Fo5TQcsu18jAZwKoB/kSSAuQBmS2OHclcROdnHvpUU\nUCeeWWpg5DnfEpEaGKV6HwK4Icb63WHkez1B8ikAd8BILewytz3O/Gw7gLNhdEJtg9HB+EaMTcXb\nx3YAI2A4pddFpBpGJ2VbAD+xr+r46t9hRNSVAFbBOHYnj8HoDH2Z5A7b+3fD6Dt4SUSqzO8eY37W\nGcCTAKoArAawBI0dohcBaAUjp74Dxg9D53iHZ3t9gXmc22E4tIesD0SkFYxOvjfjbGsBgHNg5J4v\nAPAjkvUu+3kYxnX6HMALMM6BsRIZAfBDGD+QnwH4GobDPCTOfv1gt+N/YPxwfQjgA/P1b22ff2ke\nyxbT5skkPzE/uxlGTv1NU3cvwajEUjKAGD+kCVYyojarV38uyXvSapUCM1p/nOSobNuiRGOmln5K\nMpl+jJxDRMYCeJik35SSkgESOnERKYXR1BoOoA7A8wCuJrkh/eYpipJt1ImHGy/plIEA3qIxMKUe\nRs/4j9JrlqIoiuIFL058FYDRZsdLGxjVA90TfEdRlGYCyVc1Cg8vCevESX4sIr8HsAhG7/9KGIME\nFEVRlCzjqWMz6gsivwXwBcm/Ot7XyXQURVF8QjKlgVGeSgxt9Z89AJyJ6IEYdmNC/Td9+vSs26B2\nqp1qp9pp/QWB12H3/xZjCtJaGKVV1YHsXVEURUkJT06cZMJRboqiKErmaVEjNsvLy7NtgifUzmBR\nO4NF7QwXvjs2Y25IhEFtS1EUpSUgImAmOjYVRVGUcKJOXFEUJYdRJ64oipLDqBNXFEXJYdSJK4qi\n5DDqxBVFUXIYdeKKoig5jDpxk5oaYPly47+iBI3qS0kX6sRh3FijRwNjxhj/9UZTgkT1paQTdeIA\nVq0CVq8G6uqANWuM10r4mDFjRrZNSArVl5JOdNg9GiOlNWuAkhJg6VKgqCjbVilOzCHK2TbDN6ov\nJRZBDLtXJ25SU2NESKWleoOFlVx14oDqS3FHnbjSoshlJ64obugEWBlCKwvCwfTp07NtQtpQjSnJ\n4ikSF5FbAFwI4wHJHwG4lOQBxzrNMhK38plWU1jzmUrQqMZaLhmJxEWkJ4ArAQwlORjG04DOTWWn\nuYRWFijpRjWmpIKXdEo1gAMAviMiBQDaANiSVqtCRFmZER0VFhqVBaWl2bZIaW6oxpRU8JpOuRLA\nnQC+BfASyUku6zTLdAqglQVK+lGNtUyCSKckfFCyiPQGcD2AngCqADwpIueTfMS5rn0wRnl5ebN5\nxl1RETBiRLatUJozqrGWQUVFBSoqKgLdZsJIXEQmAjiJ5JXm8iQAx5K8xrFes43ElXAwY8aMnB21\nqShuZKROXESGAPgHgOEA9gN4EMA7JP/sWE+duJJWtE5caW5kpDqF5AcA/g7gPQAfABAAc1LZqaIo\nihIMOmJTyRk0EleaGzpiU1EUpYWjTlxRFCWHUSeu5AzNee4URUkWzYkriqJkCc2JK4qitHDUiSuK\nouQw6sQVRVFyGHXiiqIoOYw6cSVn0HlTFKUpWp2i5Aw6YlNpbmh1So6jz1VU0o1qrPmjTjxLWM9V\nHDPG+K83mRI0qrGWgTrxLKHPVVTSjWqsZaBOPEvocxWVdKMaaxkkfDybkh6KioClS/W5in7QuVP8\noRprGWh1iqIoSpbISHWKiPQTkZUissL8XyUiP09lp4qiKEow+IrERSQPwGYYD0r+wvGZRuKKoig+\nyEad+IkA1jsduKIoipId/DrxcwA8mg5D0o0OelDSiepLyRaeq1NEpBDABABTY61jn9uivLwc5eXl\nKZgWHNagB6uXfulS7anPRWbMmBHK+VNUX4pXKioqUFFREeg2PefERWQCgJ+SHBfj89DmxJcvN0at\n1dUZNbOvvQaMGJFtqxS/hHXuFNWXkiyZzomfhxxNpQQ96EGbzoqddAyqUY0pXvEUiYtIGwAbAfQm\n6SqrMEfigHEzBDHoQZvO2SOskTgQnL6sbanGWgZBROI62Mcn2nTOHmF24kGiGms56FS0WUDno1DS\njWpM8YPOneITnY8ie7SUuVNUY4ofNJ2iKIqSJTSdoiiK0sJRJ64oipLDqBNXFEXJYdSJK4qi5DDq\nxJWcIYzzpihKttHqFCVnaCmDfZSWg1anKIqitHDUiSuKouQw6sQVRVFyGHXiiqIoOYw6cSVnaClz\npyiKH7Q6RVEUJUtodYqiKEoLx5MTF5FDReQJEVkrIqtF5Nh0G6YoiqIkxut84ncDeI7k2SJSAKBN\nGm1SFEVRPJIwJy4ihwBYSfK7CdbTnLiiKIoPMpUTPxLANhF5UERWiMgcETk4lZ0qSjLo3CmK0hQv\nkfjRAN4EcBzJd0VkNoAqktMd69FeAlZeXo7y8vLgLQ6Qmhpg1SrjmYb6CKzwk4tzp6jGFDsVFRWo\nqKhoWJ45c2b6n3YvIocDWE6yt7k8CsDNJE9zrJdT6ZSaGmD06MbnGC5dqjdZ2Mk1J64aUxKRkXQK\nya0AvhCRfuZbJwBYk8pOw8CqVcbNVVcHrFljvE6Wmhpg+XLjv6JYqMaUTOC1TvznAP4pIu8DGALg\n9vSZlBnKyozoqLAQKCkxXieDFW2NGWP815tMsVCNKZmgRY/YrKkxoqMePYCNGxvzln7ymMuXGzdX\nXZ1xs772GjBiRGbsb2nkWjoFUI0p8dERmylSVGRER+PHN0Y5W7b4i3qCiraUxOTi3CmqMSXdtOhI\nHGga5UyZAsya5S/qsaKt0tL0dlxppUNukisaU31lniAi8RbvxLdsAcaOBT7/HCgoAA4cAFq1Aurr\njajHT0VBOm8CrXTIXXJBY6qv7KDplBSpqTGauZ9/DnTubNxckYgRId13n/+bK52dT0FWOiiZI1c0\npvrKXVq0E7cLd+tWoHdvo3lbWgqcc46/SCTdN4HmRXOTXNGY6it3adHpFCuyWbPGEO5zzwGbNiWX\nd3RuKx3N0Uzl3pXgyCWNqb4yj+bEAyBI4epNkF5mzJiRk/OnqMaUWKgTV1oUuVgnrijx0I7NNKND\nnZV0oxpTUkWdeAz8VALojagkg1eNqb6UeKgTj4HXSgCd10JJFi8aU30piVAnHgOvJVdaX6skixeN\nqb6URKgTj0FRkVHC9dpr8Uu5tL42c+Ti3Cnx8KIx1ZeSCK1OCQAt+1LSieqr+aIlhoqiKDlMEE68\nwOOOPgdQBSACoJbkMansVFEURQkGrznxCIBykkPD7sCDKsfSsi7FjSB1oRpTgsCrExcf62aNoGq7\ntaxLccOvLlRjSibw6pgJYJGIvCMiV6bToFQIqrZby7rCSbbnTfGjC9WYkim8OvGRJIcBGA/gZyIy\nKo02JU1Qtd1a1hVOZs6cmdX9+9GFakzJFL6rU0RkOoAaknc63qe9jre8vBzl5eVB2OgLZzmW25NQ\nvEzpmc2yLn1MljthmADLTRe5pjHVV/aoqKhARUVFw/LMmTPTX2IoIm0A5JHcLSLfAfASgJkkX3Ks\nF7oSw3iPnApr7a0+Jis2YXDiTnJNY6qvcJGpWQwPB/C6iKwE8CaAp50OPKzEa9IWFRkPpw2bgDVX\nmlvkmsZUX82PhE6c5GckjzLLCweRvMPrxrNdQpWLecegbM72uc8U2T7OXNNYkPZm+9wrJiQD+TM2\n1Uh1NTlkCFlQYPyvrmZWqK4mly/P3v7jUV1NLlvW1LZUbQ7LuQ+a6dOnRy2H5ThzTWNB2BuWc5/r\nmH4zNd+b6gYaNuRw4suWGRcYIAsLDdE4ieXEWgLpvAm8nPvmgGosPqqx8BOEE0/bAJ5EzTYvgx2a\nc3NNn1yeOqqx+KjGWgip/gpYf3BE4mT8ZluiX/Lm3lyzjq+wMD3HF+YmfpCoxmKjGgs/CCASz9os\nhonqaJcvNyKoujrj1/6114ye/uZEGEvQmhOqMdVY2Mn5qWjjCczLYAlFSYRqTAkzOe/EE5GuKIIk\nVq9ejaeffhofffQRdu7ciYMPPhh9+/bF8ccfjxNOOAEFBU1n6Y030k1HwaWfGTNmBD5/Sro0VldX\nh6VLl2LRokVYt24d9uzZg3bt2qGsrAynnHIKhg4dCpHoezeRhlRjzY9m78SDpKYG+PDDCNav/zfu\nuut27NixA2eccQaOPvpoFBcXY8+ePVi7di2ee+45bNmyBf/93/+Na665Bq1atWr4fryReToKLv2E\nccSmnZoa4J139mD58r/gz3++E126dMGpp56K0tJSFBUVYceOHVixYgXmz5+P4uJiTJs2DaeffjpE\nJKGGVGPNkyCceFo7NjNNvLrrfv0+IvB9Hnzw9/j4408zEonE3M67777L8ePHc8CAAXzjjTdIRneS\nFRSQixc3rq/lVpkh7Brr1espAj146KFn8Y033o+5jbq6Oj711FMcOHAgx48fz82bNzfR0Jw50ftQ\njTVPEOY68UwTq9IgEonwF7/4E4GOBP7KgoJ6T/XEkUiETz75JDt16sRZs2axqirCsjLjjAFkWRlZ\nWWl8p7IyvVUAbva1RMKqsW+//ZannXYFgT4EXolZCeO8fvv37+fMmTPZuXNnLliwqEFDBx1E5ucb\nGlu0yPhOJipNWrq+soE6cZPqavL++5tGKjU1NTz33HM5ePBRHDDgExYWkn36GE7X+f1YpWafffYZ\nhw0bxsmTJ/OFF+qYn9+4jz59Gr9TWdm03CqoG6O5l8J5JYwa+/TTT3nUUUfxxz8+l2Vl1a4aS3T9\nXnnlFXbu3Jl/+cuDnDOHDRoDjNfWd9xK+oLQmOore6gTZ6MA8/ONCMaKVD755EsOHTqUl156Kb/9\n9ltWVho3l/2msEhcT1zN448/nqef/mMOGrSv4Ua1O/R01iBrU9ogbBpbvHg5Dz/8cN57772MRCIx\nNebl+n388cfs2bMnb7vtfxt0YznyeKNRg9CY6it7qBNn01z13LnkihX/Ye/evfmb3/ymIfe9bFmj\n0y0oiBaql6bqvn37eMYZZ/C0087k0qW1CVMoQd4Y8exrSc1g59wpmcJNY48//jSLi4v57LPPRq3n\npjGvqZAvvviC/fr1429/+0cuXmykU+J9JyiNJbKvJWks06gTZ7QAe/cm//jHt3n44Z05d+7cqHUW\nLCBbtzaO+KCD3FMqiUaf7du3j+PGjeOkSZNYX18f9ztB5zBjNaW1GZx+nBq79tp5PPzwznzzzTej\n1omnMa+jGzdt2sRevXrx/vvvT/idIDUWa1+qsfTSIpy4PQqIFRFUVho3F/A8gWL26LGgYR27CBM1\nT+Pt22LPnj0cPXo0f/azn8WtcLG+b90Y6YhmtBkcDF41duSREQK3ETiSffuuS5vGPvnkE3bt2pWP\nPfaYp++rxnKXZunEnTeUdXOUlRl/bhHBsmWkyEMEOhF4o6EEcNkyo3fffnO55cRj2RErAqmqquKg\nQYM4a9YsT8ezaFFs21Mh3RULzZVkNLZ0aR1FfkpgCIEtadfYBx98wI4dO/L111/3dCxWek81lltk\n1InDeIDECgALY3ye8gE5RW2/OQoK3DsSq6oivPzyO1hQ0IPAGgJk//5GZG6/MQsLjf+LF3sTYqII\nZOPGjezatSvnz5/v6Xj8Rmixthdrbmir3LE53GTpzMEmo7GtW7/lqFE/Yps2/0VgFwGyVy+ypCS9\nGnv++efZuXNnfvrppwmPJVFHu5/z46axxYsbyx1znTDl+DPtxK8H8I90OnGnqBcvbowC7DeKFRHs\n3FnH4uJrCQxi//6buXAhuXChlVqJ3o7f2da8RCDvvPMOi4uL+c477yQ8Hj8RWjx73CKt5pS3TPex\n+NXYxo07+J3vjKLIOSwp2cfHHjP0lZeXGY3dd9997N+/P7dv357wWPr0SS1ijnXuVV/pI2NOHEA3\nAIsAlGciErcL0Znzs6LOJUv2csyYswmMJbCzIfqwVwgAhrCTvVBeOqPmz5/PLl26cOPGjXGPx0+E\n5ka8qK055S3jHUsQ1SleNGZFnWvWbOKRR5ZS5HoC9Q0jKe0/zKk6Ay8amzJlCsvLy7l///64x+I2\nVsEPsc59S9FXNsikE38CwFEAxqbTiZOJRV1dTZaUbCUwmoccchZLSvY2uSHtTUxnFUo6mDVrFgcN\nGsSqqipXe4OYczle1Bb2vKWf5mu8Y8mExqz95+e/w8LCbvzVr/7YxFEG9cPslbq6Op5++um85JJL\nmnSmB6Uva1tu576l6CsbZMSJAzgVwJ/M1+UwnnafNifuxN5xs2wZed997xHoQeBXzM+v5z33NL2R\nvPwQxLvofnNmkUiEkydP5imnnMLa2lp/B+iDRM4njLnxZJqvsY4zHRpzdnLefz8p8jCBYublzeec\nOU0j3FT15XUdO7t37+bQoUP5u9/9zt8B+iTWsTlbK2HRWJD6ygaZcuK3A9gEYAOALwHsBvB3l/U4\nffr0hr8lS5akfICNUZFRd5uX93fm5xeze/cnWFBgvOe3Oet20Z03cjI5swMHDvCkk07i5MnXZk3g\nYcv3kcE2X4N24vbzVVZGlpTsp8gUivRmXt5HadFXrHW8sHnzZnbv3p1///sTqjGTsKVHErFkyZIo\nP5nxEsNk0ynJ/nI3XqBdBM4nMJAFBR9w8eLo3KSfi+ccfXfPPdFlZfZqBb+i2LRpJ1u3LmFe3r1Z\nEXiQDw4OKtoKsvkaS2Op64vMy/uYwDACp7GgYBunTk1dX4WFRke7cyh+Ko5n6dIVzM8vZn7+W6ox\nhi894peccOKp/HJXV5Pf/e4SAr2Yn/9TFhTsaZL79nvx7PZYs8U5qwyS3e7dd5MiGwgcwfz8Z6OG\nXWcickp0TrxeC7/XzEt6Kojmq5vGUtXX4MH1zMv7C/Pzi9mly19YUBBpkvtORl9WztxeKWUNxU9W\nu5WV5E03kSILCXRhQcHnDdvLVGQeRo2FKT3il5wY7ONW0hVrPmb7+1u3buWkSZPYrVt3/uEPT8ec\nJTBW/s6eR3f73FllkGi2uHi2VlczappaYBnz8zty2bIPMt78jGe71wjQT6SYqeOrriYvu2x6k+27\n2Rrrhne+v3LlSg4fPoJlZcfxnXfWNDl38c5ldbXRanOrnba+5xwEZK+U8quxykoj6ABIEVLkLh50\nUBk//rgq4+mN5qixTP4Q2skJJ+6MTAYONBzmwIFN50ouKCAHDdrL3/3uLnbs2JE33HADa2pqfJ0U\nt0g7XsRg2bVwobfBDG6CcqsH/+UvH2WPHj24cOGW0OTsvEaAznMT77xkIicZ7yZ2HlNlpWFzPI2V\nlHzFq666hh07duTcuXNZX1/v2x7n3PKxHL1lW58+5Lp13tMMzuO9/357kEDefHOEl19+NUeMGMf8\n/NpQ6Mtuey5pLNOBlp2ccOJk4y/3ggXRQrQc7IIFZF7eAQJzCHTn6NET+NFHHyV1UpwONVHuzqro\niDXIwXnTxYr83G7q3/72tywtHcTS0h2hydl5bXpWV7NhJr144vbSvPYSFccj8VTBjce0aFE8je0g\ncAuB9pw48Tp+/fXXiXeewB5rH7EcSzyNxToHbsdrj8StybVqa2t50kmnsF27i1hQUB8KfZG5p7Fs\ndo6G2om7nUDnDWbcADvZvv2dBI4kcAJ7917uSYjxLpx1wb1WF8RyzLEce6xa2sWLo8sdI5EIp0yZ\nwmOOOY4vv1zjSdRhKd0ivYs7Xlor3jn0kw+Ndc69aWw9Dz10CoEOBK5gv34bU9aYl0jcTqyRol5a\nF/aUyty50WMf9uzZw+OOG8Wzz76WVVWJJ2QLk77IcGjMa+shHYTOiVsiiRfZWk3dwsI1zMv7KfPy\n2lHkfALLGyYVSkSiC2Q51AULvA3GcLuI8cTlNdIgDUd+2WWX8cQTT+S+ffuSPqZskKq4gxoBaEXY\n9pRXvJu3rIzMy4uwoGAxRSYwL6+YIjcR+DxQjS1caFQ3eRlQ5jyXiaqg/Ghs165dHDp0KG+99dak\njydbhEVjlZVsGBOQSULnxC2RxJqMp76+ng899AwPPvhkAoezXbtf87XXKn1fRC/Na7+Cdd40znxm\nKhe3rq6OEydO5Lhx47hnz56kjilb+HEmbt+NFUF7vebWutZYgUSloLt37+Ydd/yVrVqVEihlx473\n8/3394ROY/ZzkCgn7IWtW7dy4MCB/NWvfuU6RXJY9UWGR2OaEwdiTsbzxRe7OHv2bH73u9/lgAFH\nMy/v7wT2RaUunA40UclavAsUlGDdHreVbJO0traWF154IceOHctqly/HO6YwNoO9Eq8Z7OXGjc4/\nT49ZCvrRR5/xxhtvZIcOHTh69OnMy3uZQCTUGrNajM6ccLLXe+vWrRwyZAinTJniOjw/mbxyLhCk\nxlp8Ttze8fLee+Rtt63jpZdew3bt2vGcc87hG2+8waqqSMxfzkWLjBSIl7m3412goByinzymF+rr\n63nllVfy2GOPdZ2Vzrqp7VFZWJvB8QjSIdgjcQBRHXubN0f4i19U8JRTzmT79u05ZcoUrl+/3vX6\n21N9Xud3T0Zjfo89aI1t376dw4cP59VXX826ujrX43GW3qrGNCfe4MSNGy3CvLzFbNPmhwQ6slOn\naVy79osow92ioug66+hfxVgXLN6FdLsB/YrVbx7TC5FIhDfeeCP79u3LdevWue4vVvli2JrBbng5\nx16iYKc77uE1AAAdmUlEQVSTmTOn0Ynn5+/jTTc9xIMOOopAf3bt+mdu2VLTZBvO1EWsAV5BaSzZ\nFEvQGquqquLxxx/PU089tUmrTzXm/nkqKZ1UCJ0T79p1HoFBbNWqhHl5cwjs8SQKt7JAK4URr5PU\n7w2TjFjdnEEQv9hz5sxhp06d+Morr8S1L5tRQjKkmku2UlhW34rVF1FZaQg+L+83LCjozGHDTmJe\n3nO0poiNdy3d9JUOjSXrDNOhsQMHDvCqq67ioEGDoqZJVo01Ps4xP99bZVE6CZ0Tb9t2HO+440Vu\n3hydMkk0s54zEh84sLGqJMg5joMQaxC/2Fbq6A9/eIUdO3binXfeyUgkEreJHvZhxfZ0RaxUhvN6\nOitFqqsNx213tH36kK+++jEPOmgQAfCQQ67gyy9/FHWuEnUM2te1Okft08gGWUUThDMM4npXVpJ/\n/WuE06ffySOOOIKLFi2Ka2NL0ph9KgSAXLQowpdffpkXXnghDxw4kNFjCp0Ttw+rt4bJex1IU13d\ntM7aet+tSsQesfm5YVIRaxB5OOcPVr9+Gzhs2HCedtpp3LZtW07cTGT0uXBGPvYpEtw+i1VjvWxZ\ndKoDeJ8iZ7GoqCNFLqalsTlzGrcda7BIrOZyrOkbgtJYqtcvCI05BwY9/vjL7NKlC2+99VbW1tbm\nhMbcrl9QGmt8KlOEwEL263cM+/fvz3nz5qkT793biKLtol+0qGm5od9mqrNKxP7DkKkHPySTvnHD\nbbTfa6/t55QpU9i9e/eGiCnMOM9FvDyuWzTrpgmysZkLvEvgdAKdWVDwR65cuZtDhpB5edPZqpVx\nEwapL2vfzUVjziH6c+eSX331FU888USOGjUq5jM7w4LfvL1fjfXsWU/g/xEYSmAwu3d/krt2+Zt6\nIShC58TdZgR0+0X020x1XpBkp6FNhaA6f5yRuD1KeP7559m9e3dOnjzZ9SlBYSFWRUWsaiC3pq/T\nQRrN3DcJnEqgK4G7CXzLgoLGCLpbt8bzZjWTg9AX2bw05jZEnzSqo2bNmsUOHTpw9uzZvueMyRR+\n8/ZeNbZzZx27d3+cQBmNaYef8tSnkk5C58SdHUbOJ4lbuSk/uUM3p5fsNKGpEGTnT3V17NF+u3bt\n4hVXXMEePXrwxRdfTM3oNBHrpolXjuf8zJ6q6N59FYcP/yGB7gTuI7A3SkdOx2zlyoPQl7V+c9NY\nZaWhrwULmm5n3bp1HDlyJEeNGtWkQioMJJO3j6Ux46HWEfbs+SyPPHIQgWMIPGumUlJ7eHkQhM6J\nW51MVl47Vq6RdD/pVvPY3kkVq5MiG3m9oPbppdn8wgsvsGfPnjzvvPP45ZdfprZDxs61JpuD9XIu\n4m3byH9vJnA5gY5s334W+/ff59q5bW3LqSX7e87OzVj2qcYM6urqOHv2bHbo0IEzZszg3r17U9sh\ng9WY1/MQb9vGPDpvEygn0J/Tpj3F0tJITI1lg9A5cbe5B9xGPboRK80QZHQSBqqrjZylW77Oye7d\nuzl16lQWFxfz3nvvbTJ4w88+E5XQBTH029pmvPlzSKO1cc01xmyCwM0EdjY4T7fObQu3+S2qq73N\nhGetqxqLZtOmTTzzzDPZp08fvvTSSyntM0wa+/TTTzlq1DkEuhC4n0AtFy5s1Eu2nbdFpp6x2RrA\nWwBWAlgN4PYY67mezGQmibeaOda62YiI4pFK9GqJzs/zG1evXs2xY8fy6KOP5ttvv+3bXi8ldFYU\n6qcKw3kOnHlIZ9781Vf38447ZrNTp04cP/4SApui0iP2Gz/etlVj8b+XjMaeeeYZHnnkkZw4cSIr\nk+jFDYvGnn32a06efC3bt2/Piy/+DYHdDft2lhqGYZqBjEXiANqY//MBvAlgpMs6rjeR1ygnVpQU\nNrykQmLhFPrUqd6rHiKRCB966CF27tyZP/nJT7hp007PIoyXYxwyxH0Eo5ftOc+BMy1hzZ9TWlrP\nbt0eJXAki4rGc/nyD6M63woKjGkaEm1bZLpqLAGpaGzPnj385S9/yeLiYs6ePZs7dtTmjMZKSnaz\nU6f/IdCBHTpcw/Xrt7KykmzdunGfVvo/lfMbNNl4xmYbAG8DKHH5LO7kOn4miQ9LU8dJomZqol93\nu9DdnjrkJTrYsWMHL7/8ahYUdGZe3jwOHlwfd3/2Wu5YOWIrHeE1nRAr6rKcpDUSrrKSvPfel9m3\n79EEvkfglYb13ToqLXtjbVs1lhmNrV27lmPHHs+DDhrM/PxXY2rCua1saGzjxlrecstcdujQhSJn\nE/gkaY1lg0xG4nlmOqUawP/GWCdUzdGgSdRMjffr7nSmc+akVttsdAq+S2AERY7hAw+8FdderxUa\nsW7AWE1at6jLyk336fMBTzrpFPbu3ZsPPvgYBw+ub1LNYh+dmZ/f6FhjRdWqscxp7I03IszLe5xA\nd4qcxwULNse01U96JCiN5edH2KPHAvbvP5Bjx47lK6+85VoxZdeYVa7qteWWCbIRiR9iplPGunzG\n6dOnN/wtWbIk3cefUZy/3nPnRl/8eJGDU+xuIvITHTRus56dOv2Nhx9+BC+77DJu3bo1oT1+SOQ0\nnDek8ePyGYGLCHTi9dffzf3795N075RsHNwTnd6IdbPHewRgcyCcGtvNdu1+yXbtOvCOO+5oeLBJ\nkOMm/GvsdQKjCJTwj398pmHaXTeNrVtHtmpl2GmvmffacguaJUuWRPnJrFSnAPgVgF+4vJ/2E5BN\nEv16Oz+35ouJNZrRKSIv0YE9YrFX/ZSV7eLPfnY9i4uLeffddzcMrU412og16MKtOb5161ZeffXP\nmZ/fnnl5v2ZZWZWnCDLWyDq3/ajG0qsx5zm3a2zAgE84btwP2bdvXz7//POBRbN+NPbhhx/ylFNO\nY2FhD+bnP8jBg+sSasw+lYMVicc63myQqeqUYgCHmq8PBvAagBNc1svEMWcV+00RqwnonC+mrMx7\nLtD+/UTVGW437urVq3nCCSewrKyML774IquqIinP4+F0GnYbKivJF17Yzptv/jXbt2/Pa6+9lp9+\n+lWTfcZy1G77sEeRzhtSNZY+jbmdczcH+8wzz7BPnz6cMGEC33vv40Amg4ulMask8b331nHSpEns\n1KkT77rrLn799V7PGouXlglD52amnPggACvMnPgHAG6IsV4mjjkUJBKA27B0r2KPF1F4GeoeiUT4\nxBNPsF+/fhw9ejQrKipSPlZ7dGTZkJ9fyUMO+QWBdjzssEv54YcbYn4/UUWIM2KM1VSfPn16SseS\nS2RaY36Guu/bt4933HEHi4uLedFFF6U8F0ssjQErCJzNvLxi/uIXM2JORZFIY7HSMmHo3AzdYJ+W\ngtf5jP02Naurybvv9h5RxMvr1dbW8qGHHmKvXr05bNh/8ZFH5rO2ttbXcTojwaqqCHv1WkLgfALt\nCFxHq9Y71sOH400LGm+/Yel4yhbp0Jh1PRcs8O6w42ls165dnDZtJg89tAPPO+8i32MY3Foa33yz\nj926PULgv2jMoTOLQE3UOAInuawxdeJZwmv+2k9T0xlNAEanX6KpCmJtyxrRNnjwAebl/YNt2hzH\n7t17cNq0aXz77bddH6Zrv6msY8zPr2ffvm/xhhumsW/fvuzZs4Qiswlsj7LVmgfGrflvVaz4qcv2\ne/6aG0FrrPF6NtZOOzXmd3vWY+7y83fwiCN+z549e3H48OG88847uWFD7JaZc7TloEF7+cQTz3Ly\n5MksLi7msGEnUuRxAvujqpfuucd9xGcuaywIJy7GdlJHROhnWzU1wKpVQFkZUFQUiAkZpaYGWL0a\nKC0Nxv7ly4HRo4H6+sb38vKAQYOApUu976OmxtjO6tVAr17AZ58Z2ywsBObOXYmPP/4X5s+fj6qq\nKhx77LEYOnQounfvjrZtO2HqVMHGjXXo0uUrjBq1EY8/vhLk2wAOx6RJE/DTn56BkpJjMWaMYM0a\nw74DB4BWrYx9FBQAtbXGNbVsrqkBvv99YO1aYOBAYNmyzF1v1Vgjy5cDY8YAdXXR76eqMWt7hYXA\nkiX1qKlZhH//+99YuHAh2rVrh+HDh2PQoEHo0qULDjqoPW66idi4cR/at9+M7ds/A/kugPcxZMhQ\nXHDBBPz4xz9Gx469MXo0sGYNIGLso7AQ2L/f2FdZWbSOsqmxVBERkJSUNpLqr4D1Bx+ReFg6FTKF\nl15wt0jcah5aD0HwgrMZbo1oc57nDRs28F//+henTZvGiy++mMcddwpFTiFwGkWu4CWX/Ia9ej3F\ngoLKmKVfVkmXfUCFvebbzZ5M5R5VY00/t85H69bRIyhT0Zh1vZ3nuK6ujh988AHnzZvHKVOm8Lzz\nzuOxx/7A1NgZFLmGHTr8gfn5L7OsbFfMfhJLY3Z78/KiUyZhyW8nA3I1nZLLJ90vfpxJdbUxRW2v\nXo2Cbd3a28ROzjSIvbffawrGT87d7Xv2Znq2J5dSjbmvt3ixMXuflVbJz/c2x0osjdlnLfVqZzLa\ndAY49rRJWPLbyZCzTjyXT3oinBGRX2cSXQGSeCY6txs4kfN13pCJhuYnOkbrewsWxO6UDSL36Kc6\nRTXmjtvcKom+61djbvoIQmP33BO/VDXb+e1kyFknTubuSY9HLLH7cSbOKMdZ/5vqj4S9E6ikxNsU\nrnbiTS0ca3RcUKjGUteYWzTs1gJLVmN2+/r0MTThJ61l7zB1fsc+oVXr1pl5ZF66yWkn3hyJNyza\nb6WKfcCH/XWqPxLGRPmNf/EifbdoyG0uCrfP7NO8BoVqLBiNOdcNUmPOfHm3bolLJZ1pGvv37d9x\natdLKWHYUSceMtLdhA/iBl6wIPpG6NXL3V63mznWrHBO25yfBYVqLPwaq6xsbI1ZnZCxOtedGrOP\nQnbrMFUn7v6Xl1JpixJFUZFRqvXaa01LtmpqjDKvmhpv23Jbv6zMKDcrKAB69gR69Gj8jB6rO9u0\niV6+555Ge4HGfa5a1VhCtmYN8Nhjxj5LS41yrz59gFdfbTxGyza3z5TgCEpjsdYtKwMGDADy84H+\n/Y1rauFFYxs3RpfJ9u5taMGLxtauNfZdWGjY8eKL0cd47LHG+wUFxv9jjklsT4sg1V8B6w8aJcUk\nXvWAM2WRaH37A4bd5jPxsm23gRHOfdpzpfZ5qeNVFKQ7B60ai00szfjVl5s+/OjX3lJwe65uIo0l\nqnZpbv0c0HRKyHC7Y+hvCtF467t9Zq/Ttg+hduscssxzc8SPPkqKNN2O27zUAZ8ez7SkuVNi4kNj\nQejL+XAFLxpbtKipI66uJn//+/RqLFV9ZQN14mEiTrgSK4+ZyLnHyyPGqi6wR+b2ziFrwiy3aGrd\nusabCyD79YuOrOxPUkn2BokXzSke8amxIPTl1rHpV2PV1UZtuj2fPXBgdNWJvXWZjDZyVV/qxMNE\ngjost2Zgoptp+XKyurJpeOGsJrBHPrFG09k7jQoKjIm2rE3edFP0DXb99dF2JDMvhc/To3jBp8aC\n0pez9eZXY/Y5va2/e+6JttEqSUy2bDBX9aVOPEwkWTYQ72ZKFF64fWw3w55ftJyx26g3eyQu0vhA\nWTJ+NOen6Zrk6VHsJHESoxx7rAS2D305zfCiscrK6EjcXuPd0vWlTjxseOl18dPTlCC88FsOZp84\n31nLvW6dMXrP7sDtpnl5aEMQp0dJQJAaS1Jf8cyIpbHqamNKCWu2S6dZLVVf6sRzDb/OOkF4ES/V\nEes+TvSAhlhme3logxIC/GjMg75iVZpYnwehsZasr4w4cQDdALwCYDWAjwD8PMZ6mTjm3CbZnqYY\n4YV1wzg7HRvv4wiH9NltpGps31m82N+kRbGCukw3XceOHZuZHeUyfjWWIHyNNc1COjWWq6mRZMiU\nE+8M4CjzdVsA6wAMcFkvE8ec2yTprGMR63413o8Y72Mfl/e5MKk7wTWoM++46srqjDddVWMeyHWN\nVWZPX9kgK+kUAE+hhT4oORACTNzFC66G9NnNQuzjEKxkdUG7pNqkTW7gxbu9jyhJA6oxj+SsxiLG\nj4HXkWvNgIw7cQC9AHwOoK3LZ+k/YqUJse7X6spqLu9zoXFz+W2T2qLtqBt40ZveR5SkAdVYdsiY\nxvrsZnX+YbFTQblWBO6BIJy458eziUhbABUAbiO5wOVzTp8+vWG5vLwc5eXlnratJIGXZ495fL5X\n1KZge/ZWaSlqnluK1ZuKjE1Yn61ZA5SUGBNbrFrV+NyvwkJjkowRI9JyyOajrNKybcWFTGusRw2K\nxjv0VVQU/Wy5NGss3VRUVKCioqJheebMmWAmHs8GoADACwCui7NOen+ylEYCjEyabMot2nZ+weuI\nkoBRjWWQbGnMLexvxj2dyOAshv8HYA3Ju1P6xVCCwTn92+rV/r5vm8KuyabENh1hSUn0NHaAERmN\nGNEYdcWbVi9gxo4dm7ZtKw6ypTGnvqz3MqSxnCSRlwcwEkA9gPcBrASwAsA4l/Uy8cOlkKlFJo6w\nqEne26pAaQmlAUpsVGMZAZnMiSdCRBjUthQPeMxFNsElv1hTOiKpTSnNHNVY2jH7eVLKiasTb2nU\nuHRO6l2lBIlqzDPqxJXkSDbCUhSvqMY8oU5cURQlhwnCieszNpWcQccdKEpTNBJXcgYd7KM0NzQS\nVxRFaeGoE1cURclh1IkriqLkMOrEFUVRchh14krOoHOnKEpTtDpFURQlS2h1iqIoSgtHnbiiKEoO\no05cURQlh1EnriiKksOoE1dyBp07RVGakrA6RUTmAfghgK0kB8dZT6tTlLSic6cozY1MVac8COAH\nqexEURRFSQ8JnTjJ1wHszIAtiqIoik80J64oipLDFAS5sRkzZjS8Li8v144oRVEUGxUVFaioqAh0\nm56G3YtITwBPa8emkk3Ky8sDvwEUJZtk7BmbItILhhMfFGcddeKKoig+yEh1iog8AmAZgH4isklE\nLk1lh4qiKEpw6CyGiqIoWUJnMVQURWnhqBNXFEXJYdSJKzmDlqwqSlM0J67kDDp3itLc0Jy4oihK\nC0eduKIoSg6jTlxRFCWHUSeuKIqSw6gTV3KGsWPHZtsERQkdWp2iKIqSJbQ6RVEUpYWjTlxRFCWH\nCdSJazpFCRqSiEQi2L9/P3bt2pVtc5RmSCQSyWnflZZIvK6uDtdffz0ikUg6Nq80Y7Zt24avv/66\nYXnSpEl48skn8e677+LSSxtnQV61ahX+8pe/ZMNEJcf55JNPUFtb27Dcs2dPbN26NYsWpUagTlzE\nyM/X1taipKQEeXnG5nfs2IFHH300yF0pzQh7FHTfffdh8eLFDcsPPvggJk6ciJEjR2LnzsbndRcV\nFaFr164NyytWrMDy5cszY7CSc9g1duONN2L9+vUNyxs2bEDnzp2zYVYgeH2yzzgAs2E4/Xkkf++y\nTszqlPXr12P+/Pm44YYbAAC7du1Cq1at0KZNm1RsV5oBCxcuxNNPP425c+cmXDfe3CkvvfQSvv32\nW5xxxhkAgK+++godO3ZEfn5+oPYqucfUqVPRv3//qJZcWAiiOgUk4/7BcNyfAugJoBDA+wAGuKxH\nrzz44IOcOnWq5/WDYsmSJRnfZzI0Zzs3b97Mm2++uWG5pqaGNTU1nr7rR2OXX345n3nmGZLN+3xm\ng7DbuWTJEj7wwAMNdn7zzTesra3NrlExMDWd0A/H+/OSTjkGwCckN5KsBfAYgNNT+eG45JJLcPvt\ntzcsX3fddXj++edT2aQncuUhu83Jzrq6OsyaNashgi4uLsZRRx3V8Hnbtm3Rtm3bwG174IEHMH78\neADAkiVLMGbMmNDnPZvTdc8kmzdvxj/+8Y+G5a5du2LAgAENdhYXF6OgoCBL1qUfL068K4AvbMub\nzfdSwsqfA8Att9yC4447rmF5woQJeO+99xqW33zzTdTU1KS6SyVNfPXVV1EdRdddd11DJUlBQQFq\namqwd+9eAEDr1q1x7rnnZsQuu8bmzp2LTp06AQD279+Pvn37oq6uDoBRnfDyyy9nxCbFP5FIBFu2\nbGlY3rp1K2666aaGZRHB9u3bG5b79u2LkSNHZtTGbBKKOvHOnTvjsMMOa1ieN28eSktLo5a/+eab\nhuXhw4fj448/blg+66yzojoqLr/8cmzcuLFh+corr4xanjx5MjZt2tSwfPXVV0ct/+QnP4n7+eTJ\nk6O2d9VVV0UtX3HFFVHLTnsuu+wyfP755w3Ll156adTy/Pnzo5YvuuiiqOULL7wwavmCCy7AZ599\n1rB8/vnnRy2fd955UcvnnHNO1PLEiROxYcOGhuWzzz47avmss86KWj7++OPxn//8J+r79vMzZsyY\nhk5tAJgxY0ZW+z9EBP37929w6q1atcKSJUsaorOamhr86U9/alh/27Zt6NWrV8Pyzp07ccIJJzQs\n79q1C2eddVbU8sSJE6OW7T9UVVVVOP/886M+v+CCC6KWnZ/bv79r1y6cc845Uctnn312XHt+9KMf\nRS2feeaZUcdj9R1Yy6effnrU8mmnnRa1fOqppzYs7927t6GVY30+bty4qOUf/OAHDcs7duzAySef\nHLV80kknRS3bz++2bdswdOjQqGX7/g855BAMGzasYblr16647rrr0FJJ2LEpIiMAzCA5zlyeCiOP\n83vHerlbaKkoipIlmGLHphcnng9gHYATAHwJ4G0A55Fcm8qOFUVRlNRJmO0nWS8i1wB4CY0lhurA\nFUVRQkBgsxgqiqIomSfljk0RGSciH4vIf0Tk5iCMSsGWeSKyVUQ+tL3XTkReEpF1IvKiiBxq++wW\nEflERNaKyMnuWw3cxm4i8oqIrBaRj0Tk5yG1s7WIvCUiK01bbw+jnbZ954nIChFZGFY7ReRzEfnA\nPKdvh9jOQ0XkCXO/q0Xk2LDZKSL9zPO4wvxfJSI/D6Gdt5jn8EMR+aeItArcxlSKzOFxIFCm/gCM\nAnAUgA9t7/0ewE3m65sB3GG+LgGwEkZKqZd5HJIBGzsDOMp83RZGf8OAsNlp7ruN+T8fwJsARobR\nTnP/1wP4B4CFYbzu5r43AGjneC+Mdv4NwKXm6wIAh4bRTpu9eQC2AOgeJjth+MUNAFqZy48DuDho\nG1M1cgSA523LUwHcnMkLGOPE2Z34xwAON193BvCxm60AngdwbBbsfQrAiWG2E0AbGB3aJWG0E0A3\nAIsAlKPRiYfRzs8AdHC8Fyo7ARwCYL3L+6Gy02HbyQCWhs1OAO1Me9qZjnlhOu71VNMpaRkIFDCd\nSG4FAJJfAehkvu+0vRIZtl1EesFoObwJ46KGyk4zRbESwFcAKkiuCaOdAO4CcCMAewdPGO0kgEUi\n8o6IXBFSO48EsE1EHjRTFXNEpE0I7bRzDoBHzNehsZPkTgCzAGwy91dFcnHQNoZisE+GCUVProi0\nBfAkgOtI7kZTu7JuJ8kIyaEwIt3RIlKOkNkpIqcC2EryfQDx6m2zfj4BjCQ5DMB4AD8TkdEI2fmE\nETEOA/Bn09Y9MCLEsNkJABCRQgATADxhvhUaO0WkN4w0X08AXQB8R0QucLEpJRtTdeKVAHrYlruZ\n74WJrSJyOACISGcA1mTVlTByaBYZs11ECmA48IdJLgirnRYkqwE8B+B7IbRzJIAJIrIBwKMAjheR\nhwF8FTI7QfJL8/83MNJoxyB853MzgC9Ivmsu/xuGUw+bnRanAHiP5DZzOUx2fg/AGyR3kKwHMB/A\n94O2MVUn/g6APiLSU0RaATgXRt4nmwiiI7KFAC4xX18MYIHt/XPN3uIjAfSBkffNBP8HYA3Ju8Nq\np4gUW73mInIwgJNgdLqEyk6S00j2INkbhv5eITkJwNNhslNE2pitL4jId2DkcT9C+M7nVgBfiEg/\n860TAKwOm502zoPx420RJjvXARghIgeJiMA4l2sCtzGA5P0409hPAEzNZIeGiy2PwOil3g8jD3Up\njE6FxaaNLwE4zLb+LTB6gNcCODlDNo4EUA+jkmclgBXmOWwfMjsHmbatBPABgBvM90Nlp8PmsWjs\n2AyVnTByzdY1/8i6V8Jmp7nfITACtPcB/D8Y1SlhtLMNgG8AFNneC5WdMPpqVgP4EMBDMKr4ArVR\nB/soiqLkMC2xY1NRFKXZoE5cURQlh1EnriiKksOoE1cURclh1IkriqLkMOrEFUVRchh14oqiKDmM\nOnFFUZQc5v8DxDNmvBTWwd4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe056564828>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Si at Hope and Gravesend\n",
    "dataG=session.query(Profs.Year, Profs.Month, Profs.Day, Profs.value).\\\n",
    "            filter(and_(\n",
    "                Profs.value<2.5, # remove single high outlier\n",
    "                or_(Profs.variable_name=='SILICA DISSOLVED',\n",
    "                      Profs.variable_name=='SILICON DISSOLVED'),\n",
    "                Profs.station_name.like('%Gravesend%'))).all()\n",
    "YDG, valG = plotYMDV(dataG,'r','temp')\n",
    "\n",
    "dataH=session.query(Profs.Year, Profs.Month, Profs.Day, Profs.value).\\\n",
    "            filter(and_(\n",
    "                or_(Profs.variable_name=='SILICA DISSOLVED',\n",
    "                      Profs.variable_name=='SILICON DISSOLVED'),\n",
    "                Profs.station_name.like('%Hope%'))).all()\n",
    "YDH, valH = plotYMDV(dataH,'b','Si, red=Gravesend, blue=Hope')\n",
    "\n",
    "units=session.query(Profs.unit_code).\\\n",
    "            filter(and_(\n",
    "                or_(Profs.variable_name=='SILICA DISSOLVED',\n",
    "                      Profs.variable_name=='SILICON DISSOLVED'),\n",
    "                or_(Profs.station_name.like('%Gravesend%'),\n",
    "                    Profs.station_name.like('%Hope%')))).group_by(Profs.unit_code).all()\n",
    "for row in units:\n",
    "    print('units:',row)\n",
    "    \n",
    "YD=np.concatenate((YDG,YDH))\n",
    "val=np.concatenate((valG, valH))\n",
    "sidict=gsmooth(YD,val,60)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe0559c5ac8>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVcAAAEACAYAAAAHujVXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHT5JREFUeJzt3X20XFWZ5/HvLwlJB5qXIAg2ad4E2+BSeQ2MUSlU0rwM\nYVp7EGfQRkZkJjZmtTROgjq56RcDjjTdM6OL1SIO2JF0VF7SYzckEQoXNCQRAk2TaEclSMeQEUzD\n0EiSmzzzx9mB8t46t06l6tStl99nrVqp2rXP2WffuvfJrufsc7YiAjMza68J430AZmb9yMHVzKwE\nDq5mZiVwcDUzK4GDq5lZCRxczcxKUCi4Stok6XFJ6yStGfHeVZJ2Szq4pmyBpI2SNkia3e6DNjPr\ndpMK1tsNVCJiW22hpOnA2cDTNWUzgIuAGcB0YJWk48MTas2sx0maB3wsvfxKRPyPvLpF0wLKqXsD\ncPWIsguBpRExHBGbgI3AzILtmJl1JUlvAf4TcCpwIvBvJR2bV79ocA1gpaS1ki5PDc0BnomIJ0bU\nPQJ4pub15lRmZtbLZgCrI2J7ROwCvge8P69y0bTArIjYIulQYIWkHwDXkKUEzMwGwT8CfyJpGrAd\nOA9Ym1e5UHCNiC3p359LuhM4EzgaeFySyHKrj0qaSTZSPbJm8+mp7FdIcg7WzAqLCLWy/UFSvFC8\n+taIOHxE+z+QdB2wEngJWAfsytuBGp1nkrQvMCEiXpK0H7ACWBQRK2rqPAWcHBHbJJ0ALAFOJ0sH\nrARGndCS1LfnuIaGhhgaGhrvwyhFv/atX/sF/dE3SS0HV0nxJwXrfpbGwVzSn5KlRm+s936Rketh\nwB1ppDkJWFIbWJMgO+lFRKyXtAxYD+wE5vZtFDWznrJPi9tLOjR9gz8S+B3gjLy6DYNrRDxFdmZs\nrDrHjni9GFhc7HDNzDqj6EmmMXw7zenfM3B8scS2bKRKpTLeh1Cafu1bv/YL+rtvzZra4vYR8e6i\ndRvmXMvSzzlXA2lR7nsRCzt4JONL+oem6ke8raQj6W3tyrl+pWDdy2n9BJpHrmY2MDoZ8Bxcra6x\nRp7dvO/x1I4Reb2fTdH9tuPn2u/fKlo9odUMB9cB16+BbjyU9bP0Z9Q+HrmaWen6fZRaj0euVgqP\ngAbTIAbRPA6uZtaQg2bzWp2K1QwH1z7kEWp/cRBtH+dczcxK4LSAGZ2dgtRpHo2OD49cbeC0Emya\n2bbMQOyA2f08crXCenHUNp4cAAebR65mZiXwyNXMrASeijXg/FXfrBxdN3KVtAl4AdgN7IyImZK+\nAFxAtlDXj4GP7rlxrKQFwGXAMDCvzsoFab97f5OKfuAgatZZrY4mU2y7hGztrCfI4t6OunWL3FNV\n0k+AUyJiW03Z+4B7I2K3pGuBiIgFNWtonUa2OOEqctbQgqG96V9D3RigHUhf042fj3W3dt3P9fmC\n0fV1w6Pv5yrpKOA+4M0RsUPSXwPfiYhb6+1jQtHjGlk3IlZFxO708mGyQAowB1gaEcMRsQnYCMws\n2I6ZWWkmTSr2yPEisAPYT9IkYF/gZ7ltFTymAFZK2gX8ZcSoG3pfBtyWnh8BPFTz3uZU1jHNjBLb\nPYryCNWse+0zce+3TatbXw/8FHgZWBERq/LqFw2usyJii6RDyYLshoh4AEDSZ8jysLeNvYt67qt5\nfjRwTPO7aJGDYXv0+lf9gw8+mG3btjWuaKWbNm0at99+O9Vqte37HmNU2pCkY4E/AI4iOwf1LUn/\nISK+UbetIjuNiC3p359LuoPsa/4Dki4FzgPeU1N9M/CbNa+np7I6zirSvI2TXg+Yzdi2bRte0607\nSKJSqfzKwoqLFrVnELTPlPrl1R3Zo4FTgQcj4hcAkm4H3gHUDa4NT2hJ2heYEBEvSdoPWAEsIsvB\nXg+8OyKer6m/54TW6WTpgJV0+ISWNWeQgmiedMJkvA/DqP9ZtOuEVvxGwbo/q3tC6+3AX5GdrN8O\nfA1YGxFfqrePIiPXw4A7smDIJGBJRKyQtBGYTJYmAHg4IuZGxHpJy4D1vLa2t39ru4QDqQ20FtIC\nEfG4pFuBR8imYq0D/jKv/rgure2Ra+c5uNbnkWv3KHXkemzBuj/x0tpmZsW1MFugWQ6ufcCjUbOC\nOhjxHFx7iIOoWYtyZguUwcG1SzmQmpXAI1czsxI4uA4Oj1DNOsgntPqTA6nZOPPItbc5iJp1KQfX\n7uSgadbjHFzNzErgqVid49Go2QDxyLU1DphmVpdnCxTjIGpmTfHIdTQHUjNr2SAHVwdRMyuN0wJm\nZiXoYMQrtLS2pE2SHpe0TtKaVDZN0gpJP5R0j6QDa+ovkLRR0gZJs8s6eDOzpvxawUcdkt6UYuCj\n6d8XJH0yr6lCwRXYDVQi4qSImJnK5gOrIuK3gHuBBekATgAuAmYA5wJfVloHxsxsXE0s+KgjIv4p\nxcCTgVOAfwXuyGuq6CBZjA7EFwJnpue3AFWygDsHWBoRw8CmtNbWTGB1wbbMrMstGmO8tLCbl8tp\nX1rgfcCPI+KZvApFR65BthDhWkkfS2WHRcRWgIh4Fnh9Kj8CqG1wcyozMxtfkwo+GvsgcFujpoqY\nFRFbJB0KrJD0Q7KAW2sv/ru6r+b50cAxze/CzPpOtVqlWq22f8c5X/mrT2ePIiTtQ/YNff6Y9Zpd\n8VLSQuAl4GNkeditkg4H7ouIGZLmAxER16X6dwMLI2L1iP3UXf3VU7FsPHj11+ZIi3Lfa/VvuNTV\nX/9bwbp/lL/6q6Q5wNyIOGesfTQcuUraF5gQES9J2g+YDSwClgOXAtcBvwfclTZZDiyRdANZOuA4\nYE29fTuQmllHtSfn+iEapASKNnUYcEc20mQSsCQiVkj6PrBM0mXA02QzBIiI9ZKWAeuBnWQR3kMC\nMxt/Ld4VKw023wd8vFHdhsE1Ip4CTqxT/ovUSL1tFgOLGx6pmVkntThyjYiXgUM70JSZWQ8Z5HsL\nmJmVxvcWMDMrgUeuZmYlcHA1MyuB0wJmZiXIueNVGRxczWxwOC1gZlYCpwXMzErgkauZWQkcXM3M\nSuC0gJlZCTxbwMysBB65mpmVwDlXM7MSOLiamZWggxGv6OqvSJogaZ2k5en1TElrUtkaSafW1F0g\naaOkDZJml3HgZmZNm1jwkUPSgZK+mWLbk5JOz6vbTByfBzwJHJBeXwd8Ni35ci7w34GzJJ1AtuTL\nDGA6sErS8V7qxczGXesj178A/jYi/r2kScC+eRULjVwlTQfOA26qKd4CHJieHwRsTs/nAEsjYjgi\nNgEbgZlNHb6ZWRmmFHzUIekA4F0R8TWAFONezGuqaBy/Abia14IpZGt2PyjpekDAO1L5EcBDNfU2\npzIzs/HV2sj1GOA5SV8D3g58H5gXEb/cq6YknQ9sjYjHJFVq3voqcGVE3Cnpd4GbgbObOdKhoaFX\nn1cqFSqVSm5dMxsc1WqVarXa/h3nRLzq32ePAlufDHwiIr4v6c/JBpkL61VWo1SopM8DlwDDwFRg\nf+AOYE5EHFhT718i4iBJ84GIiOtS+d3AwohYPWK/TsNa15CEfx+LkxblvhdRN9Y0se/Rn0UqU4v7\njfi/Beu+nlHtSToMeCgijk2v3wn814i4oN4+GuZcI+KaiDgy7fBi4N6I+DDwI0lnpkbeS5ZbBVgO\nXCxpsqRjgOOANcW6ZGZWnphY7FF324itwDOS3pSK3gusz2urlQzEFcCXJE0GXgE+ng5gvaRlqdGd\nwFwPUc2sG+xqfbbAJ4ElkvYBfgJ8NK9iw7RAWZwWsG7itEBzejUt8Mq/Fqv7a/uNTgs0y1domdnA\n2D5lcsGaO1puy8HVzAbGromduy2Wg6uZDYxdHbznoIOrmQ2MYQdXM7P229XBkOfgamYDw2kBM7MS\nOLiamZVgO0WnYrXOwdXMBoZzrmZmJXBawMysBA6uZmYl8DxXM7MSOOdqZlYCpwXMzEqwo4NTsQqt\n/gogaYKkRyUtrym7Mq3f/YSka2vKF0jamN6b3e6DNjPbG8NMLPRoh2ZGrvPIVhc4AEDSWcAFwFsj\nYljSIal8BnARMAOYDqySdLzvjG1m463VnKukTcALwG5gZ0TMzKtbaOQqaTpwHnBTTfF/Bq6NiGGA\niHgulV8ILE1rem8iW1sr9wDMzDplFxMLPcawG6hExEljBVYonha4AbgaqB19vgl4t6SHJd0n6ZRU\nfgTwTE29zanMzGxctSG4ioJxs2ElSecDWyPisbTjPSYB0yLiDODTwDeLNGhmNl7akHMNYKWktZIu\nH6tikQTELGCOpPOAqcD+km4lG53eDhARayXtkvQ6spHqkTXbT09lowwNDb36vFKpUKlUChyOmfW7\narVKtVpt+353MKVu+frqz1lffa7ueyPMiogtkg4lC7IbIuKBehWbWv1V0pnAVRExR9IVwG9ExMK0\njvfKiDhK0gnAEuB0snTASmDUCS2v/mrdxKu/NqdXV3/9q/hAobqX6NsN25O0EPh/EfFn9d5v5dTZ\nzcDNkp4AtgMfAYiI9ZKWkc0s2AnMdRQ1s27QyjQrSfsCEyLiJUn7AbOB3P9lmgquEXE/cH96vhP4\ncE69xcDiZvZtZla2FqdiHQbcISnIYueSiFiRV9lXaJnZwGjl8teIeAo4sWh9B1czGxi+t4CZWQkc\nXM3MSrA9ZypWGRxczWxgeORqZlYCB1czsxJ4mRczsxJ4mRczsxI4LWBmVgIHVzOzEmzv4BpaDq5m\nNjCcczUzK4HTAmZmJXBwNTMrgee5mpmVoJM516KrvyJpgqRHJS0fUX6VpN2SDq4pWyBpo6QNkma3\n84DNzPZWG1Z/zY2FIzUTxueRLd1yQE0j04GzgadrymYAFwEzyBYnXCVp1BpaZmadtqM9U7FGxcJ6\nCo1cUxA9D7hpxFs3AFePKLsQWBoRwxGxCdgIzCzSjplZmVpdWnuMWDhK0ZHrniB6YE0jFwLPRMQT\n0q8skngE8FDN682pzMxsXLUh5zoqFuZpOHKVdD6wNSIeA5TKpgILgNbW0DUz66BWcq51YuGYS28X\nCeOzgDmSzgOmAvsDtwJHA48rG7ZOBx6VNJNspHpkzfbTU9koQ0NDrz6vVCpUKpUCh2Nm/a5arVKt\nVtu+37zA+Vz1SZ6vPtlo81GxUNKtEfGRepXVzHkmSWcCV0XEnBHlTwEnR8Q2SScAS4DTydIBK4FR\nJ7Qk+RyXdQ1J+PexOGlR7nsRrX2hrfdZpLIxR4oF9hvnxrcL1f07fWDM9vJiYa12TfoK0hA5ItZL\nWkZ2Nm0nMNdR1My6QdfeWyAi7gfur1N+7IjXi4HFrR2amVl7tWkqVm4srOUrtMxsYPjyVzOzEnRt\nWsDMrJf5rlhmZiVwcDUzK4GDq5lZCbYzpWNtObia2cDwyNXMrAQOrmZmJfA81y6Udy11q9dRd9JY\n14MX1Uv9NRvJ81zbrB1BpZl993MA6of/ZGxwOS3QgjIDaa8Yj5/BoP0nY73JwbWAbg6inRzd+edg\nVtz2He25cUsRPRtc26XoH3o3BzEzK2bXsHOuv8KBrb84hWDjZdfwgKQFygqaZfyh1ttns8dfNKj4\nPxOzcnRlcJU0AXiEbMXXOZK+AFwAbAd+DHw0Il5MdRcAlwHDwLyIWNH2I2f8RztlBdxB5NGsdcLw\nzr0PrpKmAN8DJqfHXRFxTV79Zkau84AngQPS6xXA/IjYLelastVgF6Q1tC4CZpAtTrhK0qg1tJrV\njX9o/RIYPXq2QbF7195/WY+I7ZLOioiXJU0EHpQ0KyIerFe/UEuSpgPnAX8KfCo1tKqmysPAB9Lz\nOcDSiBgGNknaCMwEVhfrQPcF0V7U6s/RAdf6UotpgYh4OT2dAkwAtuXVLRrGbwCuBg7Mef8y4Lb0\n/AjgoZr3NqeynuAA0j08lcva7pXWTjPVpEffCNwYEevz6jZsSdL5wNaIeExShbTKa837nwF2RsRt\n9bYfy8KFr2UKKpUKlUql2V0U0i8Bc7yDSl77nf75Oj/b/6rVKtVqtf07Hm5t84jYDZwk6QBghaQz\n02KFo6hRKlTS54FL0mFNBfYHbo+Ij0i6FLgceE9EbE/152fHENel13cDCyNi9Yj9Fk7D9ktwbFU/\nBJDu/SyH8ArwxY31Obb6eypp1GeRypSzSdH9Bo/nfMZrq/D96muvb1zUsD1JnwNejojr677fzC+U\npDOBq9JsgXOA64F3R8TzNXVOAJYAp5OlA1YCo05o1Quu3fKH18l8Yz8EzFZ1x+c+lB57b5A+y54N\nro8UjHenjG5P0iFk39JfkDQVuAdYFBHfrbeLVhIQ/5NsOsJKSQAPR8TciFgvaRmwHtgJzG11pkCZ\nBukPolv55Jl1zM6Wtn4DcIuygDcB+HpeYIUmg2vKLdyfnh8/Rr3FwOJG+/MfkPWTZn6f/Z/6ONm1\n95tGxBPAyUXr98Tlr9C5X8YyA77/oIpr5mfVi/9J9+Ix94UWT2g1o2eCay9xEO2sbpnFYD3glc41\nNa7BdbyDkP/4+lu35HJ957Uu0sGRa1OzBdracBNTsco7Bi97Ypl6Z6htfJQ6W+Cugp/xha23NxBp\nAY8IzAxwzrUVDqRmlqu1qVhN6Zn7uY78+u0gamZNa2EqVrN6ZuTaDcHU+VWzHue0QOc4YJoNkEGZ\nimVm1lEeubafR6hm5uBakAOmmTVlUIKrg6OZddSgTMUyM+soT8UyMyuBZwuYmZWggznXCUUrSpog\n6VFJy9PraZJWSPqhpHskHVhTd4GkjZI2SJpdxoGbmTVtZ8FHHZKmS7pX0pOSnpD0ybGaKhxcgXlk\nS7fsMR9YFRG/BdwLLEgHcAJwETADOBf4cloWwcxsfO0q+KhvGPhURLwF+DfAJyS9Oa9yoeAqaTpw\nHnBTTfGFwC3p+S3Av0vP5wBLI2I4IjYBG4GZRdoxMyvVcMFHHRHxbEQ8lp6/BGwgW4S1rqIj1xuA\nq4HamyEeFhFb9zQKvD6VHwE8U1Nv81gHYGbWMS0E11qSjgZOBFbn1Wl4QkvS+cDWiHhMUmWMqk3f\naXhoaOjV55VKhUplrN2blcvZq+4wbdo0qtUq1Wq1/TvPm+f6sypsKdaepF8HvgXMSyPY+vUa3X1d\n0ueBS8ji+VRgf+AO4FSgEhFbJR0O3BcRMyTNByIirkvb3w0sjIjVI/Y77isRmFlvaNtKBB8uGHO+\nXr89SZOA/wP8XUT8xVi7aJgWiIhrIuLIiDgWuBi4NyI+DPwNcGmq9nvAXen5cuBiSZMlHQMcB6wp\n1iMzsxK1nha4GVjfKLBCa/NcrwWWSboMeJpshgARsV7SMrKZBTuBuR6imllXaOHyV0mzgP8IPCFp\nHVkq9JqIuLtu/UFeoNDMekPb0gK/UzDm3NHjCxTWW13AN3Mxs9IMyl2x6mllXS0zszENcnBthgOx\nmTXFtxw0MyvB9s415eBqZoPDaQEzsxI4LdAa51fNrC6vRGBmVoJBTgt41GlmpRnk4GpmVhrnXM3M\nStDBqVi+t4CZdb223VvgDQVjzpYev7eAmVlHOS1gZlYCT8UyMyuBZwuYmZWgg8G14TIvkqZIWi1p\nnaQn05paSJopaU0qXyPp1JptFkjaKGmDpNlldsDMrLCdBR91SPqqpK2S/qFIU4VmC0jaNyJeljQR\neBD4Q+CPgcURsULSucCnI+IsSScAS4DTgOnAKuD4kVMDPFvAzIpq22yBwotUj25P0juBl4BbI+Jt\njfbQcOQKEBEvp6dT0ja/ALYAB6Xyg4DN6fkcYGlEDEfEJmAjMLNIO2Zm3SoiHgC2Fa1fKOcqaQLw\nCPBG4Ma0COF84EFJXwQEvCNVPwJ4qGbzzanMzGxgFAquEbEbOEnSAcA9kirAZ4ArI+JOSb9LtuTs\n2c00nu1mj6OBY5rZ3PchMOtT1WqVarXayRbTo32avkJL0ueAXwKfi4gDa8r/JSIOSiPaiIjrUvnd\nwMKIWD1iPwFDrR5/XQ66Zv2lfTnXHQVrT67bnqSjgL8pknNtOHKVdAiwMyJekDSVbHS6CPiRpDMj\n4n5J7yXLrQIsB5ZIuoEsHXAcsKZgj9rCq8qaWX0tz8VSejRUJC3wBuAWSSI7mfX1iPiupCuAL0ma\nDLwCfBwg5WOXAevJJjXM9bQAM+sOe3/9q6RvABXgdZJ+SvaN/Gu59cfzxi1lpQWa4RGtWfdrX1rg\n2YK1D/eNW8zMiuvcnVsGPrg6P2s2SDp3/evAB9d66gVccNA1630euXYlj3LNep1HrmZmJfDI1cys\nBL/sWEsOrmY2QJwW6ErOr5r1OqcFzMxK4JGrmVkJBnjkWvSrd95cVDOzfB65diXPczXrdQMycnVg\nMrPO8lSsruT/DMx63YCMXM3MOquLcq6SpgDfAyanx10RcU1670pgLtkRfyci5qfyBcBlqXxeRKxo\n94HnjSJ9osvM8rU2cpV0DvDnZAsHfHXPclb1NFxaOyK2A2dFxEnA24D3SJqVFim8AHhrRLwV+GJq\nfAZwETADOBf4clrFYGB0dmG1zurXvvVrv6C/+9a84YKP0dIq2P8L+G3gLcCHJL05r6Wiq7++nJ5O\nIQvI24CFwLURMZzqPJfqXAgsTeWbJG0EZgKr6YB6I9qxRrNl5FGr1SqVSqXt++0G/dq3fu0X9Hff\nmtfSyHUmsDEingaQtJQs3v2gXuVCwTVF7EeANwI3pnWy3gS8W9LnyU7B/WFEPEK2KOFDNZtvTmXj\nxieizCzTUs71COCZmtf/TBZw6yo6ct0NnCTpAOCelBKYBEyLiDMknQZ8Ezh2b4/azKx8nZuK1fQC\nhZI+R3aE7wGui4j7U/lG4AzgcoCIuDaV3022SuLqEfvxirBmVlgbFijcBBxVsPrWiDh8xPZnAEMR\ncU56PT87rPontYrMFjgE2BkRL0iaCpwNLAJeJAuw96cUweSIeF7ScmCJpD8jG0YfB6wZud9Wf1Bm\nZs2IiKNb3MVa4DhJRwFbgIuBD+VVLpIWeANwSzrjPwH4ekR8V9L3gJslPQFsBz4CkPKxy4D1ZNnj\nuTFe63ebmbVJROyS9PvACl6birUhr37TaQEzM2us4TzXvSFpiqTVktZJejLNKEDSQkn/LOnR9Din\nZpsFkjZK2iBpdhnH1Q55fUvvXZmO/wlJ19aU93TfJC2t+cyekvRozTZd37cx+jVT0ppUvkbSqTXb\ndH2/YMy+vV3S30t6XNJdkn69Zpue6FvPi4hSHsC+6d+JwMPALLK5sZ+qU3cGsI4sTXE08CPSqLob\nHzl9q5B9XZiU3jukX/o24v0vAp/ttb7V6dc7gfuA2an8XOC+9PyEXunXGH1bA7wzlV8K/FEv9q2X\nH6WMXCH3wgOAeieyXr3wICI2AXsuPOhKOX37LzS4qKKH+1brIuAb6XnP9K1Ov35BdlLioFR+ENmc\nbIA59Ei/ILdvx0fEA6l8FfCB9Lyn+tbLSguukiZIWgc8C1QjYn166/clPSbpJkkHprKRk3PH/cKD\nseT0bc9FFQ9Luk/SKal6P/Rtz3vvAp6NiJ+kop7pW06/5gPXS/op8AVgQareM/2C3L49KWlOqnIR\nMD0976m+9bIyR667I7sfwXSyoHMm8GXg2Ig4kewX4fqy2i/TiL69a+RFFcCnyS6q6Dk5n9seHwJu\nG58ja03OZ/ZV4MqIOBL4A+DmcTzEvZbzmV0GfELSWmA/YMd4HuMgKi247hERLwLfAU6NiJ9HSvwA\nX+G1ryObgd+s2Ww6r31F61qpb38LnEo2Grg9la8Fdkl6HVk/jqzZrJf69h2yviFpIvB+4K9rqvXc\n5zbiM5sZEXem8m8Bp6VqPdcvGPW39k8R8dsRcRqwFPhxqtaTfetFZc0WOGTPV/6aCw8ek1R7xcP7\ngX9Mz5cDF0uaLOkYci486AY5fVsH3El2UQW1F1WQ9e2DPdy3x9LbZwMbIuJnNZv0xOc2xmf2oz0j\nc0nvJcs/Qo/0C8b8Wzs0lU0APgvcmDbpmb71urJulp134cGtkk4EdgObgCug5y486OeLKur2Lb33\nQUakBHqob3mf2RXAlyRNBl4BPg491S/I79snJX0CCOD2iPjf0HN962m+iMDMrASl51zNzAaRg6uZ\nWQkcXM3MSuDgamZWAgdXM7MSOLiamZXAwdXMrAQOrmZmJfj//R9wD9Z3fYwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fe058e3a780>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# find fraser location\n",
    "mesh=nc.Dataset('/ocean/eolson/MEOPAR/NEMO-forcing/grid/mesh_mask201702.nc')\n",
    "tmask=mesh.variables['tmask'][0,0,:,:]\n",
    "f=nc.Dataset('/results/forcing/rivers/R201702DFraCElse_y2016m01d23.nc') # example for dims\n",
    "plt.pcolormesh(np.ma.masked_where(tmask==0,f.variables['rorunoff'][0,:,:]))\n",
    "plt.colorbar()\n",
    "plt.xlim(350,399)\n",
    "plt.ylim(360,540)\n",
    "plt.plot((380,398,398,380,380),(400,400,520,520,400),'k-')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "masked_array(data = [--],\n",
       "             mask = [ True],\n",
       "       fill_value = 9.96921e+36)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f.variables['time_counter'][:]"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# create daily files, use data to initialize fields\n",
    "year=2016\n",
    "fname='/data/eolson/MEOPAR/NEMO-forcing-new/rivers-climatology/bio/R201809_2_bioClim_'\n",
    "for yearday in range(1,367):\n",
    "    t2=dt.date(year, 1, 1) + dt.timedelta(days = yearday - 1)\n",
    "    #datestr='y'+t2.strftime('%Y')+'m'+t2.strftime('%m')+'d'+t2.strftime('%d')\n",
    "    datestr='m'+t2.strftime('%m')+'d'+t2.strftime('%d')\n",
    "    new=nc.Dataset(fname+datestr+'.nc','w')\n",
    "    #Copy dimensions\n",
    "    for dname, the_dim in f.dimensions.items():\n",
    "        #print (dname, len(the_dim) if not the_dim.isunlimited() else None)\n",
    "        new.createDimension(dname, len(the_dim) if not the_dim.isunlimited() else None)\n",
    "    #new_tc=new.createVariable('time_counter',float,('time_counter'),zlib=True)\n",
    "    #new_tc[:]=f.variables['time_counter'][:]\n",
    "    \n",
    "    new_run=new.createVariable('no3',float,('time_counter', 'y', 'x'),zlib=True)\n",
    "    new_run[:,:,:]=const_no # other river value\n",
    "    new_run[:,400:520,380:]=no3dict[yearday]/mwN*1000.0 # convert to muM; Fraser value\n",
    "\n",
    "    new_run=new.createVariable('sil',float,('time_counter', 'y', 'x'),zlib=True)\n",
    "    new_run[:,:,:]=const_si # other river value\n",
    "    new_run[:,400:520,380:]=sidict[yearday]/mwSiO2*1000.0 # convert to muM; Fraser value\n",
    "    \n",
    "    new.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create constants file\n",
    "f=nc.Dataset('/results/forcing/rivers/R201702DFraCElse_y2016m01d23.nc') # example for dims\n",
    "year=2016\n",
    "fnameC='/data/eolson/MEOPAR/NEMO-forcing-new/rivers-climatology/bio/R201812_bioConst'\n",
    "new=nc.Dataset(fnameC+'.nc','w')\n",
    "#Copy dimensions\n",
    "for dname, the_dim in f.dimensions.items():\n",
    "    #print (dname, len(the_dim) if not the_dim.isunlimited() else None)\n",
    "    new.createDimension(dname, len(the_dim) if not the_dim.isunlimited() else None)\n",
    "\n",
    "#new_tc=new.createVariable('time_counter',float,('time_counter'),zlib=True)\n",
    "#new_tc[:]=f.variables['time_counter'][:]\n",
    "\n",
    "new_run=new.createVariable('no3',float,('time_counter', 'y', 'x'),zlib=True)\n",
    "new_run[:,:,:]=const_no # other rivers\n",
    "new_run[:,400:520,380:]=np.mean([no3dict[yearday]/mwN*1000.0 for yearday in range(1,367)])\n",
    "\n",
    "new_run=new.createVariable('nh4',float,('time_counter', 'y', 'x'),zlib=True)\n",
    "new_run[:,:,:]=const_nh # other river value\n",
    "new_run[:,400:520,380:]=nh4_cst/mwN*1000.0 # convert to muM; Fraser value\n",
    "\n",
    "new_run=new.createVariable('sil',float,('time_counter', 'y', 'x'),zlib=True)\n",
    "new_run[:]=const_si # other river value\n",
    "new_run[:,400:520,380:]=np.mean([sidict[yearday]/mwSiO2*1000.0 for yearday in range(1,367)])\n",
    "\n",
    "new_run=new.createVariable('dia',float,('time_counter', 'y', 'x'),zlib=True)\n",
    "new_run[:]=dia\n",
    "\n",
    "new_run=new.createVariable('phy',float,('time_counter', 'y', 'x'),zlib=True)\n",
    "new_run[:]=phy\n",
    "\n",
    "new_run=new.createVariable('mes',float,('time_counter', 'y', 'x'),zlib=True)\n",
    "new_run[:]=mes\n",
    "\n",
    "new_run=new.createVariable('zoo',float,('time_counter', 'y', 'x'),zlib=True)\n",
    "new_run[:]=zoo\n",
    "\n",
    "new_run=new.createVariable('don',float,('time_counter', 'y', 'x'),zlib=True)\n",
    "new_run[:]=don\n",
    "\n",
    "new_run=new.createVariable('pon',float,('time_counter', 'y', 'x'),zlib=True)\n",
    "new_run[:]=pon\n",
    "\n",
    "new_run=new.createVariable('bsi',float,('time_counter', 'y', 'x'),zlib=True)\n",
    "new_run[:]=bsi\n",
    "\n",
    "new_run=new.createVariable('tur',float,('time_counter', 'y', 'x'),zlib=True)\n",
    "new_run[:]=const_tur\n",
    "\n",
    "new.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "session.close()\n",
    "engine.dispose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "f.close()\n",
    "mesh.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 [('no3', 6.2599999999999998, 6.4369529952319278), ('nh4', 4.0899999999999999, 4.4271332178136902), ('sil', 55.689999999999998, 58.655887856134882), ('dia', 0.0, 0.0), ('phy', 0.0, 0.0), ('mes', 0.0, 0.0), ('zoo', 0.0, 0.0), ('don', 0.0, 0.0), ('pon', 0.0, 0.0), ('bsi', 0.0, 0.0), ('tur', 1.0, 1.0)]\n"
     ]
    }
   ],
   "source": [
    "new=nc.Dataset(fnameC+'.nc')\n",
    "print(yearday,[(ikey, np.min(new.variables[ikey]), np.max(new.variables[ikey])) for ikey in new.variables.keys()])\n",
    "new.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}