Xhistogram with weights

xhistogram_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Area weighted temperature Xhistogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Collecting git+https://github.com/xgcm/xhistogram.git@2f1e8a6d7f4e1f5b56f3c36f7a8f1306675148a0\n",
      "  Cloning https://github.com/xgcm/xhistogram.git (to revision 2f1e8a6d7f4e1f5b56f3c36f7a8f1306675148a0) to /tmp/pip-req-build-taungjpe\n",
      "  Running command git clone -q https://github.com/xgcm/xhistogram.git /tmp/pip-req-build-taungjpe\n",
      "Requirement already satisfied (use --upgrade to upgrade): xhistogram==0+untagged.39.g2f1e8a6 from git+https://github.com/xgcm/xhistogram.git@2f1e8a6d7f4e1f5b56f3c36f7a8f1306675148a0 in /srv/conda/envs/notebook/lib/python3.7/site-packages\n",
      "Requirement already satisfied: xarray in /srv/conda/envs/notebook/lib/python3.7/site-packages (from xhistogram==0+untagged.39.g2f1e8a6) (0.12.1)\n",
      "Requirement already satisfied: dask in /srv/conda/envs/notebook/lib/python3.7/site-packages (from xhistogram==0+untagged.39.g2f1e8a6) (1.2.2)\n",
      "Requirement already satisfied: numpy>=1.16 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from xhistogram==0+untagged.39.g2f1e8a6) (1.16.4)\n",
      "Requirement already satisfied: pandas>=0.19.2 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from xarray->xhistogram==0+untagged.39.g2f1e8a6) (0.24.2)\n",
      "Requirement already satisfied: python-dateutil>=2.5.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from pandas>=0.19.2->xarray->xhistogram==0+untagged.39.g2f1e8a6) (2.8.0)\n",
      "Requirement already satisfied: pytz>=2011k in /srv/conda/envs/notebook/lib/python3.7/site-packages (from pandas>=0.19.2->xarray->xhistogram==0+untagged.39.g2f1e8a6) (2019.1)\n",
      "Requirement already satisfied: six>=1.5 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from python-dateutil>=2.5.0->pandas>=0.19.2->xarray->xhistogram==0+untagged.39.g2f1e8a6) (1.12.0)\n",
      "Building wheels for collected packages: xhistogram\n",
      "  Building wheel for xhistogram (setup.py) ... \u001b[?25ldone\n",
      "\u001b[?25h  Stored in directory: /tmp/pip-ephem-wheel-cache-iyaot4hb/wheels/0a/e3/93/4033fa1488deacc8f0c18b294039b9503ffe980c63c8f6ab8d\n",
      "Successfully built xhistogram\n",
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "%pip install git+https://github.com/xgcm/xhistogram.git@2f1e8a6d7f4e1f5b56f3c36f7a8f1306675148a0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from dask.distributed import Client\n",
    "import xarray as xr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table style=\"border: 2px solid white;\">\n",
       "<tr>\n",
       "<td style=\"vertical-align: top; border: 0px solid white\">\n",
       "<h3>Client</h3>\n",
       "<ul>\n",
       "  <li><b>Scheduler: </b>tcp://127.0.0.1:37931\n",
       "  <li><b>Dashboard: </b><a href='/user/dask-dask-examples-rtte9r97/proxy/8787/status' target='_blank'>/user/dask-dask-examples-rtte9r97/proxy/8787/status</a>\n",
       "</ul>\n",
       "</td>\n",
       "<td style=\"vertical-align: top; border: 0px solid white\">\n",
       "<h3>Cluster</h3>\n",
       "<ul>\n",
       "  <li><b>Workers: </b>2</li>\n",
       "  <li><b>Cores: </b>4</li>\n",
       "  <li><b>Memory: </b>2.00 GB</li>\n",
       "</ul>\n",
       "</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Client: scheduler='tcp://127.0.0.1:37931' processes=2 cores=4>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "client = Client(n_workers=2, threads_per_worker=2, memory_limit='1GB')\n",
    "client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.Dataset>\n",
       "Dimensions:  (lat: 25, lon: 53, time: 2920)\n",
       "Coordinates:\n",
       "  * lat      (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0\n",
       "  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0\n",
       "  * time     (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00\n",
       "Data variables:\n",
       "    air      (time, lat, lon) float32 dask.array<shape=(2920, 25, 53), chunksize=(292, 13, 27)>\n",
       "Attributes:\n",
       "    Conventions:  COARDS\n",
       "    title:        4x daily NMC reanalysis (1948)\n",
       "    description:  Data is from NMC initialized reanalysis\\n(4x/day).  These a...\n",
       "    platform:     Model\n",
       "    references:   http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly..."
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds = xr.tutorial.open_dataset('air_temperature',\n",
    "                              chunks={'lat': 13, 'lon': 27, 'time': 292})\n",
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<xarray.DataArray 'temperature' (time: 2920, lat: 25, lon: 53)>\n",
       "dask.array<shape=(2920, 25, 53), dtype=float32, chunksize=(292, 13, 27)>\n",
       "Coordinates:\n",
       "  * lat      (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0\n",
       "  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0\n",
       "  * time     (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00\n",
       "Attributes:\n",
       "    long_name:     4xDaily Air temperature at sigma level 995\n",
       "    units:         degK\n",
       "    precision:     2\n",
       "    GRIB_id:       11\n",
       "    GRIB_name:     TMP\n",
       "    var_desc:      Air temperature\n",
       "    dataset:       NMC Reanalysis\n",
       "    level_desc:    Surface\n",
       "    statistic:     Individual Obs\n",
       "    parent_stat:   Other\n",
       "    actual_range:  [185.16 322.1 ]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "temperature = ds.air.rename(\"temperature\")\n",
    "temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from xhistogram.xarray import histogram as xhistogram\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "cos_lat = xr.ufuncs.cos(\n",
    "    xr.ufuncs.deg2rad(temperature.coords[\"lat\"]))\n",
    "cos_lat = cos_lat.where(temperature)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "temperature_bins = np.linspace(temperature.min(),\n",
    "                               temperature.max(), 20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_hist_area_weighted = xhistogram(\n",
    "    temperature,\n",
    "    bins=(temperature_bins, ),\n",
    "    dim=(\"time\", \"lat\", \"lon\"),\n",
    "    weights=cos_lat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_hist_flat = xhistogram(\n",
    "    temperature,\n",
    "    bins=[temperature_bins, ],\n",
    "    dim=(\"time\", \"lat\", \"lon\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_hist_area_weighted /= temp_hist_area_weighted.sum()\n",
    "temp_hist_flat /= temp_hist_flat.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_hist_area_weighted = temp_hist_area_weighted.persist()\n",
    "temp_hist_flat = temp_hist_flat.persist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAESCAYAAAAIfCk9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeXhU5fnw8e+dyTKB7CEJOwn7DkIgoCKIC7hrXaqiFZci1t2ftlrbqm2t1Wp9xdbijgvVKoqi4q7gAkHCvsgatrAlJGTfM8/7x5mEIWSZgUwmmbk/1zVXZs45z5l7DmHuPOfZxBiDUkop1ZAgXweglFKq7dIkoZRSqlGaJJRSSjVKk4RSSqlGaZJQSinVKE0SSimlGhXs6wBaUqdOnUxycrKvw1BKqXZlxYoVh4wxCQ3t86skkZycTEZGhq/DUEqpdkVEdjW2T283KaWUapQmCaWUUo3SJKGUUqpRftUm0ZCqqiqysrIoLy/3dShKuc1ut9O9e3dCQkJ8HYoKcH6fJLKysoiMjCQ5ORkR8XU4SjXLGENubi5ZWVmkpKT4OhwV4Pz+dlN5eTnx8fGaIFS7ISLEx8dr7Ve1CX6fJABNEKrd0d/Z5pVX1fg6hIAQEElCqdYwa9YsBg0aRLdu3bjtttuaPHbRokUsWbKklSLzPz9sPcTQhz7ns/UHfB2K39MkoQDIyMjgjjvuaPKYnTt3MnTo0Ab3zZkzh3379nn0nk2drz167rnnWLhwIY8++mizx2qSODFf/XyQaofhzrdXsWJXnq/D8WuaJNqImhrfVp1TU1OZNWvWcZc/niThT2bOnElmZiYXXnghhw8frtv+0UcfkZaWxkknncSZZ57JwYMH2blzJ7Nnz+bpp59m5MiRfP/99z6MvH1Kz8xlRPdousaEc+NrGWzPKfZ1SH7L73s3uXrkow1s3FfYoucc3DWKhy4Y0uxxF198MXv27KG8vJw777yTGTNmEBERwT333MPnn3/OU089RXh4OPfccw/FxcV06tSJOXPm0KVLF1588UVeeOEFKisr6du3L2+88QYdOnQ45j1qamro168f27dvp6CggLi4OBYtWsRpp53GhAkTePXVV+nSpQu3334769ato7q6mocffpiLLrqIRYsW8eSTT/Lxxx+Tk5PD1VdfTW5uLmPGjOGzzz5jxYoVde/x61//miVLltCtWzc+/PBDPvnkEzIyMpg2bRrh4eEsXbqUjRs3NvhZVqxYwQ033ECHDh049dRTW/Tfos6n98OBdS17zs7D4Jy/N7p79uzZfPbZZ3z77bd8/PHHddtPPfVU0tPTERFeeuklnnjiCZ566ilmzpxJREQE9957b8vGGQAOl1Sy6UAR957dnwtGdOXS/yzhuld+4v3fnExipN3X4fkdrUm0kldeeYUVK1aQkZHBrFmzyM3NpaSkhKFDh7Js2TLS0tK4/fbbmTdvXt0X6YMPPgjAL37xC5YvX86aNWsYNGgQL7/8coPvYbPZ6N+/Pxs3buSHH35g9OjRfP/991RUVJCVlUXfvn159NFHmTx5MsuXL+fbb7/lvvvuo6Sk5KjzPPLII0yePJmVK1dyySWXsHv37rp9W7du5dZbb2XDhg3ExMTw3nvvcdlll5GamsrcuXNZvXo1wcHBjX6W66+/nlmzZrF06VIvXem2JSsriylTpjBs2DD+8Y9/sGHDBl+H1O4t25ELwLje8fSK78jL140ht7iSG+Ysp7ii2sfR+Z+Aqkm48xe/t8yaNYv58+cDsGfPHrZu3YrNZuPSSy8FYPPmzaxfv56zzjoLsP5i79KlCwDr16/nD3/4A/n5+RQXFzNlypRG32fChAl899137NixgwceeIAXX3yRiRMnMmbMGAC++OILFixYwJNPPglYXYRdkwDADz/8UBfr1KlTiY2NrduXkpLCyJEjARg9ejQ7d+48JobGPktBQQH5+flMnDgRgGuvvZZPP/3Ug6vopib+4m9tt99+O/fccw8XXnghixYt4uGHH/Z1SO1eemYe4SE2hnePAWBEjxiemzaKm17P4Na5K3npulRCbPr3b0sJqCThK4sWLeKrr75i6dKldOjQgUmTJlFeXo7dbsdmswHWAKohQ4Y0+Bf29OnT+eCDDxgxYgRz5sxh0aJFjb7XhAkTmD17Nvv27ePPf/4z//jHP+puOdW+z3vvvceAAQOOKnfw4MG658aYRs8fFhZW99xms1FWVnbMMY19lvz8/IDr2llQUEC3bt0AeO211+q2R0ZGUljYsrc+A8XS7bmkJscSGnwkEZw+MJFHLx7K/e+v4/fvr+OJy4YH3O+at2i6bQUFBQXExsbSoUMHNm3aRHp6+jHHDBgwgJycnLov1qqqqrpbE0VFRXTp0oWqqirmzp3b5HulpaWxZMkSgoKCsNvtjBw5kueff54JEyYAMGXKFJ599tm6RLBq1apjznHqqafyzjvvAFbNw7UhtjGRkZEUFRU1+VliYmKIjo7mhx9+AGj2s/iDhx9+mMsvv5wJEybQqVOnuu0XXHAB8+fP14ZrD+UWV7D5YBHjescfs+/KsT2544x+vLsii6e/2uqD6PyT1iRawdSpU5k9ezbDhw9nwIABjBs37phjQkNDmTdvHnfccQcFBQVUV1dz1113MWTIEP7yl7+QlpZGr169GDZsWN2XcUPCwsLo0aNH3XtMmDCBt956i2HDhgHwxz/+kbvuuovhw4djjCE5OfmohlaAhx56iKuuuor//e9/TJw4kS5duhAZGUlxceM9SKZPn87MmTPrGq4b+yyvvvpqXcN1U7fN2qPaW2/Tp09n+vTpAFx00UVcdNFFxxzbv39/1q5d24rR+YefdljdXRtKEgB3n9mPAwVlzPp6K12i7Vw1tmdrhueXpKlbC+1Namqqqb/o0M8//8ygQYN8FFH7VFFRgc1mIzg4mKVLl3LLLbewevVqX4cVcPR391h/+nA981Zkseahsxttd6iqcXDTaxn8sO0QL/5qNJMHJrVylO2PiKwwxqQ2tE9vN6lj7N69mzFjxjBixAjuuOMOXnzxRV+HpBRgjY9ITY5rsmE6xBbEc9NGMahLJLfOXcWaPfmtGKH/0dtN7dSjjz7Ku+++e9S2yy+/vK6r6Yno169fg20VSvnSoeIKthws5uKTujV7bMewYF6ZPoZfPLeEG+Ys5/3fnEyv+I6tEKX/0STRTj344IMtkhCUai+WZTbdHlFfYqSd124Yy2XOwXbv3XIy8RFhzRdURwmI203+1O6iAoP+zh4rPTOXjqE2hnWLdrtMn4QIXroulf0F5dz4WgZllTpzrKf8PknY7XZyc3P1P51qN2oXHbLbdYoJV0vdaI9oyOhecTxz5Umsycrn9rdWUl3j8FKE/snvbzd1796drKwscnJyfB2KUm6rXb5UWXKKKtiWXcxlo4/vmkwd2plHLhzCnz7cwEMLNvDXi4fqYDs3+X2SCAkJ0SUglWrn0jOPzNd0vH41Ppl9+eXMXrydrjHh3Hp635YKz695/XaTiEwVkc0isk1E7m9g/zQRWet8LBGREe6WVUoFhvTMXCLCghnaNeqEzvPbKQO4eGRX/vH5Zt5bkdVC0fk3ryYJEbEB/wbOAQYDV4nI4HqH7QAmGmOGA38BXvCgrFIqAKRn5jImOZbgE5y4LyhIeOKyEZzcJ57fvbeW77fqbejmeLsmMRbYZozJNMZUAm8DR81RYIxZYoypnRwoHejublmllP/LLixne07JCd1qchUaHMTsa0fTNzGC37y5ktJKnV68Kd5OEt2APS6vs5zbGnMjUDt3tKdllVJ+KN05X9P4Pi2TJACi7CH85vS+FFVUs/fwsTMZqyO8nSQa6j7QYF9UETkdK0n8zpOyIjJDRDJEJEN7MCnlf5ZuzyUyLJjBXU6sPaK+pEhrYN3BwooWPa+/8XaSyAJ6uLzuDhyzELKIDAdeAi4yxuR6UtYY84IxJtUYk5qQkNBigSul2oZlmbmMSYk74faI+pKirHEoBwvLW/S8/sbbSWI50E9EUkQkFLgSWOB6gIj0BN4HrjXGbPGkrFLKvx0sLCfzUAnjW6g9wlVilLMmUaRJoileHSdhjKkWkduAzwEb8IoxZoOIzHTunw38CYgHnnMObql21gwaLOvNeJVSbUtLjI9oTIfQYCLDgsnW201N8vpgOmPMQmBhvW2zXZ7fBNzkblmlVOBIz8wl0h7M4BMcH9GYxKgwsrUm0SS/n7tJKdV+pWfmkZYShy3IO1NoJEXZteG6GZoklFJt0oGCcnYcarnxEQ1JjAzThutmaJJQSrVJ3myPqJUUZSe7qEJniW6CJgmlVJu0dHsuUfZgBrXw+AhXiVF2KqsdFJRVee092jtNEkqpNil9Ry5pveO91h4BkBSlA+qao0lCKdXm7MsvY1duqVdvNYG1xCnogLqmaJJQSrU5R9oj4rz6PrU1iewirUk0RpOEUqrNSc/MJTo8hEGdvdceAVqTcIcmCaVUm7M0M5e0lDiCvNgeARAeaiPKHky2JolGaZJQSrUpWYdL2ZNX1qJTgzclUQfUNUmThFKqTVmWaa0f4e1G61pJOjVHkzRJKKXalKWZucR2CGFAUmSrvF9SpNYkmqJJQinVpqRn5pKWEu/19ohaiVF2sovKddR1IzRJKKXajD15pWQdLvN611dXiZFhVNUYDpfqqOuGaJJQSrUZteMjxvfp1GrvWbtCnbZLNEyThFKqzUjPzCOuYyj9EiNa7T11ao6maZJQSrUJxhjSM3MZ19v74yNc6VrXTXM7SYhIkoi8LCKfOl8PFpEbvReaUiqQZB0uY29+mXtdXx0O2L8GWqCxOSHSOTWHJokGeVKTmIO13nRX5+stwF0tHZBSKjAt3e7B+hHLX4LnT4Nlz5/w+9pDbESHh+j8TY3wJEl0Msa8AzgAjDHVQI1XolJKBZz0zFzi3WmPcDhg2X+s51/8AbIyTvi9k6J0hbrGeJIkSkQkHjAAIjIOKPBKVEqpgHKkPSIekWbaI7Z+AXmZcN5TENUF3p0OpXkn9P661nXjPEkS9wALgD4i8iPwOnC7V6JSSgWU3Xml7CsoZ5w78zWlPweRXWHUdXD5HCg+CB/cYtUwjlNCZJi2STTCrSQhIkGAHZgInAzcDAwxxqz1YmxKqQBRNz6iuUF0BzfCjsUw9iawhUC30TDlb7DlM1jyzHG/f1KUnZziChwOHXVdn1tJwhjjAJ4yxlQbYzYYY9YbY3R4olKqRSzdnkuniDD6JDTTHrFsNgTbYfT1R7aNuQmGXAJf/wV2/nhc759UN+q68rjK+zNPbjd9ISKXSrM3DJVSyn1We0Qe43rHNd0eUZoHa/8Hw6+ADi41DhG4YBbEJsO8G6A4x+MYjoyV0HaJ+jxtk3gXqBCRQhEpEpFCL8WllAoQu3JLOVBY3nzX1xWvQnU5pM08dp89Cq54Hcrz4b0bweFZx8vE2lHXOjXHMdxOEsaYSGNMkDEm1BgT5Xzt3bUFlVJ+b2ndfE1NJImaKvjpJUiZCElDGj6m81A490mrzWLxEx7FULuMaY7WJI4R7O6BInJaQ9uNMd+1XDhKqUCTnplLQmQYvTt1bPygnxdA0T44/59Nn+yka2DXElj8OPRMgz6T3YqhriahPZyO4XaSAO5zeW4HxgIrAPf+FZRSqh5jDEu3uzE+In02xKZAvylNn1AEznsS9q2C934NM7+HqK5NlwHCgm3EdgjR200N8OR20wUuj7OAocBB74WmlPJ3Ow6VkF1Uwfim2iOyVkDWT5B2MwS58ZUV2tFqn6gqsxqya9zriKkD6hp2IrPAZmElCqWUOi7pdetZNzE+Ytl/IDQSRk5z/8QJ/eGCZ2D3UvjmL+4ViQzT+Zsa4EmbxLM4p+TASi4jgTXeCEopFRiWZuaSFBVGSmPtEYX7YcN8GPNrqweTJ4ZfDruXwI/PQM/xMOCcJg9PirKzLfuQZ+8RADxpk3CdRasaeMsYc3wjV5RSAa92vqaT+zTRHpHxstWdNW3G8b3JlMdg7wqYPxNu/g5iezV6aFKUVZNwOEyrrmfR1nlyuynGGPOa8zHXGPOjiNzptciUUn4t81AJOUUVjY+PqCqHjFeh/1SI6318bxJih8tfs9adeHc6VDd+Oykpyk6Nw5BboqOuXXmSJK5rYNv0FopDKRVgatePaLTRev08KD0E4xoYPOeJuBS4+N+wbyV88cdGD0uM1G6wDWn2dpOIXAVcDaSIyAKXXZFArrcCU0r5t/TMXDpH2ekV3+HYncZA+n8gcbA1gO5EDboAxt0K6f+GXuOtuZ7qSXROzZGjjddHcadNYgmwH+gEPOWyvQhodhZYEZkKPAPYgJeMMX+vt38g8CowCnjQGPOky76dzvepAaqNMaluxKuUauNq52ua0K9Tw+0RO3+Ag+utOZlaarq4sx6BrOXw4e2QNAw69T1qt6513bBmk4QxZhewCxjv6clFxAb8GzgLq8vschFZYIzZ6HJYHnAHcHEjpzndGKNdDpTyI9tzijlUXNF419dlsyE8zprMr6XYQuDyV2H2BHj3OrjpKwgJr9udEFF7u0lrEq7cbpMQkXEislxEikWkUkRq3JjgbyywzRiTaYypBN4GLnI9wBiTbYxZDujU40oFiKV14yMaaI/I2wGbPoHR04/6Em8R0d3hFy9YtZRPf3vUrtDgIOI6huqo63o8abj+F3AVsBUIB24Cnm2mTDdgj8vrLOc2dxmsKcpXiEiDfeBEZIaIZIhIRk6O51MEK6VaX/r2XLpG2+kZ10B7xE8vggRZ60R4Q7+zYMK9sPJ1WP3WUbsSI8PI1prEUTwacW2M2QbYjDE1xphXgdObKdLQzURPln46xRgzCjgHuLWhSQaNMS8YY1KNMakJCQkenFop5QtNrmddUQSr3oAhF0O0J39PemjSA5A8AT65x3pPp6QoO9lakziKJ0miVERCgdUi8oSI3A00MW0jYNUceri87g7sc/cNjTH7nD+zgflYt6+UUu3Yz/uLyC2pbHg969VvQUUhpN3i3SBswda6FFWlcGhr3eakqDBtuK7HkyRxrfP424ASrC//S5spsxzoJyIpzgRzJbCgmTIAiEhHEYmsfQ6cDaz3IF6lVBv0RvpOwoKDmDww8egdDofVYN1tNPQY4/1A4lKsn4d31G1KjLSTU1RBja51XcetaTmcvZQeNcZcA5QDj7hTzhhTLSK3AZ9jdYF9xRizQURmOvfPFpHOWFN+RAEOEbkLGIzV5Xa+szoaDPzXGPOZR59OKdWmZBeV896KvVwxpjudnL2J6mz7CvK2wy9eap1gYpOtn3lHkkRSVBgOA7klFXULEQU6t5KEMaZGRBJEJNTZS8ltxpiFwMJ622a7PD+AdRuqvkJghCfvpZRq2+b8uJMqh4ObTm1gmo305yCiMwy+6Nh93hDaESKSjq5JOMdKZBdqkqjlyQR/O4EfnaOuS2o3GmOaWSpKKaWguKKaN9N3cc7QziTXn/U1exNkfguT/wDBoa0XVGwK5O2se+k6oG5ot+jWi6MN8yRJ7HM+grCm5FBKKbe9/dNuCsurufm0PsfuXDYbbGEw+vrWDSouBXYcWYH5yPxN2g22lttJwhjzCFiNyMaYkuaOV0qpWlU1Dl75YQdpKXGM6BFz9M7SPFjztrX+Q8dOrRtYbIr13lXlEGInwZkktBvsEZ6MuB4vIhuBn52vR4jIc16LTCnlNz5eu499BeXMnNhALWLla1Bd5v1urw2JSwEM5O8CIMQWRKeIUK1JuPCkC+z/A6bgnPnVGLMGOGZwm1JKuTLG8PziTPonRTBpQL0BrzXV8NNL1sC2zj5YDTnW2Q027+husNk6VqKOpyOu99TbVNOCsSil/NDiLTlsOlDEjNP6HDvCetNHUJgF43xQi4CGx0pEhen8TS48SRJ7RORkwIhIqIjci/PWk1JKNeaF7zLpHGXnwhFdj92ZPhtielmrz/lCh3gIjTx6rESkXedvcuFJkpgJ3Io1Qd9eYKTztVJKNWhtVj5Ltudyw6nJhAbX+7rZuxL2pEPazRBk802AIhCXfFRNIikqjEPFFVTXOHwTUxvjSe+mQ8A0L8ailPIzz3+XSWRYMFeN7XnszmWzITQCTrqm9QNzFZsC2UduiiRG2Z2jrivrxk0EMk96N/UWkY9EJEdEskXkQxE5ztXJlVL+bnduKZ+u28+0cb2ItIccvbPoAKx/H0ZeDXYfD1qLS7F6NzmsJlZd6/pontxu+i/wDtAF6Aq8C7zVZAmlVMB66YdMbEHC9ackH7sz4xVwVFkzsfpabArUVEKhNUF1ksvUHMqzJCHGmDeMMdXOx5t4tjaEUipA5JVU8k7GHi45qduxt2yqyq0k0W8KxDcwbqK11evhVDc1h/ZwAjxLEt+KyP0ikiwivUTkt8AnIhInIo0sVKuUCkSvL91JeZWDGafVuyNtjLXQT0kOjG8j/V7qjZXoFBGKiE7NUcuTuZt+6fx5c73tN2DVKLR9QilFWWUNry3ZyZmDEumbWG+at8WPw+q5MPF30HuibwKsL7o7BIXU1SSCbUHEdwzTAXVOnvRuSvFmIEop//Duij0cLq3i5vpTcKyaC4segxFXW8uHthVBNojpecy6EtlFWpMAD5KEc+Gh84Bk13I6VbhSqlZ1jYMXv89kVM8YUnvFHtmx7Wv46A7oPQkueMYan9CWxKXUGyth195NTp60SXwETAfisaYKr30opRQAn204wJ68sqOn4DiwDt65DjoNgCteb931ItxVu66EsfriWGtda00CPGuT6G6MGe61SJRS7VrtRH69O3XkrMFJ1saCvTD3CgiLhGnv+n5MRGPiUqCiAMoOQ4c4EiLt5JZYo66DbR5Nced3PPn0n4rI2V6LRCnVri3dnsu6vQXcNKE3tiCB8gKYezlUFFkJIrqbr0NsXL0eTklRYRgDh4o9Wq3ZL3mSJNKB+SJSJiKFIlIkIoXeCkwp1b48/10mnSJC+cWoblBdCf+7Fg5thl++7ptpwD1Rf6xE5JFlTAOdJ0niKWA80MEYE2WMiTTGRHkpLqVUO/Lz/kIWb8nh+lNSsAcHWY3UOxbDhc9Cn8m+Dq95scnWz7x6A+o0SXiUJLYC640xOspaKXWUF77LpEOojWvSelndXNe8BZN+b83N1B6EhENkl7qaRGKUc/4m7QbrUcP1fmCRiHwK1F057QKrVGDbm1/GR2v28avxyURvetsaMHfSNTDxt74OzTOxKXU1ifiOoQQJ5GhNwqMkscP5CHU+lFKKV37YgQF+02MnfHindXvp/P/X9sZCNCcuBbZ/A1ijrjtFaDdY8GzE9SMAItLRGFPivZCUUu1FQWkVb/20m5n9i+m08G5IHAyXvwa2kOYLtzWxKVC0H6rKICTcGlCnk/x5tJ7EeBHZiHPJUhEZISLPeS0ypVSb9+ayXcRUHuSu7AetMRDT3gF7O+3PUtfDaSegA+pqedJw/f+AKUAugDFmDXCaN4JSSrV95VU1vPvDBt6J/CchNeUwbR5ENbCOdXtRb6xEQqSdHK1JeJQkMMbsqbeppgVjUUq1Ix+u2MnfKv9O15q98Ms3IWmwr0M6McesKxHGoeJKqgJ8rWtPksQeETkZMCISKiL34rz1pJQKLI4aB3Ff3cPJto3Ihc+2nWm/T0R4LIRFHzNWIifAu8F6kiRmArcC3YAsYCTwG28EpZRq23bM+z1nVS9i0+A7kZFX+TqcliECcclH1SRAB9R50gV2gDFmmusGETkF+LFlQ1JKtWkrXqPPz/9hge1Mzv3FQ76OpmXFpliz1gKJzqk5An1dCU9qEs+6uU0p5Y/KC+HT+zEf38WimhHkn/53goNtvo6qZcWlQP5ucNTUjboO9BXqmq1JiMh44GQgQUTucdkVBfjZb4hS6hjGwIb58PnvoegAiyPP5/fFl/PVWD9crDI2BRxVUJBFfHRPbEES8N1g3alJhAIRWAnFdbGhQuAy74WmlPK53O3w5qUw73qISGTVlHlMz7mKy04eRIdQT+5WtxMuPZxsQUJCRJi2STR3gDFmMbBYROYYY3Y1dpyIPGuMub1Fo1NK+UZVOfz4DHz/FNhC4Zwn2Nrzl/xq9jIGdo5kxmm9fR2hd7iOleg9iURd69r9NommEoTTKQ1tFJGpIrJZRLaJyP0N7B8oIktFpMLZrdbtskopL9j+DfznZFj0Nxh4Hty2nENDpnP96ysJC7Hx8vQxRIT5YS0CrMGAttAjs8FG6lrXXv2XFhEb8G/gLKxus8tFZIExZqPLYXnAHcDFx1FWKdVSig5Y7Q7r34O43nDtfOgzmfKqGma8mE5OUQX/u3k83WLCfR2p9wTZIKbXUSvUrdx92MdB+Za3/xwYC2wzxmQCiMjbwEVA3Re9MSYbyBaR8zwtq5RqAY4aWP4SfPNXqK6ASQ/AKXdBiB1jDPfNW8vK3fk8N20UI3vE+Dpa74tLcRkrYSevpJLKagehwYG51nVLJomG5gXuBrhO5ZEFpLl5vhMpq5Ryx94V8PHdsH+NNcX3uU9CfJ+63U9/tZWP1uzjt1MHcO6wLj4MtBXFpsCupWAMiZFWN9ic4gr/rkE1oSWTxDMNbGsocbi7sp1bZUVkBjADoGfPnm6eWqkAV5YPX/8ZMl6BiCS47FUYcslRa0DMX5XFrK+3cvno7twysU8TJ/MzcSlQWQSluUctY6pJohkikgo8CPRylhPAGGOGYz2Z00CxLKCHy+vuwD4339KtssaYF4AXAFJTU3VpVaWaYgysfQe+eBBKcyFtJpz++2Om916+M4/fzVvHuN5xPHrJMKS9LSB0Ilx6OCVG9QcCe0CdJzWJucB9wDrA3WkRlwP9RCQF2AtcCbi76O2JlFVKAdRUQ952a6qJg+thx/ewNwO6pcI170GXEccU2ZVbwozXM+geG87sa0YH3r14l7ESSb2HAwT0gDpPkkSOMWaBJyc3xlSLyG3A51ijs18xxmwQkZnO/bNFpDOQgTWC2yEidwGDjTGFDZX15P2VCihl+XBwg5UMapNC9s9Q7fwrOCgEEgbC+U/DqOkQdOyXf0FpFdfPWY4BXp4+hpgOAbhScUwvQCBvB3FDQwkOErIDeF0JT5LEQyLyEvA1UJdWjTHvN1XIGHOsAxgAACAASURBVLMQWFhv22yX5wewbiW5VVapgOdwWL1vDq6HA+uP/CzYfeSYDvGQNBTG3GT97DwUOg2A4Ma/9CurHcx8cwV78kp588Y0Ujp1bIUP0waF2K3xEod3EBQkJEQG9gp1niSJ64GBQAhHbjcZoMkkoZRqAcbAD0/D5k8heyNUFlvbJQji+0GPMZB6PXQeZiWFyM5HNUI3f3rDHz9Yz9LMXJ66fARpveO99EHaidiUurESiVGBPaDOkyQxwhgzzGuRKKUa9/NH8PUj0HUUjJxm1QyShkLiIAg58V43z3+Xyf8y9nD75L5cOrrBin1giUuGrV8CkBQZxq7cUt/G40OeJIl0ERmsI56VamUVxfDZ/ZA0DG78EmwtOwb2s/X7efyzTZw3vAt3n9m/Rc/dbsWmQPFBqCwhMSqM5TvzfB2Rz3jy23YqcJ2I7MBqkziqC6xSyksWPQaFe+HyOS2eINZm5XPX/1YzonsMT10+gqCgAOrq2pS6Hk47SYq0c7i0iorqGsL8bf0MN3jyGzfVa1EopRp2cAOk/wdGXQc9xrboqffll3HjaxnEdwzjxV+lYg8JvC/ARrmMlUiKsroJZxdW0COugw+D8g2PZoF1zgRbhtVgXftQSnmDwwEf3wPhMXDmwy166uKKam6Ys5zyyhpevX4MCc7pJ5STy1iJuhXqArQbrNtJQkQuFJGtwA5gMbAT+NRLcSml1vwX9qTDWX+GDnEtdtoah+GOt1axNbuYf00bRf+kyBY7t98IjwV7jDXqunat6wDtBuvJUMq/AOOALcaYFOAM4EevRKVUoCvNgy/+CD3Hw4iWnWjgr59s5JtN2Tx84RAm9k9o0XP7FedssEnOmkSgdoP1JElUGWNygSARCTLGfAuM9FJcSgW2rx6C8gI476kGR0Yfr9eX7uTVH3dywykpXDuuV4ud1y85x0rEdgglxCYcDNAV6jxpuM4XkQjgO2CuiGQD1d4JS6kAtucnWPk6nHw7JA1pkVOWVlbz+KebeG3pLs4YmMiD5w1qkfP6tbgU+HkBQaYmoFeo8yRJXITVaH03MA2IBv7sjaCUClg11VZjdVQ3mNgyK/Yuy8zlvnlr2Z1XyvSTk/nd1IHYtKtr82JTwFENBXtIiAwL2DYJt5KEcynRD40xZ2JNyfGaV6NSKlD99DwcXAdXvAFhESd0qtLKap74bDNzluykZ1wH3p4xjnGBPt2GJ1xng42KYsehEt/G4yNuJQljTI2IlIpItDGmwNtBKRWQCvfBt3+DfmfDoAtO6FTLd+Zx37tr2JlbynXje/G7cwbSIdTbqxX7maPGSqSRnhmYo649+a0pB9aJyJdAXUo1xtzR4lEpFYg+e8C6vXHOEx5NzueqrLKGJ7/YzCs/7qBbTDj//XUaJ/fp1MKBBojILmALc9YkJlJQVkV5VU3ADTr0JEl84nwopVratq9g4wdw+h+O3Obw0Ipdedz77lp2HCrhmnE9eeCcQXQM09rDcQsKgthkyNtBQl/ngLrCCnrGB9aoa7d/g4wx2g6hlDdUlcPC+yC+L5ziecW8vKqGf365hRe/z6RrdDj/vSmNk/tq7aFFxKVY8zc517rOLirXJNEYEVnHsdNwFGCtKvdX5xgKpZSnfnga8jLhVx9CsGfTY6zcfZh7311DZk4JV6f15PfnDiJCaw8tJzYFdnxPUqS1WFMgLj7kyW/Tp0AN8F/n6yuxZoItAOYAJ9bSplQgyt1uJYmhl0HvSW4XK6+q4emvtvDid5l0jrLzxo1jmdBPR0+3uLgUqCqhs60ICMxR154kiVOMMae4vF4nIj8aY04RkWtaOjCl/J4xsPBeq/Yw5VG3i63Zk8//vbuGbdnFXDmmBw+eN4hIe4gXAw1gzh5O0eVZhNqCOBiAk/x5kiQiRCTNGLMMQETGArUduXXktVKe2jAftn9j9WaK7Nzs4RXVNTzz1VZmL95OUpSd124Yq3MveZuzE4Ec3klCZCI5erupSTcBrzin5hCgELhRRDoCj3kjOKX8Vnmh1eW183AYc1Ozh5dWVvPL59NZt7eAK1K784fzBxOltQfvi+kJiHOsRA+tSTTFGLMcGCYi0YAYY/Jddr/T4pEp5c8WPWYtj3nlfyGo+X73f/pwA+v3FfCfaaM4Z1iXVghQAdatwOjuzrESZ7M1u9jXEbU6T9aTiBaRfwJfA1+JyFPOhKGU8sT+tbBsNqReD91HN3v4+yuzmLcii9tP76sJwhdikyEvk6SowJzkz5M5iF8BioArnI9C4FVvBKWU33I44JP/g/A4OONPzR6+PaeYP3ywnrEpcdxxRr9WCFAdI86aMjwhMoyi8mrKKmt8HVGr8qRNoo8x5lKX14+IyOqWDkgpv7bqdcj6CS6eba1+1oTyqhpunbuSsOAgZl15EsG2lltXQnkgNgVKD9Gtg5UcsovK6RXf0cdBtR5PfuvKROTU2hcicgrW1OFKKXeUHIIvH4Jep8CIK5s9/K+fbGTTgSL+ecVIOkfbWyFA1SBnD6eeHAACb0CdJzWJmcDrLu0Qh4HrWj4kpfzUlw9BZbG12lwzE/h9snY/b6bv5ubTenP6wMRWClA1yDlWonPNfiAy4NolPEkShcaYESISBWCMKRSR45uJTKlAs2sprH4TTrkLEpteFW53bin3v7eWkT1iuHfKgFYKUDXKWZOIrdgLDCQ7wJYx9eR203tgJQdjTKFz27yWD0kpP1OaBwtuh+geMPG3TR5aWe3gtrdWgsCzV51EiLZD+J49GsLjsBftJjQ4iGytSRxNRAYCQ4BoEfmFy64oQG+UKtWUylJ46yrI3wXXvA+hTTd4PvHZJtZmFTD7mlH0iAus2UbbtLgU5PAOkqLC9HZTAwYA5wMxHD2JXxHwa28EpZRfqKmGedfDnmVw+auQMqHJw7/aeJCXftjBdeN7MXWojodoU2JTIOsnkiLt2nBdnzHmQ+BDERlvjFnaCjEp1f4ZAx/dCVs+g3OfhCGXNHn4vvwy7p23hiFdo3jg3KbbLJQPxKXAhvfp0ieIjQcDqybhyQ3PS0QkSkRCRORrETmks78q1Yiv/2w1VE/8HYxtusJdXePgjrdWUVXt4F9Xjwq45THbhdgUMA76heaTHWA1CU+SxNnOBuvzgSygP3CfV6JSqj1L/w/88E8YPR0mPdDs4U9/tYWMXYf52y+GkdIpcAZptSvOHk69bdkUVVRTUhE4E197kiRqp5w8F3jLGJPnhXiUat/WzYPP7oeB58N5/2x2PMT3W3N4btF2fpnag4tGdmulIJXHnGMlujsH1AVSN1hPksRHIrIJSAW+FpEEoNmbcyIyVUQ2i8g2Ebm/gf0iIrOc+9eKyCiXfTtFZJ2IrBaRDA9iVar1bfsa5s+0RlRf+nKzs7tmF5Vz9/9W0zchgocvHNJKQarjEtkZgsNJqNoPEFDdYD2ZKvx+EXkca1BdjYiUABc1VUZEbMC/gbOwblEtF5EFxpiNLoedA/RzPtKA/zh/1jrdGHPI3TiV8om9K+B/10LCAGv675Cme4fXOAx3vb2a4opq/vvrcYSHajtEmyYCsclEl2cBcDCAahLujJOYbIz5xnWMhBxdhX6/ieJjgW3GmExnubexEotrkrgIeN0YY4B0EYkRkS7GmP0efA6lfOfQNph7OXSMh2nzIDym2SLPfbuNJdtzeeLS4fRPimyFINUJi0shPDcT0JpEfROBbzgyRsI4f4rzeVNJohuwx+V1FkfXEho7phuw33n+L0TEAM8bY15wI16lWk/hfnjD2b31mvkQ1fz4hmWZuTz91RYuHtmVy1O7ezlA1WJiUwja/i32EAmoAXXujJN4yPn0FuBSINmlnGmojIuGWu3ql2nqmFOMMftEJBH4UkQ2GWO+O6qwyAxgBkDPnj2bCUepFlSWD3Mvg9JcmP4xdOrbbJHc4grueHsVveI78tdLhtWvlau2LC4FqS5jUMcybbhuxAdYtYkqoNjl0ZQsoIfL6+7APnePMcbU/swG5mPdvjqKMeYFY0yqMSY1IUEXhVetpKoc3r4acjbDlW9Ct1HNFnE4DPe+u4bDJVU8e9VJRIR5Mr+m8jlnD6ehHXK1JtGI7saYqR6efznQzzlb7F7gSuDqescsAG5ztlekAQXGmP0i0hEIMsYUOZ+fDfzZw/dXquU5auC9G2HXj1Yvpj6T3Sr28g87+HZzDn++aAhDu+nKv+2Oc6xE/5BD/BhAA+o8SRJLRGSYMWaduwWMMdUichvwOWADXjHGbBCRmc79s4GFWGMvtgGlwPXO4knAfGd1PBj4rzHmMw/iVarlGQOf3AObPoapj8Owy9wqtmr3YR7/bBNTh3Tm2nG9vByk8oroHiBB9JKDWpNwJSLrsNoIgoHrRSQTqMDZcG2MGd5UeWPMQqxE4LpttstzA9zaQLlMYIQbn0Gp1rPoMVgxB069B8bNdKvI+r0F3PzGCjpH23n8suHaDtFeBYdCdHe6OA5QUllDcUV1QNwydOcTnu/1KJRqD356ERY/DiddA2f8ya0in284wF1vryauYyivTB9DdHhI84VU2xWbQvzhvYDVDTYiIcLHAXmfO72bdrVGIEq1aRvmw8L7oP85cP4zzU63YYzhxe8zeezTTQzvHsOLvxpNYqQuv9LuxaUQuW8tYK113VuThFIBzhhruu/3Z0CPNLjsFbA1/d+mqsbBHz9Yz9vL93DesC48dcUIndnVX8SmEFJxmEhKyS4KjHYJTRJKNaTkEKx5G1a+Doc2Q8IguOotCG16tbiC0ip+898V/Lgtl9tO78s9Z/UnKEjbIPyGs4dTTzkYMFOGa5JQqpajBjK/tRLDpoXgqILuY+DCZ2Hopc0uPbort4Tr5yxnT14pT14+gstG62hqv+McK9EvOCdgejhpklAqfw+sngur3oSCPRAeB2NnwKhrIdG9VeKW78xjxusZGODNG9NI6x3v3ZiVbzhrEoPtuawLkFHXmiRUYKquhM0LrVrD9m+sbX1Oh7P/AgPOheAwt081f1UWv5u3ju6x4bw8fYwuHOTPwiKhQyf6mBy+0pqEUn4oexOsegPWvGXNuRTV3Vpi9KRpEOPZ3F/GGJ7+cguzvtnG+N7x/OeaUcR0CPVS4KrNiEuhe+5BcrQmoZSfqCi2urCuegP2LIOgEBh4Lpz0K6v20MziQA0pr6rh3nfX8PHa/VyR2p2/XjyM0GBPpkJT7VZsCknZ33GwpBxjjN8PjtQkodofRw2UF0B5vvWzLN963tjPrAyoLIZO/eHsv8LwKyHi+CeDzCmqYMYbGazanc/vpg5k5sTefv9FoVzEpRBd+S7VleUUV1QTaffvAZKaJFTrMwaqy51f4gX1HvlHvvwbTAAFUFFIk7PU20LBHmMt/mOPgSGXwEnXQo+xzQ6Ca86Wg0Vc/+pycksqmH3NKKYObX79COVnYlMQDN0lh4OFFZoklPJY4T7Y+gXsW330F77ro6ay6XMEh4M92nqEx0BkF6unkeuXf2M/Q8JPOBk0ZPGWHG6buxJ7qI13bh7P8O7Nr0Cn/JDrWImicvom+veoa00S6sQ5amDvStj6OWz5HA5Y0xYQHgcdOzm/6GMhNvnIF/9RD+eXe93rKI96F7WGN9J38fCCDfRPiuTl61LpGhPu65CUrzjHSvSS7IAYUKdJQh2fsnzY/jVs+QK2fWn1FBKbNXXFmQ9DvynWX/7t+F59jcOQnpnL/5bvYcGafUwemMgsXSxIRSRiQjrSqzowpgzX33blHmOsVdi2fm4lht1LwdRYNYS+Z0H/KdD3DOt1O+ZwGFbuPsxHa/bxyboDHCquoGOojVsm9eHeswdg0yk2lAgSm0zvymy+05qECmhV5bDze+sW0tbPIX+3tT1pGJx6l1Vb6J56XF1I2xJjDBv2FfLRmn18vHY/e/PLCAsO4oxBiVwwvCunD0zUCfrU0eJSSM5ZzbwAmORPk4Q6VuYiSJ8NOxZDVanViNx7Epx6N/Q7G6L9Y06irQeL+GjNPj5au58dh0oIDhIm9k/gvikDOHNwkt5WUo2LTaar+ZycgjJfR+J1+r9AHbF/DXz1sDVNRWQXGDkN+k+F5FMhxD/WQtidW8pHa/fx0Zp9bDpQRJDA+D7x3Hxab6YO7awjppV74lIIpYqawn2+jsTrNEkoOLwTvvkrrHvXalOY8jcYc1Ob62F0vA4UlPPxWqvGsGZPPgCje8XyyIVDOGdYZ10MSHkurjcA4cW7/X7UtSaJQFZyCL57Epa/BEHBMOH/4JQ7rW6o7YQxhoKyKvYXlLO/oIz9BeUcKChnX345BwrL2J9fzo7cEoyBod2ieOCcgZw3vAvdY5teF0KpJjm7wXZ2HKCoopooPx5Qp0kiEFWWwNLn4MdnoKrEGo086QGIalujh40x5JdaCeBAYZn1xV9QXpcQap+XVdUcVS5IICnKTpdoO4O6RHHxSd04b3gX+gTAUpOqlUT3wCHB9JKDZBeWa5JQfqKmyprkbtHfofggDDwfzngIEvq3WgiV1Q7ySio5VFxBbkklh4oqyC2pILe4kkPFtdut17nFlVTWOI4qbwsSkiLD6BxtZ1DXKCYPTKRztJ2uMeHWz+hwOkWEEmzTyfaUF9mCqYzoRq/8gxwsrKBvYqSvI/IaTRKBwBj4eQF8/WfI3QY9x8Mv37TmMmph1TUONh8sYuXufLYcKCK3pOLIl39xJQVlVQ2WCw0OIiEijPiIUBIj7QzqHEV8RBidIkLpGhNOl2g7XaLDSYgM07EKqk1wxCTTs2AP2/18QJ0mCX+380f48k+wNwMSBsJVb1s9llqooS23uIKVu/NZtfswK3cfZm1WAaWV1u2fKHswiVF24juGMqhzFJ0iQol3JoJOzgQQ39F6HREW7NeNf8r/hHTqTa/dGSz183UlNEn4q4Mb4KtHrEFwkV3hwn/BiKvAdvz/5FU1DjYfKGLl7sOs3HWYVXvy2ZVbCkBwkDCkaxRXpPbgpJ4xjOoZS/fYcP3iV34rJKEPMVJCQV420MfX4XiNJgl/k78bvn3MWnktLArOfATSbrZmRvVQTlEFK3cfZtXufGctIZ/yKquNIDEyjFE9Y5mW1pOTesYyrFu0jkpWgcXZw4m8ncB4X0biVZok/EVxDnz/JGS8AgicfBuceg90iHP7FDUOQ8bOPD5Zt59vN2ezJ88aTRpiE4Z0jebqsb2sWkKvWLpG27WWoAKbc8rwkMKdvo3DyzRJtHflBbDkX7D031BdZo2SnnS/21Nn1CaGhev2s3D9AXKKKrCHBHFavwSuG5/MST1jGdI1SmsJStUXmwxAROke38bhZZok2quqMmsQ3PdPQdlhGHwxTP4DdOrXbFGHw5Cx67CVGNbtJ7uogrDgICYPTOTcYV2YPDCRjjpvkVJNC+1IcUg8sRV7WbhuP+cM7eyXtWv9Jmhvaqph9Zuw6HEo2gd9JsMZf4KuJzVZzOEwrNh9mE/W7ufT9fs5WGglhtMHJHLecE0MSh2P4E69GZB9iAvnrmR873j+dMFgBnWJ8nVYLUq/FdoLhwM2zodvHoW87dB9DPziBUiZ0ESRhhPDpAEJnDe8K2doYlDqhNgT+zKs6FsendyPf3y9k/Nmfc/VaT35v7MGENvRPyaL1G+Its4Y2PY1fP2ItSxo4mC48i0YcE6DYx1qHIZVuw/zybr9fLruAAcKywkNDmJS/wTOG96FMwbpFNhKtZie45E1bzHtx3O4dNx0ZhVO5Pmf9vDRmv3cfWY/rhnXq92P/hdjjK9jaDGpqakmIyPD12G0nN3LrOSw60eI6QWnPwjDLjtmkZ/sonK+23KIxVty+H5rDvmlVZoYlGoNxsCO7yD9OdjyGdhCKeh7MX8/PIm3dsfQPymCP50/hFP7dfJ1pE0SkRXGmNQG92mSaIMOrIdv/mL90nVMhIm/hVHXQbBVfa2qcbBi12EWb8lh8eYcNu4vBKBTRBgT+ycwcUACpw9IINKPJx1Tqs05tA2WzYbVc6GqlNyENJ4oOIN3Cgdz5uAu/OG8QfSK7+jrKBukSaKtMwbyMmHvStjyKax/H+xR1rTdaTMhtCNZh0vrksKS7bkUV1QTHCSM7hXLxAEJTOyfwKDOUQTpvEZK+VbZYVjxGvz0AhTuJT+8B/8qOZN5NRO5esIgbj29b5trC9Qk0dYUZ8PeFS6PlVBuLYZDaASM/TXlY29n2QEHizfnsHhLNttzSgDoFhNelxRO7hOvtQWl2qqaKmtizaXPwd4MSoMieL1yEp/Yz+f6cydw8chubeaPOp8mCRGZCjwD2ICXjDF/r7dfnPvPBUqB6caYle6Ura9NJomKYti/+uiEUOAcfCM2TOJgyhJHkBc9lN3hg/i5uivfbc8nPTOXimoHocFBjOsdb91G6p9An4SOftkXWym/tucnSH8Os3EBDmP4tGYs38VfwdWXXsrIHjG+js53SUJEbMAW4CwgC1gOXGWM2ehyzLnA7VhJIg14xhiT5k7Z+nyeJGqqrIn19q3EZK2gJmsFttzNiLHmOyqwd2OnfSCbgvqzoro3S0u7sbcEHPX+CXondKxLCmkp8YSH6mhnpfxC/h7MsuepWj6H0OoiVjr6sqHHNE6+8AbiIjoQHmojLDio1f8Q9GWSGA88bIyZ4nz9AIAx5jGXY54HFhlj3nK+3gxMApKbK1vf8SaJivIS1rz7GLaaCqSmgiBHJUE1lQTVVGBzWK9tNZXYTAU2RyXBjkpsjgqCHVUEG+t1sKkkxFQShJUQ8kwkqx19WOPowxrThzWO3hwmiriOoSRGhpEUZScpqvanve5152i7rrmslL+rKKY84w1Kv/sXcRVZVBob5YRRQQjlhFJJCNUSSpWEUh0URnVQGI6gUBy2MBzBYRibHRMchgTbIcROULCdkNhupF5w83GF01SS8HbrSTfAdWKTLKzaQnPHdHOzLCIyA5gB0LNnz+MKsqqqirHbnwWgwoRQQQgVBFPh/MeqJIRSQqkS63mVRFEpIVQ7t1XZwqgihOqgMPI69qYgdgTB8ckkRdsZEGXnNGcySIgMIyxYawVKBbywCOyn3IJ9/M3krFxAzsbFmOpyqKqw5mCrriCopoLgmgrCasqxOYqw1VQQXGX9MRpiKgmlklBTRYhY67dss/WB40wSTfF2kmiozlS/6tLYMe6UxRjzAvACWDUJTwME6NgxioL/yyI4JAybzUZ4kBAh0mYalZRSfiooiITUi0lIvfi4T2FqqqgoLyWxquFVH0+Ut5NEFtDD5XV3YJ+bx4S6UbZFSFAQ0ZH+u0atUsp/iS0Ee8dovHWT2tvjxZcD/UQkRURCgSuBBfWOWQD8SizjgAJjzH43yyqllPIir9YkjDHVInIb8DlWN9ZXjDEbRGSmc/9sYCFWz6ZtWF1gr2+qrDfjVUopdTQdTKeUUgGuqd5N7Xt6QqWUUl6lSUIppVSjNEkopZRqlCYJpZRSjfKrhmsRyQF2uWzqBBzyUThtjV4Li16HI/RaHBHo16KXMSahoR1+lSTqE5GMxlrsA41eC4tehyP0Whyh16JxertJKaVUozRJKKWUapS/J4kXfB1AG6LXwqLX4Qi9FkfotWiEX7dJKKWUOjH+XpNQSil1AjRJKKWUalS7TRIi0kNEvhWRn0Vkg4jc6dz+DxHZJCJrRWS+iMS4lHlARLaJyGYRmeK76FtWY9fCZf+9ImJEpJPLNr+7Fk1dBxG53flZN4jIEy7b/e46QJP/P0aKSLqIrBaRDBEZ61LGX6+FXUR+EpE1zmvxiHN7nIh8KSJbnT9jXcr45bU4LsaYdvkAugCjnM8jgS3AYOBsINi5/XHgcefzwcAaIAxIAbYDNl9/Dm9eC+frHljTre8COvnztWjid+J04CsgzLkv0Z+vQzPX4gvgHOf2c7HWl/f3ayFAhPN5CLAMGAc8Adzv3H5/IHxXHM+j3dYkjDH7jTErnc+LgJ+BbsaYL4wx1c7D0rFWtAO4CHjbGFNhjNmBtX7F2PrnbY8auxbO3U8Dv+XopV/98lo0cR1uAf5ujKlw7st2FvHL6wBNXgsDRDkPi+bIao/+fC2MMabY+TLE+TBYn/k15/bXgNo1RP32WhyPdpskXIlIMnAS1l8Irm4APnU+7wbscdmXxZEvUr/hei1E5EJgrzFmTb3D/P5a1Pud6A9MEJFlIrJYRMY4D/P76wDHXIu7gH+IyB7gSeAB52F+fS1ExCYiq4Fs4EtjzDIgyVirYOL8meg83K+vhafafZIQkQjgPeAuY0yhy/YHgWpgbu2mBor7Vf9f12uB9dkfBP7U0KENbPOba9HA70QwEIt1i+E+4B0REfz8OkCD1+IW4G5jTA/gbuDl2kMbKO4318IYU2OMGYl1Z2GsiAxt4nC/vhaeatdJQkRCsP4DzDXGvO+y/TrgfGCacd5kxPproIdL8e4cqWq3ew1ciz5Y91PXiMhOrM+7UkQ648fXopHfiSzgfedth58AB9aEbn57HaDRa3EdUPv8XY7cRvHra1HLGJMPLAKmAgdFpAuA82ftbciAuBbuardJwvmX4MvAz8aYf7psnwr8DrjQGFPqUmQBcKWIhIlICtAP+Kk1Y/aWhq6FMWadMSbRGJNsjEnG+sUfZYw5gJ9ei8Z+J4APgMnOY/oDoVgzfvrldYAmr8U+YKLz+WRgq/O5P1+LhNpejiISDpwJbML6zNc5D7sO+ND53G+vxfEI9nUAJ+AU4FpgnfNeI8DvgVlYvRK+tP6fkG6MmWmM2SAi7wAbsW7F3GqMqfFB3N7Q4LUwxixs6GA/vhaN/U68ArwiIuuBSuA6Zw3TX68DNH4tfg08IyLBQDkwA/z6dwKsnl6viYgN6w/jd4wxH4vIUqxbjzcCu4HLwe+vhcd0Wg6llFKNare3m5RSSnmfJgmllFKN0iShlFKqUZoklFJKNUqThFJKqUZpklCtwjktwioR+diNYxc5Z99cK9aMvv8Sl9l8myi30KU/fHFzbJRFrgAACX9JREFUxzdQ/m4RKReRaJdtqSIyy42yySJytafv6Qsi8nsfvGe7uT7qaJokVGu5E2uSOXdNM8YMB4YDFRwZ6NQoY8y5zhG1x+sqYDlwics5M4wxd9Q/0DnOwFUy0Ca+BJ3jAZricZJo4PN6Kpk2cn2UZzRJKK8Tke7AecBLLtuCRWS5iExyvn5MRB6tX9YYU4k1i21PERnhPPYDEVnhXBtghss5d4rLmhnObW+IyEUur+c6Jz6sH2MfIAL4A1ayqN0+qbb2IyIPi8gLIvIF8Hq9U/wdaxLB1c4aiU2stU2WO2tEN7ucb7GIvCMiW0Tk7yIyTaz1DtY540BE5ojIbBH53nnc+c7tTZ33WxH5L7CuseskIn8Hwp1xznX+hb/e5fPeKyIPO58vEpG/ichi4E7nyOX3nO+9XEROaeA6JjtjXul8nNzQ9alfTrVhvp6rXB/+/wDmAaOBScDHLtuHYNUuzgJWAaHO7YuA1Hrn+AD4pfN5nPNnOLAeiHe+3smRNTOKnT8n8v/bO/9YL6s6jr/ejF9B4YzCMi1wMxG9YnJ1MSMt0aVttRQFp64Ny+koq41aLpgXW1nTTVP8QTK4TqzAEAe5BKsBM3N4r/HjAiGTXO3ScBU5wRzQ/fTH5zzc5359nu/3e2m37/3C57V9x3nOc55zznPueM45n3M+7wPPpPBJwJ9J541U5D8PmI8PnF6n98yJo3UG2oBO4D0Fz1e+2y3AvBQeAXTgWlqXAv/CvYBHAN3AgpTuG8D9KdwOPJfqcyYuqzKyRr4HgQm5OpS104FcmvFAV+56LtCW+zs8nLv3M+BTKfxRXPKjsh1GASNT+Eygo6h94tc8v2aW5QiagDQCfsPMOrNZQ4a5/METwBpgqvmsoTSrXPh2SZlJ6HT8Y/SPoofMbIOkhySNA64GVlrveSN5ZgFfMrMeSU/jEg0PFaRbbWb/rlLPjCuA8yTNSNcnpXoeAl62JFEt6TX8ICDwGcBncnmsMLMeYLekPcDEGvluMj//IKPudqrC8lx4OjBJOvqnGCPpfebnVWQMAxZKOh/4Dy7THjQx0UkEA83FwBckXYWPhMdIWmZmN6b7LfjI+pSyDJKNvQXYmTqa6Xin8rak9SnfajwB3IB3BLML8j8P/4Bmel/DgT0UdxIHa5R1NFvg62a2tqKsS/E1loye3HUPff9PVmrmWI18D1Zc19NOR+hrdq5Mk3/fISm/ap3kt4B9wOSU/p0qaYMmINYkggHFzO4ws9PMlWhnAb/LOghJVwNjgU8DDxTtYJLLXd8N/NXMtuIj5/3pwzcRPyOiFu34GRuY2faC+9fjJpbx6Xcq8BFJH+vHq76FHxOasRa4LdUfSR+XNLof+QFcK2lIWqc4A9jVj3yrtdPh7Hn8gz5O0lhJI3CJ/TLWAV/LLtJsoajcv6UZ0E1Atohe2T5BkxCdRNAQ0gLzj4CbzexVYCHwk1ySJyVtxW3po/EjJcHt9EPTve/jR9RWxcz24WsfS0uSzAJWVcStSvH1shU4ImlLWphdjKuIvpIWhhfR/5n7LmADfrrirWb2Tj/yrdZOPwW2SnrSzA4Dd+Gn1v0Kl9Au43agNS2Y7wBuLUjzMPBlSS/hpqZsJlLZPkGTECqwwXGPpFG4vf8CM3uz0fWpB0nt+ELvLxtdl+DEJmYSwXGNpOyAmQebpYMIgsFEzCSCIAiCUmImEQxaJC2WNOn/WF57bmvpgOYpabKkPyQHujWSxqT44ZKWpvgt+W3D6pUr2Zx+4wrybZPULemukvq8y+GwH+/SJmluCo+U9LykOyVlznmHjjXvYPASW2CDQYuZfaXRdRhAFgNzkx/HbODbuDPfVwHMrCV1Ar+WdGHaLQQuV9JRI+/7zOzegaq4pOHASqDTzBak6PMlvT5QZQaNI2YSQcORNFrSs2nk3CVpZopfL6k1hW+Wy1Osl/SYpIUpvl3SI3JJij2SLpG0RNLOtPiblfGIpA65RMWCwor0rdMUuXxGp6S1kj4s6WxJm3JpxqfdQ4XpaxRxFrAxhZ8HrknhScBvAczsDdyHpLVmI5a/x1hJ6+TiiovIOSVKulEuB7JZ0iIlzaeytk4MBX4B7Daz7x5rvYLmITqJYDDwOWCvmU02s3Px7ZtHkXQqPsr+JC7hMbHi+ZOBz+KOXGuA+3DJj5bcXv7vmVkrLhh4idyBrpDkQ/AgMMPMpgBLgB+Y2U5guKQzUtKZwIqy9DXeuQvINKSuxT2iAbYAX5RrW03A5UxOzz23NH3U50vKe6GXcSfwgpl9AliNy2kg6exU/4vNLPOOvqGOtv4OcMTMvllH2cFxQHQSwWBgGzBd0o8lTSvYhXQRsMHM/pn29T9VcX+N+Q6MbcA+M9uWzDPbcW0igOskvYJrRJ2Dj9jLOAs4F/fA3ozrOp2W7q0ArkvhmbhsRbX0ZcwG5kjqxJ3MMkmSJbhOUwdwP/Ai7hUNbmpqAaal3001ygB3VFwGYGbPAvtT/GV4B/RyqvNluMNerbZ+AZgqKeQ2ThBiTSJoOGb2qqQpwFXA3ZLWmVl+4bXWiDkva1EpeZGNyOcCF5rZ/mSGqiblIWC7mU0tuLcceEqu72RmtltSS5X0hZjZn3AdJtIH9/Mp/gg+IyLdexHYne51p3/fkqu9XsS71WgLiyuIE/C4md3RJ7JX66mMjcDj+FrJNDPbW0f5QRMTM4mg4SQTx9tmtgy4F7igIskm3ER0svxcg2sq86jBGNzz901JpwBX1ki/C/igpKmpfsMknQNgZq/hppn59IrflaYvI9uZJGkIPvN4NF2PUpLZkHQ5btrZkcxPH8jyx+Uzugoz78tGXLcKSVfipjnwdY8ZuXq8Xy5DUrOtzWwlcA/wnOo4DCpobmImEQwGWoB7JPUAh4Hb8jfNrFvSD3HpiL24LEXdjnFmtkXSH3Hz0x7g9zXSH5JvW31AfkrdUNz0k+k+Lcc/khPqTF/E9ZLmpPDT9EqGjAPWprboptekNCLFD8P1kH4DPFbH6y8Afp5MbRuAv6Q675A0D1iXOqrDwBwze6metjazRyV9CFgt6YokGRIch4QzXdAUSHqvmR1Io9tVwBIzq9RbOuGRHxh04H/ZAnusbS3fAttqZn8/1rKDwUeYm4JmoS0tsHbhBwc90+D6DFYOALeoxJmuTvrV1krOdPhZEj3V0gbNR8wkgiAIglJiJhEEQRCUEp1EEARBUEp0EkEQBEEp0UkEQRAEpUQnEQRBEJQSnUQQBEFQyn8BL05rKyhOfsMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "temp_hist_area_weighted.plot(label=\"area_weighted\");\n",
    "temp_hist_flat.plot(label=\"flat\");\n",
    "plt.gca().legend(ncol=2);"
   ]
  }
 ],
 "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": 2
}