CR.QCTest.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# must insert this statement to render the plots within the notebook\n",
    "# this is specific to the ipython notebook\n",
    "%matplotlib inline\n",
    "\n",
    "## Import statements\n",
    "# netcdf/numpy/xray/stats\n",
    "import numpy as np\n",
    "import xray\n",
    "import pandas as pd\n",
    "from netCDF4 import Dataset\n",
    "from netCDF4 import num2date, date2num\n",
    "from datetime import datetime, timedelta\n",
    "import pytz\n",
    "import kray #### Custom functions for manipulating xray structures\n",
    "\n",
    "# import subplots function for plotting\n",
    "import seaborn as sns\n",
    "import matplotlib\n",
    "from matplotlib.pyplot import subplots\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "import kgraph\n",
    "\n",
    "# Solargeometry\n",
    "import solargeo\n",
    "\n",
    "# OS interaction\n",
    "import sys\n",
    "import os\n",
    "\n",
    "## Directory Lists\n",
    "# General directories\n",
    "dir_data = '/Users/karllapo/gdrive/SnowHydrology/proj/CloudClimatology/data'\n",
    "dir_grobs = '/GroundObs'\n",
    "dir_print = '/Users/karllapo/gdrive/SnowHydrology/proj/CloudClimatology/Graphics'\n",
    "dir_bmap = '/Users/karllapo/gdrive/SnowHydrology/proj/CloudClimatology/data/basemap'\n",
    "\n",
    "# Individual data directories\n",
    "dir_cv_cimis_ucipm = '/Users/karllapo/gdrive/GroundObs/CA_CentralValley.UCDavis/RAW/CIMIS/UCDavis_IPM.Daily'\n",
    "dir_cv_cimis = '/Users/karllapo/gdrive/GroundObs/CA_CentralValley.UCDavis/RAW/CIMIS/CIMIS.Hourly'\n",
    "dir_sio_cdwr = '/Users/karllapo/gdrive/SnowHydrology/proj/CloudClimatology/data/GroundObs/YOS.SIO.Obs'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processed Site: bee\n",
      "Processed Site: dan\n",
      "Processed Site: dep\n",
      "Processed Site: for\n",
      "Processed Site: fsl\n",
      "Processed Site: gin\n",
      "Processed Site: gld\n",
      "Processed Site: hod\n",
      "Processed Site: lee\n",
      "Processed Site: mer\n",
      "Processed Site: olm\n",
      "Processed Site: oth\n",
      "Processed Site: piu\n",
      "Processed Site: pri\n",
      "Processed Site: smi\n",
      "Processed Site: smj\n",
      "Processed Site: stb\n",
      "Processed Site: sun\n",
      "Processed Site: vvw\n"
     ]
    }
   ],
   "source": [
    "###################\n",
    "## READ SIO DATA ##\n",
    "###################\n",
    " # Read supporting station information\n",
    "os.chdir(dir_data+dir_grobs)\n",
    "stdat = pd.read_csv('All_StationSummary.v2.csv',sep= ',', \\\n",
    "                index_col=0,na_values=[-9999,'NaN']) # Read the supporting information\n",
    "stdat = stdat.groupby('Network').get_group('CDWR')\n",
    "\n",
    "os.chdir(dir_sio_cdwr)\n",
    "content = os.listdir(os.getcwd())\n",
    "num_files = len([name for name in os.listdir('.') if os.path.isfile(name)])\n",
    "\n",
    "# Empty lists and dictionary literals\n",
    "stations = []\n",
    "lat = []\n",
    "lon = []\n",
    "elev = []\n",
    "network = []\n",
    "grobs_yos = {}\n",
    "\n",
    "# time zone variables\n",
    "tz_pst = pytz.timezone('US/Pacific')\n",
    "\n",
    "na_value = ['   NaN']\n",
    "\n",
    "for files in content:\n",
    "    # Only read QC formatted files\n",
    "    if files[-10:] == 'Rad.QC.txt':\n",
    "        sitename = files.split('.')[0]\n",
    "        stations.append(sitename)\n",
    "\n",
    "        # Read SW data, asign to PST, and get SW that passes QC\n",
    "        grobs_yos[sitename] = pd.read_csv(files,sep= '\\t', parse_dates=True, index_col=0, na_values=na_value)\n",
    "        grobs_yos[sitename].index = grobs_yos[sitename].index.tz_localize(pytz.utc).tz_convert(tz_pst)\n",
    "        grobs_yos[sitename]['SWdwn_QC'] = grobs_yos[sitename]['SWdwn_Wm^-2'].where(grobs_yos[sitename]['QCFlag'] == 0)\n",
    "\n",
    "        # New data frame w/ daily means\n",
    "#         grobs_yos_daily = grobs_yos['SWdwn_Wm^-2'].resample('D', how='mean').to_frame(name='SWdwn_D')\n",
    "#         grobs_yos_daily['SWdwn_D_QC'] = grobs_yos['SWdwn_QC'].resample('D', how='mean')\n",
    "#         grobs_yos_daily['SWdwn_D_proc'] = grobs_yos['SWdwn_proc'].resample('D', how='mean')\n",
    "\n",
    "#         # List containing DataFrames with daily, processed only\n",
    "#         grobs[sitename] = pd.DataFrame(grobs_yos['SWdwn_D_proc'])\n",
    "#         grobs[sitename].columns = ['SWdwn']\n",
    "\n",
    "        # Fill in elevation/lat/lon/network\n",
    "        elev.append(stdat.loc[sitename]['elevation (m)'])\n",
    "        lat.append(stdat.loc[sitename]['lat'])\n",
    "        lon.append(stdat.loc[sitename]['lon'])\n",
    "        network.append('CDWR')\n",
    "        \n",
    "        print((\"Processed Site: \"+sitename))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xray.Dataset>\n",
      "Dimensions:     (station: 19, time: 142831)\n",
      "Coordinates:\n",
      "  * station     (station) object 'bee' 'dan' 'dep' 'for' 'fsl' 'gin' 'gld' ...\n",
      "  * time        (time) datetime64[ns] 2004-08-27T01:30:00 ...\n",
      "    latitude    (station) float64 37.53 37.9 37.63 37.77 37.48 37.77 37.52 ...\n",
      "    longitude   (station) float64 118.3 119.3 119.1 119.9 118.4 119.8 118.3 ...\n",
      "    elev        (station) float64 2.768e+03 2.987e+03 2.307e+03 1.72e+03 ...\n",
      "    network     (station) |S4 'CDWR' 'CDWR' 'CDWR' 'CDWR' 'CDWR' 'CDWR' ...\n",
      "    Grouping    (station) object 'east of crest' 'west of crest' ...\n",
      "Data variables:\n",
      "    SWdwn_raw   (station, time) float64 nan nan nan nan nan nan nan nan nan ...\n",
      "    QCFlag      (station, time) float64 nan nan nan nan nan nan nan nan nan ...\n",
      "    SWdwn_proc  (station, time) float64 nan nan nan nan nan nan nan nan nan ...\n",
      "    SWdwn_QC    (station, time) float64 nan nan nan nan nan nan nan nan nan ...\n"
     ]
    }
   ],
   "source": [
    "#####################\n",
    "## convert to xray ##\n",
    "#####################\n",
    "# Concatenate \n",
    "grobs = pd.concat(grobs_yos,axis=0,keys=stations)\n",
    "grobs = pd.DataFrame(grobs)\n",
    "\n",
    "# Convert to xray\n",
    "grobs = xray.Dataset.from_dataframe(grobs)\n",
    "grobs = grobs.rename({'level_0':'station','datetime':'time','SWdwn_Wm^-2':'SWdwn_raw'})\n",
    "\n",
    "# Fill in descriptive variables\n",
    "grobs.coords['latitude'] = ('station',lat)\n",
    "grobs.coords['longitude'] = ('station',lon)\n",
    "grobs.coords['elev'] = ('station',elev)\n",
    "grobs.coords['network'] = ('station',network)\n",
    "\n",
    "# Append station classification\n",
    "os.chdir(dir_data+dir_grobs)\n",
    "# Station data excel speadsheet\n",
    "stdat = pd.read_csv('All_StationSummary.v2.csv',sep= ',', \\\n",
    "                index_col=0,na_values=[-9999,'NaN']) # Read the supporting information\n",
    "# station data -> xray structure\n",
    "stdat = xray.Dataset(coords={'station': (['station'], stdat.index), \\\n",
    "                             'Grouping': (['station'],  stdat.Grouping)})\n",
    "\n",
    "# Merge into ground observation structure\n",
    "grobs = grobs.merge(stdat,join='inner')\n",
    "\n",
    "print(grobs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#### Elevation Angle\n",
    "d = pd.to_datetime(grobs.time.values)\n",
    "el = xray.Dataset()\n",
    "for stat in grobs.station.values:\n",
    "    el[stat] = (('time'),solargeo.AVG_EL(d,\\\n",
    "                     grobs.loc[{'station':stat}].latitude.values,\\\n",
    "                     grobs.loc[{'station':stat}].longitude.values,\\\n",
    "                     8,'END'))\n",
    "el = kray.combinevars(el,el.data_vars,new_dim_name='station',combinevarname='el')\n",
    "el.coords['time'] = d\n",
    "\n",
    "# Add to xray Dataset\n",
    "grobs['el'] = (('station','time'),el)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAegAAAECCAYAAAA1j0ToAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX1wXNV5/79Xq93FWWkTk+4MWCqKrTGYvtgTVaZxCjJ4\nGup0nB9peItxRBp7JkH5EYcKHCtWwNBWYOPYJppQh6RmMt0QjNM4JWEybTMYrBSnGA2EkBCbX2Ww\nIwsSGdxK2tq769X9/bE8V2fPnnNfVntXd3efzwyDvHv3nn253/uc85znxTBN0wTDMAzDMIGiYa7f\nAMMwDMMwxbCBZhiGYZgAwgaaYRiGYQIIG2iGYRiGCSBsoBmGYRgmgLCBZhiGYZgA4migX375ZXR3\ndwMAfv3rX2PdunXo7u7Ghg0b8PbbbwMA9u/fj+uvvx4333wznn32WQDAuXPn8IUvfAHr1q3DZz/7\nWbzzzjv+fQqGYRiGqTFsDfS3vvUtfOUrX0E2mwUA3H///bj77ruRTCZx7bXX4lvf+hZOnz6NZDKJ\nffv2Ye/evdi5cycymQwef/xxXHbZZXjsscfw8Y9/HHv27KnIB2IYhmGYWsDWQLe1teHrX/86qJbJ\nrl27sGTJEgDA+fPnEY1G8Ytf/AIdHR0Ih8NoampCW1sbjh07hhdffBFdXV0AgKuuugo/+9nPfP4o\nDMMwDFM72Broa6+9FqFQyPp3IpEAALz44ot47LHH8Nd//deYmppCc3OzdUwsFsPU1BSmpqYQi8Ws\nxyYnJ/14/wzDMAxTkzR6fcGPf/xjfOMb38A3v/lNzJ8/H01NTUilUtbzqVQKzc3NBY+nUinE4/Hy\nvWuGYRiGqXE8Gegnn3wS+/fvRzKZxHvf+14AwNKlS7F7925kMhmk02mMjIzg0ksvRUdHB4aGhrB0\n6VIMDQ2hs7PT8fymacIwjNI+CcNo2DQ4hKMnzlj/nhcNYf/9a2Z9TgDYsbHL9fNOr1Edd9OWp3A2\nnbPeN/29pG2+43nmEtYy4weyloHZa6FSWgag1fOPdl6nfL3h1CxjdHQUd911F7773e9ixYoVWLBg\nAZqamgAAf/qnf4rbb78d3/ve9/DEE09genoaPT09+MhHPoJz585h8+bNGB8fRyQSwc6dO/H+97/f\n9sMAwPh4ZVzhiUQzj1VHY/XsOoR0Ji+G9pY4+rv1E0bdWAPJYQCwfS0dN3JqomAs1WN2Y6nOFY2E\nrM8QdAMNVE7LQHVdizzW7MYaSA7j+NgEyHL5qWcvWrYbS3U+Uc86A+24gm5tbcW+ffsAAM8//7zy\nmBtvvBE33nhjwWMXXHABvva1rzmdnmEqwp7ela4NrApRWAPJYc/nGB1POR+kIRoJoTURKzjHibcm\nSj4fw1QzpL251PNsUOlZh+c9aIapVvwWoXjD6Nl1CJlsznqcZsrRSMj1+xBvIgDQmohZ/267iGM6\nmPqmUnqORkKWlmeD6MVzC1cSY+aUgeSwJYQgj9Xf3Yn2lrjWnUbGdOTUhCVE0wRGTk3g+FihkdXR\ns+sQenYdsv4tzrBnswJnmEpQSS3Pdjwveha1LOrTDlnL4iQdyOuZ/j0vGip6PcEGmpkzRBGUU9gq\n4ZZjrP7uzpJm7aaZn4Xb7ZWRUU9ncpawRWNuZ9gZZq7xS8t07qDoWYxj0elZpWURcnETdt4wNtBM\n4BlIDlvRkG6OVQnX79WoOCPf07sS7S1xiEHMJMieXYdc31DkWb7476AHiDGMCq+r3iDoWdYyPT86\nnnK9opbvD261HLr33nvvncXnKDv/+7+ZiowTi0V5rDkeq2vZAvzy9bdxYTyqncmSQE//zzn88vW3\n0bVsQdHzP/3FmPX4T38xhjOTaQDAhfGo9fjL/3XaevwDFzUXnWc2n4veQ393p3XermULcN2VC/HL\n19/G/6bPY+J/sxj/73PI5UycmUzjX4+cxJoVH7DGWrPiA/jXIyeRy5mYnjatz0r/id/ZT38xhqdf\n+A2u/dM21+9xLqjUdQhU13Vfi2N50fKZyTReOvY7/NkfX6w8Zi71LI4v/idqOR6L4PGn/x/SmRxy\nORM/fO51S6+xWBTXLLsY/3rkJBpDDchkc/jhc6/juisXKrXctWwB7t37vFbLvIJm5pRS3caAenat\n21ty2nPSnV+ctcv7Srr3ID5He01ycEg6k7OOl8eh/S7VSoPGk3NBGWaumY2WgbnXs5PbnLQ8cmqi\nIGhM1utAchitiRgy2fzetWkCG7Yf1L4nOy1zFDcTaCjvMNwYwpfWftD1a7w8rkJOwxDzLnt2HcKe\n3pVFrxFdbXIENgAYBiBXHdg0OFR0nO79yGMwTDVBWgbyhTvc5kGXW8+UYWGnZ1lncpDXogXxgnsC\nodOzeJw4+XDSM6+gmcDT392p3KfxOosuNepzdDxVJEQV4qpYhALEFi2IFzwmv+doRB3NKUeUAkAD\nF+liqhC7VXal9ExR2U7H6M5N2lXpWaxPIO5bk7ZV2R52sIFmqhq3bjWvUZ/izUKMuDQMKFfPdq+n\noBARVVR2ayJWZKTJTS4z7WLCwDDVht96FvVlGHnDWYqeZVoTMQwkh63SnYaBAgNOevfqAWMXN8No\nEG8UuqpF5LYbHU+hNRGzvbmI4qQby9jpwmhUqng2Op4qyJUUiyW4Wc0zDFOIUwUyOy2rXiPrub1l\nxiCLZUgJccXstkw9G2imLqDqXvS3HZsGh5A9nysKSlFBwpVn1fKeV2siVlAJDMiLWueqll1f4p4Z\nwzAz+cRe9VwOLQMo0rNYkIiQa3cTplmscRVsoJmaQ5ztiu4vEoRce1c+3m2NXtWxYolPcdyRUxOI\nRkIFBfJ15xSFLtbt9VomkGFqAZULWxdYKa903epZdRytpkU9k5aBwtWxDLnS5fsBPe5Wy2ygmZpC\nnu2KLmKn4+X9XvFvWn3b7Vc5BX2kM7miPTDaSzaMQvHTY+J467epUzUYplZRrT7dalmOkvaiZ9W4\nMuJEmvacZQ+Z/L5pPLd1uTlIjKkLWhMxqyKQ20ANsfqXqnSfHHWqKqhvt9cUCc/caFTu60g4VBCp\najdjZ5h6wa2W5ZW3nZ5LyasW9SsGhAHqFTIFkg0kh7Gnd6WrfWg20ExNQBe+quQmiY7SpXTCBOA6\niEMcj9xhomAJ8TGaQZOre0/vSixpm69dEdAeF6VkcA40Uw+Ik1K55KYXLVNa4rxoSKsxcUw6ByHf\nA+QocJoo0GP93Z3axhd0jNf64YZpBiv0pBoahvNYczeWSkiiO8pNQwpAnV6hcmvR+ei1tCcs70OJ\nTdhlxAL74rGGARgATKhX0E771bom70GhUtcGUPvXfa2OZbdnXC4tz4uGrPQnOueG7QdhmsWaVGkU\nmFkhi7Eg8rFOiIWK5KJFOi3zCpqpGmbbwUZevcrQDFycBY+Op4oCzURBqla15D6j/2hMGdPM70EH\na4rMMJVhNnp2q+X2lnhBtyhqcEGak9vBynqmUp2j4ylrRU75017Lmpaicw4SY6oeMVrbSTSqmbb4\nuv7uTiQSzbhj1zM4PjZhGWSdy5vSPMQbjDyjphQptzNtgqO2mXpjNloGCvVMr5f1bIeoZ3G1LL6O\nXOvrtx10dJ3rcGuseQXNVA12gRxOFYhK6Skriki1v0yzaPm1sjHXidFtsQKGqUXsGmF41TI97lbP\nqsdFPVOxElWrSRG/J9G8gmaqCrdlAMVjxf0qN7nN4cZQgavLMIqLEojHe10ZE4sWxBFuDHFnKqZu\n8apnr3UKAOCh3msK9KzLRaYoazo/1S5wu9ptb4mXfC/QwQaaCRRu3Vu61wEzrmSq+iMKURSqOFum\nc+g60YhubhJsOpObVWT1yKkJbdQnw9QCpehZZ5Dloh9y5S5Rz6J7+qYtT1Vku6jcxhlgFzcTIEoN\nGhFfJxpMOwMquslon1lVqk/E1AR0kZtO3I9yuzdF0aUMU2uUome71+iMrOzyVvVfl8+jen50PGW5\n3eVOVHO1HcUGmqkp5I5Q8j4S5U3STUA0yqaZd2e7yZsUodWBKHoO8GKY2dPf3VmgxUUL4gWeLFHL\ncq0AirZ2i9hhTvSSuWlP6RdsoJnA4KWajzxrpp7LwIxRFs8jBnnJEZkipul9Vbth+0Ff3FsMU824\n1bNcmITSp8hlLeoZ0Gs5nckp6+C70TOdmwx+UOA9aCZQ6IQsdqTR1duWiwjoOsmooHaOZLAz2Zxj\noRCAV8oMY4dKzxSI+aW1HywK+gJmNCVq221RENNUB4HJhUFkgmSURdhAM4FHJWK3r/Nad5vGcdsO\njmEY95SqZQBF7mudwZW1DFRvMSB2cTNVxfGxfCCYqt62+DcwU4uX3N+6ZhOlVAViGGb2jI6nYBgz\nHi9dLX1gZsIcjYSwd/MqZeBWrWmZDTQTePq7O7Gkbb41a6bobFXVIVWRA7Gpu9zukW4EgPsuVwzD\nlAYZ4CVt860tKSqlKTegoeNVegYK96JrVctsoJmqYMfGrqJqXl7SOOhYqqULFLq9BpLD7NJmmArQ\n392JHRu7Ch7LZHOeUrJkvdaqlh0N9Msvv4zu7m4AwIkTJ7B27VqsW7cO9957L6gR1v79+3H99dfj\n5ptvxrPPPgsAOHfuHL7whS9g3bp1+OxnP4t33nnHv0/B1BS6Un5OBfLpdfLeMxUuIFT7UbUy42aY\nIKHTMlCoZ7mfsvg6lZ7rBdt2k9/61rfwwx/+ELFYDPv27cNtt92GDRs2YPny5di6dSuuuuoqLFu2\nDOvXr8eBAweQTqexdu1afP/738djjz2GVCqF22+/HT/+8Y/x0ksvob+/3/EN1XJrNR7LGV27Od1Y\nonjlWbOupB9hGHk3WWsihuNjE1UXSMLtJmeo9uu+Fseyax2pGkulZTGTQgwMU5XVpFawYjZHtVBS\nu8m2tjZ8/etft1bKr776KpYvXw4A6OrqwuHDh/HKK6+go6MD4XAYTU1NaGtrw7Fjx/Diiy+iqyvv\nxrjqqqvws5/9rJyfh2EsdJWBWhOxguIDMrSfPXKq+owzw9QispbFvGZxi0vX5pW2sWoFWwN97bXX\nIhSa+VLExXYsFsPk5CSmpqbQ3Nxc8PjU1BSmpqYQi8UKjmXqD9lV5bS/5Ka4gdN5KGqbAkx0/ZgZ\nhnHPXGh50YJ4QVQ3adlu4l1LeMqDbmiYsedTU1OIx+NoampCKjUzm0mlUmhubi54PJVKIR5Xp7gw\ntYuuoIhTFxqnqkO6PEpyWYv707UUMMIwc0UltNzf3VlwbsMofL2o5XrZh/ZkoC+//HIcOXIEV1xx\nBYaGhrBixQosXboUu3fvRiaTQTqdxsjICC699FJ0dHRgaGgIS5cuxdDQEDo73eWmJRLNzgeVCR7L\n37HCjTOr1gaj8HE373HT4BAAWBGfiUQzxk7PCFM8PzDjsn7w8ZewY2MXNg0OFRzPVJZKXoeVHq/e\nxiq3lgEUaTmRaMbCi+NW+1XThFbL2Wx9TLptg8QAYHR0FHfddRf27duHN954A3fffTey2Sza29vx\n93//9zAMA9/73vfwxBNPYHp6Gj09PfjIRz6Cc+fOYfPmzRgfH0ckEsHOnTvx/ve/3/ENVWtABI+l\nRm4dp/pb9zqaXRsGcNkl85E9n7MeEyO5KbjEa1nAaoeDxGYI2nVfi2PZ6ddOz6KWyT09djpl1ciW\nszLEVTQdL94LgOqtDKZDp2VHA11pgnIx1tNYpfZgLvVz2UV3qo4hxIhO+XXy8U61d2sBNtAzBF1j\nlRqrVC2XMpY4pp2enerhq14jGul6oKQobqb2KbUH81whdrqxo9aNM8PIBFXL/d2d2rKcuuDNegkC\nc4KbZTAVx42LW4fY6YZFzDBzjxs9R8KFneHmRUNY8Huxgl7OrOdi2EDXObMxlrMd147jY4UuMVHQ\nwEzOI8MweeZKy07jqTIp9t+/Bnfsesb6N+tZDRtopmJidnvzGEgOF7mo2y6K40trP1h3e1MM44Ug\narleUqL8gPegmYrgdn/MLqBEfl00ElLubTEM4x9etayaUN+05akCw83FhNTwCpoJBE5BLa+dPIPp\nd1fVZJSrte4uw9QbclYFpViJWgZqPzXSK7yCZiqCXdk/cUYOoKBJuwrTzP9Xa3V3GaYacCrhKfZ1\npuP2bl6lPJeoZTbOxbCBZiqGXNSAEF1dYnCY6PJidzbDBAeqcy/X0hYn2z27DlnHDiSHC/TcwFp2\nBRtopmLo9q7E9ApxNp3O5GAYeTfY2XSOc5sZJkA47UVTZHbPrkNF3q5p1rIr2EAzc4K4aiZXmCpI\nhFxgDMMEC9HbRXq20zLjHQ4SY+aEdCZXMOsW3d+i8eY9ZoYJHj27DhVMnGU97+ldWZCKxemRpcEG\nmqkIqnxIMQJbbFsntptjiqHmIQwzFwwkh5FRdJPS6TlIZUeDyLyo3tvABprxHbkzVSQc0qZViEac\nXdtqTBNWSz6GqSRynQLqNkWIz8ld5hg1lHKmgvegmYpimvmgMDk9Y3Q8xW4wG9ibwAQR0rKs5+Nj\nE5wGWQbYQDO+09/d6Rg0ks7kWMw2sDeBCQJutAzw9Vou2EAzviDnR+7pXVlU3IAiPnl1yDDBxY2W\nAdazH7CBZsqOLj9SDhrh4JHSsQssYZhyYadl0TiLRpxXz+WDg8SYiiEGmGzYfrAuhEyrCbkfLsPU\nCnYNbqqFR/tWuQpqIz1X6t7FK2im7DjV6gXqY5YdjYSsQis60c+LFpYwNQxo9/jE4+wiPxmmXLjR\nci20k7TrvEW0t8QdCyfJerbDzXG8gmZ8QXZ/MXr2bl5VVNRBBa/CmblApWWqr318bKImJttuJhlu\nvATUt56+J3FFLnb0Mgxg0YK44znZQDO+Ue95kKpiDipEF6Fdqhnlm46Op7jZAFNRZC1zSmQxtCJe\nv+0ggLwnTDTO4gQ7Eg5Zk/GszX2CXdx1jhyhWc7z1nsepGpl8WjfKqtWcTQSwtl0rmAWrfu+DAMF\nldbYxc2o8EPPKi27nXxWC6Xcp9pb4ni0b5Wl5Ug4VFBASDynaRZ/ZwPJYaQzOdvGIWyg6xjqMqPr\nRlMu6AJmYLmv9/Su9PS6SDhUE8E4jH9UQs+k5Vpwa8+W0fEUBpLD2NO7Eq2JmGcjLzYb0cEGuk6h\n2Vu5ziX3hAVgBZfQBczMtODbsP2gq+PJbUavYxgVfupZDBSrJi37nY8talm1hy2PHwnPLFLSGXft\nc9lAM4hGQtoITSfkPEnx3wBKPm+tQ5HdFPWpit6mQBKG8UI59UyIec9B1LQqE8LPVb44Hmk5Gglh\nSdt8RCOhgoAwopSJDQeJ1SkUhUl/+0m95DzboRIsQY/LKyB6vL0lb6TFiFkOEmNEKqVnitwOMjSp\nLcXjZKdT8RhVRgUZYJ0nY3Q8pdSyHbyCrmNURe69sGlwCAAK8iTlvMl6Nc6ye0v3HTi5JUVPhLia\nvpTbTTISs9HzQHIYmwaHbPOeaY87iHqm90Teg1Jzs/duXlWkXZWWVbo9PjaBE2/pJwXkEpe1vKRN\nr2U20ExJDCSHcfTEmQIDQi5uErbc1L2ecPu57SI4iZFTE1bqBkWAM0y5ILf20RNntH2cy7nH7Tez\nSQFbv614QeFWy6bproDQ+m0HcXxswpWW2cVd58jRnrPduwJm8gCZ8iK67LgfNKOiHHoWtSwW3KgG\nqmUSYZozerbTMq+g6xgxIMRrekZ/dyeWtM23LQHIMEzlKFXP5NZe0jZf2ae93usZzCWeV9DT09Po\n7+/HG2+8gYaGBvzd3/0dQqEQ+vr60NDQgMWLF2Pr1q0wDAP79+/HE088gcbGRvT09ODqq6/24SMw\nQaBaZtgMw9gjarkaDLObwK5qxbOB/o//+A+cPXsWjz/+OA4fPozdu3fj/Pnz6O3txfLly7F161Y8\n/fTTWLZsGZLJJA4cOIB0Oo21a9fiwx/+MCKRiB+fgykBMfJTfMwNbstTMgxTGUrVc7Vr2a5GfbUb\nb88G+oILLsDk5CRM08Tk5CTC4TBefvllLF++HADQ1dWF5557Dg0NDejo6EA4HEY4HEZbWxuOHTuG\nP/7jPy77h2BKx4t7ulJpWQzDlIaXCbaX44NMayJmm1I1Lxqq2tK4nvegOzo6kMlksHr1atxzzz3o\n7u6GKUxRYrEYJicnMTU1hebm5oLHp6amyvOumYojFzDo7+7EvGi+7F+1VBZiGEatZdqD3tO7suqy\nBOyMs9vI6qDieQX9j//4j+jo6MDf/M3f4K233sKtt96K8+fPW89PTU0hHo+jqakJqdTMXkYqlUI8\n7lwVKZFodjymXPBY7gk3hgr+/r+7D1kXPpegZFRU8pqv9HjVPJas5USiGQ/1XoNNg0OcgREwPBvo\ns2fPIhbLr5ji8TjOnz+PP/iDP8CRI0dwxRVXYGhoCCtWrMDSpUuxe/duZDIZpNNpjIyMYPHixY7n\nHx+f9P4pSiCRaOaxPJA9PzMLPXbyTFXv6zCVoVLXPFAbGqvUWKKWX39zAuPjk1W591wPeDbQGzZs\nwJe//GXccsstOH/+PO6880784R/+Ie6++25ks1m0t7dj9erVMAwDt956K2655RZMT0+jt7eXA8QC\nhrwP5XZfio0zwwQLlXbd6DmdybFxDjCGaQbrdlvNM9NqGkuM3GxviRc0Y9flNtdr2c6gIEek/mjn\ndXP3ZlzAK+jKjCVrGYCjnkUtRyP6KGimMui0zJXEGFeF2+u5bGclqPZ0ECYYuNHyQHK4wDgz5cdt\n0w2nY7iSWJ3S391piVO8SFSt6qqpDm+1MS8aQntL3LGtpNx6krtZMYSdllWrZ7kQSSbL2i4X86Ih\nPNq3ylHP0UjIVStZNtB1jJgeRWLe07tyDt9RfWO3mhGbvQPummww9YNOy27ynNlzUz7aLopbqWty\nFywRt6mp7OKuY9z2kOUynv5xNp2z9g/tDHQ6kwt8H15m7vDSD5pXzP5BjS+caqC7vaeyga4zZBE7\niZnd25VDvHGSsRa/e7uShkx9IurZbVlPXjGrKWcciBgLQNtTbuIDZNjFXUfIFYTsjqu2NnO1gCze\n1kSsYFXNKx9GxI2eRS0DYC+MgmgkZO0J27mlS4W2pxYtiBdMvN1MtnkFXeN4rbkrpmxwWpW/OM3Y\n5Qpt4oz8skvm+/jOmKDiRc9yEwyA95uBQt1RqWKa6JTr/CKZ7Mw2lvyck5bZQNcwcuN1u30q1Yq5\nHsVcyXQn3ThOqS+RcAg7Nnb58I6YIONVz9XWNrJSiLrzo4+A3f1Dfs5Jy2yg6wzVzFsUfr0zG+Ms\nGneaKTvmORozEdli4xHd7xGNhDjSnrFgPdsTjYSQzea0WQ+j4yns6V2JgeRw0b6xnLIG2E90opEQ\nGoyZ5hxOrzEMOGqZDXQN4yays2fXId7bLBNyXqO432cYM3tRJFjDAC69ZL5VG5l+L90+oa7CG1Mf\nOOlZXjUzea0taZuv7R9ApU5lZC3L36s8ATeM/Go83BgqqHUualmu2OYmD5oNdI1Awn2o95qCx+1u\n6FyD1z9U+8fyd22a+bQMWhXLKx/RqLcmYmyc6wiveuZVs54Tb9lHT4u6pIhrmgyJJVNF5POZJorK\nrcpaJu/Y6HjKtZ7ZQNcAojg3DQ7hS2s/aD0OqJthcPpU+dEFgjiRzuSUUbhk1HnlXF+o9GzXDAPg\nbAsdhpEvHkL5ye0tcduJTCScr6Q4m8ULGWARMuBei0FxmlWNIqdgiP/u2XWIZ9sBpL+7E+0t+VQP\nP9I9mOpElU4lPjZyaoIn2xpME3jt5Bnr38fHJqwVrohhoMCTNZvvk1bH5dAyG+gagC6G9pY4dmzs\n4r2oOUYOFCOoBKMYpU21uGllRIEqpllYS7ln1yHlXhlTe8h6ZtSo9CQ+p4K0Sa+jc+zdvMpz8CX9\nRrKeScuj4ylLy4aBkrTMBrpGoEpCmwaHrBk1XXwEXVB7eldyFxsbyGjaoTK+Tt+pat+p7aIZ4yx3\nDMtkcwXuNl1AC1N7iJXBRIMtBorRYzojVS/oVrv03Tz51cJWjiOnJnB8bEJbLERsPiIiPmYY6upt\nbRfl7xvyKpx07VXLvAddpbgpWCCn7JDR4eAwPdFICG0XxZE9n7PNiZbLbqrSo+RjaHtBfOzEWzPH\nczR9/eKkZzGORK7d3pqI1aWenerT03e2aXCoKILaKf2RipcQYkETej39Zjo9y8heTTd65xV0FaIr\n8Ud/02pO1WaO96vsSWdyOHriDEZOea+bKyLmNIvIoqQAMdVsOhIOYSA5XHAePworMHOLnZ5v2vKU\ntownrcTqOZ5E1qi4lwzkFyNHT5xBOpNT7gWr7pOqY/b0rlRuG8qPUTcr3cShNRErSNFyaqrBBrpG\nIJHTxUg3fjFYgQ1z+ZDrZAOFYiVDKt4UVG6z6XejO9OZXMHeNc34R05NcDxBHUJ6pm5n4gqbAwj1\nLFowEyUtu5nllq1A8YSXAmplLcuGNBoJYXQ8VTThfu1k4eSetCzqWfU+dLCBrkLkPSnAPs2CghUY\nb9jdCOVUCjlf8vjYhCVUEmgmmyvpd2hNxOp+n7GWUelZhxynwLhD59HSRcartExxPaRlcVKtwzTz\n+pW9YG71zAa6SpEDFFQN2+l53tt0j2iU7Wa6TvtfsnBlMZPQKSBNHJfyn+XfkWbgTm4xpvqQ9UxG\nW47yr2ctOwVuErQKpu9wSdt87OldWbL3wa2Wo5EQnvzqdUWGV4wBUul50+CQdmwOEqsBxKhOudQc\nkHf71PM+lRdE4TltCchF91VVh+TgEqCwBm8i0Yzx8UllJSi7akObBoc4BadGEd3ZiUQz7tj1jGVw\n6lnLI6cmrOp6Om2KGqTvjDSm8j6IBlOMkAdQ1L9ZpWWgWKeqY+j8KuwCy3gFXeWIbhnCTc9nRg0V\nvHd7rNPzNHMXj6UavAPJYWv2LG9RyO5O0Q0KwKqMxNQWVESI9Etpk/WoZ9q/FZFL5orHyMeKHi5V\n8RHVloIqdYrGIS2r3NNy3235tapx6DzUXEMFG+gqR7yxq/ahuWqYN1oTMVzq0GuZbgq0l2QYsPIq\n5XMRe3pXFuWyUlCfHKmtQ3fzYGoDpwpWo+Mp222VWsNNMwk7j5cYJe0mAhvIH79h+8GiLA5xu6u/\nu7NAz8C9knmzAAAgAElEQVTMomjT4FDRhFwXrd3f3emoe3ZxVzmiO6U1EcOOjV24Y9czAPQtCxk9\no+MpLLw4bpsDTbN4p++XbgBuctZV9X91r+vv7sSDj7/k+FmY6oXSf8jFrZoA1jrluH8dH5vATVue\nUj6nqoEvj0n3ASosQqVAgeIeByJ7ele60rOTlg3TDFZM4Pj4ZEXGoX2Jah5LTAegfRAaa/22g2Uf\nj7FH1QNaLJCgiiQNN4as5iZAodjFghTyjSGRaPbjI5SVSukLqC09i3uaZKDrZbLtpu+yCHWfAooz\nKVTH2lk7uZiJ6r2I5xD1TL/bQ73XFFwbKj3LFeHkjmUi7OKuUshFKparG0gO47q7nsSG7fVrnOcq\nR5SCV1QRn+LfYoMDAAWBXroVM83e63Uvsh4Q9Sw+dtOWp+rGOAP5hYboPpb3e2V9i20endzFVBPb\nbmzastJNFGSXul2zElHP8lakGDtkF8XNBroKUTXDoCph06ZzGTum/Kj6PXtBrialqwfM1AdioZJ6\n4vhYvhwuxVvIxUB0qY+j4ylXaVS6e6NYUYy0XE49l1oNkA10lSHOtMX8OyZPtUxO7EoM0uSLmpqQ\ni9ttMQum+pCjg+u1epy430t7uHI6oy4feiA5rDTgqqwMuZHQbEroyvUKZGjyIP6+bjuW8R50lY0l\n5suKEYRM+Wgw8iU4/cIwgL2bVwEovDbEoBL5t5VvALwHXUi16llElQtfr6j2i0kDG7Yf1MZ5tLfE\nrVU3HTMvGsK05OES4zrEqoDya928T9IyUJqe7bRc0gr6kUcewSc/+Ulcf/31+MEPfoATJ05g7dq1\nWLduHe69916Qzd+/fz+uv/563HzzzXj22WdLGYpBYY6dl7KAjLs2kDKlGme75uzic6aZF6+YBw1w\nIwyGIRYtKM41pr1b0Xiq9Land2VBitbZdLG7WtSavJdM5Xnt9EyQlgH4omfPBvr555/HSy+9hH37\n9iGZTOI3v/kNtm3bht7eXjz22GMwTRNPP/00xsfHkUwmsW/fPuzduxc7d+5EJpOZ9RuuN8S9DLoQ\nxHzYenWFuSWdyZWtPKKqcIKIKkhMfE50v2WyuYI8aGCmDy3dFOx+W7vAEqY6qedGGHK5W7nWPZDX\nsqwJXVlcJ8jY22nM9BDPI9c1AGZ+TypUohurrEFizz33HC677DJ8/vOfx2233YZVq1bhV7/6FZYv\nXw4A6OrqwuHDh/HKK6+go6MD4XAYTU1NaGtrw7Fjx7wOxwio8vbquTavW9yKjFbb4p6VXHBAfoww\nDPvfIpPNWXvKFPFNjI6nivbc7NqCDiSHuZJYjRKsDUd/kCchFI8hrnrt6gyQ0ZP3kFWV96jevQyd\nX6dZp/uq2KWKyvwSx8cmrBayZOQ3bD+ojfS207LnQiXvvPMO3nzzTTzyyCP4zW9+g9tuuw3iNnYs\nFsPk5CSmpqbQ3Nxc8PjU1JTX4eoeVQELgjvblBeascv1zMUbRTqTU+79u+lqQx4QcdbfYMAxUtdN\nz1qm+qmH9EjKW5Y1JRYBsctnFh+XDbi8eJFzyVWvI3e2XEvfjZ7l89G+OU2wdd4QL3r2bKDnz5+P\n9vZ2NDY2YuHChYhGo/jd735nPT81NYV4PI6mpiakUjOzilQqhXjcuXRbJYNfqmWsf35gjeUGoYi/\nm7Y8VXeVhfxGFO6Stny5T9Xs9vWxCUQjIeuY106esfat7QLM5N9r4cXFdbVptt920YxWdmzsKvj9\nH+q9pipc3JUOZKsWPctsGhyqOY/IvKi6xrRpAmOni1292WwODz7+EsZOp5DN5iwduE0ze31sokh3\nDz7+knW/pGIgqqpi2WzOen7T4FCBnu2Q9XyBVFf7skvmW40w9t9feA+nv5207DmK+9lnn8U//dM/\n4dFHH8Vvf/tbdHd3o729HZ/5zGdwxRVX4J577sGKFSuwfPlyfOYzn8H3v/99pNNp3HTTTXjyyScR\niURsz1/tkZiVGEu3ombKh9sIebmykKrEnwy5yKny0A1ffgqZbA6LFqj30OTIfa4kVky16pm17B6v\nVcZErVDHOLlDlXhuWbu67lX0nHgckNdzuDGEYyfPIBIOWR3rZGQ921US87yCvvrqq/HCCy/ghhtu\nwPT0NLZu3YqWlhbcfffdyGazaG9vx+rVq2EYBm699VbccsstmJ6eRm9vr6NxZpxhQfvPvGioqPWc\nDpXYWxMx5Y1ALA1IgX+iJ4T2rtidXR84NceoB5zaRzq9VqU/OqeoIzcpbLrfIxoJFaVeiQFs5Bof\nOTWBedGQtYW1YftB7aTbLZwHXUVjcZ5keZBr7oqPi/V0S5kMiTcNuTY3zejF86pc4qpoVDEydCA5\njHBjKPD9oHkFbU8t69mu7rVOf/Lr7Qy3nMcMFK+sxRX0DV/WbwmK59L9HrKRlrvSESrXvp2eAdhq\nmSuJBRBVb9FaFrOfeMmBzmRzRY3XvZ5XvCmZZv4mI+eti+e99JL5RcEkYhcs0TDLbSqZ6qaWUyR1\nJTkBd5knkXBxYxkReo5aP6qOFb9f1fPRSAiP9q3SuqJF0pmcUsvy+fbfv0Z7DlnPABy1zAY6YMg1\nXMXHGPdQGgaVyBTR3SDESGtAfwOlGbfOSIsGl3I3Ke9S5rWTZ5QrDdV1wFQnqgl3PW9V2RlvSqGi\nqGrSMaU0yc0zKJ0JgFLn4vcu6pI0bPc+ZOj3orRIFbrWlqXqmQ10gKAgBqKWZ9h+Q6kOPbsOKYsb\n6Ehnctg0OKTcj6IbBAlbnpXTzUS+AVHhfZU4VdGi8nnF90/5nRRBzgQb3YS71o2zbvVrGPlKX0va\n5hdV+YtGQlaxH7Gq1/GxCet81Bd75NQENmw/aB03cmqi4N4JzNwDqCOYbg9ZRSSs7nNA41NqmEg6\nk3MVdU56dqNlNtA+I5d/szvO60XEOFNKVxpKjRChGbdcDEEuZEICBvSlP+XC+XSsWA1JPLdcoKa/\nuzPw+8+1iGolXAr1MPHWfcZIOGR9h6RN0ohOq7LBFh+Xj3PCbcUxei+q/tAiKi2LTYyc9OykZc9R\n3Ix7RNe0l+hcClayK1LC+AfNgnXBLKK4dGkYgL7YP70WyOdqZs8X7n337DqETDandAW6afLOlB9R\ny5sGh/CltR909ToxG8BtZkAtoNtGEqOdCT/ClMUgtbaL4gWFh0TkGA+7+60cyU2eNNIuBW+K18ZA\nchgbth/UuvWdAj7ZQFcYOUiAoIsDmPnh68EVFmTEdCma9YrVhtZvK6z8JOZS2nXEoa2MSDiEf35g\nTUFEsHiDkH972UjwKnru0elZvvHTY/USSyJXC5Mpx31NzosW/y1rT464pt+H/k1V3OQuWWIZT5Ue\ngRlvQWsihh0buyw9i+PJK3FVBLgKdnH7iLjHIP4gqr1IMsa0t+FUyJ0pP+Rmpj0hecvBrgShSGsi\nppwxU6AYnTedyRUElQwkh7Urj+NjhXtsKjc84x+q/r06PXOAX6GLW97qMaTUQtU2kJsudE7bV+Ri\nVk1kKXCTsGuM0ZqIaaPEyfVOngFxO9Pu/q2KjVHBK2ifESvZ2CHefO0KxTPuscvFVEHH7tjYhRu+\nrI7GtKM1ESuYFcvjq7rxELJrTWyooVqNuy2ByJQPL52SVH/LQUy1jHgti5X2xEI9wEw/ZbnHs4yb\n3GlVS8n+7k5sGhyytpFIZ+Letx1ixS9Zz3bpYnZ69uI94BV0BbHr5WyXesCURil7W+lMDh+788mC\nQC9CJUg7kYvdeUSooxXlTeq2MkoJcGMqh07P8mqLgsuCVRLKP3SrX5V+TDO/VSR/N/K171V7RM+u\nQzh64ozlzRB/G3HiIGpSh6xnivQWAzx3bOxS6pk8Zl5hAx0Q7PJqmblDnDipbrC6ur5yTiaQv6HQ\njceuVq98XroRiDS4uDkx/iPuMYuPUWSvLgWo2nAyhuL1qcs+oWjsUryDbrUnvh95S+r42ITWg9Wa\niCmNqC7qW4w8N83idpd2yL2v7WADXUGc9qY4rSp42OVzEqJB1vVxjkZCRe7unl2HrD1oO8GmMzm0\nJmJ4tG+VNdt/8qvXef8wTMWgCH+i2lfPXjx8bmM1vOK0wiXIlV7qKpY0ZufBksuQUnzJTVuest1b\nptfs3bzKmsRxHnQVUE8RnkFElbPc3hK3qhnJQS6qfTUZ0QWqOoYKG4ycmlDeAFVGe0/vSlelCZm5\np9r03GBjAFsTMderPjelPJ1QeaFME1ZlQLmGgIgYea06RnbB0+qXzqsLCnPS88ipCZxNq426rh7C\nnt6VttkY3CyjwmOp0jI417k8UCCGXYqTW1QdbMRxRAMtpnnIrerEFBtdqzsgPxmwe54+G52T200W\nEqRmGU6/dZDRBVbaBWnJTWHsjndqoiHW2AaK743kbtYFlamabCxpm49jQkld8T4hdpsSF0l2WqZj\nS/2NDaNwP9uupgEb6Dkcq5qFHETEvEU5gtKv71gUsmocEqPTaspNlKp47MKL1ekjQaIeDXStesK8\nXJ+A2lAWeKGgLnNLx5DhVBlooHzeCbEmt9NvRxHn5f6N50X1TTY4zapCqPaca1HIcwkF5MhuunIa\nZ1XqlDiOvPKmwJhyks7kuJvVHKMrUFKLtQvomvYC7feKeiko02nzOmCmjr5cqW90PGX7XrxOJESc\nKomZJqwezypKXQjYpUzyHrRPiHV7Nw0OWcFh9F+1R3UGmVKC7XR1s8Xno5EQ9m5e5bgXZxfZLe6r\ntbfEMS8a0nbHUgXFuN0HZMqHXINb1DM9R12VvBiHasnaUFXmqtS4QHF3OJ1rm3REOrP7fun41kTM\nKgwFON87TDM/SVDV3967eZWlZ1UWh+pxJz2zgfYBN5WE2K3tH6VMfqh3sw5Ki9Lls7qNMAVQYKgp\nqIRyNCmvkt6T6n0wlcNJy1RNShe9b0ctx51QZzfdatMtcjU/HWI5Xkp9szO2ixbErdW52GmOAsHs\nJuxUcKg1EbO0TPcGMUhMPo9doSIdbKArwI6NXdYP6SUHjikNUdBOK2MRuxtmOpPDhu0HlQKLRkJY\ntCBupWWo+tYCxTmgcpUp1U0+k80VdMbhVLy5R9Qz/x5qSIOluPx1E1035T9Ffekm6pQjLb83VelP\n0fMlr+TFseTJG+lZnlxksjnr3qBqTSvDBtoHVBWGxIIGvHquHNRfthzo8ihVM2MxUE1OxRJr+JJL\nTG6jJ66ixTHp2hK7ATH+oasWpipQUq/ojKZbr4I4mdVt9xDiqlV3DLmsdboXW0mKRlc12aKgUzK2\n8iRB9KAAsPpc69piiv2tKc3KTsscJOYTuu42TGVxG9wiBnjMi4ZwzsW+m/gaEreuSIPYrUwMfKFW\neHIgmRwBK7rh6HVMZWBDDDzaZ18vm2oGEF5SssRrnzxVqpWlGARq57mQx3bSM0Wci4sp8X6tqvEN\nqL0Dujr+crCcWy2zga4AFFQCVE9gSK3gdvVMN4TWRAzhxpCrKGm356aqYST00fGUdZPZsbELd+x6\nxjpW1VBAxK8qTYx3KE2yHhhIDmvTBb009/Gyn0zIRp3G021d0eMU52GnZ/GctAdthzhJEIPYyMBv\nGhxSes5UunWjZXZxVxi+uQYPsUZ2qYhuazqniLiSdxpLdZPgPehgIbYNrQfo+i4FNzEgYnCkCi/B\nZmJqYyneDzEw8PjYhPK9iXoWA9RkVFshFEDnRsu8gq4A3Ls32Iju5NHxlKdGFFSIRF5JyYFqssta\ndR4SOuV/Uv6pWO0IyN9Awo3siZlLajHf2Y7ZlO90U2HMaXKs8lTIGpMnS1Qf+6Hea3DHrmcci0LR\nbyoHi8nvzU7POzZ24f/c9WSR0W5NxKzPsHfzKut4Jy2zgfYZCr1nggfNisXavV7LE5Jryy4H1s44\ni9sfhM5lKBbHqIZSn7WK6MasF8rlKXD7vYlG0GlsCuCjGA8gP6EgHV1315MIh0OuKvoBKCqOIiPq\nWb43bBocKvi3yo0trrSdtMylPn2G62xXD17rEItlAgHY1gfWnbfBmKkkpKv/LZc3bG+J29bvDQq1\nWOqzFD3PprpVtWBXu74U3Jb0JA/WbGMzqIwn4bacp1yXX9YzoG5hKeZK22mZ96B9hI2zN+Y6P9yp\nv6yMvIck7pOJn0U+L+1pUTcr2l/e07uy4By6XrRAfqbOVJaB5LBnPRtG/bSRVV27blDVDRg5NWEF\nU4r/idAq267FqwrVfUaOGqf0OrvX0Ovs9Kz6nPL1YKdldnGXkYHksFVhhqNt55ZyzujTmZyVRiI3\ngBeR0zOobrCqWIH4mK7Zu/i4fG6uxV0d+FGLPaiIevCyZ02akl3L4ipThdM9VnwPZMxV49BYchS3\nrDnVHrZqfNKtOKGTvW1utcwr6DJBLhHa92Dj7J1ybraU090GOBddoHrMJHA6XpXPqXpfcr1nYCbI\nhZ6n6kf1FqAUFDgfWo+sN6pZLUMlceXiHHbXtFgMhFbWbu6xKmPqJiVM1KKoZy/3E9Ir0ZqIFWl8\ndDzlmKbHK2iP6LrYMMHCjZhUq+z2ljjGTqdcBfaZJqzAFLpZyEbWzUpCnJmTq08MclHtbXMlscpT\nb8WGvOyduzVeXqrx6fDSlrU1EXPUs9jvWdx7pgm31wmxqGdx5U7ntSv4IlOygX777bfxiU98At/+\n9rfR0NCAvr4+NDQ0YPHixdi6dSsMw8D+/fvxxBNPoLGxET09Pbj66qtLHS4QiD8epcIA9ZdyUQ24\naf1GdbNFwY+Op7Dw4jhef9OdF0RVfEDE6ype5UYLVhhn/VIvrmqiVC+gGLjldTwyaOXwgJHRd0qb\npPu4yoNVShyRGOVNExI5dcstJbm4s9ks7rnnHsybNw+maeKBBx5Ab28vHnvsMZimiaeffhrj4+NI\nJpPYt28f9u7di507dyKTyZQyXCCh1Q25Wtil7R9OwWOqhhheqnzJ/z564oyVh+z1PZZyLSxaMFNb\nWPe+KUiGxuHUvcoxkBzG+m0H5/ptVA1y+VuC3NsU7Wz3Wl3Pda+IgVsy9Bg1thCzJGji7nXsedGZ\nlpfiPSmTLaz7LetZR0kG+sEHH8TatWuRSCQAAK+++iqWL18OAOjq6sLhw4fxyiuvoKOjA+FwGE1N\nTWhra8OxY8dKGS4wUGQfl+usHG5WwtR5pty4ib71Eq2qYl40ZNseT4z8bE3EZt2+j/GG23QbphDV\n6jkSzruS05kcpoXGE/L91ItRdOrhboe4ytU95wXDAPbfvwZA3laI8SdidyxyqctR3io8G+gDBw7g\nwgsvxJVXXvnuwCbEVOpYLIbJyUlMTU2hubm54PGpqSmvwwUO6kAito9kg+0fXg0vldGbbcqW22Cs\nTDbnWKHIDrFQfn93Z9G1JNblJkMx1+lo9US91NouB+K1q5pspjO5okhvN6tUu+vdTnelaHI2JU3l\ngFDdBF9sfOPUbtLzHvSBAwdgGAYOHz6Mo0ePoq+vD2fOzISJT01NIR6Po6mpCanUzA0ulUohHnee\n/VeyQtJsxhKTyz9255PleDsM8u3aXjt5BtMlGrzLLplv/d9NKtK8aMia0YvMJjimwYD2/YvPNRj5\n0oAiCy+OW6VhyXiLpWLDjSH88KvX4aYtxR1zgkalq52Vc7xNg0OcyvYudtezeEyDMRO8GG4MYcfG\nLmwaHMKJtyasLRldG0Y7LoiEXG3pNLzrQlcdu/DiGS3R86rPNRt3Oo1B1+FDvddYOm27KF5U8jnc\nGMI/P7DGVsueDfR3vvMd6+/u7m7cd999ePDBB3HkyBFcccUVGBoawooVK7B06VLs3r0bmUwG6XQa\nIyMjWLx4seP5q62SmFhejpk9x06emZW7mm6qbr0a5zI57N28qqx7jLq3bxhAWAggWfiuu5quQ9Gd\n2t5S2IaS3GJfWvtBjI9PWq60IFOtlcTYrV2Im8nytFkYF3H0xBncsesZqxaAm60qETHDwm28xbQJ\ny3UuG9rs+VzRZwlrSvCK71X3vsXobDoHjSHqmd67GHTqRcuzTrMyDAN9fX24++67kc1m0d7ejtWr\nV8MwDNx666245ZZbMD09jd7eXkQikdkOFxgo4o+Dw8qLH3vJdukipplPeygnus9AEZ2iy27T4JAl\nbNGlfnxsosD9pStmwpQfzsqwx236lTzJ8VI8SHccGTfdBErWFyEWNaFz6H5nN6t8UcuisRb1LH8e\nwouWuRZ3CfAMu3I82lfa6tbLjF0+tsHIr27dVoPzu9YyVSESc/AffPylIvd40KjGFbQbbbtx+dYy\nbrSlOobqXffsOoRsNgcT3ifkTgZaxZK2+coVrBjT4YcVVJ1XrChGeqbtABVcqKQEOHCkcpS6heBG\ncFS+UyV4LyVb7Y4rh/ipCpGYg8+em/LjduItGme/bu5BphTjTK8jPU8LLmQ35yTcVAKT34u891up\n2uiqzyTmXLv5HFzqswTqTZClUo5oY78M0aN9qyyjOnJqomDPetqcuRHI9Xx1n0mXMrFogX2U/7xo\nyMoGkM8tjie649g4lx+vXjFKl+F7QTF234kqatvP75D2sO00o+oe56XRhgjpWdXYIxoJWeV73W6j\n8AraA/VW6m+2BPnmJbvNM9niqmKA+6jTgeSw8vV2KVjRSAj7719juWbllTG9TmU4nKojMc6Ieva6\nZVXuWu/lZK5W9XYd3HTY9Ur3CzdpkTo3uu51VD3tod5rMD4+qZzw0ef0cq2xgXYJ7zv7Q1BchKW+\nB5oV6xq4O5130+AQvrT2g+jZdaig2lAQvpNaZrZ6DvLvM1fvzauxpQI/uq5vfuE0Til1DeR8ZtoG\npftBqZ+NDTQzp8zVzeTRvlVl2cv16q4Tb0pUVlRezbvZGrj03XxvhgkKXrU0mwI/bih18u/2NaIb\nm7bEbtryFM5lZiYbuhQtt2PwHrQLePVcfdjtF5OR3NO7sqherpfzl4pTigWVBBSr1XHluvKiqtrG\nOGOnK6+UyzgbRn7CLWqknOfXoSrXeTatTw+j97Z386oCLdt1puMVtA20mc9BOdWBW/ey2Ge5lMpG\nXo+VERu56xAjTakNHgWW0MqbcY+YouY1liQo2zBBoJTSu6W8Tj6H3etl97IuAMuP39GNlkXE68+N\nbeE8aA28aq4+aEXkdj/LD8Hq8jR1PWfdnlMl5B/tvG7W79dPgpIHLX7fpeSs+53nzuhxa+D90LLY\nNlP+/cXaBH5qmV3cGriaUPVBKRxemseXG12OZSabK1gJk8trSRvvJVeSUgxttRpnSu2pZtwGWPmh\n5UUL4touc8ToeMpq0NPeErd1V5cCG2gFXAiC8QqJVLe/TEUavPYPp5l6td9o55J6bRNL5SgZ79hp\nGch/rxu2H7QWBKPjKcfYEtKyl257bKAFBpLD1pfOVD+0Sq1Ee0YS6UByGP3dncoxS7muaPbupncs\no8dpJcRUB5VqtTo6nkLPrkO28Qriqt2Ntun6k3tF28F70O/Ce861B82CN2w/aFtgwI+9KzfnjEZC\nWHhxvKCDl+q1VJJUvAnMi4YC39EqKHvQIqxzphxQgwxRk+0tcYQbQ7YBnBRhLl6DDQbw5Fd5D5qp\nI6iASM+uQ7azVbczWS/YTQZEKBpbTKty834MA4E3zkFiIDlsrYTcGudoJIRH+1bVnVu8GnBKnyyF\nR/tWeX6N7JEZOTVR0O7W7bVjV9Og7lfQYvoF7z3XH6rVaam4WTnrWu6RmOX3QZGkAAoikRdeHOdu\nVgK6FXSpEdyqfr9MsJmtN0x+vdP5dJql53SPy8117Lxhdb2CJvGOnJrwnM/G1AbljNaPhGcK5be3\nxJWzcrtrTFUoxTTzhlnsoMZ50P5DAVZ8Tygvfu4hz3apKb9e1rO8Ik5nclYNf7eomvCcTeuvsbor\nVEIr5od6ryl4nPel6hNdg/dSESM5da0yVTeSdEbdrEN+DRfN8IZYGII1PvdU27VLetYt4HSfx2li\n51bPdWWgRXfXpsEhdmszAMp306C+zYTddaUTppMbltzd4UbeG3WDWLVJRXtL3Pea0ExwcTKQlJVR\niqfNzbW1aEHcVst17eJmmHJyfGzC2jKxW61FIyHL0IqQYbZzmdGNIuj7z0FA3MLS7Qc65a4ytY1T\n7+qRUxNYv02femsY+uJEAJQ6FxkdT9lqua4MNBUsaG+ZCbDh1TMjMht3t9tVWDqTK9hTlp+zuybF\nrjmMPU6rHrqxOt1EGUaHaeYn5qr7xuh4ynFbJZ3J2Wq5rgw0kDfSNGvW7REy9QdVAvMj7UqFbMy9\nTgzsAkuYPE6FSY6PTaBn1yFtYRmGcYOqHKlhuF/82Wm57gw0cdOWp3j1zFjM9R6kbnzOwy0daimp\nM74UqS3342aYUhCvs3LdT+rSQA8kh3kFwhRxfGxmrzIoKyqq3cuG2hsDyWErAHSumi0w9YOb+0Up\n/dzr0kAztYvdismJUntD+0U0ErK2YbiOtHsoOIw9ZEylcOq6RVXOvEaD11WaFQDbusxMbSDXyNVR\nrpxiu4pCs4HOx/m7s0NXvY1h7ChnzYFSNVxXK+j129g41zpeqj+V81rwapyD4kKvJajmtmqVsnez\n91rLTH3jNmDUzy2oullBc8Q24xdi2T6nWTfV/iaX9djpVEnxEOVuDF/tbBocslYp8uRHLiDDMG6x\n07NYJ1+ur10u6mYFzftRjF/Ie9d2q2MKRBs5NVFSdSLDyLtsuZuVHjGvmVqOlrPmOlMfpDM52xx5\n6gHvVJhIB3W8WtKm72ZV8wa6Z9chTqNgSqYUV7TdClpu8u5m9SxGf1Jq0KbBIe9vrIbZsbHL+p5E\nY6wrCMMwbhCvH9W9QFz4uY17IcPsZjuupg30hu36Em0M44bZ7lOXY3+KZuccyW2PWISIMM38JF3V\nKYxhnBD1HwmHSuobLZ9Ptkkn3tJPIj3vQWezWWzZsgVjY2PIZDLo6elBe3s7+vr60NDQgMWLF2Pr\n1q0wDAP79+/HE088gcbGRvT09ODqq6/2/IFKZSA5zAFhTCAgw6pzg7mJFhX3rRk1tM+8p3dlgdeM\nJzuEZx4AABK/SURBVOlMOchkcxhIDqO9JT7rzApxz7qs7SZ/9KMf4cILL8SOHTvwP//zP7juuutw\n+eWXo7e3F8uXL8fWrVvx9NNPY9myZUgmkzhw4ADS6TTWrl2LD3/4w4hEIrP6YG5h1xYTBMqVIiXX\n77abddcjH7vzSevvct1EGUaEerOXA7f2ybOLe/Xq1di4cSMAYHp6Go2NjXj11VexfPlyAEBXVxcO\nHz6MV155BR0dHQiHw2hqakJbWxuOHTvmdbiS4dUzU078rOTl9loVj+NKeDNs2F4cY0JlPmU4vY0J\nAk6FTQjPBvo973kPYrEYpqam8MUvfhF33HEHpqenredjsRgmJycxNTWF5ubmgsenpqa8DueZgeSw\nUrAMMxvmyk3KBsUZ1Y1O1yKQJ+5MNVFSkNibb76JT3/60/j4xz+ONWvWoKFh5jRTU1OIx+NoampC\nKjUTTZlKpRCP+9vWjUr8sQiZuaScRtU0Z1e+tB6R3ZBcx5wpFYq6Luf5vASOet6DPn36NNavX4+t\nW7fiQx/6EADg8ssvx5EjR3DFFVdgaGgIK1aswNKlS7F7925kMhmk02mMjIxg8eLFjudPJJodj9HB\ne05MELjskvk48daErRuaCo3IxzQYwLQ0wVx4cdzxfEFkNlouJw08uWFKxDTd6a/hXUPupOcLIiE8\n1HuN637ung30N77xDUxOTuLhhx/Gww8/DADo7+/HwMAAstks2tvbsXr1ahiGgVtvvRW33HILpqen\n0dvb6ypAbHx80utbslbODBMEjp0841gmcNrMRxs/+PhLOHriDIC8wPf0rgQAq2RlayKGL639YFVe\n46Vo2Q+qbWLDBIvs+RzOOWxx/eO7pWRFPVORHGCmkuXDf7MS4+OTWPB7MVd6NkwzWA5hr6KuxhsX\nU/u4SZ1qb4kj3BjC62/O5DmLebyikQZQUMrSNIEf7bzOnzdfJsploOnmRpMXGS5ExPgJFRWxg7pV\niXoWr1dRy/3dnUV2S6flmi5UwjB+Q3tKMm6mvSOnJnD0xBmrotDIqQkrl1dsmSiXBQ3WlNpfqKdz\nOpNT1tNX1dhW/R4MUypuAkSp3KeoZ7peZS17qQvPBpphSsAwZoLB3HpwopGQo/Gg/Ei5drRcCate\n9lXFz53O5Ipubroa2xwYxviN0zWmul6BvMZ9y4MOEuzeZiqBKoKa8hjdrmYNI+/yovxcXVQ2laaU\nZ+31tGoWkZsVyCsQVXU1Wq0wjEy5siEoXqS9JW57TpV98nLfqMp2k9w6jqkU5RI0Gd7WRMwyHrp9\najfG5dJL9B1w6on+7k7eg2ZcU66JLrmwWxMxx3POZhFZdQaaV81MJSnnylUu1xkJOwefiJBLrTUR\nw46NXeV7YwFGldNMQTYAippjMEylkPXsJjBUJBoJoTURQ7hR7yqvOgPNMNUMCdgw8oZ27HTKdRqQ\nHOVdizgZXlq50MSGvWnMXEL93yPhvLH1sngkPdvVC6i6Pej+7k6uqsRUPV4L79M1TwapFvtBk3fM\nKdJVDBzTBYkxTKWgFpJerkXyhg0kh221XHUraN5vqm3aW+I4Pla/5Vp1brJIOGQZdFUQWbUzkBxW\nRraqclBptQKwgWaCjer6Je+Zmwl6Va2gVXmQTG0xcmqiKHK3VolGQth//5qC2rymWZzHK3uM5JSr\nakeuoU/7zDoy2fxqpdbd/Ux10ZqIYUnb/IL0q9ZErEjPTlUGRarKQNfaqoFRUwtBgG62YWgVKBsa\n+fOTC43OaZr1kwetSqOi74O71jFBY8fGroJrlrZsRGR3uJ2Wq8ZAs2ubqRWoQw4Z5oHksKvJpzjz\nlhtqVDOya1tX4EGmXrdBmGAh69krdlquCgPNxpmpNuyMB60AKUDEzT4qFUaotQpZPbsOKb8rXUU1\nGfaqMXMN6bln1yFXeqaiRW7aTgY+SIzTKJi5xm3jC6+uebtgOErdIANEbrM9vSsxkBy2zZ2sJpz2\n02ttv52pXdKZHI6eOKPc3opGQshkc1YVsYHksLXifvDxl7TnDLyBroX9SKa6ceNKdVtb1+15TVNv\nnJxyJ6sJdlMztYbqms5kc9rCRHZFhwLt4ubVMxNUKuFqFoUu7lfXiy4GksNswJnAIq+UnQJDReMs\n6rlq86B59czUCl7LAMps2H4QixbMuNE3DQ7VdLlP0+QcZyaYGAYKtFgKVMfb6RyBXUFzYBgTZOT0\nHyfjO9uVoGywTrxV/ZNXJ0+AKsWKYeYa3eTRNIElbeomNqrcZzcT0EAa6Hpx4THBpMFwdleVMnue\nbYlacU+67aLqLubipukNr6CZoNKaiCm3udxOnNOZnKsJaCANNLu2mbnG64qXciHt9qYXLXBOq3B6\nT4aRjxivdve2G+PLEdzMXGKn1eNjE57SHlXBYSOnJtDeEteuuoEAGmh2bTNzTSlFQCLhfJ6y3ax4\n5NSEVfqP/pNLAzpRK0FTbvKXa+WzMtUH1cvWQf3dxXzm9pY49t+/xpOn7PjYRPVGcTNMJahkdzSx\nAMfoeAo7NnZpZ+K691VKShfDMO5xUy87k80XJxG9QZsGhzxNLE0TuGnLU9rnA2Wgr7vrybl+C0wd\n4nWlFo2EtMazv7vTdkVMbSbTmRzSmRw+dqf+mieX9mzfb5CopzQxpnpJZ3Lo7+4sanQhQhXE6L+R\nUxN47eSZksbSEag0q1qqL8zULmKrQ7nSVynYpVwsWhAvGAeoTA52ORENMseXMLXMtKlPqZSrAxKX\nXqLfgw6UgWaYakA2MlQonwyR216vIv3dndiw/WCBsA1jpqCB3AM66HnQmwaHkD2ff79slJlqZCA5\nXFImwaIF6p72FKdCvc/pebvIb8M0g+Mws3P3MUyQERuzq5q0y9Ase140hAW/FytaJYvnpdW5bOh+\ntPO6Mr378kNadvNdMEwtYBjAZZfMR/Z8TjsppaY3cpqhTsuB2oNmmKBTroAymhbT3pXOiNHz1ZoT\nrGpYzzC1iCqwTL5fUNcrt7CBZuoeu6Av+XFZhJReQekW9Lfd+cQ95GnJpU3nkN+TnBNs1+Q9CCxp\nm28ZZnZxM9WIqEU3UDer0fGUdV/Yu3kV2lviBVpOZ3IFmRh2WmYDzdQkXla66UxOGxktP64qnkH7\nxP3dndZeNBUlkYuX0D6U+Bi9Vxqrv7sTrYlYwWRg0YJCkdsFlgSFkVMTbJyZqkSM/yAvEBldUbuq\nxyiqm+jv7izSr6htOy2zgWZqlnJHO0cjISxaYO+upb0lKkqyp3dlwfNiX2cqVCKfk86RzuQKVtxk\nwKORUKADxADg6Anv6SYM4xd2XjIVezevKtAygCItA3k993d3WnqeFy2+59B5KG1SjCtxwlcDPT09\njXvuuQef/OQn0d3djZMnT/o5HMNYkDGbzZ6xbOBJjOSyooAPOwaSwwXBYzQrB/Iz6x0bu6xztrfE\nC56nz5HO5Ar2oLmJBFNviKvUUjTdmohh7+ZVRZ4r2Y1tGMCjfauU5xAnzrRylvW8//41Wi0Daj3b\n4WsU97//+7/jmWeewQMPPICXX34ZjzzyCP7hH/5BezxHcTPlRowidrsnSpHV5K4eHU9ZxtkNlG5F\nr6fxVKJNJJoxPj6pPIcY2S0GWvV3dyKRaHb1XuYK1jJTTlTRz+0t+paPZMQj4RAWXhzHl9Z+0HqO\ngrScJteANy0D7vVME4LWRAwP9V6jHd9XA71t2zYsXboUf/mXfwkA6OrqwtCQvjk1i5rxArmKxAtf\nLgYgGzaiZ9chZLK5oqAvEoxKZKUiilxGJ2in17KBZmodueCHPMEmI0kGN5NVx5IsaZtvGWg3RtYO\nOy0D7vQsT7ztDLSvhUqmpqbQ1NRk/TsUCmF6ehoNDbz1XS+QwXQq3kFikVescoEO+XhCFo6TkNzM\nnsuF15tAuV7LMHONrqqWfAxQbIzlVSsws5IlBpLDlpbdtDAVqUTqotN9aHQ8ZVt0yPcV9LJly/DR\nj34UALBy5UocOqTPAeNZd/khVwqtFsWVJv3ysnGkYx7tW6WspkN7oGOnUzibzllGGECBWABojagX\n1/FAchjhxpBVmcpvo+U0Cw7CWLyCZioN3Q9URlA2xKoCO+L+b4MBnE3PrCJl9684gdbdS1QrYdF4\n0/1F9oiJk34vq2ivLm75eNV3Qe9DV6jE1xV0R0cHnnnmGXz0ox/Fz3/+c1x22WV+DhcoGoyZHNcG\nozCUnkq77b9/DTYNFrr8vUTn0mvLFdGrcrXYuV9mcy4v5y31PcyGShrAoBtbZnZQnuull8zHibcm\nLMPkhnnRENoumun/TZ2P2i6KF0TKU0/h106eQTQSwv7711jP3bTlKaQzOVx6yXzrPJsGh4oi7Ze0\nzS+4l6juLw/1XoNNg0N47eQZ6/42Lxqy7mUn3pooeL+684jjhxtD+OcH1lifTXzvNKZMuDFU8Ddp\nSHevEDW28OJ4wdhu9acbUzeWfLzIwovzBtsp28HXFbRpmrj33ntx7NgxAMADDzyAhQsX+jUcwzAM\nw9QMgarFzTAMwzBMHo7WYhiGYZgAwgaaYRiGYQIIG2iGYRiGCSBsoBmGYRgmgPiaZuWW6elp3Hvv\nvXjttdcQDocxMDCASy65pKxjvPzyy/jqV7+KZDKJEydOoK+vDw0NDVi8eDG2bt0Ko0yNfrPZLLZs\n2YKxsTFkMhn09PSgvb3dl/FyuRy+8pWv4I033oBhGLjvvvsQiUR8+2wA8Pbbb+MTn/gEvv3tb6Oh\nocG3sf7qr/7KKnLz+7//+/jc5z7n21iPPPIInnnmGWSzWXzqU59CR0eHL2P94Ac/wIEDBwAA6XQa\nR48exXe/+10MDAz49ntVmkpoGaiMniupZaDyeq6UlgHWc8mYAeDf/u3fzL6+PtM0TfPnP/+52dPT\nU9bzf/Ob3zTXrFlj3nzzzaZpmubnPvc588iRI6ZpmuY999xj/uQnPynbWN///vfN+++/3zRN0/zv\n//5vc+XKleZtt93my3g/+clPzC1btpimaZrPP/+8edttt/k2lmmaZiaTMT//+c+bf/EXf2GOjIz4\n9j2eO3fO/PjHP17wmF9j/ed//qf5uc99zjRN00ylUubXvvY1X79D4r777jP3799fkbEqid9aNs3K\n6bmSWjbNyuq5Ulo2TdbzbAiEi/vFF1/EVVddBQBYtmwZfvnLX5b1/G1tbfj6178O892MsldffRXL\nly8HkK8Pfvjw4bKNtXr1amzcuBFAfjXR2Njo23h//ud/jr/9278FAJw6dQrvfe978atf/cq3z/bg\ngw9i7dq1SCQSAPz7Ho8ePYqzZ89iw4YN+PSnP42f//znvo313HPP4bLLLsPnP/953HbbbVi1apWv\n3yEAvPLKK/iv//ov3Hjjjb6PVWn81jJQOT1XUstAZfVcKS0DrOfZEAgDravZXS6uvfZahEIzlVxM\nIfX7Pe95DyYny1fW8T3veQ9isRimpqbwxS9+EXfccUfBZyn3eKFQCH19fRgYGMDHPvYx3z7bgQMH\ncOGFF+LKK68EkP8O/Rpr3rx52LBhA/bu3Yv77rsPd911V8Hz5RzrnXfewS9/+UsMDg7ivvvuw513\n3unr9QHkXXC33347AH+vxbnAby0DldNzpbUMVEbPldQywHqeDYHYg25qakIqNVPv2e+GGuK5U6kU\n4vG4zdHeefPNN3H77bdj3bp1WLNmDXbs2OHreNu2bcPp06dx4403IpPJ+DLWgQMHYBgGDh8+jKNH\nj6Kvrw9nzsyUqSvnWB/4wAfQ1tZm/f2+970Pv/71r30Za/78+Whvb0djYyMWLlyIaDSK3/3ud76M\nBQATExN44403cMUVVwDw/1qsNJXWMuDvd1hpLQP+67mSWgZYz7MhECvojo4Oqw1lJWp2X3755Thy\n5AgAYGhoCJ2d5Wu+cPr0aaxfvx6bNm3CJz7xCV/H+5d/+Rc88sgjAIALLrgADQ0N+KM/+iNfxvrO\nd76DZDKJZDKJJUuWYPv27bjyyit9GevAgQPYtm0bAOC3v/0tUqkU/uzP/syXsf7kT/4EP/3pT62x\nzp07hw996EO+XR8vvPACPvShD1n/9vNanAsqrWXAv++wkloGKqfnSmoZYD3PhkCsoD/ykY/gueee\nwyc/+UkA+ZrdfkDRdH19fbj77ruRzWbR3t6O1atXl22Mb3zjG5icnMTDDz+Mhx9+GADQ39+PgYGB\nso+3evVq9PX14VOf+hTOnz+P/v5+LFq0yLfPJmIYhm/f4w033IAvf/nLWLduHYD89fC+973Pl7Gu\nvvpqvPDCC7jhhhswPT2NrVu3oqWlxbfv8I033iiIavbzWpwLKqVlwH89V1LLwNzp2U8tA6zn2cC1\nuBmGYRgmgATCxc0wDMMwTCFsoBmGYRgmgLCBZhiGYZgAwgaaYRiGYQIIG2iGYRiGCSBsoBmGYRgm\ngLCBZhiGYZgAwgaaYRiGYQLI/wfJZujdkeFrywAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113c25e50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#### Example of raw, QC, and proc at gin flats (known shading)\n",
    "fig,axes = plt.subplots(1,2,sharey=True,figsize=(8,4))\n",
    "axes[0].plot(grobs.el.sel(station='gin').values,grobs.SWdwn_proc.sel(station='gin').values,'.')\n",
    "axes[1].plot(grobs.el.sel(station='gin').values,grobs.SWdwn_raw.sel(station='gin').values,'.')\n",
    "for ax in axes:\n",
    "    ax.set_ylim(0,1200)\n",
    "    ax.set_xlim(0,75)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xray.DataArray 'SWdwn_proc' (time: 2976, station: 18)>\n",
      "array([[False,  True, False, ..., False, False, False],\n",
      "       [False,  True, False, ..., False, False, False],\n",
      "       [False,  True, False, ..., False, False, False],\n",
      "       ..., \n",
      "       [False, False,  True, ..., False,  True, False],\n",
      "       [False, False, False, ..., False,  True, False],\n",
      "       [False, False, False, ..., False,  True, False]], dtype=bool)\n",
      "Coordinates:\n",
      "    network    (station) |S4 'CDWR' 'CDWR' 'CDWR' 'CDWR' 'CDWR' 'CDWR' ...\n",
      "    longitude  (station) float64 118.3 119.3 119.1 119.9 118.4 119.8 118.3 ...\n",
      "  * station    (station) object 'bee' 'dan' 'dep' 'for' 'fsl' 'gin' 'gld' ...\n",
      "    elev       (station) float64 2.768e+03 2.987e+03 2.307e+03 1.72e+03 ...\n",
      "  * time       (time) datetime64[ns] 2004-08-27 2004-08-28 2004-08-29 ...\n",
      "    latitude   (station) float64 37.53 37.9 37.63 37.77 37.48 37.77 37.52 ...\n",
      "    Grouping   (station) object 'east of crest' 'west of crest' ...\n"
     ]
    }
   ],
   "source": [
    "#### Further refinement:\n",
    "# Remove the 'smj' station -- excessively bad\n",
    "grobs = grobs.drop(['smj'],dim='station')\n",
    "\n",
    "# - remove months with > 10% gaps\n",
    "# - set 0 daily means to nan\n",
    "## Resample, exclude months w/ > 10% gaps\n",
    "grobs_day = grobs.resample(freq='D', dim='time', how='mean')\n",
    "grobs_day.SWdwn_proc.values[grobs_day.SWdwn_proc.values == 0] = np.nan\n",
    "numdays_permonth = (~np.isnan(grobs_day['SWdwn_proc']))\n",
    "print(numdays_permonth)\n",
    "# .resample(freq='M', dim='time', how='mean')\n",
    "# grobs_month = grobs_day.resample(freq='M', dim='time', how='mean')\n",
    "# ds.SWdwn.values[numdays_permonth.values < .9] = np.nan"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}