precip_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#precip analysis\n",
    "import glob\n",
    "import xarray as xr\n",
    "import os\n",
    "import numpy as np\n",
    "import warnings\n",
    "import sys\n",
    "import pickle\n",
    "import matplotlib.pyplot as plt\n",
    "sys.path.insert(0, '/home/disk/p/jkcm/Code')\n",
    "from Lagrangian_CSET.CSET_data_classes import CSET_Flight\n",
    "from Lagrangian_CSET import utils\n",
    "from Lagrangian_CSET import met_utils as mu\n",
    "from Lagrangian_CSET.LoopTimer import LoopTimer\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_precip_from_leg(flight, start_time, end_time):\n",
    "    precip_data_location = r'/home/disk/eos4/jkcm/Data/CSET/precip_retrievals'\n",
    "    precip_file = glob.glob(os.path.join(precip_data_location,\n",
    "                                  'RF{:02}_*.cdf'.format(flight)))\n",
    "    if len(precip_file) == 0:\n",
    "        return([np.nan])\n",
    "    with xr.open_dataset(precip_file[0]) as data:\n",
    "        data['time'].values = data.time.values.astype('datetime64')\n",
    "        precip_data = data.sel(time=slice(start_time, end_time))\n",
    "        prec = precip_data['lwf']\n",
    "        col_max = prec.max(dim='height', skipna=True)   \n",
    "    return col_max.values*3.6*24  # in mm/day, from g/m2/s originally\n",
    "    #open the right precip file\n",
    "    #chop out the right segment\n",
    "    #get the colmax precips\n",
    "    #return that whole array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = utils.read_CSET_Lookup_Table(rf_num='all', sequences='all', legs=['b'], \n",
    "                                 variables=['Date', 'ST', 'ET', 'uhsas_conc1_M', 'Lat_S', 'Lat_E', 'Lon_S', 'Lon_E'])\n",
    "start_dates = [utils.CSET_date_from_table(x['Date']['values'][i], x['ST']['values'][i])\n",
    "               for i in range(len(x['Date']['values']))]\n",
    "end_dates = [utils.CSET_date_from_table(x['Date']['values'][i], x['ET']['values'][i])\n",
    "             for i in range(len(x['Date']['values']))]\n",
    "lats = (x['Lat_S']['values'] + x['Lat_E']['values'])/2\n",
    "lons = (x['Lon_S']['values'] + x['Lon_E']['values'])/2\n",
    "times = [utils.CSET_date_from_table(x['Date']['values'][i], (x['ST']['values'][i]+x['ET']['values'][i])/2)\n",
    "               for i in range(len(x['Date']['values']))]\n",
    "Nds = x['uhsas_conc1_M']['values']\n",
    "Nds[Nds=='-'] = np.nan\n",
    "Nds = Nds.astype(float)\n",
    "locs = np.array([lats, lons, times]).T\n",
    "flights = x['rf']\n",
    "all_precips = []\n",
    "mean_precips = []\n",
    "for (start_time, end_time, rf) in zip(start_dates, end_dates, flights):\n",
    "    precips = get_precip_from_leg(rf, start_time, end_time)\n",
    "    all_precips.append(precips)\n",
    "    mean_precips.append(np.nanmean(precips))\n",
    "    \n",
    "    \n",
    "locs = np.array([lats, lons, times, all_precips]).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_loc = r'/home/disk/eos4/jkcm/Data/CSET/Python/'\n",
    "with open(os.path.join(data_loc, 'nep_scenes.pickle'), 'rb') as f:\n",
    "    nep_scenes = pickle.load(f)\n",
    "\n",
    "lats_nep = nep_scenes['lats_nep']\n",
    "lons_nep = nep_scenes['lons_nep']\n",
    "int_local_days_nep = nep_scenes['int_local_days_nep']\n",
    "scenes_nep = nep_scenes['scenes_nep']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {},
   "outputs": [],
   "source": [
    "min_dists = []\n",
    "best_scenes = []\n",
    "failed_locs = []\n",
    "precip_by_scene = {-1: [],\n",
    "                 0: [],\n",
    "                 1: [],\n",
    "                 2: [],\n",
    "                 3: [],\n",
    "                 4: []}\n",
    "for i, loc in enumerate(locs):\n",
    "    lat, lon = loc[0], loc[1]\n",
    "    dec_day = loc[2].timetuple().tm_yday + (loc[2].hour + loc[2].minute/60)/24\n",
    "    local_day = dec_day + lon/360\n",
    "    matching_days = int_local_days_nep == int(local_day)\n",
    "    dists = [np.sqrt((lat-lat2)**2 + (lon-lon2)**2) for lat2, lon2 in zip(lats_nep[matching_days], lons_nep[matching_days])]\n",
    "    if dists:\n",
    "        best_i = np.argmin(dists)\n",
    "        if dists[best_i] > 2:\n",
    "            failed_locs.append(loc)\n",
    "        else:\n",
    "            min_dists.append(dists[best_i])\n",
    "            assert len(scenes_nep[matching_days]) == len(dists)\n",
    "            scene = scenes_nep[matching_days][best_i]\n",
    "            best_scenes.append(scene)\n",
    "            precip_by_scene[scene].append(loc[3])\n",
    "    else:\n",
    "        failed_locs.append(loc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "type 0 (open): total count is 4\n",
      "type 1 (closed): total count is 4\n",
      "type 2 (homogen): total count is 0\n",
      "type 3 (disorg): total count is 19\n",
      "type 4 (too clear): total count is 0\n",
      "type -1 (NC): total count is 33\n",
      "failed locs (too far): 8\n",
      "0.20154757192328254 0.2576920221147485 0.004395898967535648\n",
      "3.4220401135117426 7.21162939659615 0.07864288016311183\n",
      "4.036152806837552 12.6696389531464 0.014658555050205357\n",
      "4.98424976257631 14.343632712552168 0.04369541852303674\n",
      "failed locs (too far): 8\n"
     ]
    }
   ],
   "source": [
    "scene_dict = {-1: 'NC',\n",
    "              0: 'open',\n",
    "              1: 'closed',\n",
    "              2: 'homogen',\n",
    "              3: 'disorg',\n",
    "              4: 'too clear'}\n",
    "\n",
    "\n",
    "\n",
    "for s, locs in precip_by_scene.items():\n",
    "    print(\"type {} ({}): total count is {}\".format(s, scene_dict[s], len(locs)))\n",
    "print(\"failed locs (too far): {}\".format(len(failed_locs)))\n",
    "\n",
    "\n",
    "\n",
    "for s, locs in precip_by_scene.items():\n",
    "    if len(locs)>0:\n",
    "        mean = np.nanmean(np.concatenate(locs))\n",
    "        nin = np.nanpercentile(np.concatenate(locs), 90)\n",
    "        ten = np.nanpercentile(np.concatenate(locs), 10)\n",
    "        print(mean, nin, ten)\n",
    "        \n",
    "#     print(\"type {} ({}): total count is {}\".format(s, scene_dict[s], num))\n",
    "print(\"failed locs (too far): {}\".format(len(failed_locs)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 209,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAEYCAYAAAD4czk4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VdW9///Xh9EwhwCiTAGjItGAgBpAkVGgKIPXAWnR67WK0iIOcC9aNAFsiy0g/enVIsrsTylOYKmgFs1XsVQggBK+AjIpIigmVBkLyfr+sU+OOclJzgk5kE3yfj4e5wFZe332XnuR5MPae+21zTmHiIiIn1Qp7waIiIgUpuQkIiK+o+QkIiK+o+QkIiK+o+QkIiK+o+QkIiK+o+QkIiK+o+QkIiK+o+QkIiK+U628GxArjRo1comJieXdDBERKWTdunUHnHONSxNTYZJTYmIia9euLe9miIhIIWa2u7QxuqwnIiK+U2FGTiKxkjh+WZnid00ZGKOWiFReGjmJiIjvKDmJiIjvKDmJSLlJT0/HzMq7GWH16NGDHj16RFV39OjR3HDDDae3QYV89dVX3HTTTdSvX5969epx44038uWXX0Yd/7e//Y3u3btTp04d6tWrR+fOnVm5cmVwe48ePTCzsJ/+/fsH6z311FOkpKSQl5cX0/OLaXIysxZm9qqZ/cvMfjCz182sZZSxrphPh1i2UUQklrZv387MmTNJS0s7Y8c8cuQIvXr14vPPP2fevHksWLCAbdu20bNnTw4fPhwxfubMmQwePJhOnTrxxhtvsHjxYm6++WaOHDkSrPPss8/yj3/8I+Qzffp0AAYNGhSsd++99/Ltt98yb968mJ5jzCZEmFktYCVwHLgDcMATwPtmluKci9xjMBeYWahsa6zaKCISazNmzKB9+/Z07tz5jB1z1qxZ7Nixgy1btpCUlARASkoKF154ITNnzuShhx4qNnbXrl088MAD/PGPf+SBBx4Ilvfr1y+kXrt27cIet0aNGgwbNixYFhcXx+23387UqVO58847y3pqQbEcOd0NtAGGOOfedM4tAQYBrYCRUe7ja+fc6kKfI5HDRMSvNm7cyNChQ0lISCAuLo6LL76Y3//+9yXGLF++nC5duhAXF0f9+vUZMmQIW7ZsCamzdetWhg4dSpMmTTjnnHNo2bIlN998MydPngw59qBBg4iPjycuLo5u3brx4YcfFjneK6+8Qtu2balZsybJycm88cYbUZ3b8ePHWbhwIcOHDw8pz8vLo27dukyaNCmkPCcnBzMr8yhj6dKlpKamBhMTQOvWrenWrRtLliwpMXb27NlUqVKFe++9t1THPHr0KIsXL+aGG26gYcOGIduGDRvG5s2b+fjjj0u1z5LEMjkNAlY7577IL3DO7QRWAYNjeBwROUt88skndOnShe3bt/PUU0+xbNkyHnroIfbs2VNszPLlyxk4cCB16tRh0aJFPPfcc2zatImrr76ar7/+Oljv+uuv5+uvv+a5555jxYoVTJkyhZo1awbvfWRmZtK1a1eys7OZNWsWr732GgkJCfTp04d169YF9/Pee+8xfPhwLrzwQl5//XXGjRvHmDFjiiTDcFavXs3Bgwe55pprQsq3bt3KoUOHuPzyy0PK169fDxAsd85x8uTJiJ/c3NyQ/WRlZXHppZcWaU9ycjKbN28usc0fffQRbdu25ZVXXuGCCy6gWrVqJCUl8b//+78lxr3++uv8+OOP3HHHHUW2dejQgXr16rF8+fIS91EasXzOKRkIl7KzgJuj3Md9ZjYOyAVWA2nOuaL/zRGRs8LYsWNJSEhg9erV1KpVC4BevXqVGDNhwgTatGnD22+/TbVq3q+oLl26cNFFFzFt2jSmT5/OgQMH2LZtG0uWLAm5/1FwBDNu3DhatmzJypUrqVGjBuBdurr00kuZPHkyb775JgBpaWm0bduWJUuWUKWK9//1Sy65hNTUVC6++OIS27p69WrMjJSUlJDyzMxMADp27BhSvn79emrWrMkll1wCQEZGBj179izxGADXXnstH3zwQfDr7Oxs4uPji9Rr2LAhOTk5Je5r79697N27l3HjxvG73/2OCy64gMWLF/PrX/+akydPMmbMmLBx8+fPp0mTJgwYMKDItipVqpCSksLq1asjnku0YpmcGgLheiUbKNqLRS0E/grsxbsUOA5YaWZ9nXMfhAsws3uAewBatoxq3oWInCFHjhxh1apVjBs3LpiYIjl8+DCZmZk8+uijwcQEP12yysjIACAhIYE2bdowfvx49u/fT48ePbjwwguD9Y8ePUpGRgaPPvooVapUCbnU16dPH1566SUAcnNzWbNmDePHjw8mJoCrrrqKaNbq3Lt3L/Xq1Qsmv3zr1q2jSZMmNGvWLKQ8MzOT5ORkqlevDkCnTp1Ys2ZNxOPUrVu3SFm4WY7OuYj7ysvL48cff2Tu3LnceOONgPcfhl27dvH73/+e+++/v8i+9+7dy3vvvceYMWNC/l0Katy4MVu3xm6KQKxXiAjXM1HNE3XOjSjw5YdmtgTYhDep4upiYp4Hngfo3Llz5H8VETljcnJyyMvLo3nz5qWKcc5x3nnnFdnWtGlTdu/2lmgzM959913S09N55JFH+P7772ndujXjxo3jvvvuIzs7m9zcXCZPnszkyZPDHisvL48DBw5w4sQJzj333CLbw5UVduzYMWrWrFmkPDMzs8ioCbyRU9euXYNf16lThw4dIk9ILpws4uPjyc7OLlIvJycn7IiqoISEBLZt20bfvn1Dyq+77jqWL1/ON998w/nnnx+ybeHCheTl5YW9pJcvLi6Oo0ePRjqVqMUyOeXgjZ4Kiyf8iKpEzrkfzWwZcFdZGyZyJpV1+SOoGEsgxcfHU6VKlZD7RNHEmBn79u0rsm3fvn0kJCQEv27Tpg3z58/HOcfGjRt55plnGDVqFImJiXTv3p0qVarwq1/9ittvvz3ssapUqUKjRo2oXr06+/fvL7J9//79tGrVqsT2JiQkFLmM5pxjw4YN3HfffSHl3377LVu2bOFXv/pVsOxUL+slJyeTlZVVpN7mzZvDzrIrKDk5Oezlt/xRV8ERZL758+fTvn172rdvX+x+s7OzadSoUYnHLo1YJqcsvPtOhbUDSr5DVzwj/GhMpEKrCOv71apVi6uvvpqFCxfy+OOPExcXFzGmdu3adOrUicWLF5Oenk7VqlUB2L17Nx9//DGjR48uEmNmdOjQgenTp/Piiy+yadMmBgwYwDXXXMPGjRvp2LFj2F+4AFWrVuWKK67g1VdfJT09PVjvn//8J7t27YqYnNq2bcuJEyfYs2dPcIS4fft2Dh48GGx7vqeffpq8vLyQkdKpXtYbNGgQY8eOZceOHbRp0wbwpoivWrWKKVOmlLivoUOH8uKLL7JixQpuuummYPmKFSto3rw5TZs2Dam/du1asrKygs84FWfnzp1ceeWVEc8lWrFMTkuBqWbWxjm3A8DMEoFuwPjS7szM6gEDgX/GsI0icgZNnTqVa6+9li5duvDwww/TvHlzduzYwYYNG3j66afDxkyePJmBAwdy/fXXM2rUKA4dOkRaWhr169fn4YcfBuDTTz9lzJgx3HrrrSQlJZGbm8vcuXOpVq1acMLF9OnT6d69O/369eOuu+7ivPPO48CBA2RmZpKbmxv8JT5x4kSuu+46hgwZwsiRI/nuu+9IS0sr8ks6nO7duwPerMT85JQ/GeKFF16gRYsWNGnShHfeeSc4fXzt2rV07NiRuLg46tate0rPR919990888wzDB48mCeeeAIz47HHHqNFixaMHPnTkzsZGRn07t2b2bNnB0eQP/vZz+jZsycjR47kwIEDtGnThldffZV33nmHOXPmFDnW/PnzqVatWpHp8gUdPHiQrVu3Mnbs2FKfS3FimZxmAb8GlpjZBLwRz2TgKwo8WGtmrYDtwCTn3KRA2VjgYuB9fpoQMRZoCvw8hm2USiAWl9UkNq644gpWrVrF448/zujRozl+/DitWrUq8WHN/v37s2zZMiZOnMgtt9xCjRo16NGjB3/4wx+C90KaNm1Ky5YtmT59Onv27OGcc87hsssu469//SudOnUCvJlya9asYeLEidx///3861//onHjxnTs2DHkGZ/8CRLp6enceOONJCUlMWPGDP70pz9FPL/ExESuvPJK3nrrreDkgszMTBo2bMiUKVMYP348hw4dYsiQISxevJjbbruNRYsWFTsjLlq1a9dm5cqVPPjgg4wYMQLnHL1792bGjBnUqVMnWM85R25ubsjSQmbGm2++ySOPPEJaWho5OTm0bduWl156qUgCOnHiBC+//DL9+/cv8R7csmXLqFGjBkOHDi3TeRVk0czuiHpn3lJFTwF98S7J/R14wDm3q0CdRGAnMNE5lx4ouwFvdHUxUB/4Ae/5qCecc59Ec+zOnTs7vWyw/PnhcpSSkz8u61UWc+fOZcyYMXzzzTfUqlUrONHg3XffLeeWnTkDBgygUaNGLFiwIOx2M1vnnCvVEDGma+s55750zv2Hc66ec66uc25IwcQUqLPLOWf5iSlQ9pZzrptzrpFzrrpzLsE5NyjaxCQiUl5GjBhBs2bNePbZZwFvRl7+6K0y2LBhA++//37M1xbUywZFKiDNGDxzqlatyuzZs8nMzGT37t18//33lSo57du3jzlz5oQspRQLSk7iK7okJ2ej1NRUUlNTgegehK1ICr4+I5b0PicREfEdJScREfEdXdaTELqsJiJ+oJGTiIj4jpKTiIj4jpKTiIj4jpKTiIj4jpKTiIj4jpKTiIj4jqaSi0hYfljEVyovjZxERMR3lJxERMR3lJxERMR3dM9JRE4L3bOSstDISUREfEfJSUREfEfJSUREfEf3nCoQve5CRCoKjZxERMR3lJxERMR3dFlPRHwpFpepNR397KWRk4iI+I6Sk4iI+I4u64mIFEOrXJQfjZxERMR3zDlX3m2Iic6dO7u1a9eWdzPKRM8piUhhFWH0ZWbrnHOdSxOjkZOIiPiOkpOIiPiOkpOIiPhOTJOTmbUws1fN7F9m9oOZvW5mLaOMPcfM/mhm35jZUTP7h5l1j2X7RETk7BCzqeRmVgtYCRwH7gAc8ATwvpmlOOcOR9jFi8BAYBywA/gVsMLMujjnNsSqnaeLJjOIiMROLJ9zuhtoA1zsnPsCwMw+BbYBI4HpxQWaWXtgOPBfzrk5gbIMIAuYBAyKYTtFRM4alfVZq1he1hsErM5PTADOuZ3AKmBwFLEngEUFYk8CrwD9zKxmDNspIiI+F8uRUzKwJEx5FnBzFLE7nXNHwsTWAJICfz9tdFlORCqis3UB3Vgmp4ZATpjybCC+DLH524sws3uAewJfHjKzLVG0sziNgANliBeP+jE21I+xoX6MAXuyzP3YqrQBsV5bL9xyExZFnJ1KrHPueeD5KPYfuQFma0v7BLMUpX6MDfVjbKgfY6M8+jGW95xyCD/CiSf8qKig7BJi87eLiEglEcvklIV376iwdsDmKGJbB6ajF479N/BF0RAREamoYpmclgKpZtYmv8DMEoFugW2RYqtTYOKEmVUDbgXecc4dj2E7ixOTy4OifowR9WNsqB9j44z3Y8xWJTez2sBG4CgwAe8e0mSgLpDinDsUqNcK2A5Mcs5NKhD/CtAP7yHcncB9wPVAV+dcZkwaKSIiZ4WYjZwCK0D0ArYCC4CX8JJMr/zEFGBA1TDHvhOYg7eqxDKgBdBfiUlEpPKpMO9zEhGRikOrkouIiO8oOYmIiO8oOYmIiO8oOYmIiO8oOYmIiO8oOYmIiO8oOYmIiO8oOYmIiO8oOYmIiO8oOYmIiO8oOYmIiO8oOYmIiO8oOYmIiO8oOYmIiO8oOYmIiO9UK+8GxEqjRo1cYmJieTdDREQKWbdu3QHnXOPSxFSY5JSYmMjatWvLuxkiIlKIme0ubYwu64mIiO9UmJFTRZA4fll5N4FdUwaWdxNERDRyEhER/1FyEhER31FyEhEpo9GjR3PDDTecseN99dVX3HTTTdSvX5969epx44038uWXX57W+P79+2NmTJgwIaT8qaeeIiUlhby8vFM6l+IoOYmIlMH27duZOXMmaWlpZ+R4R44coVevXnz++efMmzePBQsWsG3bNnr27Mnhw4dPS/zLL7/Mxo0bw2679957+fbbb5k3b16ZzqswJScRkTKYMWMG7du3p3PnzmfkeLNmzWLHjh28+eabDBkyhMGDB7N06VJ2797NzJkzYx5/8OBBHnzwQaZPnx52f3Fxcdx+++1MnTq1zOdWkJKTiPjG8uXL6dKlC3FxcdSvX58hQ4awZcuW4Pb09HTMjM8++4yePXtSq1YtzjvvPB5//PEil5U2btzIoEGDiI+PJy4ujm7duvHhhx+G1Mnf37Zt2xg4cCB16tShVatWTJo0KarLVMePH2fhwoUMHz48pDwvL4+6desyadKkkPKcnBzMrEyjjKVLl5KamkpSUlKwrHXr1nTr1o0lS5bEPP6///u/SU5O5rbbbit2n8OGDWPz5s18/PHHpTyb4ik5iYgvLF++PJggFi1axHPPPcemTZu4+uqr+frrr0PqDhkyhD59+vDmm28yfPhwJk+eHJIIMjMz6dq1K9nZ2cyaNYvXXnuNhIQE+vTpw7p164oce+jQofTq1Ss4mkhLS4sqgaxevZqDBw9yzTXXhJRv3bqVQ4cOcfnll4eUr1+/HiBY7pzj5MmTET+5ubnBfWRlZXHppZcWaUtycjKbN2+O2ObSxH/00UfMnz+fZ599tsR9dujQgXr16rF8+fKIx4+WnnMSEV+YMGECbdq04e2336ZaNe9XU5cuXbjooouYNm1ayGWlu+++m/HjxwNw3XXX8cMPPzBt2jQeeOABGjRowLhx42jZsiUrV66kRo0aAPTr149LL72UyZMn8+abb4Yc++GHH+bOO+8EoE+fPqxcuZKXX345WFac1atXY2akpKSElGdmZgLQsWPHkPL169dTs2ZNLrnkEgAyMjLo2bNnxL659tpr+eCDDwDIzs4mPj6+SJ2GDRuSk5MTcV/Rxp84cYKRI0cyduxYLr744hL3WaVKFVJSUli9enXE40dLyUlEyt3hw4fJzMzk0UcfDSYm+OlyU0ZGRkj9W265JeTrYcOG8cILL7Bp0yY6depERkYGjz76KFWqVOHkyZPBen369OGll14qcvyBA0MfPr/00kuDo5yS7N27l3r16gUTYL5169bRpEkTmjVrFlKemZlJcnIy1atXB6BTp06sWbMm4nHq1q0b8rWZFanjnIu4n9LEP/nkkxw9epTf/OY3Ue2zcePGbN26Neo2RKLkJCLlLicnB+cc5513XpFtTZs2Zffu0KXZzj333LBff/3117Ru3Zrc3FwmT57M5MmTwx4vLy+PKlV+uqvRsGHDkO01a9bk2LFjEdt97NgxatasWaQ8MzOzyKgJvJFT165dg1/XqVOHDh06RDxOwWQSHx9PdnZ2kTo5OTlhR0SFRRP/5Zdf8tvf/pYXXniB48ePc/z48WC948ePc/DgQerWrUvVqlWD5XFxcRw9ejTi8aOle04iUu7i4+MxM/bt21dk2759+0hISAgp279/f9ivmzVrRoMGDahSpQqjR49mzZo1YT8FE1NZJCQkFLmU5pxjw4YNRe43ffvtt2zZsiWkPCMjg+rVq0f89O7dOxiTnJxMVlZWkbZs3ryZdu3aRWxzNPE7duzg2LFj/OIXvyA+Pj74AZg6dSrx8fF89tlnIfHZ2dk0atQo4vGjpZGTiJS72rVr06lTJxYvXkx6enrwf+S7d+/m448/ZvTo0SH1//KXvwTvOQG88sor1KlTh0svvZTatWtzzTXXsHHjRjp27BizRBRO27ZtOXHiBHv27KF58+aA99zTwYMHQ0YVAE8//TR5eXkhI6VTuaw3aNAgxo4dy44dO2jTpg0Au3btYtWqVUyZMiXivqKJ79ChA++//36R2J49e/KLX/yCu+66K2S2H8DOnTu58sorIx4/WkpOIuILkydPZuDAgVx//fWMGjWKQ4cOkZaWRv369Xn44YdD6s6aNYu8vDyuuOIKVqxYwQsvvEB6ejoNGjQAYPr06XTv3p1+/fpx1113cd5553HgwAEyMzPJzc2N6pd4NLp37w7AJ598EkxO+ZMhXnjhBVq0aEGTJk145513grP/1q5dS8eOHYmLi6Nu3bqlfj7q7rvv5plnnmHw4ME88cQTmBmPPfYYLVq0YOTIkSF1MzIy6N27N7Nnz+b222+POr5Bgwb06NEj7PFbtWpVZNvBgwfZunUrY8eOLdW5lESX9UTEF/r378+yZcs4ePAgt9xyC/feey+XXHIJH330Eeeff35I3SVLlvDuu+8yaNAgFi5cyIQJE3jssceC2zt27MiaNWtISEjg/vvv57rrrmPMmDF89tlnwYQSC4mJiVx55ZW89dZbwbLMzEwaNmzIlClTSEtLY8SIEfz4448sXryYevXqsWjRIuLi4k75mLVr12blypVcdNFFjBgxgp///Oe0bt2alStXUqdOnZC6zjlyc3NDntkqTXy0li1bRo0aNRg6dOgpn1dhVpoZHn7WuXNnd7a/bFCvzBApWXp6OhMnTuTEiRMhs/rK09y5cxkzZgzffPMNtWrVom/fvgC8++675dyyM2fAgAE0atSIBQsWhN1uZuucc6UaImrkJCJSBiNGjKBZs2bBB1XXr19Pp06dyrlVZ86GDRt4//33Y762oD/+6+EDZR21VJQRh/pBpHSqVq3K7NmzyczMZPfu3Xz//feVKjnt27ePOXPmFJkgUVZKTiJy1khPTyc9Pb28m1FEamoqqampQOkehq0I+vfvf1r2G9PLembWwsxeNbN/mdkPZva6mbWMMtYV84n8hJqIiFQoMRs5mVktYCVwHLgDcMATwPtmluKci/yiEZgLFF6zPXbrYZxGfpjMICJSUcTyst7dQBvgYufcFwBm9imwDRgJhH8ZSKivnXOxWzlQRETOSrG8rDcIWJ2fmACcczuBVcDgGB5HREQquFgmp2RgU5jyLCDygk+e+8zsuJkdMbOVZnZN5BAREaloYpmcGgLhXiaSDUReKhcWAqOAPsA9QAKw0sx6FBdgZveY2VozW/vdd9+VvsUiIuJLsZ5KHm4OZdEXh4QLdG5EgS8/NLMleCOxJ4Cri4l5HngevBUiStdUOR30nJSIxEIsR045eKOnwuIJP6IqkXPuR2AZcEUZ2yUiImeZWCanLLz7ToW1AyK/2D48I/xoTEREKrBYJqelQKqZtckvMLNEoFtgW6mYWT1gIPDPGLVPRETOErFMTrOAXcASMxtsZoOAJcBXFHiw1sxamdlJM3u8QNlYM5tlZsPNrIeZ3YE3Bb0pMCGGbRQRkbNAzCZEOOcOm1kv4ClgAd4lub8DDzjnDhWoakBVQhPjFmBo4FMf+AEvOd3lnPskVm0UEZGzQ0xn6znnvgT+I0KdXRSaweecewt4K2yAiIhUOnqfk4iI+I6Sk4iI+I6Sk4iI+I6Sk4iI+I6Sk4iI+I6Sk4iI+I6Sk4iI+I6Sk4iI+I6Sk4iI+E6s3+ckUiZlfR8U6J1QIhWBRk4iIuI7Sk4iIuI7Sk4iIuI7uuckFU5Z71tVhHtWuncnZzuNnERExHeUnERExHdielnPzFrgvQm3L94LBd/DexPul1HEngNMBn4BNAA2AP/jnPs/sWyjSCS6LChS/mKWnMysFrASOA7cATjgCeB9M0txzh2OsIsXgYHAOGAH8CtghZl1cc5tiFU7RU63WNzv8QMlaSlPsRw53Q20AS52zn0BYGafAtuAkcD04gLNrD0wHPgv59ycQFkGkAVMAgbFsJ0iIuJzsUxOg4DV+YkJwDm308xWAYMpITkFYk8AiwrEnjSzV4DxZlbTOXc8hm0VkdNMIy8pi1hOiEgGNoUpzwLaRRG70zl3JExsDSCp7M0TEZGzRSyTU0MgJ0x5NhBfhtj87SIiUknE+iFcF6bMooizU4k1s3uAewJfHjKzLVEcqziNgANliBeP+jE2Kn0/2pMx2U2l78cYKWs/tiptQCyTUw7hRzjxhB8VFZQNtCwmNn97Ec6554Hno21gScxsrXOucyz2VZmpH2ND/Rgb6sfYKI9+jOVlvSy8e0eFtQM2RxHbOjAdvXDsv4EvioaIiEhFFcvktBRINbM2+QVmlgh0C2yLFFsduLlAbDXgVuAdzdQTEalcYpmcZgG7gCVmNtjMBgFLgK+AmfmVzKyVmZ00s8fzywIP2S4CZpjZL82sN/AK0BpIi2EbSxKTy4OifowR9WNsqB9j44z3ozkXbh7CKe7MrCWhyxf9HW/5ol0F6iQCO4GJzrn0AuVxwG/xHsZtAGzEW77og5g1UEREzgoxTU4iIiKxoFXJRUTEd5ScRETEd5ScRETEd5ScRETEd5ScRETEd5ScRETEd5ScRETEd5ScRETEd5ScRETEd5ScRETEd5ScRETEd5ScRETEd5ScRETEd5ScRETEd5ScRETEd5ScRETEd6qVdwNipVGjRi4xMbG8myEiIoWsW7fugHOucWliKkxySkxMZO3ateXdDBERKcTMdpc2Rpf1RETEdyrMyEkkX+L4ZWWK3zVlYIxaIiKnSiMnERHxHY2cJIRGHSLiB0pOElNKbhKtHj16APDBBx8AkJ6ezsSJE3HOlV+jTsHo0aPZtWsXb7311hk75ldffcWDDz7Iu+++i3OOPn36MGPGDFq2bFli3IoVK3jyySfZvHkzOTk5NG7cmK5du5Kenk67du2C9Xr06EFGRkbYffTr14/ly5cHv37qqaeYM2cOGzZsoEqV2F2MU3ISXylrcpOz1y9/+Uv69+9f3s0ole3btzNz5kw+/vjjM3bMI0eO0KtXL2rWrMm8efMwMyZMmEDPnj359NNPqV27drGx2dnZdOrUiVGjRtG4cWO+/PJLpkyZQmpqKp999hmtWrUC4Nlnn+WHH34Iif3HP/7BQw89xKBBg0LK7733Xp588knmzZvHnXfeGbPzVHISEV9o3rw5zZs3PyPHOn78ODVr1izzfmbMmEH79u3p3LlzDFoVnVmzZrFjxw62bNlCUlISACnBUhHeAAAY9ElEQVQpKVx44YXMnDmThx56qNjY2267jdtuuy2k7Morr6Rt27a8+uqrPPzwwwAho6iCx61RowbDhg0LKY+Li+P2229n6tSpMU1OmhAhIqfdK6+8Qtu2balZsybJycm88cYbReqkp6djZsGvt27dytChQ2nSpAnnnHMOLVu25Oabb+bkyZPBOsuXL6dLly7ExcVRv359hgwZwpYtW8Lud9OmTfTr1486depwyy23BLe//PLLtG3blnPOOYfLLruMpUuX0qNHj+Blx+IcP36chQsXMnz48JDyvLw86taty6RJk0LKc3JyMDPmzZsXsb9KsnTpUlJTU4OJCaB169Z069aNJUuWlHp/CQkJAFSvXr3YOkePHmXx4sXccMMNNGzYsMj2YcOGsXnz5piOIDVyEqmA/HB5NP/+4Xvvvcfw4cMZOHAg06ZN47vvvmPMmDGcOHGCiy++uNj466+/ngYNGvDcc8/RqFEjvv76a/72t7+Rl5cHeIlp4MCB9OrVi0WLFnHo0CEef/xxrr76ajZs2ECzZs1C9jd48GDuuusu/ud//id4b+Tdd9/l5z//OYMGDWLatGkcOHCABx54gGPHjnHRRReVeH6rV6/m4MGDXHPNNSHlW7du5dChQ1x++eUh5evXrwcIljvnyM3NjdSNmBlVq1YNfp2VlcXgwYOL1EtOTmbx4sUR9weQm5tLbm4uu3fvZvz48TRt2rTIiKig119/nR9//JE77rgj7PYOHTpQr149li9fTteuXaNqQyRKTiJyWqWlpdG2bVuWLFkSTAqXXHIJqampxSanAwcOsG3bNpYsWRJyj6PgKGXChAm0adOGt99+m2rVvF9lXbp04aKLLmLatGlMnz49ZJ/3338/Y8aMKdK2du3a8cYbbwRHbZdddhmdOnWKKjmZGSkpKSHlmZmZAHTs2DGkfP369dSsWZNLLrkEgIyMDHr27FniMQCuvfba4KQR8O4bxcfHF6nXsGFDcnJyIu4P4KqrrmLdunUAJCUlsXLlSpo0aVJs/fnz59OkSRMGDBgQdnuVKlVISUlh9erVUR0/GkpOInLa5ObmsmbNGsaPHx8yk+uqq66ipLUwExISaNOmDePHj2f//v306NGDCy+8MLj98OHDZGZm8uijjwYTE/x0eSvcTLOhQ4cWadvatWt55JFHQi4nduzYkdatW0c8t71791KvXj1q1KgRUr5u3TqaNGlSZOSWmZlJcnJy8PJZp06dWLNmTcTj1K1bt0hZwfbmK80sxwULFvDDDz+wY8cOpk6dSt++ffnoo4/C/pvs3buX9957jzFjxoT0dWGNGzdm69atUbchEiWnCsQPl3IqAk2Hj50DBw5w4sQJzj333CLbwpXlMzPeffdd0tPTeeSRR/j+++9p3bo148aN47777iMnJwfnHOedd16R2KZNm7J7d9Gl3ArXzW9buBFDSW3Ld+zYsbCTKjIzM4uMmsAbORW85FWnTh06dOgQ8TiFE1F8fDzZ2dlF6uXk5IQdUYWTP3q76qqrGDBgAImJiUyZMoU///nPReouXLiQvLy8Yi/p5YuLi+Po0aNRHT8amhAhIqdNo0aNqF69Ovv37y+yLVxZQW3atGH+/Pl89913rF+/nl69ejFq1Cjefvtt4uPjMTP27dtXJG7fvn3Bm/wFFf4ln9+2b7/9ttRtA290V/gymnOODRs2FLnf9O2337Jly5aQ8oyMDKpXrx7x07t375B9JScnk5WVVaQ9mzdvDjvLLpIGDRqQlJTEF198EXb7/Pnzad++Pe3bty9xP9nZ2TRq1KjUxy+OkpOInDZVq1bliiuu4NVXXw1OZAD45z//ya5du6Lah5nRoUOH4D2kTZs2Ubt2bTp16sTixYtDJhXs3r2bjz/+mGuvvTaqtnXu3JnXXnst5JLYunXr2LlzZ8T4tm3bcuLECfbs2RMs2759OwcPHgyZwADw9NNPk5eXFzJSyr+sF+kzc+bMkH0NGjSI1atXs2PHjmDZrl27WLVqVZFnkKKxf/9+Pv/8cy644IIi29auXUtWVlbEURPAzp07S5zgUlpRXdYzsxbAU0BfwID3gAecc19GiOsM3AN0B1oCB4APgQnOuZ2F6u4CWoXZzVDn3JvRtFNE/GfixIlcd911DBkyhJEjR/Ldd9+RlpZG06ZNi4359NNPGTNmDLfeeitJSUnk5uYyd+5cqlWrRq9evQCYPHkyAwcO5Prrr2fUqFEcOnSItLQ06tevH3xeJ9q2DR06lHvuuYcDBw6Qnp5O06ZNI6520L17dwA++eST4PNZ+ZMhXnjhBVq0aEGTJk145513gtPH165dS8eOHYmLi6Nu3bqn9HzU3XffzTPPPMPgwYN54oknMDMee+wxWrRowciRI4P1MjIy6N27N7Nnz+b2228HvPtuHTt2JCUlhXr16rF161aeeuopqlWrFrbP5s+fT7Vq1YpMly/s4MGDbN26lbFjx5b6fIoTMTmZWS1gJXAcuANwwBPA+2aW4pw7XEL4MCAZ+P+ALKAZ8Biw1sw6OOe+KlR/BZBeqGwLZ4FY3O/RvQqpiPr06cNLL71Eeno6N954I0lJScyYMYM//elPxcY0bdqUli1bMn36dPbs2RN8Bumvf/0rnTp1AqB///4sW7aMiRMncsstt1CjRg169OjBH/7wB84///yo2ta3b19eeuklJk6cyNChQ0lKSmLatGlMmjSJ+vXrlxibmJjIlVdeyVtvvcWNN94IeMmpYcOGTJkyhfHjx3Po0CGGDBnC4sWLue2221i0aFGRGYOlVbt2bVauXMmDDz7IiBEjcM7Ru3dvZsyYQZ06dYL18qeqFxyxpqam8pe//IVp06bx73//mxYtWtCjRw8eeeSRIpMhTpw4wcsvv0z//v0j3oNbtmwZNWrUKDLppCws0gwPMxsDTAcuds59EShrDWwD/ts5N72E2MbOue8KlbUCdgJPOOceL1C+C/jIOfeLUzmRzp07u/J82aAfkpMmRPiDH/6T4YfvBT/0w6nYs2cPSUlJ/OY3v+Gxxx4rse7cuXMZM2YM33zzDbVq1aJv376A9/xUZTJgwAAaNWrEggULwm43s3XOuVINE6O55zQIWJ2fmAACl+RWAUWfBCugcGIKlO0GvsMbRYmIlJujR49y33338dprr5GRkcGcOXPo27cvtWrV4pe//GXE+BEjRtCsWTOeffZZwJuRlz+yqyw2bNjA+++/T1paWkz3G809p2Qg3JoYWcDNpT2gmV0CNAH+b5jNN5jZEaAqsB6YovtNInK6VK1alX379vHrX/+a77//ntq1a3PNNdewePHisNPUw8XPnj2bzMxMdu/ezffff1/pktO+ffuYM2dOyHJKsRBNcmoIhHvsOBuIblJ9gJlVA/6MN3J6sdDmt4A1eJf8zgV+DbxhZiOccwuL2d89eBMuIi4VH4kfLoOIyJlVo0aNsOv8lUZqaiqpqalA6R6ErShO10ry0T6EG67Hiz6iHNkzQFdgoHMuJOE550aH7NzsDWA18HsgbHJyzj0PPA/ePadTaI+vKEGKiHiiSU45eKOnwuIJP6IKy8x+jzfKucM5906k+s65XDNbDDxpZuc5576J9lgi5ckPk2NEznbRJKcsvPtOhbUDNkdzEDP7DTAeuN85F346RzGhgT/P+lGRiIhEL5rktBSYamZtnHM7AMwsEeiGl3BKZGb34z0X9Rvn3NPRNixwf+pm4EvnXNE1SkQqMF3ilcoumuQ0C29ywhIzm4A3ipkMfAUE19UIPL+0HZjknJsUKBsGzACWAyvNLLXAfn9wzm0O1LsNb1r63wL7PRf4FdAJCH1to4icFbSArpRFxOTknDtsZr3wli9agHep7e94yxcdKlDV8KaAF3x2qn+gvH/gU1AG0CPw951408v/iHd/6wjezL3+zrkVpTslEZHYUIItP1HN1gusofcfEersotAMPufcfwL/GcX+VwO9ommLiFQOmlhSuWlVchER8R0lJxER8R0lJxER8R29pl1E5DTRfbNTp+QkIhWWnhc7e+mynoiI+I6Sk4iI+I6Sk4iI+I6Sk4iI+I6Sk4iI+I6Sk4iI+I6Sk4iI+I6Sk4iI+I6Sk4iI+I6Sk4iI+I6Sk4iI+I6Sk4iI+I6Sk4iI+I6Sk4iI+I6Sk4iI+I7e5yQi4mNlfSfV2fqyQo2cRETEd6JKTmbWwsxeNbN/mdkPZva6mbWMMvYcM/ujmX1jZkfN7B9m1j1MvSpm9oiZ7TKzY2a20cz+o7QnJCIiZ7+Il/XMrBawEjgO3AE44AngfTNLcc4djrCLF4GBwDhgB/ArYIWZdXHObShQbzIwFvgNsA4YBiw2s+udc38r3WmJiAjE5lX15XFpMJp7TncDbYCLnXNfAJjZp8A2YCQwvbhAM2sPDAf+yzk3J1CWAWQBk4BBgbImeIlpinNuaiD8fTNLAqYASk4iIpVINJf1BgGr8xMTgHNuJ7AKGBxF7AlgUYHYk8ArQD8zqxko7gfUABYWil8IXGZmraNop4iIVBDRJKdkYFOY8iygXRSxO51zR8LE1gCSCtQ7DnwRph5RHEdERCqQaC7rNQRywpRnA/FliM3fnv/nQeeci1AvhJndA9wT+PKQmW2J0J6SNAIOlCFePOrH2FA/xob6MQbsyTL3Y6vSBkT7nFPhpAFgUcRZlLHR1gttlHPPA89H0Y6IzGytc65zLPZVmakfY0P9GBvqx9goj36M5rJeDuFHLvGEHxUVlF1CbP72/D/jzaxwMipcT0REKoFoklMW3j2hwtoBm6OIbR2Yjl449t/8dI8pC6gJXBCmHlEcR0REKpBoktNSINXM2uQXmFki0C2wLVJsdeDmArHVgFuBd5xzxwPFy/GS1c8Lxf8C2BSYHXi6xeTyoKgfY0T9GBvqx9g44/1oRecgFKpgVhvYCBwFJuDdG5oM1AVSnHOHAvVaAduBSc65SQXiX8GbKj4O2AncB1wPdHXOZRaoNwV4AHgUyMRLYCOBwc65t2JxsiIicnaIOCHCOXfYzHoBTwEL8CYp/B14ID8xBRhQlaKjsTuB3+KtKtEAL9H1L5iYAn4DHALGAE2BLcAtSkwiIpVPxJGTiIjImVZhViU/Q4vTPmRmbwXqOTNLj/mJ+FAZ+/Z3ZvaOmX0f6LP/PM3N9S0za25mTwe+v44E+iMxythKuTBytH0W7c9wCce528w+N7PjZrbFzO6N5XmcKbHur7J+35nZEDNbH4jdbWYTzKxqNLEVIjkVWJy2Ld7itCOAC/HW56sdxS5exFtD8HG8+2Hf4C1O26FQvbuBJsCbMWq678Wgb0cDccBfT1sjzx5JwC14j2B8WMrYyUA68AwwAFiNtzDyz2LZQB+Kts+i/RkuwszuBmYCrwH9gcXAs2Z2X9maXi5i3V+n/H1nZv3w+nRNIPZPePMWfhfVmTjnzvoP3n2qXCCpQFlr4CTwUITY9niTPO4sUFYN757X0kJ1qxTY7oD08j53P/dtoT5LCvTZf5b3OZVjX1Yp8PdfBvojMYq4JnjLe00sVP534NPyPq/y7rPS/AyH2X814FtgXqHy2XgrIlQv7z4or/4q6/cdsB7IKFT2ON7M7KaR4ivEyIkzszgtzrm8WDb6LFGWvq2sfRZWGfqi0i6MHGWfRf0zHEYXoDFF+3YBkABcXaoGl7MY99cpf9+ZWQugQ5jYBXiPFw2I1MiKkpzOxOK0lVVZ+lZiQwsjl6wsP8P5CwwU/h6vyH17JhbkDtuvgf/YHokQC1Sc5HQmFqetrMrStxIbp7QwciVSlp/h/G2F4yty357WBbkLbQt3nOKWxAtRUZITnP7FaSsz9U/50vdoycrSP/l1KtMzNad1Qe5CdU75+7aiJKczsThtZVWWvpXY0MLIJSvLz3Bxo4CGhbZXJGdiQe6SRlcNIsQCFSc5nYnFaSursvStxIYWRi5ZWX6G8++fFP4er8h9eyYW5A7br4FnrmpFiAUqTnI6E4vTVlZl6VuJDT8sjOxnZfkZ/gfelPFwfZuNNyu1ojntC3I7577EW6ouXOwJ4O1IjYz2ZYN+Nwv4NbDEzAouTvsV3sN1QPjFaZ1zG8xsETDDzKrz0+K0rSnUsWbWGUjkp6TezsxuCvz9b2Fmv1QEp9y3gfJr8abqNg0UdTazQwDOuVfPyBn4SIHvl06BPweY2XfAd865jECdk3jP3dwF4Jz71syeAh4xsx/5aWHkXkQxnf9sF6nPSvkz/AWw2znXG8A5d8LMHsN76PZr4D28fv0vYLRz7t+n+/xiLVb9VZrvOzP7O9DKOVdwZuSjwF/NbCbwMnA53kO4f3LO7Yt4IuX90FisPkBLvKeRfwB+xFvFofDDZ4mEeXgWbwWD6cA+4BjwT6BHmGPMDcSH+ySernMr708Z+/aD4vqsvM+rnPqyuO+fDwrVmVsormrgB3s33vTeT4Gbyvt8fNRn0f4M7yoYV6B8JLA10LfbgFHlfd4+6a+ovu8CP+e7wpTfiDeCOg58ifcQbtVozkMLv4qIiO9UlHtOIiJSgSg5iYiI7yg5iYiI7yg5iYiI7yg5iYiI7yg5iYiI7yg5iYiI7yg5iYiI7yg5icgZZWaPmtkWM8szsyHl3R7xJyUnETnT/g78DPg/5d0Q8S8lJ5FSMrOnzeyt8m7HqTCzD8zsg/Jsg3Pun8657eG2mdmDZvapmel3UyWnbwCRUjCzC/AWCZ1Y3m2poP4MNAHuKO+GSPlScpIKwTzVz8ChHgA2OufWnoFjnVXMLNPMDhTzaRHNPpxzR4H5wNjT21rxOyUn8RUz+8TM/mJmk8xsu5kdC1zm6V2o3kYzm2Nm95nZ/8V7gVn/wLabzSzDzH4ws0NmtsbMBpb2GGHaVhPvZWn/f5ht6WbmzKytma0ws8Nm9qWZ3RnYPsLMPg+05/3ACKzMsRHaOywQd9zMssxsaJg6SWa2wMx2mtlRM9thZs+ZWXyBOjcF2tc+TPwHZvYPAOdcR+dco2I+X0XT5oBX8N6V1rUUMVLBKDmJbwTexnkZMBC4GngQuA3vpZivm1mjQL0awCVAv8D2x4DrgY/NbCreL7csvERyG/AxUK80xyhGKtAA+LCEOouBZcAQYB0w28x+h/cyt/HAncDFhElwZYwNYWZ9AvW24b1T54/AnwLxBZ0P7MEbEfYDJgG9gb8VqPMmsBfvcmbBY1wMXEuBl07GyAa8d4f1j/F+5WxS3i/G0kef/A+QgvdStAwKvJAM7xegA4YGvu4Y+PqDQvWGBcqLfQlftMcoJvZ/gDygRpht6YH42wuUxQMnge+BegXK7w/UbVXW2BLaugrYDFQpUHYVhV46FyauGl7SdsDlhc7vX0DtAmXTgRwgrpT/zhPwEuJxvFek7wGaFqrzId4rw8v9+1Kf8vlo5CR+kv9a6Uedc7kFyj8P/JkQ+PPywJ9jC9WbCLzlSn79e7THCOd84AdX8qu7387/i3MuB/gWWO2c+yHMsQrfhylLbJCZVQWuAF51zuUV2Oc/8d4EW7BujcBzR5+b2VG8y6P5I8OCo6zngVp4o0zM7By8SQvznXefKGrOuSecc82dczWdd8mvuSv62u7v8PpbKiklJ/GTjsBe59yqQuX5v6T2BP68PFAvOCnBzNoAF+FdGovFMcI5B+9/+yXJKfT1v4spy99frGILagRUB/aH2Va47Pd4o6KFeJc6r8S7DBhyDOfcXmAJcG+g6GagIbG/pJfvKN6rxKWSqlbeDRApoCPwdZjyW4Ej/PQ/+suBNYXq5CeXvTE6Rjjf411u87sDeCOgc8NsOxfYXeDrYXijnyfyC8ysTjH7fRb4u5l1wrv/9KFzbnNsmlxEQ7zzkEpKIyfxhcBDl+2BVoFJC/nl5wOjgGecc4cD9VLwJgwUlJ+Ukst6jBKa+TlQ3cyaR39mZ17gcuUa4KaCD7Oa2VVAYqHqtfASWUF3FrPflcD/xbvX1A3vmaTTpTWw5TTuX3xOyUn8oi1QG++y2Vwz621md+CNZD4H0gL1LgTqAJmF4ncCHwETzex+M+thZrea2Twzu6iUxyhO/nI7V57yWZ45aXjn+6aZDTSz/wT+AhS+t7McuMPMRpnZdWb2Z6CkKdx/BrrjjWpei32zwcwa4F2i1fJGlZiSk/hFx8CfP8Obrv0W8Ae8Kc29nXPHAtvzJ0OEJCfnnAP+A3gV7wHOFXjTp6vgJa7SHCMs59wu4BPghlKf3RnmnHsP+DnepIbXgXF408ULj0ZGA0uB3wKLgLoEJj0UI/+e3lznXKT7b6dqIN69tTdO0/7lLGDez7RI+TKzp/CmgEe1kkB5HSMwAvkTcJ5z7kis2na2MLO78SZBXOSc++I0HeNt4IBzbsTp2L+cHTRyEr/oSNH7SH48xgK8CRWjyt6cs4eZtTOzG/Cm6795GhNTB6AnWruw0lNyknJnZgZ04DQmp1gdIzDZ4L/wZvZVJs/i3WPaCvz6NB6nKXDn6Up+cvbQZT0REfEdjZxERMR3lJxERMR3lJxERMR3lJxERMR3lJxERMR3lJxERMR3lJxERMR3lJxERMR3/h/xEcxc9VIGwgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(nrows=3, sharex=True, figsize=(6,4))\n",
    "# ticks = [10, 200, 100, 600]\n",
    "\n",
    "for axi, scene in zip(ax, [1,0,3]):\n",
    "    prcps = np.log10(np.concatenate(precip_by_scene[scene]))\n",
    "#     axi.set_xticks(ticks)\n",
    "#     axi.set_ylabel(\"#\")\n",
    "    axi.hist(np.array(prcps)[~np.isnan(prcps)], density=True, bins=np.linspace(-2, 2, 21))#bins = np.logspace(np.log10(0.01), np.log10(30), 21), log=True)\n",
    "    axi.annotate('{} ($\\mu$={:0.2f})'.format(str(scene_dict[scene]), 10**np.nanmean(prcps)), xy=(0.99, 0.95), xycoords='axes fraction', \n",
    "                 verticalalignment='top', horizontalalignment='right', fontsize=16)\n",
    "ax[2].set_xlabel(r\"$prcp$ (mm day$^{-1}$)\")\n",
    "# ax[2].set_xticks([np.linspace(-1, 2, 6)])\n",
    "ax[2].set_xticklabels([10**i for i in ax[2].get_xticks()])\n",
    "# plt.suptitle(r\"Accum concentration\", y=1.1)\n",
    "# ax[2].set_xscale('log')\n",
    "# ax[2].set_xticks(ticks)\n",
    "# ax[2].set_xticklabels(ticks)\n",
    "\n",
    "fig.tight_layout()\n",
    "\n",
    "fig.savefig(\"/home/disk/user_www/jkcm/misc/agu/precip_pdf.png\", dpi=600, bbox_inches='tight')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 191,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "light drizzle frac: 60.715990%\n",
      "heavy drizzle frac: 0.000000%\n",
      "light drizzle frac: 17.175573%\n",
      "heavy drizzle frac: 0.000000%\n",
      "light drizzle frac: 37.646722%\n",
      "heavy drizzle frac: 0.000000%\n"
     ]
    }
   ],
   "source": [
    "for axi, scene in zip(ax, [1,0,3]):\n",
    "    prcps = np.concatenate(precip_by_scene[scene])\n",
    "    light_drizzles = np.logical_and(prcps>0.1, prcps<2)\n",
    "    heavy_drizzles = precips>2\n",
    "    \n",
    "    print('light drizzle frac: {:%}'.format(sum(light_drizzles)/sum(~np.isnan(prcps))))\n",
    "    print('heavy drizzle frac: {:%}'.format(sum(heavy_drizzles)/sum(~np.isnan(prcps))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "radarlidar_file = r'/home/disk/eos10/imccoy/CSET_RadarLidarData/RF07_COMBINED_HCR_HSRL_data_mask_version4.cdf'\n",
    "data = xr.open_dataset(radarlidar_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:                (layer_no: 20, range: 699, time: 54745)\n",
       "Coordinates:\n",
       "  * range                  (range) float64 40.0 60.0 80.0 100.0 120.0 140.0 ...\n",
       "Dimensions without coordinates: layer_no, time\n",
       "Data variables:\n",
       "    absolute_time          (time) timedelta64[ns] ...\n",
       "    lat                    (time) float32 ...\n",
       "    lon                    (time) float32 ...\n",
       "    alt_msl                (time) float32 ...\n",
       "    ant_elev_angle         (time) float32 ...\n",
       "    ant_az_angle           (time) float32 ...\n",
       "    beta                   (time, range) float32 ...\n",
       "    beta_molecular         (time, range) float32 ...\n",
       "    beta_std               (time, range) float32 ...\n",
       "    ratio_bscat            (time, range) float32 ...\n",
       "    dBZ                    (time, range) float32 ...\n",
       "    vel                    (time, range) float32 ...\n",
       "    sp_width               (time, range) float32 ...\n",
       "    SNR_HCR                (time, range) float32 ...\n",
       "    combined_mask          (time, range) float32 ...\n",
       "    hsrl_attenuation_mask  (time, range) float32 ...\n",
       "    layer_top              (layer_no, time) float32 ...\n",
       "    layer_bot              (layer_no, time) float32 ...\n",
       "    layer_count            (time) float32 ...\n",
       "    lidar_cloud_base       (time) float32 ...\n",
       "    mask_flag              (time) float32 ...\n",
       "Attributes:\n",
       "    Resolution:       HSRL/HCR gridded to 2 Hz, 20 m vertical resolution\n",
       "    Nav Data Source:  HCR data stream (corrected altitude)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fa852d74898>"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEDCAYAAADX1GjKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XucHHWZ7/HP0z23JIBJIO4qEoF19UhEUFlFQS5BBT0S1oN3XFfYDV7ZlXATUV+KICt3XUVBvB1QYVflEN3DnQCKQQ13oqAcuSq6hAmE3Gamu57zR1VP19T0pXqmZqp65vvOq16Tqf511a+reuqp37XM3RERkdmtlHcGREQkfwoGIiKiYCAiIgoGIiKCgoGIiKBgICIi5BQMzGwnM/uhmT1jZhvM7MdmtjiPvIiICNh0jzMws7nA3cAQ8CnAgdOAucDL3X3TtGZIREToyWGfy4FdgZe4+4MAZnYP8Hvgg8C5OeRJRGRWy6NkcAMw4O77JNbfDODu+09rhkREJJc2gyXAfQ3WrwV2m+a8iIgI+VQTLQTWN1g/CCxIs4EddtjBd9555yzzJCIz0O23377O3RdNZhsHHzjPnxqsptvfPUPXuPshk9lfXvIIBhA2GidZqzeY2dHA0QCLFy9mzZo1U5EvEZlBzOyRyW7jqcEqv7omXWfH8vN+v8Nk95eXPKqJ1hOWDpIW0LjEAIC7X+Tue7n7XosWTSrQA/Di09q3U+91ZHZt2fu8/exMtvPST53XUfoXfuuLmew3addzz2n62is/2P19APZ+T/PPV2QTPd+t/h5e9MWpP5+7njOxfbzppo9nnJPxHAhS/utmeZQM1hK2GyTtBvxmujJh1ZYFEQDWfHtFZvsLetrvL42+Zzp8g2ez33FabXYGzIru5bxzMEETPN1BqyvBFH2F4rw8sS/NtQecn3FOxnOcEU9XTdTN8igZrAT2NrNdayvMbGdgn+i1aWGV6dpTqNqXzXZKI8W40raMMU1ee+HXz5qSvLSz508/1fF7qr3TcAWcCqWJfT+8RTCYqvuJMQp+uGdDySCPYPAN4GHgSjM7zMyWAVcCjwEXTuWO9/xIvShaKkCgX3JSZ1U+AHOe6vALN1Wxo8U3p9rX+C/bAuO1157Ebp/s/HNPxtBI5wVgL/BELa8/LPugGvQ2/6JMy7EocDBwnKqnW7rZtH/loxHGS4HfAZcA3wMeApa6+8ap3He1v/6NK2+dyj2lE/R2/p5SJd0Xbt//FV0w2iR/8/M/1nkmAG9xB1qZ13i9jRjDlR4qczrf315HTW87RJGriX525QlNX7OJlgxaXAlsGm54p6X0MQkBnmrpZrnc/7j7o+5+uLtv5+7buvvfu/vDU73f+EWqvHV6T9zIXON17xzbiDyRYFAZSHfKRtso2vyRXfWnr3SeiTbbbfaHXRoxAofK3M6P/cYdJ361sAm8tTTN1YhZiQeDV3w4fQC1FiXl8tBkcpQ2A8W9kDpQxVMt3azAheHsxRvJ7rogu8bhdPs2ti4ID/duJ4fVJMm7sTe9+nMtt/HGvU9l64L2V7Z9Dj+bSq0UNMGGuXa8p/l2e5vNLhXASLWMt6iSaGYy14pSqfNb29KIs99bz5z4TnNise9UUE4fBZMdKvb4WL0qL6v2rpYKfiVSyaBLNesS2uoCNtW8p15NNbJNmI9k8XvL8+a23kZvmeH57f/Ag9i+pqwudgLHsjwMlUqZYGAiF+fx6w584791tI037n1q+sQGlTnd9+dhsahZaf11aqk6UP9/0D8NfzcFriZyYMQ91dLNuu/bnkKzLqEtu89NMS9DrbPB8PbVaN3YL09loPVfxPBzelN9hqDXRqugJlqH3FaLEkezPA4viA5Af+et943qtFdd94l0743qrTb8TZPGjAaCsjG8zdidvmHf01K/Py8eq6OrdtA2kyytVftjv0xH19Kp+p5mwFNWEamaqGD2/Oi5Y3oNxQV9+Z2soKf+x9mzsTy6Li558Um65acnpvsMXm8zsPLUtP4lt/t3H6gf82b1z+UtJdyN0tMTaCyZhEolPN5btw+P7wGHtB+YVZmTuCACG3eaQMv3NIvfnNZKoKnel7gYj2nPmo4/mwIHAxyqKZduNuOCQbNujQDBnBzPlsHm54cX0GpUTZKsB3/OQ+1b6spb2t+m/fJ7x9WL+R3c1f3dP57Lm16Tsiolsd0tO9RX9G5qfJznPW64t++qWBtlvfNXwgb31777nNH9HXTgFzq+Qx8e6mH/t5xJaTjM18bn9XLAIV/kkBef2PQ9lXk2ms9aw3/fxs5LNAe8qbOqrMnyoH4eOrr5SbT6j+lqOi3BYBr2MUHhCOR0Szcr8ClIZ/ma94/53TysJmnEe/I9XdU50f6j7CUvijfeePLo//f4lyZ98VNe3GtdIzvpSRO2NaTrU2mJSNazpf77yLyxO1269AwA+p92zMB7W5+HWm8jq4TbcYPylvC1p146QHWgs/o+D4zKnNJor5jqAAw9p4wFza9y1T7wqHT11JLwmGxYXN/vfoeemaq/f3XONPdRjX2kVj2EXvGhsaXnZPtVKdag7C2CSrNSeMcK3JsIjGrKZUJbN7vazNzMTkusX2BmF5vZOjPbZGbXm9numXykBro+GKz+085jV3iL0ZRT1LMmlYD60Y6y0WoUdNDkGpK23SMYDQbpP3O1H57deaB9QsbWTQPc9dV6O03yIvTs4rC+ZWCwSqkUUN7c+mtXO3+lKBhY4KMlnc3PM8pbO+v36ZUSw9uW6NtUv+JV+4yrHmx+MTcPj8cbXn86lXnOm179ObY8t/6Zb/nJiWyd3/5Cv3XBNAcDa/L/hOT0KOPGEsR+97nNo4pnNM0KBks+kf1gxO36nvuqyW4jbEC2VEunzOw9wB4N1hvhjAyHAMcAhwO9wCoze8GkPlATXR8M7ls2tlrDqi367+fdY6FWMomO+sBgiww1OTPl4cbr43X2QL300cEXtDLX2PzclF+JRIx51T/X918eHvvipueFeRiaH267XTVRz8YoCETXoKHnlEaDYNDnY0oGB7/iM+3zWjWGt7XRAXt3fm3FmAGIjVgAI9uClwwvO5tfML5rzq+/27578q+/M71dmEvRDc8eHzuP0nDzzzhuUF0wNm11oH4OS33Ng0GyXWXCpuhvc8Pwf98+2W2E4wyyLxmY2XzgPKDRl2QZsC/wD+7+A3e/OlpXAprXb05C1weDJHNveledZ2O/lxk92qUt4X9GWnVuaZLXWnVJUrzOft+3nTX6/mCk8Sne73+O70Pf0R+225gAsHX7+v6TJYNacB7arkS1Wmr7rat1I63drc59sjpaWrAqrH9xveP7NXemaOOoGL2bnL5n6iWK2y9ufZGu9kO137nh5k/iBsPzSl01AV9lbutqouS5Tk7PUk3ZvpZZDz1vXhqejO36/3rSJQOAwC3V0qEzgbXu/oMGry0D/uTuq2or3P0Z4CfAYRP+IC3MvGAQTGxkb9yeH81+6oOgBxgJvyy1i12jvvM11Wa1NU3OWLxPebxXklcaf0E3/XX9r3j/t5zJQft/gco8Tz0Ng1dtTAkgvv/bv5m40EZZCPpheEtv04tqrbG1Jwp4tTvbn19xAlYN5+TpfdbGVNe89t3tp5oubS2x5lsreHan9NGudm7euPepmMOvLjmOBz57bMv3JEeY56FWLVjth55NY899vH6/FlxHv+vB2GqaeOmtWffkV3zo3Oym7fCpmfbCOxh413QbZF8yMLN9gfcDH2mSpNUTIReb2Tadfo52ZlwwuPNrK5oWOa1Nw2VNswboyfAe6B0M/wLL0UWu1V1Vs7uk3g2N/zDj7SSVOda22B2/kFcHjPUvHiDohVKTaqhxqobHWqeDFoPQRueBCsBHSpRGjNe9Y/yF86Zrw3EDtZlZ43np2Qyb/rrMto8FzP1LfV9znmyf4VpQadVgnGSVcKls20dpqP459zm8+QX/F/9xfOrtH/zyTzdcf8jLTkm9jUbin7AnMf9WvGqsNj+Ul6J1pfq6V3z43JalCoBXLj83DBgZlpamIhg8u/mJDKqJjCqlVAuwg5mtiS1HJ7dnZr2Ek3Ke7e4PNNltqydCQsqnQnZixgWDVnpa1H3GTcUkZUGvj/6B9Q+2TgvNR332bmm8PiiPvbiP1vlWGp/iX11y3Oj/K/0lhhaEXSmbdQsdx2Hg6XTHc+0Xwzvquf8dwHCJ0nDz5zssXXrG6AUmXnLqf9rpe9bBYXhbG31Y0I03nNxgK2OlDnDx91TC78Gq6z7BvMdtdGqKTgJKK8+8bH7D9Vffd/qkthsMh1/eUnV828xQ/JFSsY4Fe3zsvPDmozZoPRh7YfZg/Lm64xsrRtPWdDIX0viM2+j3pIg6qCZaV3sIV7Rc1GBzJwFzgFYn25jAEyEnY0YGg3jjV1xvX7peKMlePq9cPrlqo9cccQ5eCidqA7jn/PBL3+pOqFkja9N59uPpg1idcIq7rd5NAUEfUHL6Nqa72FmlNGb2zFKT6qhx73OjZ5ONa2Su2bhj3+g03f3P1NPcfvEKbvv+cWzZoYT3wK0/TH8XXt5aKxmkfgtz/xLQ+2z4vru/ciwj24QX2Z9f0XzG0E4s+NWf6zPLZsi3hPm0aljyjH93gz4f7VLav9545QfPpTzkBH1hbKgNPAt6GB138Mrl5zZtdypVx/YIndQzQlJ+f/LgGMNeTrW0Y2aLgVOATwP9ZjY/akgm9nuZsATQ7ImQ0OKpkBM1M4NBk7vqcsrRuOPO6SS/p7Vqm+QgoIlMG7ztY41vc+OTiZWH66WQNGe4Z3MVL4X5GZmb8sMmCwUpYkjvlgAbMczh1h81v5jXqukaVVUEva3bWhopR+l/eelxrRPGbPvIljH7WX1Z+vemcdWDZ3U0kVxao912o8F98aqhoKdeIqsMRDcWFpaCSiOx76PXj/0d31gBQ40vcndcuGJsdeYkrialoam5FG0353mZdC0NKKVaUtgVGAAuJbyg1xaA46P/707rJ0I+OhXT/ec4W8/U8SbBoCfl7JVBYpbGOy6aXPfA0ghgMLLj2At509k9W1h1fX0+nqUHncGNN5zMgW/8N4Jl9VbzoKf+x2xz2t+ubf6rXoLecH4aTzlKzRJ3cs26vMZtWVDGqt6y2mZgsMrmvwq/lo06Atx3dudVCds81nlldHnTMGvP7PwJaZ2Ys24C9Vdt1EpBvZucrdsb5diT8Xo3GXd/pX784v//H585b7REbT72fPaub37H27++vv3JPJOgPDQ1JYMNWybfZgBMeEBZA3cBBzZYv4owQHwTeJBwjMGRZra/u98MYGbbAYcC388qM3EzsmTQbIRr6gFYGXchLI+Edd09/WMvzC3rslP8Za3/237e+LrPs2VR75jyerwBPM1EdVsXlkbz0pPyOQ/JaqGgp3XjKsDwdgaBtZxN85b/OpE13wqD77w/Z/M4ut4tnQeDVF1WJ+mGVZ/MZDvx6aZrVZGjd/YXxm5kYqe2J9FF2ar19975tRVjvn/9LcbDxJ91cOfXJ37TVJqCYPDmF2bTBuFuVL2Uamm/LX/a3W9KLtHLj0S/byQMBquBS83s3WZ2cLTOCLukZm5GBgOaBIO0A7Cy7tXQszWak6iaONytjn6KwOUl2LDr3LCKJ9abpzTioz2WLMWoay/X+5lX5qQ8Rg3GEmxZ2P7rVN4Kw/PTHeCs6ud//uNstlNU8Tv8WhvRnV9fMe6C3xMrid577tgL5dy/OL3P1n+P36j0PdP8OzSuG/EEWQAvOyHjEci9PWyz3Y7ZjDPAUi1ZcfcAeCtwHXABcAVh5eyB7v5YZjuKmRHB4DXvG9vXvDSJOYiazgk0CT//8QnhoJpEQ9zdX25x55IicPVvcIKeqDE2tun4g3uqTep740ojzraPBlC1DhqQE9kteduqIu8J57yZqioBGXunft9ZY79fa89s/n27/eIVY6rgemM10pO540/Nsp9i/qoHz8qk703YgNyTapnwPtzN3T+VWDfo7ke5+0J3n+vuB7n73ZP+QE3MiGCQbBhsdjdcDVIU4wyGt8skW2NYkP2zBQaequAG8x7dPO610ViSIhjMWRfwy0uPozxk3Pb9dA2lyQv6/N9D/7Otq3XuO+tYrFqfbqKR/d/cfnppae43Z2RTNXLvee2386p/St/LbslJrW+ySkNk3mnykJeclEmVb8YNyIXVlbm//w9/afl6b2/ji1IpUfUSr2utKQ85/Zl32iL8RnXyZU9RTXTzVScR9MGmxeMr4Ws9YcqbSm3/aFdffvyY90zEr7+zIlV1TO9GmPvn5p/t5qtOmngmZFqlrSLa45jzxjx/vJHyMJm31V39wBczGYEMUHVLtXSzrgwG7YqT28xp/FyAZANyoykf7vz6inH1qVkoVYyH39d+gNRo+haTjMUFPbB1/vh5c2oPdC9vNYafM7EeQq30TrBj28BTTv+Gbp/5XTpx978fi5eclx3fvHTg1noupYmyDJ440+EI5K7VlblPNton73z7exp3pxxXTTSNgbxnfE1OS2n/MMpbaVzPH62a+0djYDBtD6EORpFO8Nj9+jsr+MV/ph8wJjNDuxuN5KjnzGRU2gi8lGrpZt2Z+8T3qm/j2G9Rs2AwUhlbf56sFpnsSONW+prMKdRMchbJZrwE2/5hEzYy9qDU3n/3V45NNY3y3u89B6uk7yvenzLAiEB4c9Oq5tN7Jj/BZMP9ZrCNcKI6lQy6Qq3Ou6a3yZV0eGjsty1ZTTSVD1u682ud9chIO7T/jm+s4LrbPoMl5o8pDcEbX/f51Pu77fvHsf1vK2N6IrWy5tvTO0+/dK+X/+t5lIfg3nOaV79ag7mUsvDss3/MZKK6ES+nWrpZVweDvd/bePrikSZTfnqyXj3RtHD7xSvGlA52X5F9N9O0Oi1JJIOHBc51v2g8M2YzM70/vuTjni8d27JbK1DoZ0W4k9mgsyLrztzXJuRs8rSqTcN9DdcnZ19s1H2uNhsjZNOgVZvtslO9mzusVkrUyU6mZ5BIlnY/rv1NVbGvo+kGnGU56CwPhT4FzdR6yjSrqtg01CQYNJnOuZksvqC3/HRiT6ir9nX2xUpObbHw/iaPRBOZZmmqPNuWHHLkzI6SQVdOVNeu18HQUJOP1WFRtDqns/RZ+bsPnNtx+0Vy/pjrfza5h6SIZCVtVfpUttlNVrc3DqfRlcGg1uvgde88u+HTpUY2Ny4ZUO3sbvs3X8jnbmUiD1GfirERIllI+92cyEOIpoMzoecbd52uDAY1o4/sS2g67UOB7zxEZruJTE8+HRwYmcS8Q92iO8s+UTXR8LwmwWCwSYflBo/vE5GZ5c2LP57xFsc/+L7Z0s26M9xZ4mdCs6kcpuppSiJSHFc9en6m23Po+tHFaXTnJ4yqe5p1v/TexuunYu4TESmWg1/e2fiaNFQyKKromDedbrlJdVDfYHfGPhFJ75p70o+8T8PdZkXJoCuDQbs7/GZd1O7/fDEbqESkuMIG5O6eaiKNrgwGL93lr1q+ruogEcmOdf2AsjS6Mhi0U5mn+fJFZrt93n52JtsJG5C7uz0gjRkZDJo1IIvI7BH0ZHcB1wjkLmXDM//EiUhr1Yyej6ARyCIiXazTmX9b6faH3acxI4NB74aZH8VFpLWsOgC5w0jykbkz0IwMBgseUJuBiGQjrCaa+cFgRn7CX/3vJoPRRHL2+mVn5Z2FWaO8NbubwtkwAnlGBgORovrZSj1adLpkdaxrXUvTLN1sRlYTiYhkZ3ZUEykYiIi00e3PN05j5oc7EZFJCHsTlVMt7ZjZ283sR2b2iJltMbMHzOwMM9s2kW6BmV1sZuvMbJOZXW9mu0/Zh0TBQESkpdqgs4zaDI4HqsAngUOArwEfBq4zsxKAmRmwMnr9GOBwoBdYZWYvyP4ThlRNJCLSRobVRIe6+5Ox3282s0Hgu8ABwI3AMmBfYKm7rwIws9XAQ8CJwL9klZm4tiWDrIs1ZjZgZmeZ2RPR9lab2X5ZfigRkaxk2ZsoEQhqfh393DH6uQz4Uy0QRO97BvgJcNjkPk1zaaqJsi7WfBNYDnwGeCvwBHCNme056U8jIjIFAi+lWiZo/+jnb6OfS4D7GqRbCyw2s20muqNW0lQTZVasMbM9gPcCR7n7t6N1NxN+yFOj7YiIFIa7UUl/od/BzNbEfr/I3S9qltjMdiS89l3v7rX3LQQebpB8MPq5ANiYNkNptf2EGRdrlgEjwOWxdBXgMuBgM+vvKPciMiGvP0wjoTvRQTXROnffK7a0CgTbAFcCFeDI+EuMPul97Fuy/ExJEy3XTLRYswR4yN03N0jXB7xogvkRkQ787EqNhE5rKkYgm9kAYdX6rsDB7v547OVBwtJB0oLo5/oJfpSWOg4GLYo1jTIYL9akSdfoAIiI5CrLYGBmvcCPgFcDb3H3exNJ1hLeOCftBjzq7plXEUGHwSCDYs2Eiz9mdrSZrTGzNU8+2ajmSkQke1mOM4g63XwPOAg4zN1va5BsJbCjme0fe992wKHRa1MidTDIqFjTLt1gg9cAcPeLavVwixYtSpttEcnAvm+b3W0MAZZqSeGrwDuAs4FNZrZ3bKn1vFwJrAYuNbN3m9nB0ToDzpyCjwekDAYZFmvWAruY2dwG6YaBB9NmXESmz8+vmL1tDO5QCUqplhTeHP08hfCCH1/+OdyfB4Td7q8DLgCuIOzef6C7P5btp6tLM+gsy2LNSsLxB++IpesB3gVc6+5DE/kQIiJTKcNBZzu7uzVZPhtLN+juR7n7Qnef6+4HufvdU/kZ04wzqBVrTicq1sReezyqLooXa04grBY6mUSxxt3vMrPLgfOj0sZDhAPYdgGOyODziIhkqtZmMNOlKddkXaw5Evg2cBrwX8BOwCHufsekPomIyBRxt1RLN2tbMnD3ndNsyN0HgaOipVW6LcCKaBERKbzZ8DwDzVoqItKCO7OimkjBQESkJaOarqdQV1MwEBFpo9vbA9JQMBARaaE2N9FMp2AgItKKh+0GM52CgYhIG+pNJCIyy7kakEVEBFRNJCIiqDeRiMis565gICIiqGupiIgwO9oMZn4TuYhMi6UHnZF3FqaEYwRBKdXSzVQyEJFM3HjDyXlnYcrMgoKBgoGISEtqQBYREWBWFA0UDERE2lDJQERklnMgCBQMRERmNwdUMhARkdkwzkDBQESkHQUDEZHZztSALCIiqGQgIjLrObh6E4mICHrspYiIzIZqou6eZk9EZDp4yiUFM9vJzH5oZs+Y2QYz+7GZLZ6SfHdAJQMRkVYyHHRmZnOBG4Eh4B+jrZ8GrDKzl7v7pkx2NAEKBiIibWQ46Gw5sCvwEnd/EMDM7gF+D3wQODezPXVI1UQikqn93npm3lnIXmDplvaWAbfVAgGAuz8E3AocNkW5T0XBQEQydctPT8w7C5kzT7eksAS4r8H6tcBuWea5U6omEhFppYPGYWAHM1sT+/0id78o9vtCYH2D9w0CCyaUv4woGIiItGSdNCCvc/e92qRpFFpyH8igYCAi0k52DcjrCUsHSQtoXGKYNgoGIiLtBJltaS1hu0HSbsBvMtvLBKgBWUSkldo4gzRLeyuBvc1s19oKM9sZ2Cd6LTcKBiIibWTYm+gbwMPAlWZ2mJktA64EHgMunLIPkIKCgYgUxmved07eWWgso+koohHGS4HfAZcA3wMeApa6+8YpyHlqajMQkcLYOn/m35+6+6PA4XnnI0nBQEQKI+jLOweNpawC6moKBiJSGEFv3jlowEk71URXUzAQkcIoZDCAWfE8AwUDESmMan/eOWhsNlQTzfzWGhHpGl5yXvVPuc3i3FyGD7cpKpUMRKQwgt6Clg66/EKfhoKBiBSGlyHoK1ZjbQcDyrqagoGIFIb3OJW5xQoGgHoTiYhMJy97IXsUqWQgIjKdSuBF7NaiYCAiMn285MULBrOkzWBCh93MrjYzN7PTEusXmNnFZrbOzDaZ2fVmtnuD9w+Y2Vlm9oSZbTGz1Wa230Q/hIjMEKWwqqhwZkHX0o6DgZm9B9ijwXojnI/7EOAYwomYeoFVZvaCRPJvAsuBzwBvBZ4ArjGzPTvNj4jMMAVsq7Ug3dLNOgoGZjYfOA9Y0eDlZcC+wD+4+w/c/epoXQk4MbaNPYD3Ase6+zfc/QbgncCjwKkT+hQiMmOkf9ywZKnTksGZwFp3/0GD15YBf3L3VbUV7v4M8BPgsES6EeDyWLoKcBlwsJkVcciJiMxmqiaqM7N9gfcDH2mSZAlwX4P1a4HFZrZNLN1D7r65Qbo+4EVp8yQiM1BBG5AzetJZYaU67GbWS/hItrPd/YEmyRYC6xusH4x+LkiZbmGTPBxtZmvMbM2TTz6ZJtsi0m2KekFVyWDUScAc4PQWaYzGhyNZA5g23RjufpG77+Xuey1atKhVUhHpUhZYMdsMFAzAzBYDpwCfBvrNbH7UkEzs9zLhnX2ju/paiaBWGmiXbrDBayIyG8QuqLXZS/OexdRQb6KaXYEB4FLCC3ptATg++v/uhHX+Sxq8fzfg0djDntcCu5jZ3AbphoEHO/kAIjIz7HHMeRAADktOOo+RbcIiQu4XWbUZjLoLOLDBAmGAOJDwAr4S2NHM9q+90cy2Aw6NXqtZSTj+4B2xdD3Au4Br3X1ooh9GRLpX0BdWE5kDVn/qWSEusrOgmqjtdBTu/jRwU3J9OMaMR9z9puj3lcBq4FIzO4GwxHAyYSnrzNj27jKzy4Hzo4bph4APA7sAR0zu44hIt3IjLBkAQU/9TjUo55WjmC6/0KeRWScudw8IRxNfB1wAXAFUgQPd/bFE8iOBbwOnAf8F7AQc4u53ZJUfEek+tbr3IHabGvTll5+a2VBNNOGJ6tzHt/m7+yBwVLS0eu8WwlHMjUYyi8gsFAaC6LISKyWURnLLUl2XX+jT0KylIlIcAaONyDW5dzX1AjRiTwMFAxEpBqM+2ij5/7ypZCAiMj2CHurPDrDwechQjIfddHt7QBoKBiJSHNGFP3zITVgkCHoKUDSYBcGgADFXRCQsAdTaBzxWMshd2jEGUxQwzOzFZvYlM7vHzDZGDwVbGT0OoFH65WZ2v5kNmdkDZvahNPtRMBCR3O350XOp9kddNAOjVC1/SnsGAAAPK0lEQVRAaSBi5N619E2Eg3u/SziI9yPAIuCXZvaqMXk1W044qeiPCB809p/ABWb24XY7UTWRiOSuOhCNPPawVGABlKO5CEqV/Otocm4zuAz4qruP5sLMbgQeBv6V8NECtZkcTgcucfdToqSrzOz5wOfN7GJ3b9pRVyUDEcldfOoJCwADq4TrykWYoCbHaiJ3XxcPBNG6Z4DfATvGVr+WsMRwaWITlwDbEz6JsikFAxHJX21sQTQCuVY6AKgO5JmxSMHmJjKzhcDLgN/GVtcmCk0+ZGxt9HO3VttUNZGI5K5UDeeusSggeMnp2RK+lntvomJONfHvhM0Z58fW1R4NkHx4WMsHh9WoZCAi+auN8I2qiUpVozQCrzninFyzNSp9yWCH2hMZo+Xo5KbM7A1m5imWmxplxcxOBt4LfMzd41P+16LmhEKXSgYikjuLBwMHq4ZzEo1sY1iQ/215B9NRrHP3vdqk+QXw0hTbSj4nnqib6BeAT7n7txIvx0sAT8TWL0y83pCCgYjkrjTiYFGX0lj9e7XfwtdylmU1kbtvBu7vOA9m/0A4I/Q57t7oEcS1toEljA0GtbaC37TavqqJRCR3vZvDQWelYUafkl6qOkE5LCXkKudBZwBm9jbCaf8vdvfjmyRbDaxj/HNh3kdYKri11T5UMhCR3P36Oyt4yefOw6rhxd+i29TysBdibqI8p6Mws/2AHwD3AN8xs71jLw+5+50A7j5iZp8mHGT2R+B6YCnhIwWOcffhVvtRMBCRYoguuBaEASEoG1bNf6K62gjkHC0F+oFXMP7u/hFg59ov7v51M3PgOOAE4FHChuYL2u1EwUBECqM0EgaCcjWamyg+lXWO8mzEdvfPAp/tIP2FhFNSdETBQEQKwRxKlbAk0LPFR6emyH3G0BnwsPs0FAxEpBCsEvYqqswxgrLRE83AUISnjBVw0FnmitA0IyJCqVp/3nGp4vUpKQpQTZR3b6LpoJKBiBRCaQTwcNK68nC9IbkIwWA2lAwUDESkEEpD9UdeWgDmTnnEqPbmnTO6/q4/DQUDESmEvo1O0AO9G8PeO14qTiAoQrvFVFMwEJFCCMphT6LejU6pAkGPE/RZOCo5RwUYZzAtFAxEpBDMgWgE8ugAtGrYsJw7n/nRQL2JRKQQgjKUR5xSxbntB8fhJauPOchZzs9AnhYqGYhIIZiHvYgqA2H3Ias6pRHLvzfRDOg2moZKBiJSCG6w+rLjRqd+sOgRmOUhZ9+3nZVr3mp5abd0M5UMRKQQkm0DFkST1Bn8/IoTcslTPC8znYKBiBSKW1gv1DMUUBkoc+sPm03fP10ZYlY0ICsYiEgh1O6+e7aG//GSUaoU4yLc7Y3DaajNQEQKYfTxllHJ4NYfHk95uCBXYc1NJCIyPfo3hI0GHus9VIQxBhp0JiIyjW75yYnj1uX5UJlR7sXIxxRTMBCRYomVDLyU9yCDyMyPBQoGIlIs8Wce1wag5U3VRCIis50DqiYSEZle8ZJBqZJfPsaY+bFAwUBEiqs8VIyhv6omEhGZbrEL7+jYg5ypN5GIyDTr3RxeePd/8xe55aqTcs4NM2JAWRoKBiJSKEFv2IOoOlCMCRLCQWczPxooGIhIodTmIyrUTKFFyssUUTAQkULJfZbSBlQyEBHJSWEabdVmICIiMDvmJipGC42ISEJh5iWC8OE2aZZpYGbvMTM3s8ebvL7czO43syEze8DMPpRmuyoZiEgh/ezKfB91OcqL05htZvOB84A/N3l9OXAhcAZwPXAQcIGZmbt/rdW2VTIQkUJauvSMvLNQV5ySwZnA3cA1yRfMrAc4HbjE3U9x91Xu/ingO8Dnzay31YZTBwMze4uZ3WJmG81sg5mtMbOlsdcXmNnFZrbOzDaZ2fVmtnuD7QyY2Vlm9oSZbTGz1Wa2X9p8iMjscOONJ+edhboCPOnMzPYB3gd8tEmS1wKLgEsT6y8Btgf2bbX9VMHAzD4IXAncDrwNeAfwn8Dc6HUDVgKHAMcAhwO9wCoze0Fic98ElgOfAd4KPAFcY2Z7psmLiMh0syBItUzZ/sO7+ouAs9z9wSbJlkQ/70usXxv93K3VPtq2GZjZzsD5wAnufn7spXgxZRlh1Fnq7qui960GHgJOBP4lWrcH8F7gKHf/drTu5iizp0bbEREpDqeTQWc7mNma2O8XuftFGeTiJKCfsC2gmYXRz/WJ9YOJ1xtKUzI4ivBQfL1FmmXAn2qBAMDdnwF+AhyWSDcCXB5LVwEuAw42s/4U+RERmTaGY55uAda5+16xZVwgMLM3RL2B2i03RelfBJwCfMzdt7bMamhCFVZpehPtC9wPvNvMPg28EHgYOM/dvxqlWcL4ogmEd/zvN7Nt3H1jlO4hd9/cIF0f8CLqRRoRkWLItnH4F8BLU6SrXSe/DNwI3Bb1JoLwemnR70PuvoWxJYAnYtuplQgGaSFNMHh+tJwFfBL4f4RtBl8xsx53/1K0s4cbvLe28wXAxihdsggTT9e0GGNmRwNHAyxevDhFtkVEMpJhMIhuhu/v4C27Ed6EN7p2rge+BHyc+o30EsYGg1pbwW9a7SRNMCgB2wIfcPcfR+tujNoSTjazLxMWTxodreSokbTpxomKWxcB7LXXXjN/OKCIFENnbQZT4d3AQGLdJ4BXEd6Y1wafrQbWAUcQjjGoeR/hDfetrXaSJhg8BfwtcF1i/bWEvYeeF+2o0V39guhnLaINAo1u6xfEXhcRKZSp7CnUjrvfllxnZh8grB66KZZuJKrKv8DM/kgYEJYStvse4+7DrfaTpgG5WR1+7W4+iNIsaZBmN+DRqL2gtq1dzGxug3TDQLMuUyIiOUk54KwAM5u6+9eBDwPvJOzx+R7Chuevtnwj6YLBFdHPgxPrDwYed/c/E44x2NHM9q+9aGbbAYdGr9WsJBx/8I5Yuh7gXcC17j6UIj8iItPHKVwwcPcPuHtyDFfttQvd/cXu3u/uf+vuF6TZZppqov8LrAIuNLMdgD8AbwfeBBwZpVlJWF91qZmdQFgtdDJh6eHMWCbvMrPLgfOjQRQPEUaxXQjruUREiqcgcxNNpbbBwN3dzP6ecLDD5wjr9+8HjnD370dpAjN7K3A2cAFhY8dq4EB3fyyxySMJ5884DZhPOM/GIe5+RzYfSUQkW3q4TcTdNxDOh9FsTgzcfZCwoeKoNtvaAqyIFhGR4lMwEBGZ5dyhOvPriRQMRETaUclAREQUDEREZjsHZsEzkBUMRERacvCZ32agx15KrnY/7ry8syDSmhM2IKdZuphKBpKre885Nu8siLSnNgMREVEwEBGZ9YoxCd1UUzAQEWnFgRynsJ4uCgYiIu2oZCAiMttpOgoREXHwWTDOQMFARKQdjUAWERG1GYiIzHbu6k0kIiKoZCAiIo5Xq3lnYsopGIiItKIprEVEBNAU1iIvO15TTMvs5oAHnmrpZioZSEv3na0ppmWW89nxcBsFAxGRNmZDA7J5F3aZMrNngQfyzscstwOwLu9MiM5DGy9090WT2YCZXU14nNNY5+6HTGZ/eenWYLDG3ffKOx+zmc5BMeg8SFbUgCwiIgoGIiLSvcHgorwzIDoHBaHzIJnoyjYDERHJVreWDEREJENdEwzMbCcz+6GZPWNmG8zsx2a2OO98FZWZvd3MfmRmj5jZFjN7wMzOMLNtE+kWmNnFZrbOzDaZ2fVmtnuD7Q2Y2Vlm9kS0vdVmtl+DdCUzO9nMHjazrWZ2t5kd3iSPy83sfjMbivL3oeyOQHGZ2dVm5mZ2WmK9zoXkx90LvwBzgd8D9wF/DxwG3Av8P2Be3vkr4gLcBvwHcASwP/Bx4OlofSlKY8DPgMeB9wCHADcT9lt/QWJ734vevxw4CPgxsAXYM5HudGAIOB44ELgQCIC3JNItj9afHqU7Lfr9w3kfuyk+L+8BniCc5eC02HqdCy25LrlnIFUm4V+BKvCi2LpdgAqwIu/8FXEBFjVY9/7oIrQ0+v2w6PcDY2meAwwCX46t2yNKd2RsXQ/hwL+VsXXPjS4+n0vs9wbgnsR7/xv4biLdt6KLX2/ex2+Kzsl84M/RxT4ZDHQutOS6dEs10TLgNnd/sLbC3R8CbiX8I5IEd3+ywepfRz93jH4uA/7k7qti73sG+Aljj+syYAS4PJauAlwGHGxm/dHqg4E+4NLEfi8FdjezXaLfXwssapDuEmB7YN92n69LnQmsdfcfNHhN50Jy1S3BYAlhFVHSWmC3ac5LN9s/+vnb6Ger47rYzLaJpXvI3Tc3SNcHvCiWbgh4sEE6qJ+rJdHP5L6T6WYMM9uXsGT2kSZJdC4kV90SDBYC6xusHwQWTHNeupKZ7QicClzv7mui1a2OK9SPbbt0C2M/n3b3ZH/lRulosM1kuhnBzHoJ6+vPdvdmc2rpXEiuuiUYQFhPmmTTnosuFN1VXknYxnJk/CXSHdepSEeTtDPRScAcwgbaZnQuJFfdEgzW0/gOZQGN75IkYmYDwEpgV+Bgd3889vIgzY8r1I9tu3SDsZ8LzCx5wWmUjgbbXJh4vetF3Z9PAT4N9JvZfDObH71c+72MzoXkrFuCwVrqdZtxuwG/mea8dI2oeuJHwKsJuxPem0jS6rg+6u4bY+l2MbO5DdINU6+XXgv0A3/TIB3Uz1WtPjq572S6mWBXYICwgXZ9bIGwy+d6YHd0LiRn3RIMVgJ7m9mutRVmtjOwT/SaJJhZibA/+kHAYe5+W4NkK4EdzWz/2Pu2Aw5l7HFdCfQC74il6wHeBVzr7kPR6qsJL0hHJPbzPuC+qAcYwGrCbouN0g0S9hKbKe4i7LufXCAMEAcSXsB1LiRfefdtTbMA8wj/YO4l7Ga3DLgb+AOwTd75K+ICfI2oLzuwd2J5QZSmBPwCeAx4N2F3xJsILwI7JbZ3GeFd7D8TBpgfAluBVybS/Vu0fgVwQJSPADg0ke5D0frTonSnRr9/NO9jN03nJznOQOdCS65L7hlInVFYTFjlsQF4Fvg/wM5556uoC/BwdMFptHw2lm4h4QCjQWAz4aCkPRpsbw5wLuGgqa3AL4EDGqQrA58CHiHs2ngP8PYmefwg8Lso3e+Bj+R93Kbx/IwJBjoXWvJeNGupiIh0TZuBiIhMIQUDERFRMBAREQUDERFBwUBERFAwEBERFAxERAQFAxERQcFARESA/w9n2TZq1lIXYwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dbz = data.dBZ.T.values\n",
    "dbz[dbz<-2000] = np.nan\n",
    "\n",
    "# data.dBZ.where(data.dBZ<-2000) = np.nan\n",
    "plt.imshow(dbz, aspect='auto')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fa852d2ba90>]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb8AAAEBCAYAAADo06t3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXm8HFWZ//9+uvvu+36zkoSwJQMEjSyiBtEZcAOdURFlZoSfIMqIfEUcEWZGBYQZcdRRmYFxGYddEDUwgAyCGcSwREKAsIYkJOTebHffl+7z++NU3Vu3b+9d3V3ddd6vV7/q3qqnqk5t53Oe52yilMJgMBgMBj8RKHQCDAaDwWDIN0b8DAaDweA7jPgZDAaDwXcY8TMYDAaD7zDiZzAYDAbfYcTPYDAYDL7DiJ/BYDAYfIcRP4PBYDD4DiN+BoPBYPAdoUInIFsCgYCqqqoqdDIMBoOhqBgdHVVKKd86QEUvflVVVYyMjBQ6GQaDwVBUiMhYodNQSHyr+gaDwWDwL0b8DAaDweA7jPgZDAaDwXcY8TMYDAaD7zDiZzAYDAbfYcTPYDAYDL7DiJ/BYDAYfEfR9/PLhulwhFDQW/o/FY7w4At7+eAxC9h+cISDQxOcsKIlrWO80TPCum//PqFNWVC46N0r2Tc4wfYDwzy5o5flrTUcs7iBzvpK/vbty1jY6M7gAd+890V++viOlGwvee9hCEIoKAC867A2jl7c4Eo6Nu/q41M/fpLRyXBcm9UL64koeNfhrTy7q5/B8Wle6h4E4PL3Hcln1x2adTq27R/i15u7+OGj2/jL4xZxz+Y9c7ZXlQX5/CmHooCe4Qm6BsY5elEDAqxd1sxJhyZ/H3pHJtm8q4/3HNUxZ/3g+BTHfP2htNP8lqWNnHHsQv727csQkZg2f9x2kBVttXQ2VKZ9/JGJaX7/ygE+cMyCtPf1BENDcMMNsGwZLF0Kzc0wNaW31dTA5CRUV0MgoH9lZdDamvZptu0f5r3/umHm/x3Xvj/u8zAkRpRShU5DVtTU1KhMOrnvPDjCKdf/nn87+zjOOHZhDlKWGd/931f5/u9e4z/OeQsX3vIMADuv+0Baxzjzh39gy5sDWaXj7969ki+fdkRWx7BZ9tX/yXjfU45o47/OPd6VdKz79qO80TOa1TFeu+Z9lGVZYLrotmf4n+e6M9r3yM46HrzkXUntPvyjx3l2dz/Pff0vqK8sm1m/fksXF9++OaNzA/zu0nUc2lYbc9uyr/4PTdVlbP7Hv0j7uBffvpn1W7q4/+J3smphfcbpKxjXXw+XXZbePt/9LlxySVq72M/V5pFL17EizvNIhoiMKqVqMtq5BPCt5/fa/mFAf3QX3745bYHJFfsGxwHoG53K+Bh/tqhhRvwWNVZx52dPpKWmgjf7RqmrLCMQgMaqckRgcjrCxHSEiekwNRUhqsqCHH7lA0xH3CsUHdpWQzii+PoZq1ncVMV0RFEWDFAWCFBZHqCpupxwRDEdUVSGAiggohQf/feNhF1Mx1+s6uA/H9vBZacdweKmKk45op03+0ZprilHEKrKgtRUBIkoENFeeECEsckwP3hkGz99fAcRFwqLoYBQXxni/i++k/a6SvYNjvP8ngEaqspor6tgSXM1ocBsaX46oggFhItue4Zt1nubjJ09ukAYDs9NbzgSAeC2z5xAfVUZS5qq2XZgmLKgsKSpGoDhiWnKggG2HxzmrYc0EYnAfc91cdndzxFJ8jwyfW/39OvBRkYnpzPav+C88QYEg/DMM7BnD/T2Qnk5KAUjI9DVBR0d+v9IBC6/HF54Ie3TnLiiZY74tdRUuHkVvsK34vf5W/805/+xyTBV5cECpcZdnFGQ45Y2stjK1A7rqJtnWxYMEP39ZOvZRBMQ4YiFdZxyRHtcm7IYt74s6G44R0QL3EXvXjmzrqEqfkjVvg+VZUFaasvdSwfQVFM+81yWNFezpLk6rn3IujehQMC1ENfCxiqWtepC/1sPaZqzralGX6szfFkq30bOeP11OOaY2V8yfvhD6OlJ+zShgBAMCK9/6/0ZJNLgxFsVXnlkKqpEPD4Vvx6omAmY+gDPUdwVDYaYvP46HJpGfXBzM/T15S49LiAii0XkByKyUURGRUSJyLIom7UicpOIvGzZ7BKRW0VkeYzj7bSOEf37cAzb861jTojIKyJyYZw0flhENovIuIi8ISJXikhKJTXfip9fCBjtm8FL9dvmsZQQ09OwfTscdljq+zQ16dCot1kJfBzoAx6LY/MJYDXwb8D7gK8CbwE2iciSGPa/BU6K+m1wGojI+cCNwC+B04G7gBtE5HNRdqdZNk9b5/4+cCXwrVQuzrdhz5NWtLBxe/phh2Kj2D0/t/WqyG9H2kTfPg/pf+mwc6cWwHTEr7kZNm1K+1Qqv3GD/1NKdQCIyGeAWC2Z/lkpdcC5QkQeB3YA5wP/GGV/UCn1RLwTikgIuAa4WSl1hbX6URFZCFwlIj9WStkVy9cBf1BKXeCwqwWuFJHvKqX2Jro433p+0XUYpZonFHMzaC+mvVjEI9md8+CtLV5ee00vDz889X2yCHvm69EppSIp2ByIse4N4ACwKIPTngS0AbdErb8ZaAHeAWB5lWvi2JWhPcGE+Fb8giUcD3Rm0CaTcwc372OxCKghRWzxW7kysZ2TpiYYHYWJidykqYCIyFFAO/BSjM0fsuoGJ0TkiRj1fautZXRT2K3WclUiO6XUDmDUYRcX34pfa60/mgiXsManjZdEx4terSFDDhzQpaO2ttT3aW7Wy8I2egmJyCbH74LkuyTGClv+B9rz+0nU5nuBLwCnAZ8CxoFficg5DhvrxhB9Y3qjtsezs9c1x1g/B9/W+V35gaO4/aldhU5GXNzKqIu9zs9t/HI3PKTzKeGlxkhp09enPblAGr5Ek9W9pLcXOjtzk67kTCul1rp8zB8Cbwc+oJSaI0xKqS84/xeRXwFPANcyG760P9FkL0Qiu5Q+c996fjUVc3W/qD++BBgPw98U29Mvyte1t3dWzFIlQ8/Py9mUiFwLXACcp5RKOoaeUiqMbsm5WETsce2iPTyi/u9NYgfQ6NgeF9+Kn9dxKxMo9rBnnlu35YV8XpFp7ZkH+vpmxSxVnJ5fmnixgCAiV6C7OXxRKXVzOrtaS/vNtOv2VkfZ2XV4Lyays/ohVjvs4uJr8Tv7+FjdUEqLYg57up3ybPJ9cTk1uX4qSVt7Fp1P6GHssGc6eKPOzxVE5GLgauAKpdQP0tgvBHwM2OXolrAROIiuE3RyDtqbexxAKbUL2BLHbgp4INn5fVvnB/D2Q1u5/andQPHVkaSKFzw/L91bEwY2ZMzwsB7AuqUF9u2DJUv0DA5PPglnnZXesbLw/PKJiHzU+vOt1vJ9InIAOKCU2iAinwC+BzwIPCIiJzp2H1RKvWgd52zgTOB+YDfQAVxkHfdsewel1JSI/AO6U/se4GHgVOA84AtKqUnH8b8G3CciNwK3A8ehO7l/P1kfP/C5+JUqTrHxSmZvPA1D0fPQQ/CNb8Telk43B4CGBh2/9Lj4oevknNxgLTcAp6BHYBFreXqUrW0DutN7O/BtdD3dKHpkltOVUr917qSU+g8RUcClwGXALuDvlFI3RNndb4nzPwGfBvahR3e5JpULM+JX4hRz2LNUKebGVcWbchfo6tLL55/XMzVUVMwOTn3SSekdKxiExkbPhz2VUgkzEKXUp9HCk+w4T6A9uFTPeyN6iLNkdvcA96R6XCdG/EocL4Q9s8FNnXDjWK6lp0DPJdPkG88dLX6hEKxalV63hnhkML6nrwsfLuPrBi/OF6mIC+MJKWbHLxdpz/SQxXof43mZxXo9BaWrCxYscEf4IOMhzkxBxB18LX5+IFDsrl8Jko9yllfqelOlKMqetvi5RXHM7FCyGPHzKG55oi/sGXDnQAZXybU0FW+9oodFu7sbFi5073hFMKdfKWPEz8IrnandLrBv2mk+LhuvPONsSFfUis0D9DRdXe6Kn/H8CoqvxS8S8V5m6HajDJP3RVHM98OFtHvfI/Ro+kZHtVAtymSWnjjYnl8az8Tzj6+I8LX4DYxNJTcqEO4Nb5bZgbziJXntY3flvuTxmrwvdhrPl0nsbg5uil9TE4TDMDSU3n6ev1nFga/Fb/XC+kInIeeMTYXT3scr35aXWrW5nZJchyNNuNNl9loDhrjZ4MUe4syEPguCr8Vv7bJmvvb+I/U/xVFATpsiKfjnBXMvDBmzb59eujn9UAmN71mM+Fr8YP7URqXGX77FxTBNCWD8IUNG2J6fm+JXJON7liqlnfP7mNbaCp762ns80c+vWOqd8oVX6lMNabB3rx6SrKXFvWNm4PmZd8c9fO/5lSoiHuvgnmFSvPaxu6XjhXoy3rqbs3g1XTPs3Qvt7VoA3SJDz89DX3VRY8TPwmsfX3aZrNeuJkM89JUXa/uReG+Cd6/Hownbu9fdkCeYOr8C43vx81KLQnAvU/LWVXkHL7SCNFHgIiQX4ldVBeXlps6vQPhe/AyGQpAvDXb7NL4V7r17oaPD3WOKmCHOCogRPwuvfNReSYfB4MQDDnPhiER0Vwe3PT8wQ5wVEN+Ln1c/aq+mq5hxo9Vp0ZdNiv4CCkBfH0xN5Ub8mpvTEz/z/FzD9+Jn8DZue8KZFircrBvOp3cfv8GLKV2lTC76+NlkEPY0j84djPgZPEspf+O5bmhVyvcu7+RS/EzYs2AY8bPwWn+ybDD1hgaDi3jM8zO4g+/Fr1RLyCY0Mh9TJvAuni6w5drzGxrSdYqGvOJ78TP4Cy+UCUopyuA2niy07d0LlZVQn4NZYOyO7v39KZmbN8c9jPh5FE+XhPOI126DW+OU5jqTj5dKI7wZYHdwz8VDy2CIM68NzFGs5Ez8ROQUEVExfv1Rdk0i8mMROSgiIyLysIgcnat0xcMrYuPJkm+B8NK98FJa3KDELie35GJ0FxszxFnByMesDhcDTzv+n7b/EN3eej2wHPgC0AdcDjwqImuUUm/mOnGllql5EY+UKzxTwMmGVC/BvNYusncvrFyZm2ObaY0KRj7Cni8ppZ5w/DY5tp0BvAP4a6XU7UqpB611AeAreUib5yiFDDoWXsmMvdC/LdNnXPiU5x5Pvv/d3b70/ERksYj8QEQ2isioFblbFsOuUkS+LSLdIjJm2b8rhl1ARC4XkZ0iMi4iW0Tkr+Kc+3wReVlEJkTkFRG5MI7dh0Vks3W8N0TkShFJaeqNQtf5nQF0KaUetVcopQaAe4Ez85kQr31z2eTRnsxADIYEeKBMEpuJCejpgUU5mhTaFr8UPb88z425Evg4OiL3WAK7nwDnA/8IfBDoBn4rImui7K4Cvg78EHgf8ARwl4i832kkIucDNwK/BE4H7gJuEJHPRdmdZtk8bR3v+8CVwLdSubh8hD1vFZFWoB/4LfBVpdQua9tq4IUY+2wF/kZEapVSw7lMXKlWHpfqdRWaYi9XmIJRmnR36+XChbk5fmOjXqbT4CV/n/b/KaU69DnlM8BfzE+LHAt8EjhPKfUza90GdB7+TbSDg4i0A18GrlNKXW/t/qiIrASuA+637ELANcDNSqkrHHYLgatE5MdKKbtfyHXAH5RSFzjsaoErReS7Sqm9iS4ul57fAPAd4DPAqWjVfy+w0boRAM3oUkU09pvQFOvAInKBiGwSkU3T09OxTAylgouZtVdaOnohFZ71tLzGnj16mSvxC4V0FwoPhj2VUpEUzM4ApoA7HftNA3cAp4lIhbX6NKAcuCVq/1uAo0VkufX/SUBbDLubgRZ0NRkisgRYE8euDO0JJiRn4qeU2qyU+rJS6l6l1Aal1PfQLmwHuhEM6GqMWHlBwk9TKXWTUmqtUmptKJQP59WHeCCHzoX36pU8P191j257el4pQOSNri69zJX4QbEPcbYa2KGUGo1avxUtdisddhPAthh2AKscdjA/IpiSnVJqBzDqsItLXpVDKfWMiLwKvM1a1Yv2/qKxPb68FYfyHEv3NMYrKH6K7Rl69vOzxS9XdX5QyCHOQiLibIB4k1LqpjSPkSx61+xY9qv5GW0sO2IcM1U7e10sXZlDIdwmp7e3lRhxZLRq78p1fd9MajyIZzMDgy/J12fiOdHu6tKzrTcnzUszp3Ce37RSam2Wx0g1epeOHXFsU7VL6S3Ka2tPEVkLHA48aa1aDywSkXUOm3rgQ9Y23+G5j7+EcKNA4ZVjZHoOU6ZKkz17YMGC3H6YaXh+HiwUJ4ve9TqWTTI/3h/LjhjHbE7RDqDRsT0uOfP8RORWYAfwDLql53HoDux7gB9YZuuBjcAtInIZs53cBfiXXKWt1Cmlehm3ryXj+fxczvwK7UmZ1sApsnMnHHJIbs+RpufnsSe3FfiIiFRH1futAiaZrePbClQAhzK33s+um3vRYQe6Tq87RbuNtpHVD7HaYReXXHp+L6BbAv0M3cXhEuAe4ASl1EGYaU30QeB/gRuAXwFh4N1Kqd05TNs8PFiiyopS8CBL4RoKTam913llchK2bIHVq5PbZoPt+RXnw1qPbl35MXuF1V3hLOAhpdSEtfpBtBh+Kmr/c4AXrIYqoIXsYBy7XuBxAKu73JY4dlPAA8kSnjPPTyl1LXBtCna9wHnWL++Y/NVQipiCgws8/TQMD8Of/3luz9PcrIV2dBRqanJ7rjQRkY9af77VWr5PRA4AB6xW/M+KyJ3A90SkDB3t+xx6yMoZYVJK7ReR7wKXi8gQOiJ4Frob3JkOuykR+Qd0p/Y9wMOWzXnAF5RSk47kfQ24T0RuBG5HRxevBL6frI8fFKbBiyEBxVn4Kw68c2u9kxJDAp59Vi/f9rbEdtniHN/TY+KHHl3FyQ3WcgNwivX3ueiO6Vej69u2AKcrpZ6J2vcKYBj4ItAJvAJ8XCl1r9NIKfUfIqKAS4HLgF3A3ymlboiyu98S538CPg3sQ4/uck0qF2bEz6OUVMndU3m9N25sST3fUmXLFmhpyW03B5g7vueSJQlN8/0pKaWSvqlKqTHgS9YvkV0YLZBXp3DMG9FDnCWzuwddnZY2hR7bs+B4YaBjP1Ay99lTQp4+JrKQBlu36vq+XL+7ac7sUDLfUoHxvfgZvI1XMms3sxsvXJPX8k8P3JL5vPIKHHlk7s/j4ZkdShkjfiWIFzJXN3A7g/bSfcmX+BRbtxfPaPLBg3o2hyOOyP25zJx+BcGIn4WXMkbIPj2eyUQ8htc8HoNH2b5dL3M1ia0T4/kVBN+Ln9fyQpM5G9zEdGbPEFv8VqzI/blqayEYTMnz81ohvZjxvfgZDOngRhgxqyNkefpiC4MWDFv8li9PbOcGImkNcWaKM+5gxM/CZArexN2nkvnR3PbIM/HI3GzlZzLQJOzYAR0d+et3V9zTGhUlvhc/E2b0LqU8n1+x4puw2/bt+fH6bJqbjfjlGd+Ln9dwZdaA7A9hKBHcimjkr3WqR9i5M//iZxq85BUjfh4l28zGdIT1LvmZOLm4nr+nUhsOw65duZ/NwYkJe+YdI34Wvgnn+BgvzcVnyiZz8dTn190N09P5Fb8UPT/TNsE9fC9+JhPyFxnP5+duMgqG1wt5nrjPXV16uXhx/s7Z1AT9/drrTIYnblLx43vx8ypez6TSIZtLyU+I0IeYDDQ+/f162diYv3PaHd0HBvJ3Tp9jxM9jeMUTdVtyMrksr9wLt8mPnJtCQ8bYApRP8TNDnOUdI34WJquYpVRHBfGSE1mad7hEsMWvoSF/5zRDnOUd34tfKWb0XsrkvUa2z7t4bm3pvdd5oxDil6LnZ75t9/C9+BkM+aaQGZjJO1NgYEDH3Gtr83dO2/NLIexpijXuYMTPwjSs8CZeeSqu95sscIVmKUY8XGNoSAtfII/Zowl75h3fi5/XGlUYDc4dpo+UISUGB6G+Pr/nNA1e8o7vxc+reE2USwVzXz2Kl0p9Q0NQV5ffc5aX60G0jeeXN4z4GQx5Jq/ZvIc0JRU8MSxfITw/MEOc5RkjfhZFlkcYCoRbdcO5zuK9oCFFSyE8P0h5iDNPFBBKACN+HiWbPNbUbcUmm3taMvlNluLtpehkziiU+BnPL68Y8fMYbmWypZJZu53ZlshtSUqy+5b+++GXO0fhwp5mWqO8YsTP4FlKNbyT1241pXkLc4vx/HyBET8LX4RzDJ4hU103r2mOUaqwnp+HxE9Efi8iKs7vQctmWQKbxqjjVYrIt0WkW0TGRGSjiLwrxnkDInK5iOwUkXER2SIif+X29YXcPmCxUarehZfwygACbqSikFeS0ZvqjVtfPExM6Ln8CtXgZXwcxsagqiqmSZ6/pc8D0aWAk4B/BdZHrb82xrqhqP9/AnwAuAzYDlwE/FZETlJKPeuwuwr4MnAF8CfgE8BdIvJBpdT9GV7LPHwvfob84JUyRqaFHY8kP2XiXabRwiQMDuplocKeoOv94ogf5O9bUkq9OP/ccj4wCdwRtWm7UuqJeMcSkWOBTwLnKaV+Zq3bAGwFvgmcYa1rRwvfdUqp663dHxWRlcB1gGviZ8KeM5hswYuU6lMptJgW+vzReOY5D1nOSiHCni0tetnTk/9zp4CIVAEfA+5VSqUbnz0DmALutFcopabRInqaiFRYq08DyoFbova/BThaRJZnkvZY+F78PJcJxMgF0g51eCYnyQ6vPRuDDyik59fRoZf79uX/3Knxl0Ad8PMY264VkWkRGRCR9SJydNT21cAOpdRo1PqtaLFb6bCbALbFsANYlXHqozBhT4+SbWjDK2FGL+GRqkdDDDzzuhbS8+vs1Mu9e/N/7tT4G2A/8IBj3QRwI/AQcAA4Evga8EcROV4p9ZJl1wzE6sfR69huL/vV/BJ/tF3WGPEzGNLADQE1IuxhCun5pSB+Lr86IRHZ5Pj/JqXUTbEMRWQh8F7g+1a4UqdHqW7gQofpY1ZL0K3oBivn2IcgdvKjyz2p2mWNET8Lr2VIXkuP73HZlS5UK2PzXiWhkJ5fXR1UVsIrryQ0c/HNmVZKrU3R9hx0NVmskOcclFK7ReQPwNscq3uBpTHMmxzb7WWTiEiU9xdtlzWmzs8z8RaN19JjKE1MF5842OJXCM9PBN76Vnj00fyfOzl/A2xRSm1J0T7ag9sKLBeR6ii7VejWo9scdhXAoTHsAOa1QM0U34tfMeDr0rqLF2/GPDUkpZBhT4ATT4Tubk999CKyFt0QJanXZ9kvBU4GnnSsXg+UoVuL2nYh4CzgIaXUhLX6QbQYfirqsOcALyildmRyDbEwYU+DZ8mFc+IFh8eIsIexPb/a2sKcv6MDRkdheLhwAjyfvwGmgduiN4jId9BO1EZ0g5cjgMuBCPAt204p9ayI3Al8T0TKgB3A54DlOIROKbVfRL4LXC4iQ8AzaIE8FTjTzYsy4mdRStmRKyOZeKjkWYp4QIMNsRgc1MIXKFBQzNndIYb45fuztITqbOBBpVSsPhhb0SL2aXQ3iIPAI8A3lFLRlZfnAtcAVwONwBbgdKXUM1F2VwDDwBeBTuAV4ONKqXvduCYb34uflGg2lM11ecE78ip+99pK/vqHhgrT2MXG2eJz5cqYJvmsr1VKTQFtCbb/FPhpiscaA75k/RLZhdECeXXqKU0fU+dXBJR4dpM/spnPz71U5JXoS87Uo891fuuZd3xwsLDhRtvz825fv5LBE+InIktE5G5rdIBBEbnHqjTNG16J8nklHW6SzSW5fTu84NXm4xknu0wP3IaYFDxdAwPQ0FC48y9cqJddXYVLg08ouPhZTV8fQY8M8LfAXwOHoQczrcn9+XN9hswodLq8MImsRx+NKxT6+RriUGjxa22F8nLYs6dwafAJXqjzOx9YARyhlNoGICLPAa8Bn0VPn5FfJifha1/T04vU1UEwqFtgBQI61woE5v5dVgahEFRYY7MqpY9RV6e3B4MQDs/uEw7P7hMKwdTUzLHWPr2L8M4+llW8yllbdhBQCvlJN5SFoKZGj/ZeWan7A9mjwLuMyZcNOcXuTlBWpr+ZQjUuicXAACxZUrjzi2jvz4hfzvGC+J0BPGELH4BSaoeIPI5u2pp/8XvySfjOd+auq7Gc0EhE/5TSy3DYVTfpo9aPB/TEWQD8NobhJz8Jt97q2nn9gCtPqUiGN/Ns9PzGG+HCC+euCwS4OxBiWoSKH4QgFJwtNFZU6G8vENChwPFxXWAsK5v9iWhvLRDQBcOqKv0LBPT+lZWzy5oaXSitr9c25eV6vQhUV8PLL8M731mYe2OzaFFc8Sv5Bkd5xAvitxr4TYz1W3F0iMw1c14qe1T1LVvg6KN1bhWvdKqU9txAf5gi+ldWpluO2WIZDM7+bXt7U1Oz2yyP8J8ffIl7n+3iK+89jGsf2U5EhMe/vI7Q9JSeaHNsDL7wBdi+Pf61lGLFoUtk2grW7TBloVoZF/zN+OMf9fL66/X7PzkJk5Pc++QOegbG+PCahbRXh2YLluPjOuoSDsMpp+hox/T07Pdjf0MDA/pbtCeDHR7W/09M6HX2cnh4dls8jjkmL7ciLosWwTPRrf9nMZEZd/CC+CUa7Ts3cT0HMV+k/fv1sr19VsziHkB06RFmlzZ2GDQN+tsHebMxwvCSZXTXWx1uly6FoEN8Dz1Ue6cJyDazLngmacgpBatz3LUL3vEOuPTSOat//sM/8NybA5xw0cm0L2nMbRoiERgZ0WI4NaXFNRLR/9fWwnLXpozLjMWLYf16LdCmcjhneEH8IM1RvEXkAuACgPJowXEDW/xaW90/thu0tOR0wksvjftonNi5pOvVe+7+9fbCsmWFTUMgoEOf3hlBZS5Ll2oh7unxbh5UAnihprmP2HM0NRHbI0QpdZNSaq1Sam0o5I5+z8kkDhyA5mYdnvQiLS3Q36/DPyWM2yLslXBwpvU26dwO7xRfoujry1lDrZLBbnCze3dh01HieEH8tqLr/aJZhYsjeMcjZoayf78OeXqEeVllS4te9ro2u4dv8Ew42LPqlGPiiJ9HyiXewBa/XbvmbTL3yT28IH7rgRNFZIW9QkSWoUcFX1+QFHlM/OZhh0JyGPo0H9lcSnUYvLwyNaUbmyTw/DwUcS8cSTw/c4/cwQvi95/ATuA3InKmiJyBbv25G7gOqUzCAAAgAElEQVSxICk6cADa4g5nV3hsz+/gwZwc3nxbpUtBCzX9/Xppwp6JaW/XjedM2DOnFFz8lFIj6OkqXgVuBm5FT3dxqlJqOH/pcPxTQPFLKXOyxS+O51dKTpub/Zq8cl+84FUXxJPts6rwm2NV8Wu8cG8KTiCgW3wa8cspnmjRoZTaBfxVYc4elQmEw1pUvOz5pRD2LAXvLRfX4JX7kq90eKpTtF1HHcPzM6G8KJYsiVnnZ3CPgnt+nqOvTxc/C9TEOFYmMK80nOOwJ7jscXko/zVkR1bP0vb8TNgzOUuWxPT8zKfkHkb8ojlwQC+93L+mulp3oM9Vg5eczKBeGkV7vwq5K0/PiF/qLF2qhzgLh2NsLI1vqdAY8bOY8XRs8bPn1fIiIjnv6F6KZCNcbmp3IfWzoNptxC91lizRwmfm9csZvhe/eZmac2gzL9PamtOwZ6niFQ+04MkoxPmN+KVOgr5+BnfwvfjNw4PiF7P+Lceen1fCe15Jh8EF+vp0yD7GkITmOUdhRnnJOUb8LGY+vv37Z8OKXiaB+GWbkRTaKbEpuHeUK/KY0XtKVHp7k3p9JfvM02XpUr00nl/O8L34zfvW9u/XwuLVcT1tkoQ9vRLeKzXcagWb6352+X78KY2basb1TJ2GBt2ozY5EWXiqMFPk+F785uH1oc1sWlt1SToSKXRKioZs8o1SKUoUdHBvI36pI6LvVd/8sf1NudYdjPhF4xHxS+phtLZq4bOHjDKkhMk3NAXJQPv6Eo7uYogijvgZ3MGIXzQeET9nWCxmYd3uh5iDFp9eCpmaME9mePK+Gc8vPYz45RTfi9+8jN4j4peUPIzyUni8I8IGFxgY0HVZhtQw4pdTfC9+c5ic1C9bMYhfDj2/UsWN+i5PelQJ8ExyIxE9nVF9faFTUjzEFL/8PVEROUVEVIxff5Rdk4j8WEQOisiIiDwsIkfHOF6liHxbRLpFZExENorIu2LYBUTkchHZKSLjIrJFRFwf+9mInxNbSIpc/DyT4XmRDJ1Jr0SCU322XknvDCMjuuRgxC914jV4yX9KLgZOcvzeO5MWHTpbD5wOfAE9QUEZ8KiILI46zk+A84F/BD4IdAO/FZE1UXZXAV8Hfgi8D3gCuEtE3u/mRXm8PX/+UAo44L0O7nFJ4vll+4G42SrQUzMLFDFe07O0GBzUy7q6mJvNOxKDpiYdKg6HIRgsZEpeUko9EWfbGcA70FPQPQogIhvR09J9BS2ciMixwCeB85RSP7PWbQC2At+0joOItANfBq5TSl1vneNREVkJXAfc79ZF+d7zm5OheHB0l7jU1ORscOtceA1FnXGXIHl/HkNDepnE8yvIPINexW4ZOzBQ2HQk5gygyxY+AKXUAHAvcGaU3RRwp8NuGrgDOE1EKqzVpwHlwC1R57kFOFpElruVcN+L3xyKSfxEfDG+p5v+gPEtCojt+ZmwZ+rYLWML3+jlVhEJi0iPiNwmIksd21YDL8TYZyuwVERqHXY7lFKjMezKgZUOuwlgWww7gFWZXkQ0JuxpoVDFJX6gxc+ehaIEMR6od0jW/SVplNyIX/rEEL88N7gaAL4DbAAGgeOArwEbReQ4pdR+oBnYGWNfa+ZimoBhyy6Witt2zY5lv5pf7xJtlzW+F7853/T+/XrQ3WL5QNvaciZ+xkuKTbHdl4KO6OIkSZ2fIQZxPD8XC4UhEdnk+P8mpdRN9j9Kqc3AZsf2DSLyf8BT6Lq8K9HlyVgvWXQq3bbLGt+L3xzsPn4FbCoXK6+Km3+1tcGOHa6nwXhH8ym2uqh46S2YFqZY52dwYItfb29iu8yZVkqtTWcHpdQzIvIq8DZrVS+xvTF7NIM+h93SBHa9jmWTiEiU9xdtlzWmzs9JsXRwt2lvnzfwLXiotO81fHRbkrWezPsoPimGPU2rTwd2g5fC1/lF4/TOtqLr6aJZBexSSg077JaLSHUMu0lm6/i2AhXAoTHsAF7MIt1zMOJnoRSeEL+08qT2dl2iHh+PcSDXklRQ3BZyLw3d5iuShD2LzbPOC7n3/NJGRNYChwNPWqvWA4tEZJ3Dph74kLUNh10Z8DGHXQg4C3hIKTVhrX4QLYafijr1OcALSinXQl2+D3vOq/M76qiCpSVt2tr08sCB2ckvXcILzqPJDkuIoSFdn15RkdzWoKms1JP/Fkj8RORWdH+9Z4B+dIOXy4E9wA8ss/XARuAWEbkMHea8HP35/ot9LKXUsyJyJ/A9ESmzjvs5YDkOoVNK7ReR7wKXi8iQde6zgFOZ23Uia3wvfjMo5QnPLy3stO7f76r4Ge8oPiaknCGDgwlDnibcGYfm5jnil+fX7wXgbPTILdXAXuAe4J+UUgd1elRERD4IXA/cAFSixfDdSqnoaejPBa4BrgYagS3A6UqpZ6LsrkC3EP0i0Am8AnxcKXWvmxdnxM8iMDoCY2OeFL+4GYPT8zMkJasMtkjLA57R6iTiZ2PCn1FEiR/k7x4ppa4Frk3Brhc4z/olshsDvmT9EtmF0QJ5dcqJzQBT52cRLKZxPW2cnp8hJfyStcZt7VkoD2toyHRzyIQY4mdwB9+Ln51JhHos76kYxS8Hnp8JQ5U2bhcCkr4tKXp+hihaWoz45Qjfi59N6GARil9dnW5EEOX5ZStbfvGODHnEiF9mNDfnZPxegxG/GbwsfnHrbUTi9vXzkoBlWu/kdrsbN+q/XDmGH71qI36ZYYc9rRfPl+9OjjDiZ2WwoR6rzs9uRFIstLcXR4MXj6hxpoLqdvJz3aDWc5mkqfPLjOZmPcn26Ox40KYxtjsY8bMIHdyvS6aVlYVOSnq0teWkwYtnWgkaXKVgz9V4fplhj/Ji6v1cx4ifRVnPQU+GPJOSA8/PlCxLn7w+4+lp3Y3IiF/6GPHLGUb8LIIHi1T8cuT5eQU3PRU/ebO56guW0T1MYVBrPz2btDDilzN8L352FlHWc8Cz4pcwX2hv1/UBIyP5Sk7eyEUGXuydqAslElndtTSmMzJRhyiixM8UEtzD9+JnEyrmsCfMDX260iLR4DWKdtg5M5Ft5sTw/Ir0LfAcRvyAQCRMqLenOMXPbp0aFfrMLqM0n1c0RSs8XsDM5Zc5JuyZM3wvfiJC4/gwEokUp/jlcJSXUsNzzf8LQK7uQMIBv43nlznV1boFuhE/1/G9+AG0jPTrPzwqfgkzFjO+Z1r4zYGL9+rkte4zjTo/QwzM+J45wYgf0Do6oP/wqPglJE7YM1u8UrFuvLXM8JTIG88vO8wQZznBiB/QPmy9WJ2dhU1IJtTU6NCIi2FPr2ScXkmHE68UCpLhqXSaOr/scHh+XnqsxY7vxU+AlT1vooJBOPTQQicnM0q8r59bZCMIHtTh4sH2/GprC5uOYqW5Gewp1zCNr9zC9+IHcFjPLsaXrdAzJBQjUaO8eC1U6CkvxOcU5FkMDuoIRTCY1NS8KzE49FDYtk2PlGNwDSN+wLLeLiaWe9frS5ofxJjZwWtlw2LvXF5q5NV5GBoyIc9sOPZYmJiAV18tdEpKipyJn4jsFBEV4/fhGLbni8jLIjIhIq+IyIW5SlcsGseHmGpuyecp3SUnYU9vFMGNJ5AdnogCDA6m3NLTRPRicMwxevn884VNR4kRyvHxfwt8PWrdK85/ROR84EbgWuBh4D3ADSIiSql/z3H6EIH6iREG6xtzfarcYYc9lXIl9/BK/uP6fH6uHMMDYpICnhIRM6NDdhx+uH6gr7yCOvLIQqemZMi1+B1USj0Rb6OIhIBrgJuVUldYqx8VkYXAVSLyY6XUVC4TKJOTVE9N0FfMH2d7uw6LmPBSUjJtLOApMfEgCYsERvyyo6oKli2Dl1+GI88sdGpKhkLX+Z0EtAG3RK2/GWgB3pHrBISGdB+/6YYi9vzsvn5mlBdDEgriuRrxy54jj9TiZ3CNXIvfh0Rk1KrLeyJGfd9qa/lC1Pqt1nJVbpMHZQO6GfZ0fUOuT5UxSeu9cjDKi6lrM8QiI/E04pc9lvhJJFLolJQMuRS/e4EvAKcBnwLGgV+JyDkOG2vUVvqi9u2N2p4zgpbnFy5mzy9K/LIVrlIN8RlBz5ys3ok0GrwY4nDkkTA2RmPP3kKnpGRISfxE5L1xWm5G/35v76OU+oJS6r+VUo8ppe5GN2TZhG7YMnNo2zydRIvIBSKySUQ2TWfZ9yU0oMf1DHvY80tKjLCn6wKmFBSg1Om2XnlB1/Mpwp4Q/JER08E9W6yGLu17thc4IaVDqg1e/ggclYLdaLwNSqmwiNwF/LOILFBKdTPXw+t2mNseX8zRXJVSNwE3AdTU1GT1eYes0SemPCJ+GWVWORrfk+5uOPFE+M1v4PrrYetW+PWv4c//XC+fegp+8QvdeTkQ0L9IRHdoXr4c/umfsho4YE7fwF/8AoaH4bzz4Pvf1x1+L71Ub/vZz2DjRn3+yUm9tBsJrF6t05tCB+tUcEtMStW7nsf0tH4mNTWFTklxc5TOft/5wG389yeuysspReSjwNnAWqAd2AXcA3xLKTVk2SwDdsQ5RJNSqt9xvErgKuAcoBF4Fvh7pdT/RZ03APw98FmgE91D4JtKqV+6dW2QovgppUYBN2pboz09u25vNXPFz67re9GFcyakbLAIPL9kGW5lpa5TcbvO79e/hl274Ic/hFtv1Rvuvhteew3+8z/hscfgT3+C447TO4TDOlfv7YXbb4f3vQ/e4VKbpbPO0svzzoNLLtF/2+J3ySU6tNbeDmVlOsMdGIDxcb39/vt1WrLAN2LlNiMjelldndDMEx6ql2lrg+XLKZ8Yz+dZv4wWvK8BbwLHobuuvVtE3q6UcoaCrgXWR+0/FPX/T4APAJcB24GLgN+KyElKqWcddldZ574C+BPwCeAuEfmgUup+Ny4Mct/VYQarW8PHgF1KKTtwvRE4iK4TfNhhfg7a63s81+kKeqy1Z8aZbIxRXjJOQ6LgYDisl0pBV5cWo5/8ZK7N88/rjrnd3fP3z5boXHJkRAvftdfCV786d9vLL+sS8xtv2Du7n54iI+8iY4uf8fyy5+STaXrw4eR27vEhpZSzCfkGEekFfg6cAjzi2LY9Sbe2Y4FPAucppX5mrduAdoC+CZxhrWtHC991Sqnrrd0fFZGVwHWAt8VPRM4GzkQndDfQgVb5t6LdaACUUlMi8g/oTu170AJ4KnAe8AWl1GQu0uckNDDARLAMVVGZ61PlFhfFbwZ7DjFb8ECPMQi6X+G+fbBw4fz9Ojr0ct8+d9Ix6XgN9joq/MfHZwV2wYL5+x16qC5NOPYpdg/Ora4KebsPo1ZNiBG/7FmyhIa+gwQi4eS2LhAlfDZPW8tFaR7uDGAKuNNx/GkRuQP4qohUKKUm0A0ky5nf/e0W4KcislwpFS/Mmha5au25Ax0j/jbwEHoElwngdKXUHU5DpdR/AJ8DPo4eEeZs4O+UUj/KUdrmEBoaYKCyNvGEscVAW5v7/fxs0XjFMSjP5s16+cILun4vlui0tOh6NxfETykFe/bMP7+dPlv8YolwWZlOy97SaCHnZd2O+/mkGPY0pMCSJQTD0zQNRzeOzyvrrOVLUeuvFZFpERkQkfUicnTU9tXADqsKzclWtNitdNhNANti2IGL3d9y4vlZ7u+padjfiBbIvBMa6Ke/ogRKpe3t8OSTgDuhLYWaFQ2n4Nh/P/OMXsYSv2BQi7FD/DLyNOx93nxzdt2f/jT7d3e3Dr1CbPEDPUeji+JXbEWkgqfXhD3dY8kSADr6CzOYhYgsQocoH1ZKbbJWT6Dz7oeAA8CR6DrCP4rI8UopWySbmd+lDeZ3a2sG+tV8b8T17m95q/PzKqHBQQYqvd0MO6VQlz2+p9UdIZtZFGaEyhaNcUclux0CtcNZ+RAdZ92hLbp2+tIQv+zm8/Oy3zUfz6TWiJ97WAXN5hHXPL+QiGxy/H+T1ZJ+HiJSC/wGmAbOtddbrfadExE8JiIPoj21K9DtN0C/krG+wOhXNVW7rDHiN9DHQGUtRR+UaW/XwtTnYkgklQYrsTw/0PV+btX5OY/jFL/ubv2rqIDGOA2WFizQrVMtir3Or+gwdX7uYb3jtWPDbh1xWim1NpmR1UVhPbACWKeUejORvVJqt4j8AXibY3UvsDSGeZNju71ssiY2UAnssqbQY3sWnLLeHnqrPdzNIVVyMMTZPM/N7rMXcpSZOjtj7+u2+AWDWuR27dL1iXZDlq4u7fXFUzXb8yv2Ot00SXa1efNkTZ2fezTp/L923DXxS4qIlAG/BI4H3q+USnVepWgPbiuwXESiX4RVwCSzdXxbgQogeoJV17u/+V78Qn299JSS+LkkOGVjozrjOvbY2ZX2KB3HHTe7Ll4ndlv8shQdBVq82ttnQ5sLFuj/7Tq/eCFP0OI3MaH7/Rnyz7CVUZsRXrLHGh+1bjQ/4md1Nr8VPTrXmYm6MkTttxQ4GXjSsXo9UIbu7mbbhYCzgIeslp4AD6LF8FNRhz0HeMGtlp7g97DnyAjB8TH6qrw96G5K+jHH88tuHEUB6voP6n+OOw62bJlrsGYNPP30vP3m0Nmp6wqtEXSyYt8+LabV1bBjh/47Eplt7bl6dfx9bc90715XGn8EX98Gn7gUTjsNLrtsvsHrr8M3vqE90VBo7q++fn5fRDdRSo+q89prfOPFfYxPhWl79b+gsky3fL3wQvQkKrmhfnyYzz9xF3xlg34+g4Nw7716Y3Pidgr+8sszJBRivKqGuvx5fj9Ci9U1wIiInOjY9qZS6k0R+Q7aidqIbvByBHA5EAG+ZRsrpZ4VkTuB71ne5A50K//lOIROKbVfRL4LXC4iQ8AzaIE8Fd19zjX8LX4HdQbfW1Vf/FExZ9+6xuiIQfrU91nit2bN7Eo7tOhcFw9bdLLwRGcCc3v36uPZobPOTp25dnfr47/nPcnTsXcvUJNxuM++9KrHNsDvfqd/scTvrrvg5ptnQ7HT07ou1nrXWLuWnAlQfz9cpYe+OqJ1MdMRRcVEl86atm2Dl16Cf74tN+cG1m3/Exc++Utd3q+u1oNZh0Jw7rkm7OkSozV1btb5JcMeFukK6+fkG+jRXraiRezT6FL3QXTn928opV6J2udctJBejR7ebAu6+9szUXZXAMPAF5kd3uzjSql7s74iB/4WP8srGaoogQ+zpUXXi+3bh2rMXsnneH42dl3fW96S/ABzRCdL9u3T3p3daMIW+k2bdAOfePWOznR0dzPblShzQvuiOtlXRg2O0NWlPTxn30R7/aJFVuvULMWvv1/Xe+7Zo5/54YfPngPg9tv5xI42ugfGeeTSdaxoq4XPfEYP85YliQqJnUM9+o+BgbSnMDLtkFJjrLqOujyJn1JqWQo2PwV+muLxxoAvWb9EdmG0QF6dynEzxd91flZ9xGh5VYET4gKBgO5bZzV4ybZVY22/lZEdeaTOYC+5BP793/Vg0WvW6DDa1QneTbfETyktfp2ds6LX1KTr/Wxvyl6fy3RYBJ3iF6s1bHd37DrItjb9UNyokz3jDF0AWbUKjjhifnoWLJgvJp2dsH8/KpzZ6CCpvE+dwz0Ml1dlNHefffyij8DkmLGaurw2eCll/O35WeI3XAriB7rez4XMVUS05xcMQmsrjI1pr08EPvIRbTSZZOS5OaLTmnFaaseG9bk6Ouae0+ntJRK/pibdKGfvXlc84uB+h/h1denZK5x0dcXu/mGPNrNvnx7sLxseeyz2+jlDve2eu62zE8JhqgZ0Vxi3u3woFB1DPeyrbSGTpi0z4mdq/xIyVlNPXf+uQiejJDCeHzBaXunpjy7llLnYvaBuoFeLaSCgM+50c8umJi2YWXpcTUNWv8WODmiwWuVWVc0VmERhT5E5Hd2zzfSDzvFM7TCjk3ieH7jyfGKO6/jQQ/DZzyYe59S6R9W9uRsdpHO4h711mQ3AUWyDCBSKseq81vmVNEb8gJGyEvH83BS//oOJRSUZgYBOTxbiJyI0D1l9Wjs79bRGn/40XHzx3Aw+kedn7+tS2DO0fx+89a36n2jxs2e6SNTxP9vCQKxxHU87DW66SY/BWlsbe9Z0W/zshkw5oHOoh721LVkdw4Q9EzNWXUvt+Eihk1ESGPFDe34lgUt96wDqBnqSi0oyXBCdJlv8Ojq0N/mzn2mPNNWwp0vpABAVIbjfanxTXj5f/Pr7dZ/CHHp+bQMJxGvTpvjCa63PxvP7s73bkPFxPZfj0JCeQHhqSk9qHInQPtzLvrrMxG827GlIxGhNvRa/DOtuDbMY8cPy/Erhq+vogLExKsajB05Pn9r+3uw8P5gRnWxmzJjj+TlxZvLRLS7jpSPjVGiaxoaQ6WktbgsXzm/wYothDod8axtMIH7PPTfv3DPXbBUQqnsz8/xqXn2J+35+CcvP/QRccAFcdBG8/e3wtrfBCScQvOqblEem2Zeh5zczy7Vx/RIyVms1JnJzGEOfYsQPGCurKHBCXMLq6F43kOXwd0q55/lZmX2mNTpNw3267rCpae6GdPqNdXbCgQMEwtNZ1Sy1DzuEeOFC3Z+vvh7OP1+vT1TnZu83MkJlFoWT1kSeX6Jz19ZCbe2M55fufajatROAusc36BW//rVeWgMgBO69DyDzsKfl+hnpS8xwnfUduD19mQ/xt/iNjhKurPL8aMcpl4YtsaobzE786saHCU1PuSZ+Ys00kQlNQ306HYEsXtXOTi3oQ/2ZHwOH+C1YAN/6Flx6qR5v9L/+S/e5SzbDhHU/G4Yyfz4JPT87bfHo7KQmw7BneU/UfkNDc/6VLt2vMdOwpyE1huqN+LmFv7s6jI0RqfJWY5esoj5W5lo/0AtZ5EFxQ43pYjWvbxgdRM9tnB5ipyWeCD/xRGp9yixBaBw4CK2ZD2XXbk8ls2CBniV+3Tr4/Of134sXzxomE7/BPvoXxhrgPjmtgz26s3x0J3qbJOKXadizsjtGy1YHYmXGmXp+s2HPjHb3DcO2+Lk5gL1P8bf4jY4SqdD1RSXxzdnil6Xn1zzs6F6QDZZ4Ng9nnp6m4T5YHme4thNOSCsd9f090Loio3SICG1O8bNZsUIPabZ5sx5y7fDD40/f4/D8Mm3a35JM/JIM8l29K3okqRSIROi475dJzcIS4EBtU1K7WMwGX0riS8wZw8bzcw1/i58HPb+sIrBtetis+izr/FqGXfT8gJahzCvnm4d6ofNkV9LRMNCT1Q1uH+4lXFdPMLq+8aMf1b9kWOLXmEXhpHmoF1YdAyedpIdX27x5rkEiz2/Bgsw8vx/8gKo9u5OaHaxpJBwIpn98jOeXKsN11ryVRvyyxvd1fpFKb4lfLFLOD8rKoLmZ+sGerM7XPOSS52ft35Kh5yeRCI3D/a6loyFLj7htuI9wNmmxGiRlLH5iFSQ6O/UoL08+Od8murWn8+Xp7KRyeJCK6SSj80RzySUA3Hbs6QnNsunjJ6bBS0pEQmUMVdWasKcL+Fv8HJ5fyZQ4OzpcCXuGgzFaWKaLHfbMsIFH1cggZeHp7D3QqipoaNCeXxa0j/Qy3Z5FWsrKoLWVpgwLJ4HwNA2jA/p+BIP6eGutibjtKYMscZZYHq51H2fCt6nSokXtlbZDEpoNVGY+Z5/x/FKnr6bReH4u4O+w5+goKlkfsWKjo4O6ruz6ADUP9THU0ExjNi0sQTevr66mZbiPNzLYfSZ8m63nB9DZmb34DfcRbj8qu3QsWEBjhuJXM9hHQKm59+O+++Dll3Vfx9/8JvG8ebb4xRolxkkkoruo2F5keTldf3k2B8KJC0MRyfx9mR3Y2qhfMgZqjfi5gb89v9FRIpUlMJ2Rk44O6gd6Y5f8U6R5pI+hBhearFvjambq+dXZIpGt52cdo2GgJ+N+fqIU7SO9hDuyTMuCBTRlKML1/THuR0eHbnV6wgm6+4UdPowlIg7PL+H7cfHFs534w2HYt4/Jtg5G7f6wdqHIbmlrT6mUBWZsz9Tpr2kwYU8X8Lf4jY0R8Zjnl3XBt6NjVjQypHmol6FGl/prdXbOth5Nk5nM3g3PLwuPCyA0Mkz11ATTBRS/mTkWMy0MpBL23LYNfvQj/ffKlXDUURCJMNHeydOLVzO47j16WLN162DDBrjjDt3ZH9iy4LDM0uXA+H2JUQr6apuM5+cCvg972uLntVkdnAXztASxo4PqsRHKpiYyPnfTUB87DluV8f5z6OykZVsGzeuxxhcF98Ke/ZmLX0WPLmlPZ5sWS4Qz6fhfl20YuL0d5eyyEYtHHtHLww7T/ReVglWr6F33Hkb+dx/bb76bNUsa4fe/13Zr1gAwee993PBo5u8cM2HPzA/hFwZqG/VclpFIdoM/+Bx/i9/QEOHaGCPgFzNWi8KmTLsXRCI0D/exxY2wJ1ieX2Zhz/qBHqaCIcqybXhjpaNqYpTKicyGFqs8oMUvnE2DFysdoUiYmpHBtHed8fzSEj+HmoRCjDU0Jxa/HTv0cHIvvaQb1ViMv7wfiD8uaeS005n4w4NppGsuMw1ePFYI9SIDNQ1a+Hp79Xybhozwb7FBKRgcJFJq4mf3Jct0CK2+PkKRMMMuil/j6KAeLi1N6gZ6dIjHjeHnrJBfpoWCioO255d92BOs0WbSpL6/h5GK6vid6B3Eq9MbaW6lPVEYeudOWLp0jvDlA5lVP0MS9jZZ7+ATTxQ2IUWOfz2/0VEIhxmt1s2zdxwc4Z2HtRU4UbNkHP6xxC9pJj89DQ88oO9DIKDHamxthRdeAGCoMbNJSedhiU5tBt0v6vt76K1rzmBgtPjpyLTer/KA9nqm27JMzYz4pZ+Ouv6D9NY1kVz64reaHG1qo60rQWOJnTth2bK005YtdoMXo33JeerI4/WYshs2wAc/WOjkFC3+Fb+BAQBCjQ0wBHv6xwqcIE2sAntaoSC7Q3eyRib33Qcf+UjczfsWLU/9nImwRGcmZJcGtQO97K5zIeTpSEemHcwre/YzESwj0phleliaj9MAAA+USURBVCzxa8rA86sb6KW3toklWZx+tLmV1tdeiW+wcyd84ANZnCEzxNT5pcxUqByOPhqeyawu3aDxr/hZ82E1Lu6Al2BxU4l0eUi1zu+11/Ry40bdH29oSP/a2vjYL15m0dLsm68DszNNxGts8tpregSRyUmd8wUCOuSmFJ17trNlzXvdSYcd9kzk+e3eDbfcojuPT0/rOpUDB6Cjg0X3/ZL9tc3Zh2CdQ63FYnwcfvADGBzU6QiF9G9ggJUvbuLRo9+V1elHmttYMdJnqUzUtUxO6kl/l2Qjr5kxO5mtUb+UOOII+MMfCp2Kosa/4terPYBASwsQYWo682l33CTrkm9VFWOVNbMzoMfjjTegsRFOPHHepv0P97Moy2TMkMzze+ABuP9+WLVKjygTDuufUry45mTuPfGDuBLYaW0lHAjSlMjzu+46uOGG+evLyymTAL95yxmclm06amsZq6iKH/Z84AH4ylfmrw8GGa+s5umVb+XdWZx+pLmNivA0kf5+aIkKbdvzES5y7emnjOnnlxoz2cMhh+huJtPTunBkSBv/3jVL/KSlBThAxGPxlmwcjMGGZhqTeX47d+oPKNdYnl9tPM9v925df/HCC/Mu+md3bObl3dnNwTdDMMhgXVPisOfOnTqc9Ic/zHpcZWUgwgPPdXH9bZuzFz+gr6E1foOX3dYA0gcO6MJJOKwzuKoqvnb3czy1o5cY0hiX6Nd6tMlqHbh373zxSzYfYR7w2GfoSQT0txsO62e2NLPpsfyOf1t7WmHPYJMeJX06Ujpf3WB9CuL3xhv5adhQUcFgVR11fXEy+zff1HPhxVF7NzPDgYbmxA1edu+G5cv1yCXV1VBePi9dbiSnr74lvue3Z48+b0uLFt+KCt2606X+XCPNVqMu28uLPjcUxvObCXsaUmLNGvjUp0xpIQv8K3792qMIWKXfsJfFL82kTVRUcey2zTONeuYfTyX1/Ny8Gz21TbMd1qPZvTtuHVM2Q7TFor++JXGd365d8dPiYliur6Elfp1fV5f2vFy+dps5nl+sc0OBPT8Pf4de4vjjdf10PqI3JYrvxS/Y0AB4XPzSZH+HlYE3NsJTT8036OuD4eG4np/b2W5PXXP8Bi+255cHBhpa4tf5DQ3pwkIeGnv0NSTx/HIoPjOe374Y4md7nQXoOG2mNPImIrJERO4WkQERGRSRe0SkJOKs/hW/gQGorydYpqs9S0n8bjvnMu4+9Wz9z8UXzzfYtUsv89Sqr7e2KXaDl0hEZ7j5Er/6Fj2nX6yhxeyQXx7S0l/fQtXEqC6AxEpHDsVvsqaO8VA54vT8lIKeHt3ydsECV73OyekIJ3zrYX73UvzRYcBR4Cqdz7DoEZFq4BHgSOBvgb8GDgMeFZFUupt6Gv82eOnvh4YGRISAeF/8lFJJw4CRiCKsFBPhCD8943N89EPHw//7f7BihW64MTGh65DesCYYSlC3Mzw+Td/IJE015VmnvbeuidbndsMnP6m7YlRV6Qy2r0835shTHdNQbQOhSFi3LD3ySJ2GSETfm+ef10aJZkJ3ib56a/ScL35Rz4gQCulnA/rZvP/9uTu5CCNllVRefz08+qgOf+7fD1PWCDzr1mV02JGJ6Zjrf/jIa+wbnOD/+/kmdl4Xv/+gl7s6KKX4/u9e43sPv8a3PnI0Zx+/hIiCJ7f3sHZZMxGlqCxLf0ScyenIvH2VUnQNjLOocf4k26/sHWJ8KpzVtaTJ+cAK4Ail1DYAEXkOeA34LPCv+UyM20ixx9hramrUyMhI+jt+5CPw+uvw3HMs++r/AHBIi+7rJ4CthfZHacuOLUAzMhS1fWI6QnkwQNi6r2/0zB1LsrW2goBAeShAWTBA78gkZcEA9ZUh9g6OMzoZ/+Ve1FhFKCgz6RPRhfZQQGLuu/MrJ8EFF+i+Y7W1sztEIlpwrrtOi2EU9v0AaKwuo6osSDAgVIQCDIxNE1GKUECoLg/GvE92D7LpiCIcUax5+hF+9OvrCJRbAlxWptMRDuvGHI8/DsccMy8dX7rzWe7ZvIelzdXs6p29jx31FewbnKCyLEBHfSVBkZn77SweKCAgwsDYFL0jk5z0xnPcfsfX9DUfcYROg4ju3xYKaU/4ttt0uDiK+5/v5vO3znYqLg8F6KivYDqs6B4YZ3FTlb61Slk/nZFFlC5YRZRCKRiemOa4PS/zq1u+HPshr14Nd92lZ1OI4ot3bOY3z3bF3s/BirYath+Y/SaWNldb90Oxu3eML//ff3PR8EvI4sW6NW5n5+zy5JNjRgQefXk/5/7X0zP/11eGqK8qIyBCOKLmDBJhP5956WqtYTqiCAh0D4xTV1lGbYV+h6Kfb2VZEEF/b3YeZf/t/Aad7/2Chkq6B8YBZp6HUlpO7fu/f0inq64iREvtbMEuEBD6R/V7YlNbEWI4jqinSmttBbUVQXb2jHJISzVT0xG6BsY5pKV6Xt6QKgGB7ddmPxCBiIwqpeJ6cCLyO6BSKXVy1PoNAEqpzEpKHsG/nt/xx+sMELjstCN4qXuQYMCeCw2CgdmPzi4e2OWE2f/nbg+HFRPTYaorQgREi9TegXEmpiMc2lbD6wdGWN5aTU1FiIaqMobGpzlqQR3jUxGqyoMzH6bNkuYqdveOcdzSRjbv6ueEFc0zHqotMkpBWClE4HVHhtdSU64nNr377rRvzdc/tIqv3/siACetaGEqHKEiFCQcUUxHIkxMR2iuKZ/RDmd6nPcpFBCCASF0xFk8970vs2ZpjNFRnAeJ4mNrlzAZjhAMCIsaq9i4vYf3HtVOQISHXtzHsYsb6WyoZDqiCIrMaLtzGVEwNhnm4Zf2sfGQY7j6zqe48uNvS/uenLC8eU6GddySRhZapfNnd/ezrKWa5poKggEtuCIy83fASltQhIdf2sdmjuQjl/43v/rGX+pWnJOTWpAnJnRL0zj345wTD+HZ3f2MToYZGJuio76C3b1adA7vqOXVfcOsXljPstYaFjRU8vi2Hs5co0Oow+PTbNzeQ3NNOS//3d8jn3xLWtd/3NJGlrfWsOOgfsfee1QHCmbEbM+zOh1/tqie1QsaGByf4uGX9vHnqzq4//m9nLC8mda6CoLWte0d3EtHfQWHttUSjqgZ8autCLGspYbOhsqZwkRAZObdih4CdElzNRtePcDS5moO76idEb+3LWueue8B0Q2WAgG4/SndleRdh7fNfO+gv6EDgxM8NTJbJ/zOw1p54IUYdaNJKA8FmLT6Da89pImhiSn+bJGOMgmwfksXnfWVGYvfeSe7NPpSclYDv4mxfivwsXwlIlf41/MzGAwGHyMik8DzjlU3KaVuitr+r0qpr0btdzXwVaVUUTtPRZ14g8FgMGTMtFJqbRKbWN5RSQzH49/WngaDwWBIRB8Qa3qXJmtbUWPEz2AwGAyx2Iqu94tmFfBintPiOkb8DAaDwRCL9cCJIrLCXiEiy4CTrW1FjWnwYjAYDD4kha4ONcAWYAy4El3/dxVQBxyjlIoxSkPxYDw/g8FgMMxDKTUCnAq8CtwM3ArsAE4tduED4/kZDAaDL0nm+ZU6xvMzGAwGg+8oes9PRCLomHQmhIDsxi8yFArz7IoX8+y8QZVSyrcOUNGLXzaIyKYUOnkaPIh5dsWLeXYGL+Bb1TcYDAaDfzHiZzAYDAbf4Xfxuym5icGjmGdXvJhnZyg4vq7zMxgMBoM/8bvnZzAYDAYf4jvxE5ElInK3iAyIyKCI3CMiSwudrlJDRD4qIr8UkTdEZExEXhGRa0WkLsquSUR+LCIHRWRERB4WkaNjHK9SRL4tIt3W8TaKyLti2AVE5HIR2Ski4yKyRUT+Kk4azxeRl0Vkwkrfhe7dgdJCRB4UEWXN5eZcb56foThRSvnmB1QDrwEvAB8GzkRP5vg6UFPo9JXSD3gC+AXwKWAdcAnQb60PWDYCPAa8CZwNnA5sAA4Ci6OOd6u1//nAe4B70P0710TZXQNMAF8G3g3cCESA90fZnW+tv8ayu9r6/3OFvnde+1nPphs9tuPVjvXm+Zlf0f4KnoC8Xix8EQgDKx3rlqM73H6p0OkrpR/QFmPd31gZ6KnW/2da/7/bYdMA9AL/5lh3rGV3rmNdCHgFWO9Y125lnN+IOu/vgOei9t0P/DzK7qdWxl1W6PvnlR/QCOy1xC1a/MzzM7+i/fkt7HkG8IRSapu9Qim1A3gc/SEbXEIpdSDG6qet5SJreQbQpZR61LHfAHAvc5/HGcAUcKfDbhq4AzhNRCqs1acB5cAtUee9BThaRJZb/58EtMWwuxloAd6R7Pp8xL8AW5VSt8fYZp6foWjxm/itRoc8o9mKnqDRkFvWWcuXrGWi57FURGoddjuUUqMx7MqBlQ67CWBbDDuYfcb2BJ3R54628zUi8g60t/75OCbm+RmKFr+JXzPQF2N9L9CU57T4ChFZBHwTeFgptclaneh5wOwzSWbX7Fj2K6Wi++/EsiPGMaPtfIuIlKHr265XSr0Sx8w8P0PR4jfxA133EI3kPRU+wvIAfoOuWz3XuYnUnkcu7Ihja9D8PVCFblASD/P8DEWL38Svj9ilwiZil0wNWSIilcB6YAVwmlLqTcfmXuI/D5h9Jsnseh3LJhGJzixj2RHjmM1R232J1fXnCuAfgAoRaRSRRmuz/X8Q8/wMRYzfxG8rs/UFTlYBL+Y5LSWPFTr7JXA8uqn681EmiZ7HLjU7W/RWYLmIVMewm2S2jmgrUAEcGsMOZp+xXTcUfe5oO7+yAqhENyjpc/xAd0HoA47GPD9DEeM38VsPnCgiK+wVIrIMONnaZnAJEQmg+3a9BzhTKfVEDLP1wCIRWefYrx74EHOfx3qgDPiYwy4EnAU8pJSasFY/iM5MPxV1nnOAF6yWvQAb0U3iY9n1olv/+pln0X3non+gBfHdaMEyz89QvBS6r0U+f0AN+qN9Ht0U+wxgC7AdqC10+krpB/w7Vr8w4MSo32LLJgD8EdgNfALd1P336AxsSdTx7kB7HJ9BC+rdwDjwlii766z1XwJOsdIRAT4UZXehtf5qy+6b1v8XFfreefXH/H5+5vmZX9H+Cp6AvF8wLEWH4gaBIeDXwLJCp6vUfsBOK7OM9fu6w64Z3Tm5FxhFd2g+NsbxqoB/RXe4HgeeBE6JYRcErgTeQDebfw74aJw0fhZ41bJ7Dfh8oe+bl3/R4meen/kV88/M6mAwGAwG3+G3Oj+DwWAwGIz4GQwGg8F/GPEzGAwGg+8w4mcwGAwG32HEz2AwGAy+w4ifwWAwGHyHET+DwWAw+A4jfgaDwfD/bxSMODBa+Y2CUTAKRsEoGHEAAIyJrHltFJyqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(data.ant_elev_angle)\n",
    "ax2 = plt.twinx()\n",
    "ax2.plot(data.alt_msl,'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}