warming-stripes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Warming Stripes\n",
    "\n",
    "Let's create our own figure of \"[warming stripes](https://en.wikipedia.org/wiki/Warming_stripes)\" that use a red/blue color scale to represent annual air temperature deviations from some climate mean tempearture over a long timeseries. This notebook is based on, and adapts code from [\"Creating the Warming Stripes in Matplotlib\" by Maximilian Nöthe](https://matplotlib.org/matplotblog/posts/warming-stripes/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will use air temperature from [NASA's DayMet dataset](https://daymet.ornl.gov/), and access it through an API using the [ulmo](https://github.com/ulmo-dev/ulmo/) python package.\n",
    "\n",
    "Since we don't have ulmo installed on our JupyterHub by default, the first time you run this notebook, you'll need to install the ulmo package into our python envrionment.\n",
    "\n",
    "**Run the cell below once.** Uncomment it now, so that it reads `!pip install ulmo`, and run it. It will install ulmo on your JupyterHub. You shouldn't need to re-run it next time you open this notebook, so it's best to comment it out again when it completes.\n",
    "\n",
    "When it is complete and you've commented it back out, restart the kernel by going to *Kernel* > *Restart Kernel*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!pip install ulmo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# packages for working with our data\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "# page we'll use to access the data\n",
    "import ulmo\n",
    "\n",
    "# packages for plotting\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.patches import Rectangle\n",
    "from matplotlib.collections import PatchCollection\n",
    "from matplotlib.colors import ListedColormap"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get a pandas dataframe of DayMet data, the function we're using is [`ulmo.nasa.daymet.get_daymet_singlepixel()`](https://ulmo.readthedocs.io/en/latest/api.html#module-ulmo.nasa.daymet)\n",
    "\n",
    "We need to specify a longitude and latitude, the variables we want to download, what years we want (by saying years=None, it will respond with all years on record), and we specify that we want the result returned as a pandas dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# latitude and longitude for Seattle, WA\n",
    "latitude = 47.654\n",
    "longitude = -122.308"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see what data variables we can request, use the function [`ulmo.nasa.daymet.get_variables()`](https://ulmo.readthedocs.io/en/latest/api.html#ulmo.nasa.daymet.get_variables)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'tmax': 'maximum temperature',\n",
       " 'tmin': 'minimum temperature',\n",
       " 'srad': 'shortwave radiation',\n",
       " 'vp': 'vapor pressure',\n",
       " 'swe': 'snow-water equivalent',\n",
       " 'prcp': 'precipitation',\n",
       " 'dayl': 'daylength'}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ulmo.nasa.daymet.get_variables()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're just interested in temperatures, so we'll grab tmin and tmax:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "making request for latitude, longitude: 47.654, -122.308\n"
     ]
    }
   ],
   "source": [
    "df = ulmo.nasa.daymet.get_daymet_singlepixel(latitude, longitude, variables=['tmax', 'tmin'], years=None, as_dataframe=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take a look at what the dataframe looks like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>year</th>\n",
       "      <th>yday</th>\n",
       "      <th>tmax</th>\n",
       "      <th>tmin</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1980-01-01</th>\n",
       "      <td>1980</td>\n",
       "      <td>1</td>\n",
       "      <td>11.47</td>\n",
       "      <td>6.35</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1980-01-02</th>\n",
       "      <td>1980</td>\n",
       "      <td>2</td>\n",
       "      <td>8.97</td>\n",
       "      <td>4.41</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1980-01-03</th>\n",
       "      <td>1980</td>\n",
       "      <td>3</td>\n",
       "      <td>8.18</td>\n",
       "      <td>0.62</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1980-01-04</th>\n",
       "      <td>1980</td>\n",
       "      <td>4</td>\n",
       "      <td>4.99</td>\n",
       "      <td>0.72</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1980-01-05</th>\n",
       "      <td>1980</td>\n",
       "      <td>5</td>\n",
       "      <td>4.55</td>\n",
       "      <td>-0.35</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            year  yday   tmax  tmin\n",
       "1980-01-01  1980     1  11.47  6.35\n",
       "1980-01-02  1980     2   8.97  4.41\n",
       "1980-01-03  1980     3   8.18  0.62\n",
       "1980-01-04  1980     4   4.99  0.72\n",
       "1980-01-05  1980     5   4.55 -0.35"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a daily mean temperature from tmin and tmax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>year</th>\n",
       "      <th>yday</th>\n",
       "      <th>tmax</th>\n",
       "      <th>tmin</th>\n",
       "      <th>tmean</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1980-01-01</th>\n",
       "      <td>1980</td>\n",
       "      <td>1</td>\n",
       "      <td>11.47</td>\n",
       "      <td>6.35</td>\n",
       "      <td>8.910</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1980-01-02</th>\n",
       "      <td>1980</td>\n",
       "      <td>2</td>\n",
       "      <td>8.97</td>\n",
       "      <td>4.41</td>\n",
       "      <td>6.690</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1980-01-03</th>\n",
       "      <td>1980</td>\n",
       "      <td>3</td>\n",
       "      <td>8.18</td>\n",
       "      <td>0.62</td>\n",
       "      <td>4.400</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1980-01-04</th>\n",
       "      <td>1980</td>\n",
       "      <td>4</td>\n",
       "      <td>4.99</td>\n",
       "      <td>0.72</td>\n",
       "      <td>2.855</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1980-01-05</th>\n",
       "      <td>1980</td>\n",
       "      <td>5</td>\n",
       "      <td>4.55</td>\n",
       "      <td>-0.35</td>\n",
       "      <td>2.100</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            year  yday   tmax  tmin  tmean\n",
       "1980-01-01  1980     1  11.47  6.35  8.910\n",
       "1980-01-02  1980     2   8.97  4.41  6.690\n",
       "1980-01-03  1980     3   8.18  0.62  4.400\n",
       "1980-01-04  1980     4   4.99  0.72  2.855\n",
       "1980-01-05  1980     5   4.55 -0.35  2.100"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df['tmean'] = np.mean([df.tmax, df.tmin], axis=0)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Resample the dataframe to annual mean values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>year</th>\n",
       "      <th>yday</th>\n",
       "      <th>tmax</th>\n",
       "      <th>tmin</th>\n",
       "      <th>tmean</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1980-12-31</th>\n",
       "      <td>1980.0</td>\n",
       "      <td>183.0</td>\n",
       "      <td>15.788219</td>\n",
       "      <td>6.550685</td>\n",
       "      <td>11.169452</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1981-12-31</th>\n",
       "      <td>1981.0</td>\n",
       "      <td>183.0</td>\n",
       "      <td>16.461370</td>\n",
       "      <td>6.794301</td>\n",
       "      <td>11.627836</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1982-12-31</th>\n",
       "      <td>1982.0</td>\n",
       "      <td>183.0</td>\n",
       "      <td>15.615945</td>\n",
       "      <td>5.833233</td>\n",
       "      <td>10.724589</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1983-12-31</th>\n",
       "      <td>1983.0</td>\n",
       "      <td>183.0</td>\n",
       "      <td>16.008548</td>\n",
       "      <td>6.499671</td>\n",
       "      <td>11.254110</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1984-12-31</th>\n",
       "      <td>1984.0</td>\n",
       "      <td>183.0</td>\n",
       "      <td>15.613151</td>\n",
       "      <td>5.925534</td>\n",
       "      <td>10.769342</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              year   yday       tmax      tmin      tmean\n",
       "1980-12-31  1980.0  183.0  15.788219  6.550685  11.169452\n",
       "1981-12-31  1981.0  183.0  16.461370  6.794301  11.627836\n",
       "1982-12-31  1982.0  183.0  15.615945  5.833233  10.724589\n",
       "1983-12-31  1983.0  183.0  16.008548  6.499671  11.254110\n",
       "1984-12-31  1984.0  183.0  15.613151  5.925534  10.769342"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_annual = df.resample('Y').mean()\n",
    "df_annual.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now find the overall mean for the years up until 2020, which was the data available when I first published this notebook. We'll use this as our \"climate mean\" air temperature to compare against just for this example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11.365052455730037\n"
     ]
    }
   ],
   "source": [
    "climate_mean = df_annual.tmean[:'2020'].mean()\n",
    "print(climate_mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the annual anomaly between the annual mean temperature and this climate mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>year</th>\n",
       "      <th>yday</th>\n",
       "      <th>tmax</th>\n",
       "      <th>tmin</th>\n",
       "      <th>tmean</th>\n",
       "      <th>anomaly</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1980-12-31</th>\n",
       "      <td>1980.0</td>\n",
       "      <td>183.0</td>\n",
       "      <td>15.788219</td>\n",
       "      <td>6.550685</td>\n",
       "      <td>11.169452</td>\n",
       "      <td>-0.195600</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1981-12-31</th>\n",
       "      <td>1981.0</td>\n",
       "      <td>183.0</td>\n",
       "      <td>16.461370</td>\n",
       "      <td>6.794301</td>\n",
       "      <td>11.627836</td>\n",
       "      <td>0.262783</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1982-12-31</th>\n",
       "      <td>1982.0</td>\n",
       "      <td>183.0</td>\n",
       "      <td>15.615945</td>\n",
       "      <td>5.833233</td>\n",
       "      <td>10.724589</td>\n",
       "      <td>-0.640463</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1983-12-31</th>\n",
       "      <td>1983.0</td>\n",
       "      <td>183.0</td>\n",
       "      <td>16.008548</td>\n",
       "      <td>6.499671</td>\n",
       "      <td>11.254110</td>\n",
       "      <td>-0.110943</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1984-12-31</th>\n",
       "      <td>1984.0</td>\n",
       "      <td>183.0</td>\n",
       "      <td>15.613151</td>\n",
       "      <td>5.925534</td>\n",
       "      <td>10.769342</td>\n",
       "      <td>-0.595710</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "              year   yday       tmax      tmin      tmean   anomaly\n",
       "1980-12-31  1980.0  183.0  15.788219  6.550685  11.169452 -0.195600\n",
       "1981-12-31  1981.0  183.0  16.461370  6.794301  11.627836  0.262783\n",
       "1982-12-31  1982.0  183.0  15.615945  5.833233  10.724589 -0.640463\n",
       "1983-12-31  1983.0  183.0  16.008548  6.499671  11.254110 -0.110943\n",
       "1984-12-31  1984.0  183.0  15.613151  5.925534  10.769342 -0.595710"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_annual['anomaly'] = df_annual.tmean - climate_mean\n",
    "df_annual.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have our data all in order, we can go ahead and make the plot. (The following plotting code is adapted from [\"Creating the Warming Stripes in Matplotlib\" by Maximilian Nöthe](https://matplotlib.org/matplotblog/posts/warming-stripes/))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is our custom colormap from [colorbrewer2](http://colorbrewer2.org), we could also use one of the [colormaps](https://matplotlib.org/3.1.0/tutorials/colors/colormaps.html) that come with `matplotlib`, e.g. `coolwarm` or `RdBu`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap = ListedColormap([\n",
    "    '#08306b', '#08519c', '#2171b5', '#4292c6',\n",
    "    '#6baed6', '#9ecae1', '#c6dbef', '#deebf7',\n",
    "    '#fee0d2', '#fcbba1', '#fc9272', '#fb6a4a',\n",
    "    '#ef3b2c', '#cb181d', '#a50f15', '#67000d',\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we create bars for each year as a [PatchCollection](https://matplotlib.org/3.3.2/gallery/shapes_and_collections/patch_collection.html) of [Rectangles](https://matplotlib.org/3.2.1/api/_as_gen/matplotlib.patches.Rectangle.html), make the plot and save it as a jpg image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzIAAAB4CAYAAAApHse6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABLZklEQVR4nO3deXxU1fk/8M+5s2Yy2feQEAIh7BBACQRB3FBqQSluqD9Fa79K3Vpbl7q01qpVW2rrV1v12yqgFrGKC7hgEWVRNmVfwpaEkIXs+2RmMnPP749hkrlz753MJDOThDzv18tXm3vPnDkks9znnuc8h3HOOQghhBBCCCFkABH6egCEEEIIIYQQEigKZAghhBBCCCEDDgUyhBBCCCGEkAGHAhlCCCGEEELIgEOBDCGEEEIIIWTAoUCGEEIIIYQQMuBQIEMIIYQQQggZcCiQIYQQQgghhAw4FMgQQgghhBBCBhytvw23nWgMqONpJZ/43fb+lvMD6vtXs7IDar9sS7Hfbf/11CsB9W356E6/24qfLA+o79LFzwTU/pviar/bjo6LDqjv1FhjQO2HrnrM77Y75z0RUN+BvLb4h6sD6vutJf8IqP0tTd/43/jIvoD6Xjn9/oDah1Igr5dpn/8hoL55YWFA7dnC6/1uK4ycHFDf2zsyA2r/+OdH/G779LwxAfU9XXfa77Yj/uZ/WwB45+cFAbUP5P1fNyHX5/mMJx5F6s/vUjy3f+3n+MdVN4Bzrvr4O1Jj/R7L1O2f+t0WAJwvPOJ320BehwACfv9jzCS/mwb6OmdZEwIbCyGE9FM0I0MIISQszPnTkHLX/6ienzh/Hi57sP8E8IQQQvo3CmQIIYSEnGAyYdhfl4EJXV87os2GlppaSburnvktcmYFNmNECCFkcKJAhhBCSMhl/O5xGIdlSY5VvLAM/1p8O0RR7Dym0Wpxx7tvIio5KdxDJIQQMsBQIEMIISSkoufMRvItN0uOtezchTOvvo7Cr77Buif/KDkXm56Gn/77X5LZG0IIIcQbfUsQQggJGU1MDIb95c+SY06LBSX3PwCcnYn5/OkXcGj9V5I2oy+Zgx8/+WjYxkkIIWTgoUCGEEJIyAx95vfQp6VKjpX9/mnYSk51/sw5x5s334GGsnJJuyufeAhjL780LOMkhBAy8FAgQwghJCRir5yHhEU/kRxr+mYTala+LWvbWluH/7vuVjg7OiTHb3/7/xCXMSSk4ySEEDIwUSBDCCEk6LSJich6Qbr2xdHUhJIHHlR9TNG2nVjz8G8lx8yJCbhj9XIIWr+3PSOEEDJIUCBDCCEk6LL+9Bx08fGSY6WP/w4dlWd8Pu6rF1/BnjXSTW9HFOTjJ88/FfQxEkIIGdgokCGEEBJUCddeg7gr5kqONXz2OerfX+PX41fefjdqThZLjl36wD3IWzg/aGMkhBAy8FEgQwghJGj0Q9KR+fSTkmMdtbU49dBv/O6jvakJr197CzqsVsnxW9/8Owxee9EQQggZvCiQIYQQEjTD/vJnaKOjJcdOPfgIHHX1AfVzes8+vHf/w5JjETExGPH6P8CMhl6PkxBCyMBHgQwhhJCgiDJqED37Asmx2vfeR+MXX/aovy2vv4kdb6+WHDNNGI/Uu/6nx2MkhBBy7qBAhhBCSK9pBYb4SGllMXt5BU4/8WSv+v33Xb9A5eFCybGE667tVZ+EEELODRTIEEII6bU4kxYCY5JjJQ88CGdzc6/6tbW1YeXtd0uOGbOHwTRhfK/6JYQQMvBRIEMIIaRXGACTQfp1Uv3mCjRv3hKU/ot37JJVMYub/+Og9E0IIWTgokCGEEJIr0ToBclsDHc6UbHsr0F9jh/ek5Zujl9AgQwhhAx2FMgQQgjplUiDRvJzy7fb4KirC+pz/PDeh5KfDVlDYZo0MajPQQghZGChQIYQQkiPMQAmvfSrpH7dp0F/ntN796P6+EnJsXhKLyOEkEGNAhlCCCE9ppRW1vjZFyF5Lu/0srgFV4bkeQghhAwMFMgQQgjpsXCklbl9751elpmJyLxJIXkuQggh/R8FMoQQQnokXGllbuX7D8J6QppeFkeL/gkhZNCiQIYQ0jfSsyEsXArNg38BG0V31QciWVoZ5yFLK3OrX7tO8nPcfEovI4SQwYoCGUJI2LHxMyBcvRQsIwds2CgIS38PZI7o62GRAHmnlVk7xJCllbk1fCINZAwZGYicnBfS5ySEENI/USBDCAkfQQM2ZxGEi64B03RdBDOdDpqf/gaIiOy7sZGAKKWVtdnFkD9ve+FRtB8/LjnW6/QyUzSQNhKaX/0ZSB/Wu74IIYSEDQUyhJDwMEZCuPpOCBMKFE+zpDQIN/8ivGMiPaaUVtZmc4bluRvWStfhxM2/EvAYS0B0RiBxKJjeCDZ8DDT3PQPEJARhlIQQQkKNAhlCSOglpkG4/hdgQ3ynjwl5BWAXXRWmQZHeUEorE3l4nrveO71syBBETpncs87iUsE8giAWFQvNTx8GBI2PBxFCCOkPKJAh5yZBA5YzFZoHXgC7+Oq+Hs3gNmIChGvuBYuOlxzmjg6Imz4Eb22SHBcW3g4MGxXOEZIA9VVamZv16DG0H5Oml8X3JL3MEAkWESU7zEaMg3D1kh6OjhBCSLhQIEOCKtIgIClKBzY2v+8GwRiECbMhDBsHNmIcNIt+BjahD8czaDGwaXOh+dESMJ1Bcoa3NEJ8/2Xw/VshrlgmfZRGC81PHwEi5ReYpH/oy7QyN+9F/3E/7kF6WWyK6inhkp+A5SmnQRJCCOkfKJAhQWPSC0iO0sNs0EC45Dqw0VP7ZBwsdxpYUqb0WP4lfTKWQUunh/CjWyDkXy47xSuKIb73V6CmzPXz4R8gfrFa0obFJ0O45Vc9X/dAQqov08rc6r3WyejT02A+L4DPnIgoMKPv4hLCzb8EktJ7MjxCCCFhQIEMCZr4SK3kZzblorCPgWWNh5ApT0tiY6cCOn3YxzMoRce7UslGTJSdEg/tgPjhPwBLi/T4p2+DH9svOSaMPx/s0mtCOlQSuL5OK3OzHjuG9sKjkmMBVS+LTZX8yDtscL73D8kxFmGC5mePAl4zioQQQvoHCmRIUJgNGug00pcTS0gDUoeFbQwsNRvCyCnK5wxGsFF5YRvLYBUdoYVw3S/AEqV3sbnodK2H2fgeICqkIIkinG++AN7cIDksLPh/QM74UA6ZBKg/pJW5eS/69zu9LDIWTG+UHmusAt+0DuKOrySH2ZBsCDf8vOeD1BnBsieAzZwPJNLsDiGEBJO2+yaEdC/WpPxSYuOng58pCf0A4lLBxvrOZ2cTp4Mf3Bn6scifGYhOALLHAaePAw57H4wh9MxGDcakR0oqQAEAt1ogfrHS9W/3pbkB4psvQLj3GTDBFRQzQQPNbQ/B+dx9QEtjiEZOAtEf0srcGtZ9iiEP/arzZ31qCszTzkfrDl/vcyZbG8Pt7YDFVXRCfPcVsIwRYEOGdZ4Xpl8KfvIQ+HdfBjbA5KEQRk/vDJp43mzwg9vAt38O2NoD6+scU11djd27d6O8vBxNTU1oa2uD0WiE0WhEQkIC0tLSMGTIEGRnZ0On0/X1cEkAli9fjttuu03xXFZWFkpKSsI7IHJOo0AmXLR6sJxJ4JXFQFtzX48mqKKMGug0yndB2cg88C0fh/ZL2xwLYdIcMK9yqbypBiwmqWssE/OBVQLAw5AGI2iAhHSwpEywxCFg+gjXmBprIa5+EbBbQz+GMMtOjJDcqQcAXncG4qdvAE3+7fbOj+2H+Ok70Mz/f53HWGwChCUPQnz5ifD87Yiq/pJW5mY9fgKWw0dgGjum81jc/Ct9BzJR8WBarzTThqqu/2+3wfnPZ6F56K9gEabOw8J1S+E8fRI4fbL7gekMYKOmQUjNlhxmggA2cSb4yDzwbZ+CH94J8D6KAvtARUUF/vWvf2HFihU4edKP3yMArVaL4cOHY8KECcjPz0d+fj5mz54d4pESEho33ngjVq1apXguOzsbRUVFfvWTnJyMmpqaXveTnp6OyspKxXMPP/wwnnvuOb/6+elPf4o33njDZxtBEFBSUoLMzEyf7QJFqWXhEBkLpI+EMO8WCNf9AohJ7OsRBVVshHo8zLQ6sFHK6V5BYTBByLtEdmEilh6GeHCLdCxRscDw0aEbi94Ilp4DYdJFEC68HppJF0FIz+kMYgCAxSaCjZ8RujH0kTiTFmaj9HXAiw5B/M9LfgcxnY9bvxri4R8kx4TReRDm3dDrcZLe6U9pZW4NaxXSywSVrzYmADHJkkPc2gZYpWu2UF0O8Z2/SR+q00Nzx6NAhNn3gJIyIUxfIAtiJH1FREK4+DoI194HpAz13d85gHOO//3f/0Vubi5++9vf+h3EAIDD4cCxY8fwwQcf4KGHHsKFF14YwpEOHE8++SQYY4r/zZkzJ+D+lixZotrfkiVLgj7+wcpXEF5cXIyysrJu+ygsLFQNYtz9lJeXd9vP8ePHVYMYAJg1a1a3fQCAxWLBf/7zn27biaKIlStX+tVnICiQCTWjGUjIAGNnU2XMMRAuWtTHgwqeKKMGWq/ZGLtDeoeWjZsemifX6iBMvkRWeYhXlYAf+x5obwUvk96VECYGOYgwxbgKDJx3BYRZ10IYW+CahdH4CO7Gz4Dr3va5IzNBut6AN9RA/Gw50GELvDPOIa74M3hDreQwm7cYbHRezwdJeq0/pZW5yaqXpSTDnH++cuPoRPl7s/GMYlO+ZyvEjR9JjrHEVAi3/FJ5HY5WDzZuJjSTLgIzRMhOc6c84GMpQ6G57n6wi6/rPkAawJYuXYr77rsPbW1tfT0UQvpUd7OJW7Zs8Xk+XG0YY5g5c2a3fQDABx98gJaWlu4bAlixYoVf7QJBgUwo6YxA0lDZmgGWmQvkyCs6DTQM8tmY9g4RdW0OabvE9ODfdWQChIkXgZnjJId5QxXEQ1u7ft6/XfqwScEJqnISIvGTCWnQFFwFYeQUsNhk2d9ZDYtJAIbmBmUc/UF8pA6RBq/ZmJ1f9i4NrLUZzjeek1z8MUGAsORBQOEikYRef0src7OdLILl4CHJsfj5CtXLBA0QLZ0N55ZmwGZR7Vv86E3woiPSbiZOB7vU62ZUwhAIMxZASBsh64N32CAe2Axx5bPgx/cpPo8wLh/C/3sEbOIFrlmjc8iyZcvw2muv9fUwCOkXxo4di8RE9ayczZs3d9uHP0GKP/34ajNhwgTExsZ22wfgWhPlr+PHj+Pbb7/1u70/gv6J2fll5+OO9KCg0QHJw2TrNtyEC64a8OWAlWZjGts6YO0QwZu87qYHOZ2KjZsJFu9VPrWtEeK+rwGx6+JK9A5kktKBtN4FVclmA2YNT0CM0fcCVC6K4PVnIB7dBd7qVY1rwrmz0V5GvNdsTH0V+PE9ve+46AjEj5dLDrGoWAjjZ9P+Mn2gP6aVudXL0st+JE8vi0mWfB5zzlVnYzo5HXD+64/gLU2Sw8KCW4C4FECjAxtbAM3kS8AMJtnDec1piNs+Aa8qAVobIX6xEs6PXgWvr5K1ZYYICBcuhHDDL4H04b7HNUBUV1fj97//vc82RqMR+fn5WLRoEW6++WZcddVVmD17NlJS1DcrJWQgu+CCC1TPBSuQ6W2w429aWWlpKb755hu/2roFEvj4I6iBjFbDMCHTjElDoyFccA1gjg1m9wMHE1xBjFb9QpdFxYKdd2kYBxVcDPJKZe12J6wOV54JP7RD2n5kHuBd7rSnz50zRZZ/zm0WiHu+klcEO30SvL5a+vheppeNS1HfcZ47OsCrSiAe3AJx83sQd38JfvoI+GnpfhcYNhaIjOnVOPoD12yMV5GFnV8GbQEz/2oNxP1er6W4FAhTLwdM6n8HEnz9Ma3MrcErvUyXlISo6fldBzQ6ICpe+qC2Rv9SHxvrIC5/AdzjBgkTNBDGz3bNwqTnyB7CO+wQD2113VixexU6OX0c4qplELd+Aq5Q9IMlpkOz6G6wcRe4ZvUHsHfeeUc15cRoNOKll15CfX09tm/fjvfffx9vvfUWPvroI2zatAlnzpxBeXk5Vq9ejZtuugnx8fGK/RAy0PhKLzty5Ajq6tTXlZaVlflV9e3QoUOor69XPV9eXo7i4mLV8/4GMitWrIAoBjYz/95776G9PXgFoIIayOQkmzpTTJhOD2HkecHsfoBgQHKWbI8Cbm0D9yo/yyZfCMQmYSCKitBAI0jvijdYulLK+JFd0rQgnT4oi/5TYvQQhkn3FeGODoh7NwJW5fxr7/QyoRfpZSadBkPjpKlN3GqBWHYUzj0bIG5aDfHAZvAzxZKgip8pkly0MEEAG5ePgS7TazbGYnOCn1BOn+kpceVfwGuld85ZbDKE/PlgQ8fgXFtv1B/117QyN1txCSwHDkqOSTbHjE3pXKcIAJyLQKN8VkQNL9wL8dN3JMeYIUK2Pg8AeG05xO2fgFf6qBokOsH3bIL49vMQC39QbCKkDYcwYwFYcpbf4+xvPvvsM9Vzzz77LO69915ERKiniqanp+O6667D22+/jcrKSvz73/+mimVkwPP1Guac+5xN8Wemxd3P1q1bVc93N/PjbyDTk8X7zc3NWLNmTcCPUxO0/K+0GAPiIr1mIOLTXPnsg6lefmIGmFG6aJN3WIGaUxA3rYWw+Fedi02ZRgvhwoUQP369L0baY4zJ18ZY7E7YHB63Zy0tQPFBIGdS1+PGzwA/8F2PnzcuUofsRK8gQhQhHtgEtKjfeeD7tgFzFnSNIysXiE0AGgOrpgUAuUmR0vQaZwfEbR8Dzg7fD3Q6wI/uBvNIKWPjpoN/v0GSCjeQJJh1MHndpT9db0VOsMvJtrfC+a/noHngBTCPdEym0YLlng+enAXx8Leu1xwJif6cVuZWv/ZTmCZ03eSIu3Ie0FYJ6AyuypGeWuq7f8964etXQxw+BsI45Rt03GEHP/Y9eMUJ/zttawb/77/hPLQdwoULZRvJMr0RbOKFrlnewp1Ax8Aq237s2DHVc1dccUVAfen1eixevBiLFy/u1Ziam5uxa9cuVFRUoL6+Hq2trYiNjUV8fDzGjBmDiRMnQlCreucnURRRUlKC0tJSlJeXo66uDu3t7ejo6EBkZCSioqKQmJiI8ePHY8SIEX6vrxzMnE4nDhw4gKNHj6K+vh6NjY0wGAyIj49HRkYG8vPzERU1MGbp8/LyEBUVpTpbuXnzZlx99dWK55QCmYyMDDQ0NMiKaWzZsgULFiyQtXc/h5rhw4cjPb37zXu3bNmCEyeUP+8YY5g/fz4++eQTxfPLly/HTTfd1O1z+CMogUykQYOhifIpcMYYWOpw8FOHFB51DopNAfP6wuSODqCqxLWbeUM1+N7NYFMv7jzPho4CRkwATh4I71h7IdroezbGTTy4HRrPQCYxHUjJBKpOB/ycRp2AkSkm+WaLhduBugqfj+UnDoFbWsFMXQEmmzAdfMunPh4lxwDkJnkFqZXFfl8Q8YPfAZ6BjDnGlWJWdNDHo/ov77UxbTYn6tsCuzj0W+lxOP/2KDS3PuBa5+TBPTvDT+4BLy0E0E/ync4h/TmtzK1h7TpkPPpw58+6xERwZzMQlSj53OCiE2iqVurCt7PV9NgjL4HFe5Vwrq+EePg71VnhblUUQXz3RbAJBWDTr5BVPWMpwyDEpYIf3elabzNA+CoRW1ZWhjFjxqieD6bm5mb885//xMqVK3HgwAGfqTDR0dGYN28efvnLXyI/3/9Z8+XLl2Pbtm3YuXMnCgsLYbX6F3SazWbMnz8ft956Ky6//HLVdiUlJcjOVi/p7WnTpk0+g6Ovv/4aAHDRRRf51d+KFSt8VpsqLi7GsGHD/OorUOvXr8fLL7+MTZs2+ayMpdFoMGXKFNx111246aabYDAYQjKeYNBoNCgoKMD69esVzwc6IzNnzhycOXMGGzZs6FU/bv7Oxvha6zJr1iw89NBDqoHMxo0bcfr06aDsKdPr1DKBASNTTLKN8NyYQhWXc5I5Hsx7fwLRCdSUSC50+a7/grc2StoJs64CvDdo66cYA2K8ZmPabE7YHQpXNaePg3vtIcLG9Wx9Sk6KSRY8iSf3+nf3U3SCH5RukNeT6mWZsRGI1HtV5ypXv+MoU1sJXlkiOSSMH5iL/hPNOpj00ovbsvoQ3y0uPgLns/dALD0iO8U0Wgi550M4j9bOBJvA+ndamZut5BTa9nvdEIpNBTNFS48117puLPVEWwuc//cMeKtrU2Pu7IBYuAPi7v/2PIhx4yL4/q0Q33oO/MR+2WmmN0KYMBvChAuDtt4w1HxdTD/88MM+1wIEy8svv4zMzEz86le/wr59+7rN529ubsbq1asxffp0XHXVVaitrfXZ3u22227D66+/jr179/odxABAa2srVq1ahSuuuALz5s3D6dOB3+g7F+3fvx95eXm44oorsG7dum7L+zqdTuzatQs//elPkZOTg40bN4ZppD3jK71sz549iv/e+vp6HDoknxiYNWuWYn8//PADLBZ5Vca6ujocPnxY9fn9CWS62ztm8eLFKCgowNChysWVgrmnTK8DmeykCETolStzAQAzxwJRCb19mv4tIgqIl94l5pwDNaXyHdw77OBb10oOsag4sPMuCfUogyJGYTamUWE2xoWDH/ZaqJ2b50r1CEB6rAFRXpstihUnwIvlX/Zq+L5tXuOYCETI89t9GZ0snY2pbrX5TGlTHMdBr3FkjQKiB977Qz4b4wjdbIwnuw382C44v//CVTrXC62dCb5Yk67fp5W5NXwirV7mXUmMOx2uQKY3Sk/A+Ye74Ny9AeKWD8DLjnb/mEC0t0L8fIVrrZ1SMYCULAjTF4ClDAvu84ZAcnKy6rk9e/ZgxIgRuOuuu7Bu3bqgBzVWqxVXX3017r33XjQ3yz8r/PHJJ59g4sSJOHAgPBkTX3zxBebMmYPq6h7MGJ5D3nzzTeTn52Pfvp6ttywrK8Nll12GP/3pT0EeWfD4CmScTie++06ehr9161bXtaVCX0r9ORwObNu2TXZ8y5Ytiv24+RPI+No7RqvV4pprrgFjDNdff71qH8HaU6ZXgUyiWYfkaOlFaavVAWuH9EuOpZ0bpSQV6SOARPleMagvB6ytig/hx/fKF/5PmQPEqtcW7w8EBkQrzcY41d8Q/PBO18zUWUxnCGjRf4RekG+2aGkGL9yp8giVcRzZDd7RtfieabRgKrnuSqINWgyJkaZ7FFYr/319juP4XnCr9A5JsEtTh1pilE528+J0XZhz9xurIW5f1/3szACZ6fTEAMSZtK7Pln4gwSxd+9gf08rcvDfHlGmq7t3+Rm6tTUB9hbxKYhDxqhJX0YDqU7JznbMzE/v37Ex3qVlNTU147bXXMH/+fCQmJmLYsGFYtGgRnn32WWzevBk2Ww821IXrbu9NN92Ejz/+uEeP91RZWYl58+ahtLS01335o6ioCLfccktYnqs/WrNmDe64446AZrWUiKKIhx56CP/617+CNLLgmjZtGoxG9feuUuqX0rGkpCSMHj0a+fn50Ovl33f+9uOWnJyM3Nzu97nzlVY2d+7czr1yfK1pC9aeMj0OZAxaAdnJ0rtdDpHj2BkLalqkd2ZZavY5t8kXAGRnJJ7dK0b6b+ONVYDXviHexE0fSqt6abQQZi8MyTiDJTpC69faGAlLC1AsnQr198KdwVUJz/tusHjoW0Ds5nm92azghXul/U/0P71slNdsjNXhREl9D1JJnA7wI7uk4xg7zbVZ3wCRGSf98G21Orp/HYSC6PCYnZHfGWKxyUDaSMA8cMq2CgwYn2HG6HQzWFqObAPHvhiPdxGX/phW5mYvLUXbXuW7uNxhD3gGtc/ZrRD3b1KfnUl2zc7A1D9LuV977bUBtT916hTWrFmDxx57DBdeeCFiY2OxYMECrF69Gna7/0Hj888/321VpOzsbMydOxfXX389pk+fDpNJvg+QW3l5Oe644w6/nx8AEhISMH36dPz4xz/GDTfcgJtvvhlXX301zjvvvG7Xb6xfv17xTvq5rri4GDfddFO3a5hmzpyJ6667DpdddhmGDBnis8+7777bZ5nhvqLX6zFt2jTV80qL8ZUCEPeeNEajEeeff75fj+nt/jHd7R3jGbxMnjwZo0aNUm0bjD1lehRdMAAjU03Qel3UFldbYHOIqGmWfuAwvRFI9P1iG2gSYiPx8cs/76xA5sZbG/xbSNpQBb5P+mJiWaOB4eNVHtC3BOZKK/PUanOiw8dsjJt40HtTyiFAcka3jxsSb4TZK6WsotEGNKkvIPWF7/dK6xp3HqDtvt6FhjGMTJSmoZ2obYMf/3TlcXinl0VEguVM7FlnYZYUpYfRezYm1GtjutNYDXH7WuXZGUEASxgCmOP6YGCBy0kxSV/zMcnoyxS5WJNOcvOiP6eVudV7pZd1aqzCQC0EwatKzm6sqTw7w5KGAolDAaF/bUS9cOFCTJ48ucePt1qtWLt2LW644QaMGjUK//73v7t9TG1tLZ577jnV8+PGjcPGjRtRVFSE9evX491338W2bdtQXFyMO++8U/Vx//3vf/HFF1+onjcYDLj++uvx9ttvo6KiArW1tdi2bRvWrl2LVatW4a233sKHH36IXbt2ob6+Hi+99JLiHXS3d96RlvvOyMhAcXFx53/333+/6mPz8/Mlbb3/mz59OqZPny45tmjRItX+Fi1a5LO/jIzuv8/98dhjj6nOxOh0Orz44ouoqKjA1q1bsXr1anz55ZcoLS3FW2+9pboLvc1mw29+85ugjC/YfKWX7dy5UzIjabFYsHv3blk7z8BDqb/t27ejo6NrcqG1tRV79+5VfV5/Ahlfe8dERETIKq7dcMMNqn0FY08Zvz/1PGOWzHijbM1CdbMdta2uX1aHUwRvrHbdEXU/Pm0ExLpKAFy60FLQQq9R/qLmHOjwyGHQaRgY1L/WPb+iPNt4988ByQW4TmCdG4WbjPIPFou1KzAzGlz54qv/8jOMzPJa3N/eAtSVS0fh+dReaS5853/Bc6e4qledJcy6CmJFEeBwSB7q/fWr9DtwB5YOj9+ZwKBYiMH99/RMEWFQ3zA9JkILweuCpsniUL/EYqxrU8Tyk+DN9WDRXXfG2fgC8M0fyR929qEmgwZD4qR3rSw2J07XW5Hq+RzMx0wG93ydCeCF+8BFsXMGjRlNYGPPAy/cAzgcna9LBkgu3obHm2DQSp/naGdaGfM9m8LFrlQWxlwb87U2gZedAMvo2kiPTShwLfDtfG8w2Y0CTyLnkr9dqNsKZ1/KGfHSv0mL1dG5RqqzW6V0Li4CTo9ZG882eq+7k5wDHmmAsvOSi7Wznyfu2Zna0xBGzwDzXvAfn+7q0+Yxi+bjfQFI3xtqvzKj1vVasjq6PtT1GkGxveJ7zuMjIj3WgASz9HfHBA14hBlob4Hs88RDhF6DdnvX612vEaBR+Vz15ut9r5RWxrmfn8E6HZiPGwWizdZZelyj00GjU99AuMNq7dyQkmm1YD7aNnz2BTJ/+5h0XHYrYGmS/0O988S9X2uePD4jIAi+gwbu9Oj77GeEWpqj6Ozq1/0ZAcj7dzogHtwCVlUCNjpftlcZi4wBN0YClcdd/XFA8hfpgxK/giDg/fffx4wZM3q97qOkpAQ33XQTNmzYgH/+85+qJZJfe+011TUxWVlZ+Prrr5GUJN+/LTk5Ga+++irq6+tVFzGvWLFCtWx0RUWF35t2mkwm3HvvvWhqasITTzyh2GbTpk2Sn7VaraQymNqFO+C6O+9PFTHPNmazWbWd2WwOWVUyt7KyMrz77ruq59955x3FGT5BEHDzzTcjPj4eV155peJjP/zwQ7S0tPS78sy+AhmbzYadO3d2BhbeAYmbdyDzxz/+UXLeYrHghx9+wPTprgyU7777Dg6HegaFP4GMr0X68+fPl72WFi9ejN///veK7d17yvSmFLPfgUz+iFjVc06Ro7imK+//vOwYMK9vcpY8FJqLbwRvOAPxhy87jwsX/AR/UcnxPdXQjj99U9L58+OXDEdCpPKXgd0horyx6+InPVYP/dmLjL8sGC1pW9dmx+++PNn58y9mZyHr7CaHf1nwF0nbmoYWDL24K5r/+OWfY/Z5IxXHAEMkJF8cSUMlFXM0S6UvMOf//gp861qwK27uPMai46H52R8AAMM82pbUWjt7TjTrEGWUXzwPS3SVsfv3njLYzl5YTcuMwxgfO9HvLmnubDs0wYj0OP/yrTucXBZoeHImZQDVruorbOIFkiAGAIRx+cDZDSGjy1vR3O56YyXH6DE8SXmK32TQSCqmsdThEMbNVB/D/q4vAjapAJo75HdlNHe6vkScb70Ivt1VunBIjBGX5aovUm1st6PZdvaDIC4ZmqnqJTPF4z90lR9PyoDm+l8otmPpw8FmXQW+6Ww6RHwy/t9U9bKEByqb8X1ZIwDArNfg2knqM55HqlqwvdSV6mjQCrhxsvrds+O1rdha7Eq/0QrM5xgA6WxM52eE1+scAHjJYYhru3KVhTueBFMp+sCP7Yfzb11/K81Tb4JFKafO8KZaiLu6Nt0TxhSARci/kBkTwFOygfKjXVUEU3NkF4PujH5rh4g9p7ouhsYNMctmBwHg06X5aGzvwKJ/ft957I8LRiMvQ3m8TpFjZ1FT58+jUiPl+295M0W7ApnETLBI5X4PPTce4x5Z3xnMPH3teFwzTf3vvKuoqfOGx7DECKTG+leAgwMYplBq362soesOYtp99yD9179UbXv4ih/Dss9VsOPi+5di0Z+eVm37lzk/wrFNrs3dEm++EVl/VG97/OYlaNmxE1H5XWkbTG8Ehspnu3nNKcBdNMIUDe2L6ulInp8RbMxUaC6+UbWtWLijqwiA+zPiYpW2W9eC7/nG9YOPzwgAEIv2gRftA2+ogjBuJphXpgPTaIEMV0lj3lQDNJ7dSFajA8sY7d1dWAwfPhw7duzAtddei++//777B3TjzTffRFJSEp5//nnF87424Vy6dKliEOPpxhtvVA1k1q9fD6fTCY1G/v3rbxDjyX1xqaSwsDDg/gayzz//XHUB+pgxY7pNU5w3bx5iY2PR2NgoO2e327FhwwYsXNi/0vdnzJgBrVarGlhs3ry5M7BQShEzm82SGc+ZM2dCo9HA6ZTOmm/ZsqXzteYrrSw6OhqTJk1SPe/uS23vGEB5TcyoUaOQl5enOhPU2z1lgrJwpc3m7LeLP/s7fnwPeFkAG6j1A5xzWDv6b558qFW19GwBandYcu/rqYdTh1NEU1+sjekhxhiQnDXw1utFRHffJgw4536lkvYHpx76DSyHj3TO4pxzOqwQj/c+KAiXYcOGYceOHVi+fDkmTux9Gu2yZctw/Phx2fH29nbs2LFD4REujzzyiGt/Ox//+brYbWhoQFVVlc+xtbS04KOPPsL999+PefPmITc3F0lJSYiMjIQgCJLnuuyyy1T7cTgcPa62NhD5WnNx5MiRbv9ugiAoBjGeffQ33oGIN8+gQymQKSgokATVUVFRyMvLk7XzfKyvQMa7PyW+1rTExsZi3rx5iud8Lfp37ynTU4z7qsHmYcfJRuQq3D08VduOyiabZIZeYMB5pz4DG1cAwaNEJG+qcc3GeKWWPdiqXMVKKbXMrNfgZ+dnKgZOaqllL31XImunllq24o+vy/r1Ti377NV7MCOva38c3lTdtS6Ge43CM1XlU2nOa2fFm/gUCIt/BeaRosRPFeJU5DDFf5v3v89tyynX2hF/UstyY12zNN2lluk0DJOGRktSrVqsTtS2dvjM3M989/Gu34WgAQQNhCtuBsse1/VvqjoNcc3fsfOyrk3szEYNxg8xS6rAtVodOFzeCn52vNNKzm6w5EdqGV9zdqpaEACtDkhMg/axVyTNHH99BDh5GBCdeGvJPySpZflD4ySbYLbZHXh/f0Xn7+2Wpk3+p5YV7u9KGwHApl0GwWNzVG5rh/jGU2dfFwz/LlDPgQ53atkFmQmy4h5HylvR2N4VyLi7Pe9LhTulPlLL+DGv8rXdpJax+dd4Nvb6PNGg693BIMz4MViEdEaSW5pcpdEV3hc7O7pmMfxJLXvqS9fY/Ukt++3cUbJ+NQIwfkgUTAbfXx68qti1T4nKOMb/vTyg1LJ/3tFVTcr7fa/VMIxOi0SkQToDZbE7UdXs+33PAdRNcFW88Te17J9nGgNKLftZRqLP1DLRbgfO3pGcuv1T3ylVXl9/zr/+Tr2tV2oZW+hjh3ml1LKjKiV81VLLRk1Q6Nfj8wRwpZ/p9BCmXOoqMOOZOukjtYwpzE6F0+7du/HZZ59h48aN2LFjh+JeF9359a9/LSuvG8imkT21e/duxYvPmpoaPP3003j99dd7XXHLrbS0VHXDwCeffFI1XefCCy/0GRgoWbJkiWo53FtvvTWgRdnLly/HbbfdpnguKysLJSUlsuNz5syRpdMF03333Ye//e1vIeu/p379619j2bJliufMZjMaGxvBOUdsbCza2qQFhv7whz/g8ccflxx74IEH8OKLL0qOxcXFoa6uDna7HTExMarVAJ955hk8+uijqmO1WCxITU1VLbt8++23q1aJKy0txbBhw1Rn3Z5++mk89thjiue643dqWXK0XhbENLR1uBZfexE5XDnrFScBj0CGxSQBxsiuqXy42vkq3+vpR6OTcFluAkTOUdvagTab+t02zx67698zWPIMWhT75RxTxnpt8GNtk+dau0fheVitVGd9Ffi+LWCT53QeYlmjEdFsh0WlQpDSszkUojuRuy5OlY4r9ekdkGbGR8gW+7rXRPj8rXp2dPaLWjy4DRqPQIalZAJxXdP8AnNVKfMMYkTOcbLaorywnnOA+zkjIIqA3QZUlIBXngJLy+p63rFTIR7vusjgcP0udRqG7HjpxfvRmlav3x33v4Ia55LXAD/wHfiUOWBnZwiYIQIsd/LZvXe44t9TTSjbCgyylMPmdockiAE8XlP+lKT1bGPvZobL+7yv37f3Roc1peBpI8B0XeNnphjwmBSgSX5nVe1Xo3bcM4DpHK5T+T2r1MeIZJMsiKltscOoE6SpbKZoVzl3lXF4BjGdY/BzTb7n+16vYRibbpaV13aKHPVtfrzvPfvt6ABXyOlW4uzogNPPttzhAPeR4y1/QACzSN29Ft1EMYDKiWc/I/x5X3h+RvjTv+gAbA6g1o+7mYH8HkJsypQpmDJlCh5//HE4HA7s2bMH27Ztw+bNm7Fx40Y0NPiu+gnI15AArmAi1JRmSQ4cOIArr7wy6JtZeqcInctC/bfrr7Nbs2fPVg1kWltbsWfPHjidTlkQAyivZ5k9e7YskGloaMDBgwfR2Njos6R5d+tjfO0dA/iedRk6dCgKCgpUyy2vWLGix4GM3zkWWYnSPQ3sDhEnq7u5i1JfAW6TViPo6Z4yU4dE47Jc18aBAmNIMutUiwSE0tSxQ2HQdwV0nHPAFvjdJG98x5fgbdI3Wnykrk/qFek0DJnxRkwdFo2kaOmapFabM6ALYYlTheAtjZJDbFxXjvDQBPnmqqfrrKrBXE/xfV5V1FTKMOckREKn6XqLiJzjWE0vd+/21NIAlEjzoPvjnjK5SWYYdNKPij6vVOYvLgLVp1wbIXpgscl9XrJ2SJx8cX+bzYmT1RbUtXpd1IdhrAatgHEZ8iDG7hBR2WQfMGllZGDRarU4//zzcd999+H9999HdXU11qxZg5ycHJ+PUyqp62eCSa94P0dbWxsWLlwY9CBmsAn13y4cr42euOCCC+T7EHrYvHmzYlqZXq9X3Kdp1qxZiv1t2bLFZ1qZwWBQLN/sqbtZucsuu8xn+p+vPWN6s6eM34GMd3rSiSpL919snIOfKZIcYmkjVBqri9AJWDQxRdoPY0iK6maRbAgUTPH6cO2wBmeDtQ4b+Na1kkM6DUN8ZPjKaRp1AoYnRWDKsGhkxBuh1UhfHq7ZmF7cIeL87GxDFzZqSudGm2mx8opYSjN+vSXu9wpkUjOBFPnC6NHJ0pSkUw3taO8I7h0y8aB0916WkulXaepw0TCGiWnSNRpNlo7O4gwDgsMO1JbKv8gSMvps08lYkxaZ8dJZLodTxNHKNogcskCGabRni4mERoROwPgMM4w6aRBj7XDiYFkrBTEkbLRaLRYuXIgNGzYgOlp9fVhTU5PsWHcL+b/88ktwznv135w5cyR9vvbaazh58qTyEwK46KKL8PHHH6O8vBx2u13S19dff+37lzGI+Prb3Xjjjb3+uwVjv5JQiI+Px7hx41TPb9myRTGQOf/88xU31ExISMDYsWNlx9UCou76c+tu75hg6OnfqEerXssbrGjy80KGV3oFMsZIIC5VpbWyq8clI1qhYpBeK7h2wA6jmZO9AjFr72dj3Pix3eDl0g/E6AitqwJbCGefoowajEqLxOSsaKTEGBTX1ABAXZuj57MxZ/HDOyULcJneiORoA0Z4rb8QRY4TVcH73UqUHgdvrJUc8p6VSY0yIDZCGigfrVafUu2xU4XgLdI0Cja+IPjP00Ojks2I1EvfYwNmNsaTtc21E7sHJghAUhagCe9niFEnYGRKpOSuGeeuzYTdFQRtDhHc7lVb3xSaRf+RBg3GZZg7qzy6WexOHCpr7RwTIeGUlZXlszytUpCTkpKiWpYZQEg2mfzggw9Uz11yySX46quvsGDBAqSnp0PntbartrZW5ZGDT3p6uuq57du3q547F/h6nW/dulVxpsJXGpjSuc2bN+O7775TaN19f4DvvWOCpad7ygQcyLRaHThdF8CFTGsDuNduyoGkl41IiMDM7DjV8zERGhi04UnAEgSG6ZO8FhLagphqBED8Zg24V56/QSsgPVaPOJM2qKlm8ZE6jM8wY3xGFOJ9lICtb3WthWqxBmE2orUROCWtHjIs0QijV+rSqbr20FVG4xx8v3RmSPAKZEYnS0v4NrZ3oDIU1co4B/feMDR3MqBSkjxcTDoNzs+MxdQh0pSmRktHcF4HfaG1HrylTnKIaXWuYCZMe2wIzFVyWet1Y6K0TuHmUJvXHecQpJdFGTUYm26WpFACQJvNgUPlrX6vXyREyR//+EefF0/d8TUjk5iYKDtmMpl8psesXLlScS8Of+zdu1dxnYVS9TS3RYsW+Uwb+vLLL1XPdcdXdSlf6yDC1V+gLrzwQtVzRUVF2LhxY4/6bWpqwr59+1TPl5SU+EyHCvUsBOA7kKmtrUVdXZ3suK/AQ6m/yspKn+tbugtkfO0dEyzuPWUCFVAg4xQ5jldZAt4bmVdKZxlYsn93QbUCww15aZJjNocoSRFxp5iF4zJkXE46YqO89jgJciCD+jPgWz4G90pXY4wh1qTFkDi97KI/EAatgDHJZkweGoVRaZGyjU3dRM5R1WTDnlPNOHqmDTZH8C5oRO8Ld68P+qZ2B840+bEwthf4PumdOTZ8DBDtCpgjtAKyYuWL/EM2lsM7wD0WdTKdHmz0eSF7Pl+iDFoUDIvHNRPTMT41WpZeWDYQZ2M81VeAW6V/S2YwAfHq+/AEU06K8uJ+xRRKi/SiiWl1gEF5j6WeiInQYky6WRZUtbQ7cKi8DQ4KYkgvffXVV5g5cyYuvvhifPLJJwEFEfX19T4v9KdMUa52+qMf/Uj1MSdPngx4l/fdu3dj0aJFmDJlCurr62XnfZX8VVrH47Zr165epTv52tzxxIkTAQdsvvo7fPhwQH31hFrZXrelS5cqXtCrqaurw9NPP43s7Gx8+OGHvR1eSPkKZJQIgoCZM9X30PMVFKr1V1CgngnS3d4xwdST90RAV8RF1ZYe3SXnZ4ql6URaHVjSUB+PcLl0ZALSoqXrJtYerkaD194VOo0QlrUkBXnStLLislppSdkg4fu/hfj+K7ArpHToNALSYvRINGtVS8J602sYchIjMTc3CTfkDcH0rHgY9cp3XxxOEWX1VuwuaUZRTYhmRU4VgrfK85sBV7B8MlQpZR748QPg7dIglE10LZwbmWSG4PHLdThFnKgNXSADSwtQdFA6ljAv+o+L0OHC4Qn4yYQ0jEoySyrVudW12gfubIynmlLwDmmgzMxxQLT8Dm8wpftY3K/IYXPtSO8pSLMy8ZE6jE6PlP2dGy0dOFzRCidtDEaC6Ouvv8ZVV12F9PR03HPPPVi3bp1iUOC2Y8cOXHLJJT5Tr+bOnat4/M4770RkpPp6smXLlmHx4sWqF2aiKOLAgQNYtmwZJk2ahKlTp2LNmjWqi8WVZobcXnnlFXz11VeSY5xz/Oc//8Hll1/e49khAEhNVU/Rr62txe23344tW7bg5MmTKCkp6fxP7Xfqq7/9+/fjgQcewI4dO1BUVCTpz1cgF4isrCwsWrRI9fyxY8eQn5+PdevWqaY4nTlzBu+//z4WLVqE9PR0PPHEE35VwOtraWlpGDHC//XjEyZMQEyM+ndBeno6hg/3P/Np4sSJPvvzFVwYjUYUFhaiuLjY7/8eeugh1f56sqeM31f/p2rbUetdScdfdqsrPz2xaxEzSxshKwTgKdmsx+WjEqRjaGjHppMNmJoeC5NeI5mZiI7QwmIX0R7CjRpnTpG+0L7dcxLDJgbvDqnEmRKUN9oRa9IiNkIjm7WIMmoRodegrrVDsaqXVmDIjI1AdrwJGTERihemnqwdIiobrahutod+c1Mugh/eATZN/kVUUtsenpx8pwP80Pdg53XduWATp4NZgFFJ0rSyonpLyFNsxIPboBnZtaMuS0hFSnUVqlpDO6WfbNZjYloMMmPVF72LnKO2xY6imsBzV/sl0QnUlICnjpDs3YTYVMS1t8lulATDtKxYDPWxuF+VpUmaZmiKBhoqezWWq6amIzfVJPtMqW/rwLEzbf2pSi85x9TW1uKVV17BK6+8AsYYcnJykJ6ejsTERBiNxs4ysaWlpT77SUlJUd3pPSUlBQ8++CCefPJJ1ce/++67WL16NUaPHo1Ro0YhMjISTU1NqK2txaFDh3ym4HjLz8/HRx99pHjOYrHg0ksvxcSJE5GbmwtRFPH99993++/zx3nn+Z61f/vtt/H222/LjqvtCdNdfy+++KKsrC8A/O53v/P5uw7Es88+i7Vr18JuV87IOHnyJObPn4/ExERMmTIFiYmJsNlsqKurQ3FxMU6dOhWUcfSF2bNn+ywa4d3WnzZFRerX2J58pZVZLBb85z//UT1/xRVXYNSoUX49j9utt96KF154QfGcKIpYuXJlQKWY/Q5kels9ilcWgXkEMohPdaVJqJQuviEvVZK37RQ5Vu2p7Exrq2npwJA4vWRheqJZh/JGW8guxL0X+n+35yRumqiwaVkQNVocaLM5kWjWyVLKtAJDSrQebTYn6to6IDAgI8YVvGTGRsjy3pW0nq0MJiv3GmL88E7w8y/t3EMFcN0Nrm4ObUqZZAz7tgGegUxuHrJLm2D22giwsDqEszFuZSfAG2rAPPbVGZVsDlkgkx5txMS0aKRFq6/FcYgcx2tacfBMMzJ83OEckDpsQO1p8KSszgt6xhhGpkbiQFkL2oNY8js9xohH546UL+6vsnQftFuagNiuio1MqwfXRwDehQD8dNXUdCxbPEkWxNS22HGiB2nDhPQU5xzHjx/3ucZEzbPPPutz1uWxxx7Drl278Omnn/p8/iNHjvR6x/fbbrtNNZBx279/P/bv3y87Pnbs2B6nbeXk5GDSpEk+138EYtasWUhKSgrLXjxqcnNz8eabb+Lmm2/2WS65tra2V+uL+qPZs2fjzTff9Kttd+tZ3P35m6blq7/u9o655pprVM+pGTt2LEaPHo3CwkLF84HuKdPzxRYB4jWnJekcjDGwVOWpr+lDY5CbJP2Q+uZkPcqaui7qHB4btLlpNQwJ5tCUZM5KT0B6cqzk2Ld7/Iuee6vDyVHZZEdtawdEhSgt0qBBRqwBN+Rl4JKRSRjutf+JvD8RtS12HCpvxYGy1rAHMQCAlgbwo7slYwpHSpknfvh7cEfXv53pdJieFS9pU9NmQ50lHMEVB/cqxTwszgSj1vdblAHQMFepbrNegwSTDunRRmTHmzA62YxJ6dGYlhmLWdkJuGxkEq4ck4JrJ6bj8lHJqkGM3SniQGUz3t9Xju2lDWi1nwPpZEraW4DGM5JDGoFhVGqk32mb3REY8NjlI2Vr0UrrrGjyZ+anwwbe4RXM9jC9LMGsxx8WjZekTQJAVZOtR2sfCekLjzzyCG6//XafbbRaLd577z1ceeWVIR/PggULMH/+/IAfN3XqVPzpT3/q1XP/4Q9/6NXjPRkMBvz2t78NWn89deONN+LVV1+FwWDovvE5xJ/gxO2CCy7otk0g6258PbevYMhgMPTotQ/AZxphoHvKhK/uqOgEryoBy8jtPMTSh4Ofkq4NMOs1WDghWXKsrs2OT4/I7xK0WJ0w6QWYPNZ7mA0aWGxOtAV5E8WCydKgq6ahBcdK5DuDh1KL1QmL3TU7Y/Ja4yIIDAYfV19OkaOsqR3F9RacbmxHToz6wr5w4V9/ALG2EpXjLkN1sz38FZKs7eDH9oONndp5yOAVOBwNx2zMWbzwe/AZP3It6obrovqqcWkQOYfAXHs5McY6/7/A5IUSesPa4cThqhYcqW4ZPNWqmmvBdUbXGpmzIvQaZMYbcSqQ6owqrpuSjtEp0lRF1cX9aixNQIzHZ6IpWhaA+ePeuTkwewVUFQ3WoPw7CVFy6aWX4siRI6ioqOi+cTcSEhLw5z//GUuWLPGrvclkwrp16/D3v/8dv/nNb3q1s3t0dDSuvvpqJCQkKJ5ftWoVrr76amzYsMGv/goKCvDJJ5/gwIEDPR4TAMyfPx+vvvoqfvGLX8Bq7f37+J577kF1dTWeeeaZkJfa9eV//ud/MGPGDCxZsgS7d+/u/gEqGGO44IILcNFFF6m28bXGx2g0Ijs7W/V8MI0YMQJDhgxBeXm5z3Y5OTlIS0vz2SbQ/tTWR3W3d8zcuXN9Vhb0ZdGiRXjmmWdUzy9fvtxnQQNPYZuRART2lImMBaKlHww/mZAi27di9b4zqhdWta0dsoWpCWYd/MiqCsjMydKNMLft8S/3MNicIlDV7ErB6m5BrshdwcuWojqs2luGjSdqUVxv6fVeMEHjsIPv+QalddbQlVruhnf1Mk82h4ii+jDOElkt4CekaQImvQZmgxYmvRZGnQYGrQCdRoBGYEELYtrsDuwobcB/9ldgX2Xz4Ali3OrKwb1SXNNiDYgyqpcj9UdWfARuzc+UHGu3+1jcr8YiLYzBdAZAF1h57mGJJiyeIS2wUtNipyCGhNQjjzyC8vJy/PDDD3j++eexcOFCn/uFeNNoNJg2bRpeeuklHD9+3O8gxtPPf/5zlJeX4+9//zvy8/Nle7koYYxh7NixuPvuu/Hhhx+iqqoKK1asUK3sFRkZifXr1+Pll19GVlaWar8jR47En//8Z2zatEk1KArUnXfeiRMnTuCZZ57B3LlzkZGR4TPtrjtPPfUUCgsL8dhjj2HOnDlIS0tDRET4Nw6eMGECfvjhB2zcuBHXXXcd4uPju38QXMUXFi1ahJdffhmnTp3C5s2bfVbx8rUp6VNPPeXz7xls/szKBDJz09v+Vq5c6TOg7UlamdvkyZN9FiR47733YLH4910Z3p3gmqrBLS1gpq4PA5Y2Aih3XTiNTorEtKHStIkfyppxuEq9xLFTdO2CnRzdVQ1IIzAkmnWoag5eypT3jMx3e8OTVqamzS6ivcGG+Eid5IKLc44zLTYU11tQ0uBHDv4gxw/sABbfo3juRG34qzfxA9uAMJVebrZ2YH9lM07WdbPg/JzHgdoy8PSczjVbjDGMSDZh/+mWHv1uBAY8dGkO9B53VDh3bfIacH92K7jDDqb1qHhmigGa/A9CfvWjUZJ0U1HkKK07R4o3kH5vypQpknLJDQ0NnYuzq6urYbFYYLFYoNFoEB0djZiYGIwcORITJkwIykW02WzG0qVLsXTpUlitVuzevRvFxcVobGxEU1MTNBoNoqKikJCQgJEjR3YWAAiEIAi4++67sXTpUuzbtw+7d+9GbW0ttFot0tPTkZubi6lTp0oeM2fOHJ9rQfw1ZMgQPProo3j00Ud73RfgCriefvrpHj9+yZIlPQo6lVx00UWdMyrHjh3DgQMHUF9fj8bGRtjtdpjNZsTExCArKwujR4/2a7bCk9os2rRp0/DAAw/0evyBWLVqFVatWtVv+nv88cfx+OOPB2083vwtbtCd8AYycO0pw0bkdf7MUrKhYUUQGMP1edLpLYvdiff3d59C0WYX0WpzwuyxP4NJr0GUQUSLrff5/QmxkRgzXPrm+HZ33wYyACBy14xUi9WBSIMGR6qbcaqhHZaOc3RNQyg01YMXF4Jlj5adKgzh3jGqzpRAPLILwhj1Td18EUUOm1OEzSHC6nDC5hC9/nPC5hTRZneirs1O6yLcHDagsQqI63qfR+g1GJoQgZLawC/4lVLKKhttaO3p51FbExDTVQgCkdFAk3+prZOGxuBKr/24KptssAdxbyhCAhEXF4e4uDjVvWBCyWg0oqCgwOe+Gb0hCAImT56MyZMnh6T/wSw3Nxe5ubndN/STw+HA5s2bZcf1ej3eeOMNn5uEkv4j/IHMmSLAM5DRGzAu1YysuAgkee2x8PGhar8DkbrWDhi1gmRzt3izNijlmGfkSWdj2tpt2Hs0sDrXoWRzcNgcDhwJ43qOc4m4fzs0XoFMRbMVzdbgl+H1B9+wGs7932Lz1JsgcleKoPt/eefPXHbO7hDRMbinVnqnuRYtxiRERXR9LKbFGgLePycrTp5SVtrQjsr6XlSgs0gDGaYzgusMrupr3XhkvvS13Wixo7wh9Dt1E0JIf7Z9+3a0tsqvm5544gmMGzeuD0ZEeiKsa2QAAO2t4A3SO4lXjErEpSOl+aIn6yz4rqTR727dsxOeBMaQFKVDb1cSeK+P2bm/BA5K2Tpn8H3bZccKq/3fRyD4OFB9GiUN7ShtbEdZkxUVzVacabGhqtWGmjY76iwdaGjvQJPVgRabA212JwUxQXCi2iKrDJiTYvK7ipnAgAcvHSFJKXOKHH/acKJ36Xv2dkmFPQB+VS+7eGwy8kdIP1v/vuEkbXhJCBn0vDcrBYC8vDw88sgjfTAa0lPhD2TgSi/zNDROumGjw2vPGH+1d4hobpfeRTfqBFwy0r9FYmoK8vrX+hgSZFWnIR7+ofPHOosdpY20fmAwsnaIKPVaAG/UuVLM/HHd5HSMSZUuCH5/byUOnwnCbKnXon+YfFeL0QgMD/1YulFZeX07Vm4duJvGEUJIsHivj9FqtXjjjTeg1YY9WYn0Qt8EMlWnwJ3qaTsbjtXhTEvP9u6ob3OgwymdLblyTBJSo/Qqj/DNZNRj8hhptZ9w7R9Dwkdc/mc4172NveVN2HCshnY3H8Qqm2yyGyJpsQZER/j+csuKi8Ct06UpZacb2rF8e5DSUL2rl+kjAK3659qi84cg1yuoWvb5UdhpNpkQMsi1tbVhx44dkmMPP/wwrW0agPokkIGzA7xG+cu9utWO9Udre9w1B1DT0iGpBKLTCJibm9ij/s4bnwWdrmvBl8PhxM79xT0eH+mn2prBP1+FPRVNVCyB4GS1RZZ+NSI5QjXFTCmlTOQcL2w4AbszSIGDzQLu9C+9zKgT8IsrpItiD5c34+Pdvd/PgxBCBrrIyEjY7XZwzjv/602lNtJ3+iaQgTy9zO3dvZW9zvW3OTia2qUXo5OHRMGsD7wCxcwpIyQ/7ztahrb2cOz0TgjpK64UM2l6oVGnQVaicorZtUopZXuClFLmyeK1qZ9Ketlts7ORGiPda+a5tYU000gIIeSc0meBDOorZZvQ7ShtxLGa4GxA2NTugOg1KzM9KzbgfgrypIHMd5RWRsigcKbJLksxS42Rp5gNjYvAknx5StmbwUop8+SdXmYwARrpBn9xkTrcdbF0Xd/Wo7XYeqznM92EEEJIf9R3gQzn4Ce7djGvbrVjzYHqoHUvcqDNJk3puCA7NqAKZhqNgPyJ2ZJjtD6GkMHjRJU8xSwn2QR3lXfXxpcjoNeGMKXMk7VNvr7QK73snstyEBUhDW6eX1cY/LEQQgghfaxPSzPwiuNwttbjLZaHI9VtaLMHd21Ci9Uh2fU+MVKPMSmROFzV5tfjJ43KQFSkND1j296ioI6RENJ/2RwiTtW1Y3iSqfOYQScgKzECRTXtuCZPnlL2QbCqlKmxNANRHpUYPdLLMuMjcFNBlqT5xz+U41C5V0oaIYQQcg7ouxkZt+Y6fF/WHPQgBnBvFCm9KzorO87vx3uXXT5WUoXq+r7cX4QQEm5VTXY0WaSL7FNiDEiN0eM2ryplZaFKKfPknV5mjOxcD/PrH42SzA7ZHE4s+/xYaMdDCCGE9JG+D2RCrMVr0f+4VDPiTTqV1lIzp0g3wqS0MkIGp5PV7bIUs+wkkzyl7KuTspsnQWdtAxeln2tzJ6RgQkYM5k9Jlxx/e2spyuppTyRCCCHnpnM+kGm1OdHuUU5XYAwzh8X69dgZ3hthUiBDyKBkc4g4Ves7IFiztxKHKsMxY8tl1cvmTUrFw/Olm182t3fglQ0nwjAeQgghpG+c84EMB7CzVJqKMSMrFlq1DSHOyhmajJQEaWlTmpEhZPCqaraj0SvFzK2soR1vhDqlzJNXeln+iAQUjJTulfWPr06qjpcQQgg5F5zzgQwAbClukPwcbdRiUnqUSmuXgsnS2ZjKmiYUl1H5UkIGM6WNMkXO8adwpJR5am+VpZd5qmhox/ItJeEbDyGEENIHBkUgc6bFjuM10kpl3S36nzlZuj6G0soIIXYHR4lXitmavZU4GJaUMk8caFd/zr9+cRy2jjAGVoQQQkgfGBSBDCCflclJNCEtyqDa3ntG5ru9FMgQQoDqZjuOV7WhrtWO//v2FF7deqpvBuKVXuZ2tLIFa74vC/NgCCGEkPAbNIHMvooWNFulG8nNGh6r2DYlIQo5Q5Mlx77dTYEMIcSltqUDx85Y8O7uCvDum4dGewu4KJ91eX5dIcQ+GxQhhBASPoMmkHFy4LuSRsmx8zNjoNfIF/0XTB4h+bm5tR0HjpeHcniEEBIYLk8v23aiDt8cqemjARFCCCHhNWgCGQD4tqQBIu+6VRmh0+D8zBhZu5legcz2fcUQ6RYnIaS/aaoG565ZmYY2O5784FAfD4gQQggJH21fDyCcGtodOHimFRPTuiqWzRoeh2+9Zmpm5EkDGVofQwjplzqsQPlR3PNVO3aXNOJMk7WvR0QIIYSEzaCakQGALUXSRf8ZMUZkx0d0/hwVacSkURmSNrR/DCGk33I68Nm+MxTEEEIIGXQGXSBTWN2Gmla75JhnKeZpE4ZBo+n6tdg7HPj+YB9VJSKEEEIIIYQoGnSBDAewtUQ6KzN5SBTMeg0AYOYUaVrZ7sOlsNpod2xCCCGEEEL6k0EXyADA9lNN6HB2lS3VaQRMz3It+i/wXh9DaWWEEEIIIYT0O4MykGmzO7G7vFly7ILsOGgFhmkThkmO0/oYQgghhBBC+p9BGcgA8kX/iZF6zM1NQIRRLzm+bW9ROIdFCCGEEEII8cOgDWRKGqw43Sit8nP5qETJz4dOVKCh2RLOYRFCCCGEEEL8MGgDGQDYWiydldEITPIzpZURQgghhBDSPw3qQGbX6Sa0dzhVz2+jQIYQQgghhJB+aVAHMnYnx87SJtXzNCNDCCGEEEJI/zSoAxkA2OKVXuZ2urIep88onyOEEEIIIYT0rUEfyJxpseN4TZvsOM3GEEIIIYQQ0n8N+kAGUJ6V+W4vBTKEEEIIIYT0V4xzzvt6EIQQQgghhBASCJqRIYQQQgghhAw4FMgQQgghhBBCBhwKZAghhBBCCCEDDgUyhBBCCCGEkAGHAhlCCCGEEELIgEOBDCGEEEIIIWTAoUCGEEIIIYQQMuBQIEMIIYQQQggZcCiQIYQQQgghhAw4FMgQQgghhBBCBpz/D7Pz+G9pzz47AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 500x100 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define the shape of each bar\n",
    "rect_ll_y = df_annual.anomaly.min() # rectangle lower left y coordinate, minimum anomaly value\n",
    "rect_height = np.abs(df_annual.anomaly.max()-df_annual.anomaly.min()) # rectangle height, range between min and max anomaly values\n",
    "year_start = int(df_annual.year.min()) # year to start the plot x axis\n",
    "year_end = int(df_annual.year.max() + 1) # year to end the plot x axis\n",
    "\n",
    "# create a collection with a rectangle for each year\n",
    "col = PatchCollection([\n",
    "    Rectangle((x, rect_ll_y), 1, rect_height)\n",
    "    for x in range(year_start, year_end)\n",
    "])\n",
    "\n",
    "# Create the figure, assign the data, colormap and color limits and add it to the figure axes\n",
    "fig = plt.figure(figsize=(5, 1))\n",
    "\n",
    "# set up the axes\n",
    "ax = fig.add_axes([0, 0, 1, 1])\n",
    "ax.set_axis_off()\n",
    "\n",
    "# set data, colormap and color limits\n",
    "col.set_array(df_annual.anomaly) # use the anomaly data for the colormap\n",
    "col.set_cmap(cmap) # apply our custom red/blue colormap colors\n",
    "col.set_clim(-rect_height/2, rect_height/2) # set the limits of our colormap\n",
    "ax.add_collection(col)\n",
    "\n",
    "# plot anomaly graph\n",
    "df_annual.plot(x='year', y='anomaly', linestyle='-',lw=3,color='w',ax=ax, legend=False)\n",
    "# plot horizontal line at zero anomaly\n",
    "ax.axhline(0, linestyle='--', color='w')\n",
    "# plot a text label\n",
    "ax.text(df.year.max()+3,-.4,'Seattle, WA', fontsize=30, fontweight='bold', color='k')\n",
    "\n",
    "# Make sure the axes limits are correct and save the figure.\n",
    "ax.set_ylim(-rect_height/2, rect_height/2) # set y axis limits to rectanlge height centered at zero\n",
    "ax.set_xlim(year_start, year_end); # set x axes limits to start and end year\n",
    "\n",
    "# save the figure\n",
    "fig.savefig('warming-stripes-seattle.svg', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ulmoenv",
   "language": "python",
   "name": "ulmoenv"
  },
  "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.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}