analysis-brouwershavensegat.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "url = 'http://opendap.deltares.nl/thredds/dodsC/opendap/rijkswaterstaat/waterbase/27_Waterhoogte_in_cm_t.o.v._normaal_amsterdams_peil_in_oppervlaktewater/nc/id1-BROUWHVSGT08.nc'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "import xarray\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import statsmodels.regression.linear_model \n",
    "import statsmodels.api as sm\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# open dataset\n",
    "ds = xarray.open_dataset(url)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get first location (only 1 location in file)\n",
    "ssh = ds.sea_surface_height[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert to dataframe\n",
    "df = pd.DataFrame(ssh.to_pandas(), columns=['ssh'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "             ssh\n",
      "time            \n",
      "1978-12-31     2\n",
      "1979-12-31  8784\n",
      "1980-12-31  8760\n",
      "              ssh\n",
      "time             \n",
      "2015-12-31  52704\n",
      "2016-12-31  52560\n",
      "2017-12-31  35943\n"
     ]
    }
   ],
   "source": [
    "# check number of measurements per year\n",
    "# this data is up to date until 2018-09\n",
    "# after that waterbase was is no longer available.\n",
    "# the updates based on DDL are not yet available\n",
    "counts = df.resample('A', label='left').count()\n",
    "print(counts.head(n=3))\n",
    "print(counts.tail(n=3))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# drop incomplete years\n",
    "df = df[df.index > datetime.datetime(1980, 1, 1)]\n",
    "df = df[df.index < datetime.datetime(2018, 1, 1)]\n"
   ]
  },
  {
   "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>ssh</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>time</th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2017-12-31 23:19:59.999998976</th>\n",
       "      <td>1.80</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017-12-31 23:30:00.000003328</th>\n",
       "      <td>1.80</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017-12-31 23:39:59.999997696</th>\n",
       "      <td>1.87</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017-12-31 23:50:00.000002304</th>\n",
       "      <td>1.89</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2017-12-31 23:59:59.999996672</th>\n",
       "      <td>1.88</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                ssh\n",
       "time                               \n",
       "2017-12-31 23:19:59.999998976  1.80\n",
       "2017-12-31 23:30:00.000003328  1.80\n",
       "2017-12-31 23:39:59.999997696  1.87\n",
       "2017-12-31 23:50:00.000002304  1.89\n",
       "2017-12-31 23:59:59.999996672  1.88"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.tail()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Int64Index([2015], dtype='int64', name='time')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# double check date calculation (different timezones)\n",
    "# calculation is done while in UK, double check if back in NL\n",
    "test_df = df[\n",
    "    np.logical_and(\n",
    "        df.index >= datetime.datetime(2015, 1, 1),\n",
    "        df.index < datetime.datetime(2016, 1, 1)\n",
    "    )\n",
    "]\n",
    "# this should only show 2015\n",
    "test_df.resample('A').mean().index.year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define some extra variables\n",
    "annual_mean_df = df.resample('A', label='right').mean()\n",
    "annual_max_df = df.resample('A', label='right').max()\n",
    "year = [x.year for x in annual_mean_df.index]\n",
    "annual_mean_df['year'] = year\n",
    "year = [x.year for x in annual_max_df.index]\n",
    "annual_max_df['year'] = year\n",
    "# m to mm, sealevelmonitor is in mm\n",
    "annual_mean_df['height'] = annual_mean_df['ssh'] * 1000\n",
    "annual_max_df['height'] = annual_max_df['ssh'] * 1000\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgEAAAE/CAYAAADffhAhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmYZHdd7/H3hwkNEpYkJIQYMkxAQAEhhDbQstgySACRTYkssj+GPICiIF4i4kVzc0e5IIriEmQJyi5LcgWEMJcGhQ4wAwlkAbIQIDGQQAiJLOnM5Hv/qNNQM8x0V0/qdC3n/XqeevrU75w69f11VVd9+7edVBWSJKl7bjLqACRJ0miYBEiS1FEmAZIkdZRJgCRJHWUSIElSR5kESJLUUSYBkoYuyXySS0cdh6SVmQRIUyDJQpLvJLnZqGORNDlMAqQJl2QT8CCggEePNBhJE8UkQJp8TwPOBN4EPL1/R5I3JXltkvcnuTbJp5LcuW9/JTkhyQVJrm6OTbPv5Un+pe/YTc3x+zX3n5nk/Oa8Fyd5zqABN+d5bvO81yY5Kcmdk3wyyTVJ3plkpu/4RyU5q4nxk0nu1bfvJUkuas5zXpLH9e17RpL/TPLKpqXkK0kesZZfrjTNTAKkyfc04C3N7dgkh+62/4nAnwIHAhcCJ++2/1HALwD3Ao4Djh3wea9oHntr4JnAq5McvYa4jwXuC9wf+EPgFOC3gCOAewJPAkhyH+ANwHOA2wL/CJze1/VxEb2WkNs09fyXJIf1Pc/9gC8BBwOvAF6/nOhIXWcSIE2wJA8E7gi8s6q20/tCfPJuh723qj5dVTvoJQpH7bb/z6vq6qr6GvDRPezfo6p6f1VdVD0fAz5M78t4UK+oqmuq6lzgHODDVXVxVX0X+CBwn+a444F/rKpPVdXOqjoVuI5e8kBVvauq/quqbqiqdwAXAMf0Pc9Xq+p1VbUTOBU4DNg9UZI6ySRAmmxPp/fl+a3m/lvZrUsA+Ebf9veBW65x/x4leUSSM5NcleRq4JH0/tse1Df7tn+wh/vLcdwReFHTFXB181xHAD/dxPG0vq6Cq+m1IvTH8aP6VdX3m82B6ihNu/1GHYCkfZPkp+g1329IsvxFdzPggCT3rqqzb+RTfA+4Rd/92/c9982Ad9Prijitqq5P8j6gjWb2rwMnV9Xu3RgkuSPwOmAzsFhVO5Oc1VIc0tSxJUCaXI8FdgJ3p9eEfxTwc8B/0PtyvrHOAh6cZGOS2wAn9u2boZdwXAnsaAbbPWwIz7knrwNOSHK/9Oyf5FeT3ArYn96siCuhN1iRXkuApAGYBEiT6+nAG6vqa1X1jeUb8LfAU5ZH8e+rqjoDeAfweWA78G99+64Ffhd4J/AdeuMQTr8xz7dCHNuA36ZXr+/QG9z4jGbfecCrgEV63Qk/D3yijTikaZSqGnUMkiRpBGwJkCSpo0wCJEnqKJMASZI6yiRAkqSOMgmQJKmjpn6xoIMPPrg2bdo06jAkSVo327dv/1ZVHbLacVOfBGzatIlt27aNOgxJktZNkq8OcpzdAZIkdZRJgCRJHWUSIElSR5kESJLUUSYBkiR1lEmAJEkdZRIgSVJHmQRIktRRJgGSJI3I4uIiW7ZsYXFxcSTPP/UrBkqSNI4WFxfZvHkzS0tLzMzMsHXrVubm5tY1BlsCJEkagYWFBZaWlti5cydLS0ssLCysewwmAZIkjcD8/DwzMzNs2LCBmZkZ5ufn1z0GuwMkSRqBubk5tm7dysLCAvPz8+veFQAmAZIkjczc3NxIvvyX2R0gSVJHmQRIktRRJgGSJHWUSYAkSR1lEiBJUkeN5eyAJHcD3tFXdCfgT4ADgN8GrmzK/6iqPrDO4UmSNBXGMgmoqi8BRwEk2QBcBrwXeCbw6qp65QjDkyRpKkxCd8Bm4KKq+uqoA5EkaZpMQhLwROBtffefn+TzSd6Q5MBRBSVJ0qQb6yQgyQzwaOBdTdHfA3em11VwOfCqvTzu+CTbkmy78sor93SIJEmdN9ZJAPAI4LNV9U2AqvpmVe2sqhuA1wHH7OlBVXVKVc1W1ewhhxyyjuFKkjQ5xj0JeBJ9XQFJDuvb9zjgnHWPSJKkKTGWswMAkuwP/ArwnL7iVyQ5Cijgkt32SZKkNRjbJKCqvgfcdreyp44oHEmSps64dwdIkqSWmARIktRRJgGSJHWUSYAkSR1lEiBJUkeZBEiS1FEmAZIkdZRJgCRJHWUSIElSR5kESJLUUSYBkiR1lEmAJEkdZRIgSVJHmQRIktRRJgGSJA3R4uIiW7ZsYXFxcdShrGq/UQcgSdK0WFxcZPPmzSwtLTEzM8PWrVuZm5sbdVh7ZUuAJElDsrCwwNLSEjt37mRpaYmFhYVRh7QikwBJkoZkfn6emZkZNmzYwMzMDPPz86MOaUV2B0iSNCRzc3Ns3bqVhYUF5ufnx7orAEwCJEkaqrm5ubH/8l9md4AkTYlJGpWu8WBLgCRNgUkbla7xYEuAJE2BSRuVrvFgEiBJU2DSRqVrPNgdIElTYNJGpWs8jHUSkOQS4FpgJ7CjqmaTHAS8A9gEXAIcV1XfGVWMkjQuJmlUusbDJHQH/HJVHVVVs839lwBbq+ouwNbmviRJWqNJSAJ29xjg1Gb7VOCxI4xFkqSJNe5JQAEfTrI9yfFN2aFVdXmz/Q3g0NGEJknSZBvrMQHAA6vqsiS3A85I8sX+nVVVSWr3BzUJw/EAGzduXJ9IJUmaMGPdElBVlzU/rwDeCxwDfDPJYQDNzyv28LhTqmq2qmYPOeSQ9QxZkqSJMbZJQJL9k9xqeRt4GHAOcDrw9OawpwOnjSZCSVKXTOOyzOPcHXAo8N4k0IvzrVX170k+A7wzybOBrwLHjTBGSVIHTOuyzGObBFTVxcC991D+bWDz+kckSeqqPS3LPA1JwNh2B0iSNC6mdVnmsW0JkCRpXEzrsswmAZIkDWAal2W2O0CSpI4yCZAkqaNMAiRJ6iiTAEmSOsokQJKkjjIJkCR11jQuBbwWThGUJHXStC4FvBYrtgQkOWiA2wHrFawkScOyp6WAu2a1loD/am5Z4ZgNwMahRSRJat3i4uLUrX63VstLAS+3BEzLUsBrsVoScH5V3WelA5J8bojxqMP8UJLWh83gPdO6FPBarJYEDPIb6d5vTUPnh5K0fqb1inj7YhqXAl6LFccEVNUPl7eTHJjkXkmOXr7tfoy0r+yb64a2RmJ3fYT3Wk3rFfG0dgPNDkhyEvAM4CKgmuICHtJOWOoa++amX1utPbYirZ3N4Fo26BTB44A7V9VSm8Gou/xQmn5tNUHbtL1vut4Mrp5Bk4BzgAOAK1qMRR3nh9J0a6u1Z1xakRzYqkmUqlr9oGQWOI1eMnDdcnlVPbq90IZjdna2tm3bNuowJNHeF+Wov4DtktC4SbK9qmZXO27QloBTgb8AvgDccGMCk9RdbbX2jLoVyS4JTapBk4DvV9VrWo1kCo36vxNJ62NcuiSktRo0CfiPJFuA09m1O+CzrUQ1BWwelLrDga2aVIMmAcurBt6/r8wpgiuweVDqllF3SUj7YqAkoKp+ue1Apo3Ng+PDbpl2+fuVJtegiwUdADwN2NT/mKr63XbCmnw2D44Hu2Xa5e9XmmwrLhvc5wP0EoAvANv7bq1IckSSjyY5L8m5SV7QlL88yWVJzmpuj2wrhmGYm5vjxBNP9ENxhFyOuF3+fqXJNuiYgJtX1QtbjWRXO4AXVdVnk9wK2J7kjGbfq6vqlesYiyaY3TLt8vcrTbZBk4B/TvLbwL+x6+yAq9oIqqouBy5vtq9Ncj5weBvPpelmt0y7/P1Kk23QFQOfB5wMXE3fBYSq6k4txrb83JuAjwP3BF5I70JG1wDb6LUWfGelx7tioCSpawZdMXDQMQEvAn6mqjZV1ZHNbT0SgFsC7wZ+r6quAf4euDNwFL2Wglft5XHHJ9mWZNuVV17ZdpiSJE2kQZOAC4HvtxnI7pLclF4C8Jaqeg9AVX2zqnZW1Q3A64Bj9vTYqjqlqmaravaQQw5Zv6AlSZogg44J+B5wVpKPsuuYgFamCCYJ8Hrg/Kr6y77yw5rxAgCPo3dBI0mStA8GTQLe19zWywOApwJfSHJWU/ZHwJOSHEVvXMIlwHPWMSZJGgoXWNK4GHTFwFPbDmS35/tPIHvY9YH1jEOShs0FljROVhwTkOSU1U4wyDGSpB4XWNI4Wa0l4LFJfrjC/gBeV0CSBuQCSxonqyUBLx7gHP8xjEAkqQtcYEnjZKDFgiaZiwVJkrpm2IsFSZKkKWMSIEljbHFxkS1btrC4uDjqUDSFBpoimOQJVfWu1cokScPjdEK1bdCWgBMHLJMkDYnTCdW2FVsCkjwCeCRweJLX9O26NbCjzcAkqQ2TtFqf0wnVttW6A/6L3iV7Hw1s7yu/Fvj9toKSVjNJH+QaH5PWvO50QrVtxSSgqs4Gzk7y1qq6fp1iklY0aR/kGh97al4f9/fO3Nzc2MeoyTXomIBjkpyR5MtJLk7ylSQXtxqZtBf2k2pfLTevb9iwweZ1icGvIvh6es3/24Gd7YUjrc5+Uu0rm9elXQ20YmCST1XV/dYhnqFzxcDp5JgASdq7QVcMXDEJSHJ0s3kcsAF4D3Dd8v6q+uyNjLN1JgGSpK4ZNAlYrTvgVbvd7z9hAQ9Za2CSJGk8rDY7wMsES5I0pQZdNviFeyj+LrC9qs4abkiSJGk9DDpFcBY4ATi8uT0HeDjwuiR/2FJsY8mLeUiSpsWgUwTvABxdVf8NkOR/Au8HHkxv2uAr2glvvLhIjSRpmgzaEnA7+mYFANcDh1bVD3Yrn2ouUiOtH1vdpPYN2hLwFuBTSU5r7v8a8NYk+wPntRLZGGpzkRrnvUs/ZqubtD4GSgKq6qQkHwQe0BSdUFXLk++f0kpkY6it1cb8wJN2NYlr/EuTaLVLCd+6qq5JchBwcXNb3ndQVV3VdoDjpo2LefiBJ+3KpaGl9bFaS8BbgUfRG/xXQHb7eadWo+sIP/CkXbnGv7Q+Brp2wLhJ8nDgr+ktZfxPVfXnezt2UpYNdkyANH78u9SkGtaywcsnC72+/yOb8QEbgdtX1advZJxrlmQD8FrgV4BLgc8kOb2qJnqAotcM177wS6o9jtVRFww6O+DvgBvoXSvgJOBa4N3AL7QU10qOAS6sqosBkrwdeAwdmqUggV9SbXOsjrpg0HUC7ldVzwN+CFBV3wFmWotqZYcDX++7f2lT9iNJjk+yLcm2K6+8cl2Dk9aL61a0a3mszoYNGxyro6k1aEvA9U0zfAEkOYRey8BYqqpTgFOgNyZgxOFIrXBAabscnKguGDQJeA3wXuB2SU4GfgP449aiWtllwBF99+/QlEmd4pdU+xyro2k38OyAJD8LbKY3PXBrVZ3fZmArxLEf8OUmlsuAzwBPrqpz93T8pMwOWAsHg0mSVjLU2QGNC4Brlh+TZGNVfW0f49tnVbUjyfOBD9GbIviGvSUA08jBYJKkYRl0iuDvAP8T+Cawkx8vFnSv9kLbu6r6APCBUTz3qDliWZI0LIO2BLwAuFtVfbvNYLQ6B4NJmhZ2bY7eoEnA14HvthmIBuNgMEnTwK7N8bDaBYRe2GxeDCwkeT9w3fL+qvrLFmPTXjhiWdKks2tzPKzWEnCr5ufXmtsMo1skSJI0Jdrs2rSbYXArJgFV9aeDnCTJ31TV7wwnJEnStGura9NuhrVZyxTBlTxgSOeROs//YjRu2npPttG1aTfD2gwrCZA0BP4Xo3Ezae9JZ1CtzaAXEJK0DrwokMbNpL0nl7sZTjrppLFPWMbBsFoCMqTzSJ3mfzEaN2t9T45Dd5YzqAa3piQgyS2q6vt72PXXQ4pH6jTXgdC4Wct7ctK6DjT4ssG/CPwTcEtgY5J7A8+pqucCVNWbWotQ6hj/i9G4GfQ96aC8yTPomIBXA8cC3waoqrOBB7cVlCRp8ix3HWzYsMHurAkxcHdAVX092aXrf+fww5EkTao2u7PGYazBNBr42gFNl0AluSm9Cwqd315Y0vjzQ0n6SW10ZznWoD2DJgEn0Bv8dzhwGfBh4HltBSWNOz+UpPXjWIP2DDQmoKq+VVVPqapDq+p2VfVbXlZYXTZpc6elSeZYg/YMlAQkOTXJAX33D0zyhvbCksabH0rS+nEBoPakqlY/KPlcVd1ntbJxNDs7W9u2bRt1GJpCjgmQNK6SbK+q2dWOG3RMwE2SHFhV32lOftAaHitNJefzS5p0g36RvwpYTPIueksE/wZwcmtRSZKk1g2UBFTVm5NsB365KXp8VZ3XXliSJKlta1ks6NwkVwI3B0iysaq+1lpkkiSpVYPODnh0kguArwAfAy4BPthiXJIkqWWDXjvgJOD+wJer6khgM3Bma1FJkqTWDZoEXN8sDnSTJDepqo8Cq0492BdJ/k+SLyb5fJL3Lq9PkGRTkh8kOau5/UMbzy9JUlcMOibg6iS3BD4OvCXJFcD3WorpDODEqtqR5C+AE4H/0ey7qKqOaul5JUnqlBVbApLcrNl8DPAD4PeBfwcuAn6tjYCq6sNVtaO5eyZwhzaeR5KkrlutO2Cx+fkPVbWzqnZU1alV9Zp1unbAs9h1AOKRST6X5GNJHrQOzy9J0tRarTtgJsmTgV9M8vjdd1bVe/blSZN8BLj9Hna9tKpOa455KbADeEuz73JgY1V9O8l9gfcluUdVXbOH8x8PHA+wcePGfQlRkqSpt1oScALwFOAAfrL5v4B9SgKq6qEr7U/yDOBRwOZqLm5QVdcB1zXb25NcBNwV+IkLA1TVKcAp0Lt2wL7EqOFwfX1JGl8rJgFV9Z9JPglcWlXrskxwkocDfwj8UlV9v6/8EOCqqtqZ5E7AXYCL1yMm7ZvFxUU2b97M0tISMzMzXv1LksbMqlMEq+oGetcKWC9/C9wKOGO3qYAPBj6f5CzgX4ETquqqdYxLa7SwsMDS0hI7d+5kaWmJhYWFUYc0MouLi2zZsoXFxcXVD5akdTLoFMGtSX4deM9y83xbqupn9lL+buDdbT63hmt+fp6ZmZkftQTMz8+POqSRsEVE0rgaNAl4DvBCYEeSH9K7kmBV1a1bi0wTb25ujq1bt3Z+TMCeWkS6+ruQNF4GvYrgrdoORNNpbm6u8194tohIGlcDJQFJHryn8qr6+HDDkaaPLSKSxtWg3QEv7tu+OXAMsB14yNAjkqaQLSKSxtGg3QG7rBGQ5Ajgr1qJSJIkrYtBryK4u0uBnxtmIJIkaX0NOibgb+itEAi9xOEo4LNtBSVJkto36JiA/qV5dwBvq6pPtBCPJElaJ4OOCTh1eTvJgcARrUUkSZLWxUBjApIsJLl1koPodQO8Lsmr2w1NkiS1adCBgbdpLtn7eODNVXU/YHN7YWmcuQ6+JE2HQccE7JfkMOA44KUtxqMx5zr4kjQ9Bm0J+DPgQ8CFVfWZ5lK+F7QXlsaVVwaUpOkxUBJQVe+qqntV1XOb+xdX1a8v709yYlsBarwsr4O/YcMG18GXpAk3aHfAap4AbBnSuTTGXAdfkqbHsJKADOk8mgCugy9J02Fflw3eXa1+iCRJGifDSgJsCZAkacIMKwl415DOI0mS1smgFxC6OfBs4B7AzZfLq+pZzc//3Up0utEWFxcdxCdJ2qNBBwb+M/BF4Fh6awY8BTi/raA0HC7sI0layaDdAT9TVS8DvtdcTOhXgfu1F5aGwYV9JEkrGTQJuL75eXWSewK3AW7XTkgaFhf2kSStZNDugFOaSwi/DDgduCXwJ61FpaFwYR9J0kpSNd1T/GdnZ2vbtm2jDkOSpHWTZHtVza523EDdAUkOTfL6JB9s7t89ybNvbJB7ea6XJ7ksyVnN7ZF9+05McmGSLyU5to3n7yovDyxJ3TNod8CbgDfy48sIfxl4B/D6FmICeHVVvbK/IMndgSfSm6b408BHkty1qna2FENnOItAkrpp0IGBB1fVO4EbAKpqB7DeX76PAd5eVddV1VeAC4Fj1jmGqeQsAknqpkGTgO8luS3NNQKS3B/4bmtRwfOTfD7JG5oBiQCHA1/vO+bSpkw3krMIJKmbBu0OeCG9WQF3TvIJ4BDgN/b1SZN8BLj9Hna9FPh74CR6CcdJwKuAZ63x/McDxwNs3LhxX8PsDGcRSFI3DTw7IMl+wN3oXSzoS1V1/SoPudGSbAL+rarumeREgKra0uz7EPDyqlpxJJuzAyRJXTPs2QFPAH6qqs4FHgu8I8nRNzLGvT3XYX13Hwec02yfDjwxyc2SHAncBfh0GzFIktQFg3YHvKyq3pXkgcBm4JX0mu3bWDr4FUmOotcdcAnwHICqOjfJO4HzgB3A85wZIEnSvhs0CVj+sv1V4HVV9f4k/6uNgKrqqSvsOxk4uY3nlSSpawadHXBZkn8EfhP4QJKbreGxkiaMi0dJ3TBoS8BxwMOBV1bV1U2//YvbC0vSqLh4lNQdA/03X1Xfr6r3VNUFzf3Lq+rD7YYmaRRcPErqDpv0Je3CxaOk7hi0O0BSR7h4lNQdJgGSfsLc3Jxf/lIH2B0gSVJHmQRIktRRJgGSJHWUSYAkSR1lEiBJUkeZBEh9XC5XUpc4RVBquFyupK6xJUBquFyupK4xCZAaLpcrqWvsDpAaLpcrqWtMAqQ+LpcrqUvsDpAkqaNMAiRJ6iiTAEmSOsokQJKkjjIJkCSpo0wCJEnqKJMASZI6yiRAkqSOGrvFgpK8A7hbc/cA4OqqOirJJuB84EvNvjOr6oT1j1CSpOkwdklAVf3m8naSVwHf7dt9UVUdtf5RSZI0fcYuCViWJMBxwENGHYskSdNonMcEPAj4ZlVd0Fd2ZJLPJflYkgeNKjBJkqbBSFoCknwEuP0edr20qk5rtp8EvK1v3+XAxqr6dpL7Au9Lco+qumYP5z8eOB5g48aNww1ekqQpMZIkoKoeutL+JPsBjwfu2/eY64Drmu3tSS4C7gps28P5TwFOAZidna3hRS5J0vQY1+6AhwJfrKpLlwuSHJJkQ7N9J+AuwMUjik+SpIk3rgMDn8iuXQEADwb+LMn1wA3ACVV11bpHJknSlBjLJKCqnrGHsncD717/aCRJmk7j2h0gSZJaZhIgSVJHmQRIktRRJgGSJHWUSYAkSR1lEiBJUkeZBEiS1FEmAZIkdZRJgCRJHWUSIElSR5kESJLUUSYBkiR1lEmAJEkdZRIgSVJHmQRIktRRJgGSJHWUSYAkSR1lEiBJUkeZBEiS1FEmAZIkdZRJgCRJHWUSIElSR5kESJLUUSYBkiR1lEmA1BGLi4ts2bKFxcXFUYciaUyMLAlI8oQk5ya5IcnsbvtOTHJhki8lObav/OFN2YVJXrL+UUuTaXFxkc2bN/Oyl72MzZs3mwhIAkbbEnAO8Hjg4/2FSe4OPBG4B/Bw4O+SbEiyAXgt8Ajg7sCTmmMlrWJhYYGlpSV27tzJ0tISCwsLow5J0hjYb1RPXFXnAyTZfddjgLdX1XXAV5JcCBzT7Luwqi5uHvf25tjz1idiaXLNz88zMzPD0tISMzMzzM/PjzokSWNgZEnACg4Hzuy7f2lTBvD13crvt6cTJDkeOB5g48aNLYQoTZa5uTm2bt3KwsIC8/PzzM3NjTokSWOg1SQgyUeA2+9h10ur6rS2nreqTgFOAZidna22nkeaJHNzc375S9pFq0lAVT10Hx52GXBE3/07NGWsUC5JktZoHKcIng48McnNkhwJ3AX4NPAZ4C5JjkwyQ2/w4OkjjFOSpIk2sjEBSR4H/A1wCPD+JGdV1bFVdW6Sd9Ib8LcDeF5V7Wwe83zgQ8AG4A1Vde6IwpckaeKlarq7zGdnZ2vbtm2jDkOSpHWTZHtVza523Dh2B0iSpHVgEiBJUkeZBEiS1FEmAZIkdZRJgCRJHTX1swOSXAl8dcinPRj41pDPOQ6mtV5g3SaVdZtM1m307lhVh6x20NQnAW1Ism2QqReTZlrrBdZtUlm3yWTdJofdAZIkdZRJgCRJHWUSsG9OGXUALZnWeoF1m1TWbTJZtwnhmABJkjrKlgBJkjrKJABI8oYkVyQ5p6/s3kkWk3whyf9Ncuum/KZJTm3Kz09yYt9jHp7kS0kuTPKSUdRld0Os2yVN+VlJxuKKTGus20ySNzblZyeZ73vMfZvyC5O8JklGUJ1dDLFuC8178qzmdrsRVOdHkhyR5KNJzktybpIXNOUHJTkjyQXNzwOb8jSvyYVJPp/k6L5zPb05/oIkTx9VnfriGWbddva9ZiO/ZPo+1O1nm/fqdUn+YLdzjdXn5JDrNnafk6uqqs7fgAcDRwPn9JV9BvilZvtZwEnN9pOBtzfbtwAuATbRu7zxRcCdgBngbODu01C35v4lwMGjrs+NqNvzgDc227cDtgM3ae5/Grg/EOCDwCOmqG4LwOyo69NXh8OAo5vtWwFfBu4OvAJ4SVP+EuAvmu1HNq9JmtfoU035QcDFzc8Dm+0Dp6Fuzb7/HvVrdSPrdjvgF4CTgT/oO8/YfU4Oq27NvksYs8/J1W62BABV9XHgqt2K7wp8vNk+A/j15cOB/ZPsB/wUsARcAxwDXFhVF1fVEvB24DFtx76aIdVtLK2xbncH/l/zuCuAq4HZJIcBt66qM6v3V/xm4LFtx76aYdRtHcJcs6q6vKo+22xfC5wPHE7vb+XU5rBT+fFr8BjgzdVzJnBA85odC5xRVVdV1Xfo/T4evo5V+QlDrNvYWWvdquqKqvoMcP1upxq7z8kh1m0imQTs3bn8+M35BOCIZvtfge8BlwNfA15ZVVfRe9N8ve/xlzZl42itdYNegvDhJNuTHL+ewa7R3up2NvDoJPslORK4b7PvcHqv1bJJfN32Vrdlb2yaJ182Dl0dy5JsAu4DfAo4tKoub3Z9Azi02d7b39VY/73dyLoB3DzJtiRnJhl5UtpvwLrtzTS8biuZlM/JHzEJ2LtnAc9Nsp1eE9FSU34MsBP4aeBI4EVJ7jSaEPfZvtTtgVV1NPAI4HlJHrzOMQ9qb3V7A70PnG3AXwGfpFfXSbIvdXtKVf088KDm9tR1jXgvktwSeDfwe1W1S2tT0yIzsdOWhlS3O1ZvVbonA3+V5M7Dj3TtfN1WNSmfkz+L+oACAAACTklEQVRiErAXVfXFqnpYVd0XeBu9fizo/VH+e1Vd3zS9foJe0+tl7Prf1x2asrGzD3Wjqi5rfl4BvJdewjB29la3qtpRVb9fVUdV1WOAA+j1/V1G77VaNnGv2wp163/drgXeyhi8bkluSu/D9i1V9Z6m+JvLTeHNzyua8r39XY3l39uQ6tb/ul1Mb1zHfVoPfhVrrNveTMPrtleT8jnZzyRgL9KMok5yE+CPgX9odn0NeEizb396A3q+SG/Q1l2SHJlkBngiMPJRvXuy1rol2T/JrfrKHwacs/t5x8He6pbkFk3sJPkVYEdVndc0912T5P5NU/nTgNNGE/3K1lq3pnvg4Kb8psCjGPHr1vyOXw+cX1V/2bfrdGB5hP/T+fFrcDrwtGYk/f2B7zav2YeAhyU5sBm1/bCmbGSGVbemTjdrznkw8ADgvHWpxF7sQ932Zuw+J4dVt0n6nNzFMEcZTuqN3n9Vl9Mb6HEp8GzgBfT+m/oy8Of8eGGlWwLvotc/ex7w4r7zPLI5/iLgpaOu17DqRm8k79nN7dwJrdsm4Ev0Bv18hF5z6/J5Zun9sV4E/O3yYya9bsD+9GYKfL553f4a2DDiej2QXrPq54GzmtsjgdsCW4ELmjoc1Bwf4LXNa/MF+mY60OseubC5PXMMXrOh1A34xeb+2c3PZ09g3W7fvG+voTdQ9VJ6A3BhzD4nh1U3xvRzcrWbKwZKktRRdgdIktRRJgGSJHWUSYAkSR1lEiBJUkeZBEiS1FEmAZIkdZRJgCRJHWUSIElSR/1/5tpKi9E0TggAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot measured data (labels are a bit off)\n",
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "ax.plot(annual_mean_df.year, annual_mean_df.height, 'k.')\n",
    "ax.set_ylabel(\"{} [{}]\".format(ssh.standard_name, 'mm'))\n",
    "ax.set_title('Annual mean');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf8AAAE/CAYAAABIGcoKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcXFWd9/HPlybNDklMwJAEwuYoLiy2QLtgS2R1gcdRH8SRgAs6ooOKziNujEbMuILMuOEQBQURFJ2oMJhpaR20WDoYwCQsARESogSTsI40Cb/nj3sKbppeqpO6Xcv9vl+vevW9526/U1Vdv7rnnntKEYGZmZmVxxaNDsDMzMzGl5O/mZlZyTj5m5mZlYyTv5mZWck4+ZuZmZWMk7+ZmVnJOPmbjYGkHkkrGh3HUCS9W9I5ddjPv0j6fj1iakWbU39JL5N0h6RHJB1X79iGON4SST2buG1I2jtNf1nSP9Y1OGtqTv7W1CT1SVoraatGx9LMJHUCnwC+2OhY6k3SSZKuaXQcNfoM8O8RsX1E/LTog0XE8yOirw67+hLwsfQ+shJw8remJWkW8AoggNc3NJjmdyxwa0SsbHQgJbc7sKTog0jasp77i4hVwK34/6w0nPytmZ0IXAt8F5iTXyDpu5K+JukXkh6WdJ2kvXLLQ9J7UhPsurSu0rKNmnUlzUrrb5nmT5a0LO33LknvrjXgtJ/3puM+LGmupL0k/U7SQ5IurZ5dSZok6eeSVqfWjZ9LmpGWTZa0QtLr0vz2kpZLOnGYQx8N/DoXx5WS3jcotpskvSFNf1XSvSmmRZJeMUx9nnGZQ9Ldkl6dpreQ9FFJd0r6a6rf5LRsa0nfT+XrJN0gaZdhjlPdx8OSlkr6P6n8ecA3ge7UlL5umO13knS+pFWSVkr6rKSOtOwkSddI+lJ6nv8o6ejctntI+nU69kJgyjDPcXX9d6XXYo2kBZJ2TeV3AnsCP0uxPqO1Kj13Z6Q6rpX0HUlb55a/VtLi9Hz9TtKLBm37/yTdDDwqactBr8VWks6RdF96nJOPQdJH0vNzn6S3D1G1PuA1I9Xd2oeTvzWzE4GL0uPIIRLH8cCngUnAcuCsQctfC7wEeBHwZuDIGo97f9p2R+Bk4GxJB44h7iOBFwOHAP8MnAf8AzATeAHwlrTeFsB3yM4WdwP+F/h3gIhYA7wd+LaknYGzgcURceEwx3whcFtu/ge54yBp33ScX6SiG4D9gcnAxcBl+SQ0Bu8HjgNeCewKrAW+lpbNAXZK9X4W8J5Ux6HcSdbKsxPZa/p9SdMiYlnarpKa0icOs/13gfXA3sABwBHAO3PLDyZ7fqYAXwDOr34ZJKv/orRsLoO+aOZJOgyYR/Z+mgb8CbgEICL2Au4BXpdifXyY3byV7D2yF/Acsss1SDoAmA+8m+z5+hawYNCXiLeQJeiJEbF+0H4/Tvae2x/YDzgot++jgA8DhwP7AK8eIq5laTsrg4jww4+mewAvB54ApqT5W4EP5pZ/F/iP3PwxZM3e1fkAXp6bvxT4aJr+F+D7uWWz0vpbDhPLT4HT0nQPsGKEuAN4WW5+EfD/cvNfBs4ZZtv9gbWDyv4NuAVYCTxrhOPeARyVm98BeBTYPc2fBcwfYfu1wH6Dn5+h6gvcDbw6TS8DZueWTUuv25ZkX15+B7xoE17/xcCxafok4JoR1t0FeBzYJlf2FuDq3PbLc8u2Ta/Ts8m+dK0Htsstvzj//hh0rPOBL+Tmt0/1nTX4uRlm+7uB9wx6396Zpr8BzB20/m3AK3Pbvn2E1+JO4JjcsiOBu9P0fOBfc8uek56DvXNlhwN3jfW18qM1Hz7zt2Y1B/hlRDyQ5i/mmWdkf85NP0b2QTyW5UOSdLSka1Oz7jqyD+gRm4IH+Utu+n+HmN8+HWdbSd+S9CdJDwG/ASZWm6uT88haC74bEX8d4ZhryRI+ABHxMNlZ/vGp6C1kLSjVOn5Y2aWNB1MddxpjHat2B36SmqnXkX0Z2ECWkL8HXAVckpqavyBpwlA7kXRirrl7XapzrfHsDkwAVuW2/xawc26dp94LEfFYmtye1FoREY/m1v3TCMfaNb88Ih4B/gpMrzFWgHsHHWvXXD1Or9Yh1WNmbvngbUeMbdC+dx3iuIPtAAx5WcXaj5O/NR1J25A1q75S0p8l/Rn4ILCfpHo0Sz5KdvZX9ezcsbcCfkzW+3mXyJqZrwBE/Z0O/B1wcETsCBxaDSPF0kGW/C8E3qt0W9YwbiY7m8v7AfAWSd3A1sDVab+vILsc8WZgUqrjgwxdx42eqxTT1Nzye4GjI2Ji7rF1RKyMiCci4tMRsS/wUrJLKc/osyBpd+DbwPvIWjcmAn/IxTPaT4/eS3bmPyUXw44R8fxRtgNYBUyStF2ubLcR1r+PLElXY9+OrIl+LB0tZw461n1p+l7grEHP5bYR8YPc+iM9FxvFNmjfq4Y47mDPA26qpQLW+pz8rRkdR3b2uC9ZU/j+ZB9M/8MQyWMTLAYOlbSbpJ2AM3LLOoGtgNXA+tQx7Ig6HHMoO5C1BKxLneTOHLT8Y2Qf9m8nu4XvwkGtAnlXkF13H1y2O9ntZz+MiCdzx11PVsctJX2KrH/DUG4Htpb0mnTW/gmy56fqm8BZKYEjaaqkY9P0qyS9MMX8EFnz+JM803apnqvTdieTnflX/QWYoWFuQ4usp/ovgS9L2lFZJ8S9JA1+Poba9k9AP/BpSZ2SXg68boRNfgCcLGn/9EXxc8B1EXH3aMfKOVXSjPSafxz4YSr/NvAeSQcrs1163ncYflfPiO0T6TWYAnwKqHZsvRQ4SdK+krblme81yN4/V46hHtbCnPytGc0BvhMR90TEn6sPss5wb9Vm3uYUEQvJPnBvJrsm//PcsoeBfyL7sFwLnAAs2JzjjeAcYBvgAbK7Gv6rukDSi4EPASdGxAbg82QJ8qPD7OtnwHOrPc8BIutwdjlZ566Lc+telY51O1nz798Ypjk5Ih4E3gv8B9nZ7aNAvvf/V8men19KejjV4+C07NnAj8gS/zKyuxG+N8QxlpL1haiQJfoXAr/NrfIrstvn/izpgcHbJyeSfXFbSva6/Yis/0EtTkgxryFLisN1qiQi/hv4JFnr0CqyTnvHD7f+MC4m+7JyF9l1+s+mffcD7yJ7n68l68R60hj2+1myLzI3k/UTuTG37yvJ3m+/Svv9VX5DSdPIvmwXPjaBNQdFjNaiZmatQNIpwL4R8YFGx2JDk3Q38M70JaJpSPoyWcfDrzc6FhsfTv5mZuOkWZO/lY+b/c3MzErGZ/5mZmYl4zN/MzOzknHyNzMzK5m6/jJUM5kyZUrMmjWr0WGYmZmNm0WLFj0QEVNHW69tk/+sWbPo7+9vdBhmZmbjRtJIw1M/xc3+ZmZmJePkb2ZmVjJO/mZmZiXj5G9mZlYyTv5mZmYl4+RvZmZWMk7+ZmZmJePkb2ZmVjJO/mZmVjqVSoV58+ZRqVQaHUpDtO0If2ZmZkOpVCrMnj2bgYEBOjs76e3tpbu7u9FhjSuf+ZuZWan09fUxMDDAhg0bGBgYoK+vr9EhjTsnfzMzK5Wenh46Ozvp6Oigs7OTnp6eRoc07tzsb2ZmpdLd3U1vby99fX309PSUrskfnPzNzKyEuru7S5n0q9zsb2ZmVjJO/mZmZiXj5G9mZlYyTv5mZmYl4+RvZmZWMk7+ZmZmJePkb2ZmbaHs4/WPhe/zNzOzlufx+sfGZ/5mZtbyPF7/2BSa/CVtLel6STdJWiLp06lcks6SdLukZZL+KVd+rqTlkm6WdGBuX3Mk3ZEec4qM28zMWovH6x+bopv9HwcOi4hHJE0ArpF0JfA8YCbw3Ih4UtLOaf2jgX3S42DgG8DBkiYDZwJdQACLJC2IiLUFx29mZi3A4/WPTaHJPyICeCTNTkiPAP4ROCEinkzr3Z/WORa4MG13raSJkqYBPcDCiFgDIGkhcBTwgyLjNzOz1lH28frHovBr/pI6JC0G7idL4NcBewH/V1K/pCsl7ZNWnw7cm9t8RSobrtzMzMzGqPDkHxEbImJ/YAZwkKQXAFsBf4uILuDbwPx6HEvSKekLRf/q1avrsUszM7O2M269/SNiHXA1WXP9CuDytOgnwIvS9EqyvgBVM1LZcOWDj3FeRHRFRNfUqVPrWwEzM7M2UXRv/6mSJqbpbYDDgVuBnwKvSqu9Erg9TS8ATky9/g8BHoyIVcBVwBGSJkmaBByRyszMzGyMiu7tPw24QFIH2ReNSyPi55KuAS6S9EGyDoHvTOtfARwDLAceA04GiIg1kuYCN6T1PlPt/GdmZmZjo6xjffvp6uqK/v7+RodhZmY2biQtSv3pRuQR/szMzErGyd/MzKxknPzNzMxKxsnfzMysZJz8zczMSsbJ38zMrGSc/M3MzErGyd/MzKxknPzNzMxKxsnfzMysZJz8zczMSsbJ38zMrGSc/M3MzErGyd/MzKxknPzNzMxKxsnfzMysZJz8zcxaXKVSYd68eVQqlUaHYi1iy0YHYGZmm65SqTB79mwGBgbo7Oykt7eX7u7uRodlTc5n/mZmLayvr4+BgQE2bNjAwMAAfX19jQ7JWoCTv5lZC+vp6aGzs5OOjg46Ozvp6elpdEjWAtzsb2bWwrq7u+nt7aWvr4+enh43+VtNnPzNzFpcd3e3k76NiZv9zczMSsbJ38zMrGSc/M3MzErGyd/MzKxknPzNzMxKxsnfzMysZJz8zczMSqbQ5C9pa0nXS7pJ0hJJnx60/FxJj+Tmt5L0Q0nLJV0naVZu2Rmp/DZJRxYZt5mZWTsr+sz/ceCwiNgP2B84StIhAJK6gEmD1n8HsDYi9gbOBj6f1t0XOB54PnAU8HVJHQXHbmZm1pYKTf6RqZ7ZT0iPSIn7i8A/D9rkWOCCNP0jYLYkpfJLIuLxiPgjsBw4qMjYzczM2lXh1/wldUhaDNwPLIyI64D3AQsiYtWg1acD9wJExHrgQeBZ+fJkRSozMzOzMSp8bP+I2ADsL2ki8BNJhwJvAnrqfSxJpwCnAOy222713r2ZmVlbGLfe/hGxDrgaeBWwN7Bc0t3AtpKWp9VWAjMBJG0J7AT8NV+ezEhlg49xXkR0RUTX1KlTi6qKmZlZSyu6t//UdMaPpG2Aw4FFEfHsiJgVEbOAx1IHP4AFwJw0/UbgVxERqfz4dDfAHsA+wPVFxm5mZtauim72nwZckDr4bQFcGhE/H2H984HvpZaANWQ9/ImIJZIuBZYC64FT0+UEs7ZRqVT8m+xmNi6UnVi3n66urujv7290GGY1qVQqzJ49m4GBATo7O+nt7fUXADMbM0mLIqJrtPU8wp9ZE+jr62NgYIANGzYwMDBAX19fo0Myszbm5G/WBHp6eujs7KSjo4POzk56enoaHZKZtbERr/lLmlzDPp5MPfnNbBN1d3fT29vra/5mNi5G6/B3X3pohHU6AN9Ub7aZuru7nfTNbFyMlvyXRcQBI60g6fd1jMfMzMwKNto1/1pOQ3yqYmZm1kJGPPOPiL9VpyVNIhtlb8vc8hvz65iZmVnzq2mQH0lzgZOAO4HqwAABHFZMWGZmZlaUWkf4ezOwV0QMFBmMmZmZFa/W+/z/AEwsMhAzMzMbH7We+c8Dfi/pD8Dj1cKIeH0hUZmZmVlhak3+FwCfB24BniwuHDMzMytarcn/sYg4t9BIzMzMbFzUmvz/R9I8YAEbN/vfWEhUZnXgn8g1Mxtarcm/OsrfIbky3+pnTcs/kWtmNryakn9EvKroQMzqaaifyHXyNzPL1DrIz0TgRGAWG4/w90/FhGW2eao/kVs98/dP5JqZPa3WZv8rgGtxb39rEf6JXDOz4dWa/LeOiA8VGolZnfkncs3MhlbrCH/fk/QuSdMkTa4+Co3MzMzMClHrmf8A8EXg42z8wz57FhGUmZmZFafW5H86sHdEPFBkMGZmZla8Wpv9lwOPFRmImZmZjY9az/wfBRZLupqNR/jzrX5mZmYtptbk/9P0MDMzsxZX6wh/FxQdiJmZmY2PEa/5SzpvtB3Uso6ZmZk1j9HO/I+T9LcRlgvwuP9mZmYtZLTk/5Ea9vE/wy2QtDXwG2CrdKwfRcSZki4CuoAngOuBd0fEE5IEfBU4huzugpOqPxssaQ7wibTrz/pShJmZ2aYZMfnXIcE+DhwWEY9ImgBcI+lK4CLgH9I6FwPvBL4BHA3skx4Hp7KD02iCZ5J9YQhgkaQFEbF2M+MzMzMrnVrv898kkXkkzU5Ij4iIK9KyIDvzn5HWORa4MC26FpgoaRpwJLAwItakhL8QOKrI2M3MzNpVockfQFKHpMXA/WQJ/LrcsgnA24D/SkXTgXtzm69IZcOVm5mZ2RjVlPwlvamWsqFExIaI2J/s7P4gSS/ILf468JuIGLbfwFhIOkVSv6T+1atX12OXZmZmbafWM/8zaiwbVkSsA64mNddLOhOYCuR/KnglMDM3PyOVDVc++BjnRURXRHRNnTp1LOGZmZmVxogd/iQdTdbzfrqkc3OLdgTWj7ZzSVOBJyJinaRtgMOBz0t6J9l1/NkR8WRukwXA+yRdQtbh78GIWCXpKuBzkial9Y5gjF8+zBqhUqnQ19dHT08P3d3djQ7HzAwY/Va/+4B+4PXAolz5w8AHa9j/NOACSR1krQyXRsTPJa0H/gRUsrv7uDwiPgNcQfZlo/pDQicDRMQaSXOBG9J+PxMRa2o4vlnDVCoVZs+ezcDAAJ2dnfT29voLgJk1hdFu9bsJuEnSxRHxxFh3HhE3AwcMUT7kcVPv/1OHWTYfmD/WGMwapa+vj4GBATZs2MDAwAB9fX1O/mbWFGr9YZ+DJP0LsHvaRmS5es+iAjNrdT09PXR2dj515t/T09PokMzMgNqT//lkzfyLgA3FhWPWPrq7u+nt7fU1fzNrOrUm/wcj4spCIzFrQ93d3U76ZtZ0Ruvtf2CavFrSF4HLyYbsBaA67r6ZmZm1jtHO/L88aL4rNx3AYfUNx8zMzIo2Wm9//1yvmZlZm6npmr+kDw1R/CCwKCIW1zcks+F50Bwzs81Xa4e/rvT4WZp/LXAz8B5Jl0XEF4oIzizPg+aYmdVHrWP7zwAOjIjTI+J04MXAzsChwEkFxWa2kaEGzTEzs7GrNfnvTK6XP/AEsEtE/O+gcrPCVAfN6ejo8KA5ZmabodZm/4uA6yT9Z5p/HXCxpO2ApYVEZjaIB80xM6sPZcPp17Ci1AW8LM3+NiL6C4uqDrq6uqK/v6lDNDMzqytJiyKia7T1RhvkZ8eIeEjSZOCu9Kgum+xf1jMzM2s9ozX7X0zWs38R2aA+GvTXP+xjZmbWYkYb5Oe16e8e4xOOmZmZFa2m3v7K/IOkT6b53SQdVGxoZmZmVoRab/X7OtANnJDmHwa+VkhEZmZWmEqlwrx586hUKo0OxRqo1lv9Do6IAyX9HiAi1krqLDAuMzOrM4+SaVW1nvk/IamDrJMfkqYCTxYWlZmZ1Z1HybSqWpP/ucBPgJ0lnQVcA3yusKjMbERuurVN4VEyraqmZv+IuEjSImA22W1+x0XEskIjM7MhuenWNpVHybSqWq/5A9wBPFTdRtJuEXFPIVGZ2bCGarr1h7jVqru72+8Xqy35S3o/cCbwF2ADTw/y86LiQjOzoVSbbqtn/m66NbOxqvXM/zTg7yLir0UGYyOrVCpurjM33Y4D/69Zu6s1+d8LPFhkIDYyX+e1PDfdFsf/a1YGo/2wz4fS5F1An6RfAI9Xl0fEVwqMzXJ8nddsfPh/zcpgtDP/HdLfe9KjMz1snPk6r9n48P+alYEiYvN3Iv1bRLy/DvHUTVdXV/T39zc6jLrydUiz8eH/NWtVkhZFRNeo69Up+d8YEQdu9o7qqB2Tv5mZ2UhqTf61jvC3qUFsLel6STdJWiLp06l8D0nXSVou6YfV3wmQtFWaX56Wz8rt64xUfpukI4uM28zMrJ0VmvzJOgceFhH7AfsDR0k6BPg8cHZE7A2sBd6R1n8HsDaVn53WQ9K+wPHA84GjgK+n3xowMzOzMapX8tdQhZF5JM1OSI8ADgN+lMovAI5L08emedLy2ZKUyi+JiMcj4o/AcuCgOsVuZmZWKmNK/pK2HWbRV0fYpkPSYuB+YCFwJ7AuItanVVYA09P0dLIxBUjLHwSelS8fYhszMzMbg5qSv6SXSloK3Jrm95P09eryiPjucNtGxIaI2B+YQXa2/tzNinjkOE+R1C+pf/Xq1UUdxszMrKXVeuZ/NnAk8FeAiLgJOHQsB4qIdcDVQDcwUVJ1jIEZwMo0vRKYCZCW75SO+VT5ENvkj3FeRHRFRNfUqVPHEp6ZmVlp1NzsHxH3DiraMNo2kqZKmpimtwEOB5aRfQl4Y1ptDvCfaXpBmict/1Vk9yIuAI5PdwPsAewDXF9r7GZmZva0msf2l/RSICRNIPuhn2U1bDcNuCD1zN8CuDQifp4uIVwi6bPA74Hz0/rnA9+TtBxYQ9bDn4hYIulSYCmwHjg1Ikb98mFmZmbPVNMgP5KmkHXqezVZz/5fAqc186/8eZAfMzMrm1oH+anpzD8iHgDeutlRmZmZWcPV2tv/guq1+zQ/SdL84sIyMzOzotTa4e9Fqbc+ABGxFjigmJDMzMysSLUm/y0kTarOSJpM7Z0FzczMrInUmsC/DFQkXUbW4e+NwFmFRWVmZmaFqbXD34WSFgGvSkVviIilxYVlZmZmRam56T7da78a2BpA0m4RcU9hkZlZXVQqFfr6+ujp6aG7u7vR4bQdP7/WimpK/pJeT9b0vyvZD/TsTjbIz/OLC83MNlelUmH27NkMDAzQ2dlJb2+vE1Qd+fm1VlVrh7+5wCHA7RGxBzAbuLawqMysLvr6+hgYGGDDhg0MDAzQ19fX6JDaip9fa1W1Jv8n0mh+W0jaIiKuBkYdQcjMGqunp4fOzk46Ojro7Oykp6en0SG1FT+/1qpqvea/TtL2wG+AiyTdDzxaXFhmVg/d3d309vb6mnRB/PxaqxpxbH9JW0XE45K2A/5GdpvfW8l+avcij+1vZmbWPOo1tn8FOBD4ZkS8LZVdsLnBmZmZWeOMlvw7JZ0AvFTSGwYvjIjLiwnLzMzMijJa8n8PWTP/ROB1g5YF4ORvZmbWYkZM/hFxjaTfASsiwsP5mpmZtYFRb/WLiCfJxvI3MzOzNlDrff69kv5ekgqNxqwFVCoV5s2bR6VSaXQoZmabpNb7/N8NfAhYL6l6y19ExI6FRWbWhDycq5m1g5rO/CNih4jYIiI6I2LHNO/Eb6Xj4Vxtc7jVyJpFrT/sc+hQ5RHxm/qGY9bcqsO5Vs/8PZyr1cqtRtZMam32/0huemvgIGARcFjdIzJrYh7O1TbVUK1Gfv9Yo9SU/CNio3v8Jc0EzikkIrMm193d7Q9tGzO3GlkzqfXMf7AVwPPqGYiZWTtzq5E1k1qv+f8b2Yh+kHUS3B+4saigrDlVKhV/cJltBrcaWbOo9cw///N464EfRMRvC4jHmpQ7K5mZtY9ar/k/9Ut+kiYBMwuLyJqSOyuZmbWPmu7zl9QnaUdJk8ma+78t6exiQ7NmUu2s1NHR4c5KZmYtrtbhfXeKiIeANwAXRsTBwOzRNpI0U9LVkpZKWiLptFS+v6RrJS2W1C/poFQuSedKWi7pZkkH5vY1R9Id6TFn7FW1zVHtrDR37lw3+ZvZZvFgR41X6zX/LSVNA94MfHwM+18PnB4RN0raAVgkaSHwBeDTEXGlpGPSfA9wNLBPehwMfAM4OLU4nAl0kXU8XCRpQUSsHUMstpncWcnMNpf7DzWHWs/8PwNcBSyPiBsk7QncMdpGEbEqIm5M0w8Dy4DpZAm8OjzwTsB9afpYspaFiIhrgYnpS8eRwMKIWJMS/kLgqBpjNzOzJuEhsptDrR3+LgMuy83fBfx9dV7SGRExb6R9SJoFHABcB3wAuErSl8i+gLw0rTYduDe32YpUNly5mZm1EA921BxqPfMfzZtGWihpe+DHwAdS34F/BD4YETOBDwLn1yMISaekPgT9q1evrscuzcysjtx/qDls6gh/g2nYBdIEssR/UURcnornAKel6cuA/0jTK9n4NsIZqWwlWZ+AfHnf4GNFxHnAeQBdXV0xeLmZmTWe+w81Xr3O/IdMtJJEdla/LCK+klt0H/DKNH0YT/cfWACcmHr9HwI8GBGryPobHCFpUhpn4IhUZmZmZmNU9Jn/y4C3AbdIWpzKPga8C/iqpC2BvwGnpGVXAMcAy4HHgJMBImKNpLnADWm9z0TEmjrFbmZmVir1Sv6XDVUYEdcw/BeDFw+xfgCnDrOv+cD8TQ3QzBrHvwsxdn7OrEi1/rDP1sA7gOcDW1fLI+Lt6e/nConOzFqe7+seOz9nVrRar/l/D3g22f32vybrcPdwUUGZWfvwfd1j5+fMilZr8t87Ij4JPJp+5Oc1ZCPwmZmNyL8LMXZ+zqxotV7zfyL9XSfpBcCfgZ2LCcnM2kn1vm5fv66dnzMrmrI+dqOsJL2T7F79FwHfAbYHPhUR3yw2vE3X1dUV/f39jQ7DzMxs3EhaFBFdo61X6/C+1UF4fg3suTmBmZmZWWPVdM1f0i6Szpd0ZZrfV9I7ig3NzMzMilBrh7/vko2ot2uav53sx3nMzMysxdSa/KdExKXAkwARsR7YUFhUZmZmVphak/+jkp5FGsO/Ou5+YVGZmZlZYWq91e9DZD+6s5ek3wJTgTcWFpWZmVkba/TwzbX29r9R0iuBvyMbq/+2iHhilM3MzMxskGYYvrnW3v5vAraJiCXAccAPJR1YaGRWGpVKhXnz5lGpVBodiplZ4Zph+OZam/0/GRGXSXo5MBv4EvANPMSvbaZm+AZsZjaeqsM3Vz/3GjF8c60d/qo9+18DfDsifgF0FhNSc/LZaTGa4Ruwmdl4qg7fPHfu3Iad8NR65r9S0reAw4HPS9qK2r84tDyfnRanGb4Bm5mNt+7u7obmkVoT+JvJBvk5MiLWAZOBjxQWVZPx2WlxmuEbsJlZ2dTa2/8x4PLc/CpgVVFBNRufnRar0d+AzczKptZm/1Lzz2uamVk7cfKvUaudnTZ6AAkzs3bRjp+nTv5tyB0Uzczqo10/T0ty6L/ZAAALpklEQVTTY79MmqGDom+NNLN20Ayfp0XwmX8banQHxXb9pmxm5dPoz9OiOPm3oUZ3UBzqm7KTf/tpx+ugZoM1+vO0KE7+baqRHRTb9ZuyPc2tO1Ymrdbhuxa+5m9154F72l+7Xge1jbnvTvvymb8Voh2/KdvT3LrT/ty6096c/M1szNr1Oqg9zX132puTv5ltErfutDe37rS3Qq/5S5op6WpJSyUtkXRabtn7Jd2ayr+QKz9D0nJJt0k6Mld+VCpbLumjRcZtZlZ27rvT3oo+818PnB4RN0raAVgkaSGwC3AssF9EPC5pZwBJ+wLHA88HdgX+W9Jz0r6+RvaTwiuAGyQtiIilBcdvZlZabt1pX4Um//yv/0XEw5KWAdOBdwH/GhGPp2X3p02OBS5J5X+UtBw4KC1bHhF3AUi6JK3r5G9mZjZG43arn6RZwAHAdcBzgFdIuk7SryW9JK02Hbg3t9mKVDZc+eBjnCKpX1L/6tWr61+JAvhWGjMzG2/j0uFP0vbAj4EPRMRDkrYEJgOHAC8BLpW05+YeJyLOA84D6Orqis3dX9F8K42ZmTVC4Wf+kiaQJf6LIuLyVLwCuDwy1wNPAlOAlcDM3OYzUtlw5S3NA6U0D7fAmFmZFHrmL0nA+cCyiPhKbtFPgVcBV6cOfZ3AA8AC4GJJXyHr8LcPcD0gYB9Je5Al/eOBE4qMfTz4Vprm4BYYMyubopv9Xwa8DbhF0uJU9jFgPjBf0h+AAWBORASwRNKlZB351gOnRsQGAEnvA64COoD5EbGk4NgL54FSmoMHM2ld/nEhs02jLOe2n66urujv7290GNYCfObfmvy6mT2TpEUR0TXaev5hnwL4+nFr8WAmrcl9ZqzZtNJnv4f3rTOfjbQmD2bSetxnxppJq332+8y/znw2YjY+3GJjzaTVPvt95l9nPhsxGz9usbFm0Wqf/U7+deYe/GZm5dNqn/3u7W9mZtYm3NvfzMzMhuTkb2ZmVjJO/mZmZiXj5G9mTaWVBkqxp/l1ay3u7W9mTaPVBkqxjF+31uMzfzNrGkUOlOIz0+K02gA35jN/M2siRQ2U4jPTYrXaADfm5G9mTaSogVL8s83FarUBbszJ38yaTBFD9vrMtHgearm1OPmbWdvzmanZxpz8zawUfGZq9jT39jczMysZJ38zM7OScfI3MzMrGSd/MzOzknHyNzMzKxknfzMza1oelrkYvtXPzMyakodlLo7P/M3MrCn5B4OK4+RvZmZNqTosc0dHh4dlrjM3+5uZWVPysMzFcfI3M7Om5WGZi1Fos7+kmZKulrRU0hJJpw1afrqkkDQlzUvSuZKWS7pZ0oG5dedIuiM95hQZt5mZWTsr+sx/PXB6RNwoaQdgkaSFEbFU0kzgCOCe3PpHA/ukx8HAN4CDJU0GzgS6gEj7WRARawuO38zMrO0UeuYfEasi4sY0/TCwDJieFp8N/DNZMq86FrgwMtcCEyVNA44EFkbEmpTwFwJHFRl7WfgeWjOz8hm3a/6SZgEHANdJOhZYGRE3ScqvNh24Nze/IpUNV26bwffQmpmV07jc6idpe+DHwAfILgV8DPhUAcc5RVK/pP7Vq1fXe/dtx/fQmpmVU+HJX9IEssR/UURcDuwF7AHcJOluYAZwo6RnAyuBmbnNZ6Sy4co3EhHnRURXRHRNnTq1iOq0Fd9Da2ZWToU2+ytr0z8fWBYRXwGIiFuAnXPr3A10RcQDkhYA75N0CVmHvwcjYpWkq4DPSZqUNjsCOKPI2MvA99CamZVT0df8Xwa8DbhF0uJU9rGIuGKY9a8AjgGWA48BJwNExBpJc4Eb0nqfiYg1xYVdHr6H1vIqlYq/DJqVQKHJPyKuATTKOrNy0wGcOsx684H59YzPzJ7mDqBm5eGx/c0McAdQszJx8jczwB1AzcrEY/ubGeAOoGZl4uRvZk9xB1CzcnCzv5mZWck4+ZuZmZWMk7+ZmVnJOPmbmZmVjJO/mZlZyTj5m5mZlYyTv5mZWck4+ZuZmZWMk7+ZmVnJKPshvfYjaTXwpzrvdgrwQJ332QzatV7gurUq1601uW6Nt3tETB1tpbZN/kWQ1B8RXY2Oo97atV7gurUq1601uW6tw83+ZmZmJePkb2ZmVjJO/mNzXqMDKEi71gtct1blurUm161F+Jq/mZlZyfjM38zMrGRKnfwlzZd0v6Q/5Mr2k1SRdIukn0naMZVPkHRBKl8m6YzcNkdJuk3SckkfbURdBqtj3e5O5Ysl9TeiLoONsW6dkr6Tym+S1JPb5sWpfLmkcyWpAdXZSB3r1pfek4vTY+cGVOcpkmZKulrSUklLJJ2WyidLWijpjvR3UipXek2WS7pZ0oG5fc1J698haU6j6pSLp55125B7zRY0qk65eMZat+em9+rjkj48aF9N9TlZ57o13efkqCKitA/gUOBA4A+5shuAV6bptwNz0/QJwCVpelvgbmAW0AHcCewJdAI3Afu2Q93S/N3AlEbXZzPqdirwnTS9M7AI2CLNXw8cAgi4Eji6jerWB3Q1uj65OkwDDkzTOwC3A/sCXwA+mso/Cnw+TR+TXhOl1+i6VD4ZuCv9nZSmJ7VD3dKyRxr9Wm1m3XYGXgKcBXw4t5+m+5ysV93Ssrtpss/J0R6lPvOPiN8AawYVPwf4TZpeCPx9dXVgO0lbAtsAA8BDwEHA8oi4KyIGgEuAY4uOfTR1qltTGmPd9gV+lba7H1gHdEmaBuwYEddG9t97IXBc0bGPph51G4cwxywiVkXEjWn6YWAZMJ3sf+WCtNoFPP0aHAtcGJlrgYnpNTsSWBgRayJiLdnzcdQ4VuUZ6li3pjPWukXE/RFxA/DEoF013edkHevWkkqd/IexhKfflG8CZqbpHwGPAquAe4AvRcQasjfLvbntV6SyZjTWukH2xeCXkhZJOmU8gx2j4ep2E/B6SVtK2gN4cVo2ney1qmrF1224ulV9JzVDfrIZLmlUSZoFHABcB+wSEavSoj8Du6Tp4f6vmvr/bTPrBrC1pH5J10pq+JfRvBrrNpx2eN1G0iqfk09x8n+mtwPvlbSIrCloIJUfBGwAdgX2AE6XtGdjQtxkm1K3l0fEgcDRwKmSDh3nmGs1XN3mk33Q9APnAL8jq2sr2ZS6vTUiXgi8Ij3eNq4RD0PS9sCPgQ9ExEatS6kFpmVvP6pT3XaPbBS5E4BzJO1V/0jHzq/bqFrlc/IpTv6DRMStEXFERLwY+AHZdSrI/hn/KyKeSE2svyVrYl3JxmdbM1JZ09mEuhERK9Pf+4GfkH1RaDrD1S0i1kfEByNi/4g4FphIdm1vJdlrVdVyr9sIdcu/bg8DF9MEr5ukCWQfshdFxOWp+C/VJu/09/5UPtz/VVP+v9WpbvnX7S6yfhsHFB78KMZYt+G0w+s2rFb5nMxz8h9EqVe0pC2ATwDfTIvuAQ5Ly7Yj66hzK1lnrH0k7SGpEzgeaHgv3aGMtW6StpO0Q678COAPg/fbDIarm6RtU+xIOhxYHxFLU7PeQ5IOSU3iJwL/2ZjoRzbWuqXLAFNS+QTgtTT4dUvP8fnAsoj4Sm7RAqDaY38OT78GC4ATU8/4Q4AH02t2FXCEpEmpF/YRqaxh6lW3VKet0j6nAC8Dlo5LJYaxCXUbTtN9Ttarbq30ObmRevYebLUH2VnUKrIOHCuAdwCnkZ093Q78K08PhLQ9cBnZ9delwEdy+zkmrX8n8PFG16tedSPrmXtTeixp0brNAm4j68zz32TNqtX9dJH9k94J/Ht1m1avG7AdWc//m9Pr9lWgo8H1ejlZ8+nNwOL0OAZ4FtAL3JHqMDmtL+Br6bW5hdydC2SXQZanx8lN8JrVpW7AS9P8TenvO1qwbs9O79uHyDqgriDrWAtN9jlZr7rRpJ+Toz08wp+ZmVnJuNnfzMysZJz8zczMSsbJ38zMrGSc/M3MzErGyd/MzKxknPzNzMxKxsnfzMysZJz8zczMSub/A0lq5S3BDCAYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot annual maxima\n",
    "fig, ax = plt.subplots(figsize=(8, 5))\n",
    "ax.plot(annual_max_df['year'], annual_max_df['height'], 'k.')\n",
    "ax.set_ylabel(\"{} [{}]\".format(ssh.standard_name, 'mm'))\n",
    "ax.set_title('Annual max (values at end of period)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the statistical model\n",
    "def linear_model(df, with_ar=True):\n",
    "    y = df['height']\n",
    "    X = np.c_[\n",
    "        df['year']-1970, \n",
    "        np.cos(2*np.pi*(df['year']-1970)/18.613),\n",
    "        np.sin(2*np.pi*(df['year']-1970)/18.613)\n",
    "    ]\n",
    "    month = np.mod(df['year'], 1) * 12.0\n",
    "    names = ['Constant', 'Trend', 'Nodal U', 'Nodal V']\n",
    "    X = sm.add_constant(X)\n",
    "    if with_ar:\n",
    "        model = sm.GLSAR(y, X, missing='drop', rho=1)\n",
    "    else:\n",
    "        model = sm.OLS(y, X, missing='drop')\n",
    "    fit = model.fit(cov_type='HC0')\n",
    "    return fit, names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>GLSAR Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>         <td>height</td>      <th>  R-squared:         </th> <td>   0.317</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                  <td>GLSAR</td>      <th>  Adj. R-squared:    </th> <td>   0.255</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   5.949</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Fri, 10 May 2019</td> <th>  Prob (F-statistic):</th>  <td>0.00233</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>15:58:08</td>     <th>  Log-Likelihood:    </th> <td> -179.32</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>    37</td>      <th>  AIC:               </th> <td>   366.6</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>    33</td>      <th>  BIC:               </th> <td>   373.1</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     3</td>      <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>         <td>HC0</td>       <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>        <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Constant</th> <td>  -54.8593</td> <td>   17.873</td> <td>   -3.069</td> <td> 0.002</td> <td>  -89.889</td> <td>  -19.829</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Trend</th>    <td>    2.0208</td> <td>    0.522</td> <td>    3.874</td> <td> 0.000</td> <td>    0.998</td> <td>    3.043</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Nodal U</th>  <td>   -1.5794</td> <td>    8.248</td> <td>   -0.191</td> <td> 0.848</td> <td>  -17.745</td> <td>   14.586</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Nodal V</th>  <td>  -19.1110</td> <td>    6.761</td> <td>   -2.827</td> <td> 0.005</td> <td>  -32.361</td> <td>   -5.861</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 1.283</td> <th>  Durbin-Watson:     </th> <td>   1.341</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.526</td> <th>  Jarque-Bera (JB):  </th> <td>   0.876</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.377</td> <th>  Prob(JB):          </th> <td>   0.645</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.967</td> <th>  Cond. No.          </th> <td>    98.6</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors are heteroscedasticity robust (HC0)"
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                           GLSAR Regression Results                           \n",
       "==============================================================================\n",
       "Dep. Variable:                 height   R-squared:                       0.317\n",
       "Model:                          GLSAR   Adj. R-squared:                  0.255\n",
       "Method:                 Least Squares   F-statistic:                     5.949\n",
       "Date:                Fri, 10 May 2019   Prob (F-statistic):            0.00233\n",
       "Time:                        15:58:08   Log-Likelihood:                -179.32\n",
       "No. Observations:                  37   AIC:                             366.6\n",
       "Df Residuals:                      33   BIC:                             373.1\n",
       "Df Model:                           3                                         \n",
       "Covariance Type:                  HC0                                         \n",
       "==============================================================================\n",
       "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Constant     -54.8593     17.873     -3.069      0.002     -89.889     -19.829\n",
       "Trend          2.0208      0.522      3.874      0.000       0.998       3.043\n",
       "Nodal U       -1.5794      8.248     -0.191      0.848     -17.745      14.586\n",
       "Nodal V      -19.1110      6.761     -2.827      0.005     -32.361      -5.861\n",
       "==============================================================================\n",
       "Omnibus:                        1.283   Durbin-Watson:                   1.341\n",
       "Prob(Omnibus):                  0.526   Jarque-Bera (JB):                0.876\n",
       "Skew:                          -0.377   Prob(JB):                        0.645\n",
       "Kurtosis:                       2.967   Cond. No.                         98.6\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors are heteroscedasticity robust (HC0)\n",
       "\"\"\""
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# print the fit properties\n",
    "fit, names = linear_model(annual_mean_df)\n",
    "prediction = fit.get_prediction()\n",
    "fit.summary(xname=names)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnEAAAE/CAYAAAApN5W5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXl8VPXV/z/fmclkX0jY10QQwbAEiOCIaJCKoNaqj4qKFaSC1Fbt409autvW56mtPj4qtXUpgmiUoLhStQoyPCADEiDsOwlhTci+zT7f3x8nNzPZJ2SWeyfn/Xrd18zcuXPvd+72/dxzzvccIaUEwzAMwzAMoy104W4AwzAMwzAM03VYxDEMwzAMw2gQFnEMwzAMwzAahEUcwzAMwzCMBmERxzAMwzAMo0FYxDEMwzAMw2gQFnEMwzBdRAiRI4Q4E+52MAzTs2ERxzCM6hFCmIUQlUKI6HC3hWEYRi2wiGMYRtUIIdIBTAMgAdwW1sYwDMOoCBZxDMOonQcBbAOwEsA83y+EECuFEK8IIf4lhKgVQmwXQgz3+V4KIRYLIY4JIaoalxWN3z0thHjHZ9n0xuUNjZ8fEkIcalzvSSHEI/42uHE9jzZut1YI8SchxHAhxFYhRI0QYo0Qwuiz/K1CiILGNm4VQozz+W6pEOJE43oOCiHu8PluvhBiixDi+UZLZaEQYnZXdi7DMNqFRRzDMGrnQQC5jdNNQoh+Lb6/F8AfAPQCcBzAf7X4/lYAVwEYB+AeADf5ud3Sxt8mAXgIwP8KISZ2od03AZgE4GoAPwfwOoAHAAwBMAbAfQAghJgA4E0AjwBIA/AagE99XMcnQJbI5Mb/+Y4QYoDPdqYAOAKgN4C/AliuCFWGYSIbFnEMw6gWIcS1AIYBWCOl3AkSNPe3WOwjKeV3UkoXSOhltfj+WSlllZSyGMDGNr5vEynlv6SUJySxCcBXIDHlL3+VUtZIKQ8A2A/gKynlSSllNYAvAExoXG4RgNeklNullG4p5VsA7CDxBynl+1LKc1JKj5QyD8AxAJN9tnNKSvmGlNIN4C0AAwC0FLoMw0QgLOIYhlEz80Dip6zx87to4VIFcMHnfQOAhC5+3yZCiNlCiG1CiAohRBWAm0HWLn8p8XlvbeOz0o5hAP5foyu1qnFbQwAMbGzHgz6u1iqQFc+3HU3/T0rZ0PjWr//IMIy2MYS7AQzDMG0hhIgFuT/1QghFqEQDSBFCjJdS7unmJuoBxPl87u+z7WgAa0Gu3E+klE4hxMcAguGmPA3gv6SULd3AEEIMA/AGgBkALFJKtxCiIEjtYBhGY7AljmEYtXI7ADeAK0Eu0CwAowFsBomr7lIA4DohxFAhRDKAX/p8ZwQJxosAXI2DBWYGYJtt8QaAxUKIKYKIF0LcIoRIBBAPGpV7EaDBFiBLHMMwDIs4hmFUyzwAK6SUxVLKC8oE4G8A5iqjSC8VKeXXAPIA7AWwE8A6n+9qATwOYA2ASlAc3qfd2V4H7cgHsBD0vypBgzPmN353EMD/ALCA3LFjAXwbjHYwDKM9hJQy3G1gGIZhGIZhughb4hiGYRiGYTQIiziGYRiGYRgNwiKOYRiGYRhGg7CIYxiGYRiG0SAs4hiGYRiGYTRIxCf77d27t0xPTw/6durr6xEfHx/07WgF3h+t4X3SHN4freF90hzeH63hfdKaSNsnO3fuLJNS9vFn2YgXcenp6cjPzw/6dsxmM3JycoK+Ha3A+6M1vE+aw/ujNbxPmsP7ozW8T1oTaftECHHK32XZncowDMMwDKNBWMQxDMMwDMNoEBZxDMMwDMMwGiTiY+Lawul04syZM7DZbAFbZ3JyMg4dOhSw9WkdLe2PmJgYDB48GFFRUeFuCsMwDMP4TY8UcWfOnEFiYiLS09MhhAjIOmtra5GYmBiQdUUCWtkfUkqUl5fjzJkzyMjICHdzGIZhGMZveqQ71WazIS0tLWACjtEuQgikpaUF1CrLMAzDMKGgR4o4ACzgmCb4XGAYhmG0SI8VceGmqKgIY8aM8Xv5V199FatWrepwmZUrV+KnP/1pm9/993//d5faxzAMwzCMumERpxEWL16MBx988JJ/zyKOYRiGYSILFnFhxO12Y+HChcjMzMTMmTNhtVpx4sQJzJo1C5MmTcK0adNw+PBhAMDTTz+N559/HgCwY8cOjBs3DllZWViyZEkzi965c+cwa9YsXH755fj5z38OAFi6dCmsViuysrIwd+7c0P9RhmEYhmECDos4P7FYLPjzn/8Mi8USsHUeO3YMP/nJT3DgwAGkpKRg7dq1WLRoEZYtW4adO3fi+eefx6OPPtrqdw899BBee+01FBQUQK/XN/uuoKAAeXl52LdvH/Ly8nD69Gk8++yziI2NRUFBAXJzcwPWfoZhGIbpKbjdQG1tuFvRnB6ZYqSrWCwWzJgxAw6HA0ajERs2bIDJZOr2ejMyMpCVlQUAmDRpEoqKirB161bcfffdTcvY7fZmv6mqqkJtbW3T9u+//36sW7eu6fsZM2YgOTkZAHDllVfi1KlTGDJkSLfbyjAMwzA9FZsNOHeOhJyasmexiPMDs9kMh8MBt9sNh8MBs9kcEBEXHR3d9F6v16OkpAQpKSkoKCgI2DpdLle32sgwDMMwPRUpgaoqoKQEiIkB1NalsjvVD3JycmA0GqHX62E0GpGTkxOU7SQlJSEjIwPvv/8+AEpEu2fPnmbLpKSkIDExEdu3bwcArF692q91R0VFwel0BrbBDMMwDBOhuFxkfSspARISAKMx3C1qDYs4PzCZTNiwYQP+9Kc/BcyV2h65ublYvnw5xo8fj8zMTHzyySetllm+fDkWLlyIrKws1NfXN7lPO2LRokUYN24cD2xgGIZhmE6wWoFTp8iNmpQE6FSqltid6icmkymg4i09PR379+9v+vzUU081vf/yyy9bLf/00083vc/MzMTevXsBAM8++yyys7MBAPPnz8f8+fOblvONlfvLX/6Cv/zlL4FqPsMwDMNEHFIClZVAaSkQGwuovaQ2izgN8q9//Qt//vOf4XK5MGzYMKxcuTLcTWIYhmEYTeNyARcuAPX1NHhBC8V8WMRpkDlz5mDOnDnhbgbDMAzDRAQNDRT/JoS6Rp92Bos4hmEYhmF6JFIC5eVAWRkQFwcYNKaKNNZchmEYhmGY7uN0kvu0oUE77tOWsIhjGIZhGKZHUVcHnD9Po0615D5tSdgHzQoh3hRClAoh9vvMSxVCfC2EONb42qtxvhBCvCyEOC6E2CuEmBi+ljMMwzAMoyU8HuDiReDMGSA6mkagapmwizgAKwHMajFvKYANUsrLAWxo/AwAswFc3jgtAvCPELUx4Lz88ssYPXo05s6di08//RTPPvssAODjjz/GwYMHm5ZbuXIlzp0716V1FxUVYcyYMW3Of/fdd7vX8EvYLsMwDMOEG7sdKC6mFCKJidqLf2uLsIs4KeX/AahoMfsHAN5qfP8WgNt95q+SxDYAKUKIAaFpaWD5+9//jq+//hq5ubm47bbbsHQp6dRAiLj26EjEcXkuhmEYJlKpqaHkvVJS9QUtxr+1hVp1aD8p5fnG9xcA9Gt8PwjAaZ/lzjTOOw8NsXjxYpw8eRKzZ8/GggUL0KtXL+Tn5+P+++/Hp59+ik2bNuGZZ57Bfffdh/z8fMydOxexsbGwWCw4ePAgnnzySdTV1aF3795YuXIlBgwYgJ07d2LBggUAgJkzZ7a53aVLl+LQoUPIysrCvHnz0KtXL3z44Yeoq6uD2+3Gpk2b8Nxzz2HNmjWw2+2444478Ic//AFFRUWYPXs2rr32WmzduhWDBg3CJ598gtjYWL+2yzAMwzDhwO2mxL3V1STe1Fp54VJRq4hrQkophRCyK78RQiwCuVvRr18/mM3mZt8nJyejtrY2YG0EALfb7fc6n3vuOXzxxRf47LPPkJaWhtzcXDgcDowdOxazZ8/GrFmzcPvtZHxct24dnnnmGUycOBFWqxWPPvooVq9ejd69e2Pt2rX4+c9/jr///e+YN28enn/+eUydOhW/+c1v4PF4WrXnd7/7HV5++eWm2qy5ubnYuXMntm7ditTU1CYr4IYNGyClxJw5c/Dll19i8ODBOHbsGP75z3/ihRdewLx58/DOO+/g3nvvbXe7XdkfasBms7U6TwJNXV1d0LehJXh/tIb3SXN4f7SG90lr2tsnUgIOB70PlHjzeCifnFpQq4grEUIMkFKeb3SXljbOPwtgiM9ygxvnNUNK+TqA1wEgOztbtixYf+jQISQ2Dkf52c+AgoLuN9jtdkGvp92ZlQW8+GLHywshkJCQgMTERMTExMBoNCIxMRFRUVGIjY1tap9er0d8fDwSExOxf/9+HDp0CHfccUfjNt0YMGAA3G43ampqMGsWhRb+6Ec/woYNG5rWoRAXFweDwdA0PyYmBjNnzsSwYcMAAFu2bMHGjRtx3XXXAaAL4+zZsxg1ahQyMjIwdepUAMCUKVNQUlLS4XZra2tbbV/NxMTEYMKECUHdhtlsRstzsSfD+6M1vE+aw/ujNbxPWtNynyilsy5eBGJiAls6q64OGDkycOvrLmoVcZ8CmAfg2cbXT3zm/1QIsRrAFADVPm7XiEdKiczMTFgslmbzq6qqLnmd8fHxzdb/y1/+Eo888kizZYqKihAdHd30Wa/Xw2q1XvI2GYZhGCYYKLnfrNbIin1rj7CLOCHEewByAPQWQpwB8HuQeFsjhPgRgFMA7mlc/HMANwM4DqABwEPd3X5nFjN/qa21BsTypFix2vp8xRVX4OLFi7BYLDCZTHA6nTh69CgyMzORkpKCLVu24Nprr0Vubq5f627JTTfdhN/+9reYO3cuEhIScPbsWUR18AiTkpLi13YZhmEYJtjU15OrU6cjAdcTCLuIk1Le185XM9pYVgL4SXBbFF7uvfdeLFy4EC+//DI++OADzJ8/H4sXL24a2PDBBx/g8ccfR3V1NVwuF372s58hMzMTK1aswIIFCyCEaHeAwbhx46DX6zF+/HjMnz8fvXr1avb9zJkzcejQIZhMJgBAQkIC3nnnHej1+nbb6892GYbRJlJSYLjLRZPdTjFGHg9ZOHytHL4xR219p9MBRiOldYiKirwAcya8XLxI5bO0WDqrOwjSRZFLdna2zM/Pbzbv0KFDGD16dEC3o7UYsGCjtf0RjHOiJRzL0hzeH60J1z7xeLxCzekEbDaaFMEGkCjT6QC9nt77dh0tu5G2uhUpvesCqKONjqZOtz1xx+dIa3ifNMduBzZtMiM9PQfx8cF3n4YiJk4IsVNKme3Psj1IrzIMwzDKiL2GBprsdhJuCkKQoDIYSGAFq1P0eKgdViu9l9K7bSWTvrKM0RicNjDaxeMBqqrIAgf0HPdpS1jEMQzDRDhSklirq6Okp04nCaaoKJpiYkLfJp2OxFpLfMWd0wkUFlIbk5KA+Hj6DbtiezY2Gw1ecDh6xuCFjmARxzAME4F4PNTZ1dVRolOPh1yh0dHhEW3+4ivulOLkbjdZXSoaa/skJJCoi44ObPoIRt14PHQOlJXRse+p1jdfWMQxDMNECG43WbBqaki8ASTcYmO1bb3S68m1C5BVURGnUlJnnpxM/zE6umdbZSIZmw04f56ss4mJfJwVWMQxDMNoGLebUitUV1OMG0DWqVAEeYcDIZpbEl0usswoMXWJiTTFxJD4Y7SNr/UtJoatby1hEccwDKNBHA6yuFVUeC1SGhoQHjCUQRgAdfiKoFVcscnJ1PlHoqCNdKxWsr65XGx9aw8NG9gZXxIaH0/OnTuHu+66q8NlX3zxRTQoj+wAbr755m5VfVB49dVXsWrVqg6XKSgowOeff97tbXVGUVERxowZE/TtMEyoUdxKhYUUJxYfTx0cj+Ak4RYbS/sjLo4EXXEx7auKCm8dTUbduN006vTUKbKm9vTBCx3BljjQRW6zdW8ddXWiycwbEwMMHdr9drnd7g4T7bbFwIED8cEHH3S4zIsvvogHHngAcY1BJoESVYsXL+50mYKCAuTn5+Pmm2/2e70ulwuGnpS9kWFaICVZJcrKyGUaFdUzrW5dQQgSdACJgvJyEgbR0UBqKolfdreqj4YGGnnqdrP1zR/YEgcScAkJgZs6E4RFRUUYNWoU5s6di9GjR+Ouu+5qsoylp6fjF7/4BSZOnIj3338fJ06cwKxZszBp0iRMmzYNhw8fBgAUFhbCZDJh7Nix+M1vftNs3YoFyu1246mnnsKYMWMwbtw4LFu2DC+//DLOnTuH6dOnY/r06U3bLCsrAwC88MILGDNmDMaMGYMXG2uSFRUVYfTo0Vi4cCEyMzMxc+bMNmunPv3003j++ecBkHXvF7/4BSZPnoyRI0di8+bNcDgc+N3vfoe8vDxkZWUhLy8P9fX1WLBgASZPnowJEybgk0+oTO7KlStx22234YYbbsCMGTNw77334l//+lfTtubPn48PPvgARUVFmDZtGiZOnIiJEydi69atXT7+DKNWPB5ymRYWAqdP0+ekJK84YfxDr/daLIUgkXDiBJVoamhonoSYCQ8uF1BSQkYV5XixgOscNm+EiSNHjmD58uWYOnUqFixYgL///e946qmnAABpaWnYtWsXAGDGjBl49dVXcfnll2P79u149NFH8c033+CJJ57Aj3/8Yzz44IN45ZVX2tzG66+/jqKiIhQUFMBgMKCiogKpqal44YUXsHHjRvTu3bvZ8jt37sSKFSuwfft2SCkxZcoUXH/99ejVqxeOHTuG9957D2+88QbuuecerF27Fg888ECH/9HlcuG7777D559/jj/84Q9Yv349/vjHPyI/Px9/+9vfAAC/+tWvcMMNN+DNN99EVVUVJk+ejO9973sAgF27dmHv3r1ITU3FRx99hDVr1uCWW26Bw+HAhg0b8I9//ANSSnz99deIiYnBsWPHcN9996FlhQ6G0RouF1BbS9Yjt5tEm5rTgmgJJTeeMsr19GkSDSkpJPLayl3HBA+3m2IYy8q8cYws3vyHLXFhYsiQIZg6dSoA4IEHHsCWLVuavpszZw4AoK6uDlu3bsXdd9+NrKwsPPLIIzh//jwA4Ntvv8V991HZ2R/+8IdtbmP9+vV45JFHmlyRqampHbZpy5YtuOOOOxAfH4+EhATceeed2Lx5MwAgIyMDWVlZAIBJkyahqKio0/945513drr8V199hWeffRZZWVnIycmBzWZDcXExAODGG29savPs2bOxceNG2O12fPHFF7juuusQGxsLp9OJhQsXYuzYsbj77rtx8ODBTtvFMGpFSurMTp4k119MDHVqHE0QeJRRrspI1qoqoKiIppoaEhdM8FAqLhQW0jkfFxfcCiGBID/fgtde+zMsFku4m9IE3xrChGhxpvp+jo+PBwB4PB6kpKSgoKDAr3UEk2ifx1O9Xt+mO7W93+j1erhcrjaXkVJi7dq1uOKKK5rN3759e9N+AICYmBjk5OTg3//+N/Ly8nDvvfcCAP73f/8X/fr1w549e+DxeBDD5gpGg7hcQGUlBd4rgxXU3JlFGjqdNw+d00nuVoBc1zy6NbBISVbmixdJKMfFaSOHYX6+BXPmzIDD4cA//mHEhg0bYDKZwt0stsSFi+Li4iY1/+677+Laa69ttUxSUhIyMjLw/vvvAyDBs2fPHgDA1KlTsXr1agBAbm5um9u48cYb8dprrzUJqIrGdOeJiYmora1ttfy0adPw8ccfo6GhAfX19fjoo48wbdq0bv7T5rTc9k033YRly5ZBNlbM3r17d7u/nTNnDlasWIHNmzdj1qxZAIDq6moMGDAAOp0Ob7/9Ntz8+MxoCI+HXEmFhd60GGq3RkQ6UVEU2xwfT/FyxcVkGVVENnNpSEmjhQsLaXS10Uj7WQsCDgAsFjOcTgc8HjccDgfMZnO4mwSARVzYuOKKK/DKK69g9OjRqKysxI9//OM2l8vNzcXy5csxfvx4ZGZmNgX+v/TSS3jllVcwduxYnD17ts3fPvzwwxg6dCjGjRuH8ePH49133wUALFq0CLNmzWoa2KAwceJEzJ8/H5MnT8aUKVPw8MMPY8KECQH818D06dNx8ODBpoENv/3tb+F0OjFu3DhkZmbit7/9bbu/nTlzJjZt2oTvfe97MDbmU3j00Ufx1ltvYfz48Th8+HAz6x3DqJn6enLdXbhAMW+KJYhRB77uVqOR4hMLC0nU1dXxYIiu0NBA6ULOnKH4w8RE7Y0MNplyEBVlhF6vh9FoRE5OTribBAAQigUkUsnOzpYtA90PHTqE0aNHN30OTIqRuqZcbZ2lGCkqKsKtt96K/fv3d2+jKqa2thaJGsqB0PKcCAZms1k1F74a6Kn7w24HSktJxMXGNq/9eeCAGZmZOWFrm9pQ4/5wOGgSglytSUmhHXSipevGZiO3aX097aNg5TIM1XmSn2/Bpk1m3HVXTlBdqUKInVLKbH+W5Zg4BCanW22t5LxNDMO0i8tF1pyqKhJuSUnhbhFzKRiNNHk8FNtVWUnHs1cvsqby6FZ6UCkvp/0TSed6drYJo0aZMHJkuFvihUVcGEhPT49oKxzDMF6UuLeLFyn+h7PPRwZKdQjAW7/V4yHRkpJCMXVGY8851i4XWdwqK0nEGQyckDoUsIhjGIYJAkogd2kpdXBaGYXHdB3f+q2+1SGiosjlGh9PFrpIE3RuN1USqaqic12InlvDN1ywiGMYhgkwNhuJt4YGEm+c+abnoFQbAEjkVFaSlS5SBJ1SAq6mhiYpyeLIwi08sIhjGIYJEB4PFVovL6eOLVJigZhLQ6/3jjr2FXQGAwm6uDhvBQk1IyW5SOvq6D94PPQfOJ9h+GERxzAMEwCUwt0uF8e9Ma3xFXRKtYLycjpPlPyAShyd0Rj+FBxuNyU+tlpJuDmd1KbYWA4LUBN8KDTOkiVLkJmZiSVLluDVV1/FqlWrWi1TVFSEMWPGhKF1wDXXXNPpMi+++CIaGhqC3pb58+fjgw8+CPp2mJ6F2928cDcLOKYzFNGWmEjnS0yMN/VMcTFw4gQlGL5wgVyWNltw8tJJSQ8dNhuNJC0ro1qyx49TG4qLva7gxESO61QjbIkD8OGHdBPuDna7oWloeb9+QGPZ0KDz+uuvo6KiAvpwP7a1w9atWztd5sUXX8QDDzyAuC5kO3W73ar9z0zPobbWe+9g1ylzqeh0FCfnm57E7Sbrbm0tiS2AxJSSqkan805CtD0BXkufItgcDhKMViu9900Vq9d7rW38IKINWMSBbsKDB3dvHVarbBpufuZM58uvWrUKzz//PIQQGDduHN5++20UFRVhwYIFKCsrQ58+fbBixQoMHToU8+fPR1JSEvLz83HhwgX89a9/xV133YXbbrsNdXV1mDRpEn75y1/i0KFDSEhIwFNPPYWdO3diwYIFAKjSgYLb7cbSpUthNptht9vxk5/8BI888gjMZjOefvpp9O7dG/v378ekSZPwzjvvQAiBHTt24IknnkB9fT2io6OxYcMGxMXFtbmeliQkJKCurq7d9S9btgznzp3D9OnT0bt3b2zcuBFfffUVfv/738Nut2P48OFYsWIFEhISkJ6ejjlz5uDrr7/GPffcgw8//BDfffcdALI2fv/738e+ffvwxz/+EZ999hmsViuuueYavPbaayGtM8tEPk4nWU1qa8k6wQXqmUCjCCpfXC4Sd1VVJL6U6VLWbTBwibdIgA2jYeDAgQN45pln8M0332DPnj146aWXAACPPfYY5s2bh71792Lu3Ll4/PHHm35z/vx5bNmyBevWrcPSpUsBAJ9++iliY2NRUFCAOXPmNNvGQw89hGXLljXVWlVYvnw5kpOTsWPHDuzYsQNvvPEGCgsLAVDd0hdffBEHDx7EyZMn8e2338LhcGDOnDl46aWXsGfPHqxfvx6xsbEdrqc92lr/448/joEDB2Ljxo3YuHEjysrK8Mwzz2D9+vXYtWsXsrOz8cILLzStIy0tDbt27cLSpUvhcDiatpmXl9e0D376059ix44d2L9/P6xWK9atW3cph4lhWiElubcKC8kFlZTEAo4JHQaDNz+dEkOXkND1SakUwgJO+7CICwPffPMN7r77bvTu3RsAkJqaCgCwWCy4//77AQA//OEPsWXLlqbf3H777dDpdLjyyitR0onvt6qqClVVVbjuuuua1qXw1VdfYdWqVcjKysKUKVNQXl6OY8eOAQAmT56MwYMHQ6fTISsrC0VFRThy5AgGDBiAq666CgCQlJQEg8HQ4Xrao631t2Tbtm04ePAgpk6diqysLLz11ls4depU0/e+YvWee+5BXl4egOYibuPGjZgyZQrGjh2Lb775BgcOHOiwXQzjDw4HWdnPnaNOULG8MwzDhAt+htQI0T7BEt2pdyulxLJly3DTTTc1m282m5ttQ6/Xw+VydXk9HeHP+qWUuPHGG/Hee++1uQ7fAvdz5szB3XffjTvvvBNCCFx++eWw2Wx49NFHkZ+fjyFDhuDpp5+GrbuFcZkejZTkvrp4kSwhHPvGMIxaYEtcGLjhhhvw/vvvo7y8HABQUVEBgEZyrl69GgCQm5uLadOmXdL6U1JSkJKS0mTJy83Nbfrupptuwj/+8Q84nU4AwNGjR1FfX9/uuq644gqcP38eO3bsAECF7V0uV5fX0xGJiYmora0FAFx99dX49ttvcfz4cQBAfX09jh492ubvhg8fDr1ejz/96U9NVjhFsPXu3Rt1dXU8GpXpFnY7cOoUxb9x0l6GYdQGW+LCQGZmJn7961/j+uuvh16vx4QJE7By5UosW7YMDz30EJ577rmmgQ2XyooVK7BgwQIIIZoNbHj44YdRVFSEiRMnQkqJPn364OOPP253PUajEXl5eXjsscdgtVoRGxuL9evXd3k9HbFo0SLMmjWrKTZu5cqVuO+++2C32wEAzzzzDEa2U3F4zpw5WLJkSVNsXEpKChYuXIgxY8agf//+TW5ghukKSh6v0lIuI8QwjHoR3XHNaYHs7GyZn5/fbN6hQ4cwevTops+BSTFiRXQ0BcmEMsWIWqmtrUWihnq+ludEMDCbzcjJyQnqNrSEWveH3Q6cP0+voc75duCAGZmZOaHboMrh/dEa3ietCeU4gsZ2AAAgAElEQVQ+qasD2rEpBAwhxE4pZbY/y7IlDoERXLW1Ln5aZxgNo8S+lZZyLUiGYbQBiziGYXo8DgdZ32w2StvAWekZhtECfKtiGKbHIiXVhSwspCSqiYks4BiGaZuzZ4FNm8Ldiub0WEuclJKz+DMAupeyhdEuDgfFwjY0sPWNYZj28XiAAweA/HyqmqEmeuRtKyYmBuXl5dx5M5BSory8HDGcO6LHICVQXU3WN6eTrW8Mw7SPzQZs3kwCrl+/cLemNT3SEjd48GCcOXMGFy9eDNg6bTYbCwEftLQ/YmJiMLi7xXMZTeB0AhcusPWNYZjOqawEzGa6XwwcGO7WtE2PFHFRUVHIyMgI6DrNZjMmTJgQ0HVqGd4fjJpQap6WlFDxbx55yjBMRxQWAlu2UHm9Pn1onhqddz1SxDEM03NQYt/q6ynvG1vfGIZpD7cb2L0b2LcP6N2bkn2rGRZxDMNEJIr17cIFrnnKMEznNDSQ9e38eWDAAG088LGIYxgm4rDbyfpmtXLsG8MwzcnPt8BiMcNkykF2tgkAcPEi8M03ZIkbMCDMDewCLOIYhokYfKsuREVx7BvDMM3Jz7dgzpwZcDodiIoyIi9vA5KSTLBYyFofHx/uFnaNDkWcECLVj3V4pJRVAWoPE2bcbsqD43SSNcNmo88ej3cZ3+DOlvOV7+x24ORJ6khjYqiMkcFAk15PE6fpYwKJzUauU7udrW8Mw7SNxWKG0+mA2+0G4EBenhlXXmlC377UX2mNzixx5xqnjrpbPYChAWsRE3SkbC7WbDaa7PbmokynI9GlvCq0FF++n5X3Oh0JN5eL4pI8Hq/Ak9K7TkXgKSJPr6cLiQUe4y8eD1BRAZSVURAyW98YhmkPkykHUVFGAA7odEYkJ+doJv6tLToTcYeklB3miRBC7A5ge5gg4HaTSKuvpxghh6O5NU2xjMXGBlY8KUKuLaSkztdmo2BSt7v57xISaDIaWdQx7WO1UhCyy0Xijc8ThmE6IjvbhL/9bQPWrjXj8stzMGGCKdxN6hadiTh//p2290CE4nRSB1dTQyJJSq87My4u/J2dEF7x2BKPh9pcXU3L+Yq66GhtmryZwOJ2k/WtvJysuQkJ4W4Ro3XaCnZnIguPBzh4ELhwwYTvf9+E2Nhwt6j7dCjipJQ25b0QoheAIb6/kVLu8l2GCR9SkoWtoYECux0OEkBGo/Y6OJ2OOmal4IMi6mpq6H9GRVHMU3w8i7qeSEMDxb6x9Y0JFG0Fu7OQiyzsdmDbNkri269f8xAhLePX3xBC/AnAfAAnACiOOAnghuA0i/EHxR1ZV0cCx+0mARRpcUGKqFPweOg/V1d7RV1CAv3nmBju1CMVh4Pi3mpqyPWvkapujAZoGexusZhZxEUQSvms+npKHxJJfYS/WvQeAMOllI5gNobpHF+rVF2d100aE6PdwMyuotOhmRnc4wFqa+lC1euB5GQSdSzoIgO3m45teTkn7WWCg2+we1SUESZTTribxASI4mLg//6PjBt9+4a7NYHHXxG3H0AKgNIgtqUVQogiALUA3ABcUsrsxrQneQDSARQBuEdKWRnKdoUDp5OESkUFiRbFpcgipbmo83hI4FZUkKBLSfHG0vG+0hZS0jlfWkrHNSGBjyETHLKzTcjL28AxcRGE2w3s3UsltPr0UX/5rEvFXxH3ZwC7hRD7AdiVmVLK24LSquZMl1KW+XxeCmCDlPJZIcTSxs+/CEE7woLNRjFuNTXUgcXG9hyL26XQUtBVV3stOIqFjgWd+rFaqeKCzUYPK20NgGGYQJKdbWLxFiFYrcC33wJnzwIDB0Z2n+mviHsLwF8A7APg6WTZYPMDADmN798CYEaEiTjFZVpWRsGYBgNb3S6FloKuqooFndpxOum8r66mY8OuU4ZhukJ5ObBxI8XQaql81qXir4hrkFK+HNSWtI0E8JUQQgJ4TUr5OoB+Usrzjd9fANAvDO0KCi6X12XqckXeAIVwotNRahXAK+gqKkjQKS7X9nLaMcFHOSYXL3LcG8Mwl8aJE2SBS0gAevcOd2tCg5C+WV/bW0iIF0Bu1E/R3J26K3hNA4QQg6SUZ4UQfQF8DeAxAJ9KKVN8lqmUUvZq8btFABYBQL9+/SatXr06mM0EANTV1SHhEnN5KBUUlIS3kWD6tdnqEBOjjdwmSpUKIbwVKoJBd86RSETZHx4PWeCAyDj3u4OWrptQwPujNbxPWmO11sHjSYDNRvfwYHpXQmHhmz59+k4pZbY/y/priVOqNlztMy/oKUaklGcbX0uFEB8BmAygRAgxQEp5XggxAG0Mtmi02L0OANnZ2TInJyeYzQQAmM1mdGU7UpLLtLycXpUao5Hi2jtwwIzMzJxwN6NLKPVipSQLaHJyYGMQu3qORDobN5px2WU5sFppP0dK3qbuoMXrJpjw/mgN75PmbNliwSef5CIpaQGuusoU1AdBKWm06333BW8bXcWv26aUcnqwG9ISIUQ8AJ2Usrbx/UwAfwRZA+cBeLbx9ZNQt627WK004s5qJeHGriN1EBVFk5QUUF9bSwH1SUk0cfxc91EeXi5epCdat5tDBhiGuTS++MKCH/94BtxuOwyGXKSlbcCIET1rcIq/yX5TADwISuvhW7Hh8eA0CwDFun0kqNc0AHhXSvmlEGIHgDVCiB8BOAXKYacJbDYK2q6vp/grFm/qRAhvxQjfHHRRUUCvXjTIhOPnuobbTXkNy8q88Z5KYmomspGSrNweD103Pd1lznQfl4vSh+TlmeFyOSClBy6XA4cPm1nEtcPnALYhhKNTpZQnAYxvY345gBmhaEOgcDjIbVpdTTcxtjxoB98Rrm43iZDSUhIfvXrRdyzo2sfppPO+spI689hYRES9wp6My0X3NIeDjq/yareTldVqpamhgR5cbTY69lLS9RQVRYOMYmPpgUh5jY6ma8l3YpiW1NYCmzfTvTg7OwebNhnhctlhMBgxalROuJsXcvwVcTFSyieD2pIIxOmkzquykkfcRQJ6PXU2AB1bJQktC7rW2Gx03tfUeIUwW2C0h5JkvLqazvcLF+i9ElagiDOA5kVF0XViMNBrXBzd93zDEFwueiCqraURyU4nzVPwXXd6OvDpp8DQoVTvslcvLrfWkykuJgFnMAD9+wP9+5uwZMkGFBa+iYyMBUG3wlksdB2oCX9F3NtCiIUA1qH56NSKoLRK4/iWCdLrOdN8JKLEzwHUAfkKuuTknuly9R2sY7XSjZatztrB5SLRXVNDMYvnz5PIUoiOJjHer1/37mcGA03+uNKVGNUDB8h9JiWlBBoyhDpx5eGJiWycTmDXLuDgQUod4nvujBhhwoQJdtTWBk/ANTQAb78NbNsGDB9O56Fa+nR/RZwDwHMAfg0alYrG18uC0SitolQIKGusL8HirWegdEoAdYSKy9VopA4n0gWdy0VxnuXldLPl/IbqR0q6VykWtpISr8sb8Lo8uyvYAoGvC14ZdHT4MLB/P31OTvaKutRUFnWRRlUV1T6trqbUHqG26B87Brz+OuUVveMOYOzY8F8Tvvgr4v4fgBEtyl8xjUhJAq6wkKxwcXHsOuqptCXoLl6kTtHtJguV0ajtMlIeD8VBNTSQ1cZubz4YhFEndjudj+fOAUVFdC4CXsHWt6+6Oqe2UEoP+go1mw04coSsdQA9QAwZQi7Y3r35XqxVpAROngS2bqX7Sr8Qp/V3u4F168idn5YG/PKXZIUrLg5tOzrDXxF3HEBDMBuiZZQgX613zkxgaSnoXC7vDSA6miy1sbH0Xu050pxO6vTr6miSks51HqijXpQqGKWlJNpKS+m4KSPjk5PD3cLA0PLhwW4Hjh4l11tcHDB6NIk6jknWDnY7kJ9Px7FPn9B7MsrKyPp27BhgMgE//CHdq/2ojRBy/O066gEUCCE2onlMXDBTjGgKIVjAMe2jVIJQBI/LRe6Bigrv9/HxXtdrVFR4rSIeD1k4FGub0+kNXOc6vurFaiW39unT9MBgs9GxSkhQh2s0FERHe2Om7HaKpcrPp/9/xRVUEJ1T26iXigpg0yYK0Rg4MPTn7PbtwFtvkWBbtIhEnJrxV8R93DgxDBMAfK10AImm+npvILlS61Wp6arTNZ8CcWNTwgDcbu+rEtxutdL3BgNtn92k6kRKOmfOnydr28WLNF+JS+zVq8OfRzzR0STepKTra/Nmun4uu4ymvn3Z3aoWduywYN06MwyGHIwebULfvqHdvtUK5OZS7dXhw0nAhboNl4K/FRveCnZDGKYno9M1dwtJSS56xQWmoIyKUlI4KCJLEYU6Hc3X6UiYKZOSxsHhoFcl+WpbGI0kHhn1UlMDnDlDsWA1Nd5R8P379wxrW1dRrJEJCfSwUlxMrrKYGGDUKGDYMBqExIQHs9mC+fNnwOVywGAw4uc/34CkpNCZwE6eBF57jR6CbruNJq141joUcUKI16WUi7q7DMMwXUOI5m6hlijizOHwJlNVRJki9HyHwSvWO0XkxcSwBUJrNDTQoIQjRyhmR6ejOK9gF+OONPR6GsUK0PWzbx9QUECDIEaNovi5SB5NriacThpp/NZbSuUFN9zu0FVe8HiAzz8HPv6YRPzSpcDIkUHfbEDpzBJ3uxDC1sH3AkDI66oyTHfJz7fAYjHDZMpBdrbKgx7aQHGrqn1ABNM9bDZylR49Sol2hSA3KQu3wGA0el1mdXXAli10TV15JbnU/BkMofV7Sbi4cIGS59bWAhMn5sBsNjZZ4kJReaG8HHjjDXoomjwZmDePQli0RmddwBI/1rE5EA1hmFCRn2/BnDkz4HQ6EBVlRF7eBr75MqpBcaMfO0YDFAAaTMKu0uCiuFtdLm9y4aFDaXRre7FzfC/pOlYrsHs3iafk5OaVFw4fNmPUqJygW+F27KDBCy4X8KMfAVOnavfa6lDEcSwcE4lYLGY4nQ643W4ADlgsZr7xMgHjUiwzUpKLtK4OWLPGm2/Sn8D748ctIev8egIGA+13KSlG6tQpEhvjxrV2tfK9xH88Hhp8s307nd8tE/eOGGEK+vnb0ECDF7ZupYEtixaFPv9coGFnDNPjMJlyEBVlBEBPzyZTTribxEQIXbXMOJ3A2bMUl1VeTi68tDT/g6qPH7fguee821uyZAMLuQAhBMVJpaRQ5//ttyTwRo8GRowgtzbfS/yjqgr47js611uWzQoVR46Q+7SyEvjBD4Bbb42McJQI+AsM0zWys03Iy9vAcSxMwPHXMlNfTyPi9u8n92lSEuXEUkYd+8vhw7Q9Kd1wuUIXEN5VtG4tjIujqaWr9cor+V7SEcrAhV27aDDVoEHhacNHHwFffkmJg3/1K3pYihT8EnFCiLullO93No9htEJ2tolvuEzA6cgyIyVZ244cAY4fJ1dSaiolUL5URo2i7YUyILyrRJK1UHG1ejwUt1hcDCQlmXDzzaawCBQ1U1JCbsuaGhJP4bB6nT1LqUNOnwauvx64997Iy3np7279JYCWgq2teQzDMD2Wtqy8LpfXZVpW5k1AG4gULyNGhDYg/FLQirWwK+h03kTKiqtVpyM364gR5BLXaqB8d6moAA4dooE54UqB4/EA69cD779P5bKeeALIygp9O0JBZ3niZgO4GcAgIcTLPl8lAXAFs2EMwzBaRLHyNjSQ623fPir/lJhILtNAE4qAcAWPh1yKXSkLpwVrYXdQXK1uN1BYSO7DlBRgzBhg8ODIs/y0hcdDlrd9+yiXodFIo07DkYuyogJYvpxq52ZlAQ89FNl1czuzxJ0DkA/gNgA7febXAvjPYDWqJ8C5hRgmMqmtJfF29Ch9Tk31JpfVAi4XWQwvXqSOubSUppISmu9qfHw3GmmKjva+b3ueCdnZR2Gz7cGECYORnj4+vH8wSOj1ZIEDyDq3dSsJ3eHDgcsvp4D+SLPOOZ3kUt67l9ymcXFkeQvX/9y+HVi1is7R+fOB666LvH3eks5SjOwBsEcI8a6U0hmiNkU8nFuIYSKP+np6+j94kCxVaq7L6fFQstWWIq20lOL2fEuyRUfTfxk8GJgwgTpqp5Osiw5H88lupzxg1dXN5zkcg+F2D8auXcDKlRTgPmwYkJ5Or4MHR1ZResU65/FQihLFtahY52Jjw93C7tHQAJw44R2Yk5IS3gTU9fUGvPYasG0bpQ5ZuJAsgT0Bf2PiJgshngYwrPE3AoCUUl4WrIZFMpxbqOeRn2/Bxx/n4vbbo/lYawh/LOZWKwm3AwfIGhOoeLdA43bToIqdO2m0YFWV9zslJ91llwFXX03v+/al/5KU1H1rhu9AgFOnaNq9mwrSA7T+AQNI0CnT0KHkgtYyOp3XOme1UoUCoLl1To3nSntUVJCF+ehROmbdHZgTCA4fBpYvz0ZFBaUO+f73tVP3NBD4K+KWg9ynOwG4g9ecngHnFupZKJZXh8OO997LZcurRujMYm6zkSjau5c6tN691Zd3yukkgblzJ4mmujpycY4dC4wfTxaxvn2pUkEw0emUzPxU4gig0boVFV5Rd+oUdciK0AGA/v0nY8wYsgCOHKm+/dsVYmNp8nhotOSJE/R/Bg+mJMJpaSRa1eb+8413O39ePVbm+nrggw+ATZuA/v09+PWv6SGkp+HvJVEtpfwiqC3pQXCesp6FYnn1eDxwOtnyqhXas5jb7ZQipKCAOri0NHWJC7udOtydO4E9e8gCFBtLom3SJBJwanBdCkH7Li0NmDjRO7+62muxKypqwKZNcVi/nqyF48aRoBs7VrsuSSW1DECxWyUlVMkAoP84dCiJ69TU8P1HxSVeWUnCurpaPaXfpCS36erVFH96443Agw/mw+m8LrwNCxOdjU5VLq2NQojnAHwIwK58L6XcFcS2RTScp6znoFhepbSz5VVDtLSYX3VVDo4cIVek00niI9yuJIWGBhJs+fneOKWEBCA7m4TblVeqp62dkZxMIm3sWCAxcT/KynJw4ABZEgsKqAPX66lywoQJNCnpPrSGwUD/NzmZPjscZKE7fJjESq9e5Fru3z94rkubjQYlVFV54yQbGryWNiURtRooKQHefptCFzIygCefpP0TE+OBs4dG7Xf2/Pg/LT5n+7yXAG4IbHMYJvJQLK8ff/wmbr99AYt3jaAct2+/NWPYsBycOmXC0aPUmfrWzwwXUpI798svSbi53RRgfu21JN5GjoyM2KDoaLLUTZxIls/jx0lI795NHfrbb9MACUXQDR4cfmvRpWI0emPoALKI7dtHAl2nIzdmWhpZ7pTqHr6vLee1xOEg61VVFQmiCxfIxa4QG0vWQEVUqgWnE/jiC+Czz0jIPvAAMH16+F26aqCz0anTQ9UQholksrNNiI21IzOTBZxW8HiAPn1MGDTIhOpqEm9qsPhISZaITz/1jnq88UayuF12WWR3bDodidORI4E5cygn2e7dNH30EU19+pDgu+YaijXTqqADvHF0AJ2PdXWUd+7Agfb/l5T0KgQJHqORXt1usrgJQctER5NgU/sozsOHKW3I+fPAVVcB992njutQLfhbduvJNmZXA9gppSwIbJMYhmHCS3k5sGMHdRypqeFNn6AgJVlkPvuM6q726gXMnUu5sNRgGQw1QlDs2KBBVMy8uprcrbt3U7b+f/+brHLXXAOYTGSl1DI6HbnIo6L8E15SknDzeOhVGTmtFVFbWwusWQNs2ULC/MknycXONMffcNzsxumzxs+3AtgLYLEQ4n0p5V+D0TiGYZhQYrXSaNNDhyiQWw31MD0ech9+9hkF/PfuDcybB0ydqp04t1CQnEz1Ma+/nixW331H5bDWrKHyS2PGkKCbMEEdAzuCjRDqGnDjL1KScFuzhq7Hm28GbrutZxyzS8HfQzwYwEQpZR0ACCF+D+BfAK4DpR1hEccwjGZxuynWKj+fOpFwlQzyxeMha+Bnn1Ht1X79gB/9iPK4abFz7irHj1tQWJiLjIzoLpcVS0gAbriBpvPnqXqCxULF0GNiyC03dSrlagv3cWa8nDtHrtMjR+jYPPggWVPVgsulPkumv7eCvvAZlQrACaCflNIqhLC38xuGYRjVc/48leupriYrV7itW243jcBct44CzwcOBBYtovxqkTBQwR+OH7fguedmwOWyw2DIxZIlGy65PuyAAcB//Adwxx0kDrZuJXG8eTMdb8XdqvbYsEimro4GLvz73ySy588Hpk1Tl8CuqKABFlddFe6WNMdfEZcLYLsQ4pPGz98H8K4QIh7AwaC0jGEYJojU1pKb8uRJcsWFO+7N5QLWrx+AtWupbumQIcCjj9KABTV1ZqHg8GHK0SelBy6XA4cPmy9ZxCnodJSWZPRoGt24axe5Wz/7jAaJDB9O1rmrrgp+8mOGqKkh4fbNN5Tq5JpraMCKmgrW2+1UM3jIEGDKFI1a4qSUfxJCfAFgauOsxVLK/Mb3c4PSMoZhmCDgcFDM2549ZHUbODD8N+Z9+4B33gFKS69AejqNwMvKCn+7wsWoUZSjjyxxRowalRPQ9UdHk/XNZKKEthYLCbpVq4B336V9f801FEjfE1zXoaamhlLjfPMNXY+TJ9PgFDW5TqWkAU4AxVlmZND16JuSRQ10luw3SUpZI4RIBXCycVK+S5VSVgS7gQzDMIFAShoYsH07PfWroUxWZSXw3nvk3uvfH/jVr/ZixIhxPVa8KYwYYcKSJRtQWPgmMjIWdNsK1xG9elHw/OzZVCVi61ZyZ+fnk0Xu6qtJ0KWn91xRHSiqq8ltunEjuSanTKFap2pJJqxgtdK1mZFBltm4uHC3qH06u4W9CxqJuhOU3Fe0eO2BlcoYhtEaFRUklM6do5Qh4U434fGQFWLtWnKj3nEHiYjU1ArU1oa3bWphxAgTJkywo7Y2NLkVhSChlp4O3HMPJVDeuhUwmyllycCBJOauvrp5Ql6mcyorgc8/pzqnbjftw1tvDX8IQ0s8HnKdRkUBM2ZoI3F0Z8l+b218zQhNcxiGYQKH3U6d8f79lDRVDSlDCgvJbVdUBGRmAj/8IY08ZQLD8eMWHD5sxqhROZdswTMYyKWalUWF1nfsIEH3wQckvEePJkE3aRIF4jNtU15O4u3//o8E0tSpwC23qPN8b2ighz2lnJtWjqu/yX4FKPYtozE+biiA/lLK74LaOoZhmEtAcZ1u20YxN337hn9kZ0MD8OGHZIFLSgIWL6ZYILU/6WsJZVSr00n1brszqlUhPh7IyaGptJTE3NatwD//SWJ80iQqc5aZybnMFMrKSLxt3kzX4rXXknjr0yfcLWuNx0MDiWJjyRquNutgZ/gbEfJ3AB5QrdQ/AagFsBaAygbbMgzT06mqokSvZ8+qo1SWlGTJee89igm64QbgzjvVHWejZPt3u8ndq0zKPKORBEtMTPjjCn3xjmp1B2xUqy99+wK33w784AdU8mzrVjrXLBZywV15JVlxxo8Pv8s+1NhslCg7P59G/gJUTeTmmyn+VI3U1tIgi7FjgXHjtFn5xN/Lb4qUcqIQYjcASCkrhRAa/LsM0z3y8y2wWMwwmXK4kL3KcDiAgwdp1GlMjDpcp6WlVKB9/35g2DDg8ccpWDrcHDliwZ49ZgwbloP09ObnsVJ7MzqarBMJCd4annFx1NFVV5OrrKKCOm+djn4nBO37mBgSNaFOjeId1eoIyqhWBSG8NVwfeAA4etRbw3XPHlrmsstI0GVl0bkYiRZXq5X+b34+jbB2OMjKfMMNwE03qTd20OUi61tyMsXmqdFC6C/+ijinEEIPGswAIUQfkGWOYXoM+fkWzJnjddXk5W1gIacCpCSr27Zt5LLs0yf81iGnk0bhrVtHbtz776dA6XDme/N4yPJw+LAFy5fPgNtN5/Err2zAxImmZsXSuyLA7HaKG2toIKtGWRmJu6oqryAESBQmJgb32CijWrsbE9cVDAaywF15JR3nM2dIzBUUUPzc2rV0TioxdiNHhv/87A4NDfTfduyghxOXi6yO06bRSE41V8GQkgZZOBzkBh81StvHAvBfxL0M4CMAfYUQ/wXgLgC/CVqrGEaFWCzkqnG73QAcsFjMLOLCTE0NWQGKi6kjUUPW/aNHgZUrqRLEVVdRzrdwunQVcSWlMtrODLfbAY+HXI7Hj5sxe/aln8fR0TSlpjaf73bTtuvrabpwgVJ4uFzUyScnByeGbMQIU0jEW1sIQUlhhwyhep+VlWSpKiigtBpff00WzXHjaMrIoCB/tYoehfp6Eqb5+STc3G46p2+4geIBhw9X/3+wWunhYtgwarOaEgp3B3+T/eYKIXYCmAFKL3K7lPJQUFvGBBS3m9wevnEtLhc9nSvznE56QlEmp5Mmj4esCXo9PaEr7/V6eopR5ut0NOn19LvSUm/cjNGofXeCyUSuGoAsGCZTTrib1GNxOqmE0q5ddG4NGBD+88vlAj76iCxwaWnAf/4nddThwG4nl6fbTeLqmmsoRUZcHJCcnIMVK4xNFuVgncd6PVneEhPp8/DhlFqivJxSvZw4QSIHoGXi4sJ/DANNr17eQRF2O3DggNflun07LWM0kugbOtQ7DR4c3vgsqxUoKSHRvXMnhSm43XRef+979HCSkaF+4QZQu8vKqB+68cbIc213xZB4DECN8hshxFApZXFQWsVcEnY7Pflard5khdXV5Naor2//xPV1eej1zcWYcpFK6Z08HppazlNiYgDKtbRvn/e3Oh3F1iQleafYWG/8THS0+oNKs7NNyMvbwDFxYURKsnBt20aZ09WQsBcgd+7rr5NF8LrrgHvvpfM7lLhcdL07HCSIxo8ncdAywD6c57FeT4MD+val9lVXk1g4cYIsdULQfktM1IZA6ArR0cDEiTR5POR2LS72Ttu2kbUOoP8+YAAdv2HDvOIukOXAHA560L5wgY6BMl24QJZbhT59gJkzyXqlVC3QCpWVZLwYP57c3WrvYy4Ff1OMPAbg9wBKALjhTf6ektQAACAASURBVPYbpufMnovLRRdYfT11YpWVJNKqq8k64UtUlNcSlpAQ2osvKqp5LiCPh9pXWUk3CqfT2x5FRBoMJO769aNJEXtquplnZ5tYvIWJqiqyCqjJderxUCLY998n8fHYY9RJh5KGBrr+9XqKR8rIIHHb0XWjhvNYCDqOKSnAFVfQg+fFi5RH7/RpsqDExtK9IioqrE0NODqdV5gpSEn/31fYHTlC4k4hNRXo1y8LOl3z+EXf9219joqie6+vYKuoaP4An5xM993x4+m1f38SkmqwcncVm40svgMHUhqfcI9QDyb+PsM+AeAKKWV5MBvDtMbppM6rvJxubCUl1HEoF5ViwUpOVodFoj10Om/sTHu43fR0ePw41bZUftenD91I0tLof8bHa++mwlw6Viu5oQ4coPNHDbVOAeoEly8nV9P48cBDD9H5GSrcbur04+MpNql/f22LndhYr7BxOskFduYMWTntdnoYTU5W10NdIBHCa6XMzvbOr6mhe39xMbk36+okbDaar4S8KKEwTiedF+0RF0fnyciRXqGmPDSH2nIcDJScb1FR5MLuCaXS/O32TwOoDmZDGEKxVpWV0YVbWup1R8bFdf6ErWX0em8qAwWPh6yOe/d6xati5RswgJ5Mk5K0k12b8R+3m9xsO3bQNaCGhL0K27dTole3G5g/n1yooewsqqvJAhepbqKoKLq+KyqAu+6i++CxYyRiAPIsBNK1qGaSkiiRcGYmfU5M3IPa2px2l1e8HoqoczjIg5OcHHqPTCiprqa+IjOT8r71lD6hQxEnhHiy8e1JAGYhxL8A2JXvpZQvBLFtmmHbNgtWr86F3R7dZReF3U6WtosXSbSVlXljy+LjyQoVqaLNH5RYOt8btstFlskzZ7wxecnJFDitCLuevM+0jhL3tn07WRvS0tQjUurrgXfeIRfX8OHAwoWhLSHkcNA9ol8/SlnSckRoJGIwkPV14ECyyp47Bxw+TOeIXh+8Ua5axR+vRyShWGnj4sgirdbcdMGiM0tc47giFDdOxsaJacRisWD27Bmw2+1YvTrXr9xhdjvdgI4c8Vra9Ho6CXu6aPMHJXbOd4i4knRy9276Pj2dAoLT0nrOE1kkUFlJI05Pn6bjq6YSOAcPUqml6moqWH/LLaGzDEpJ4g2g+pNaSOkQDGJj6b8PH07HobiYBF1FBQn9lBT1WGuZ4GK30/3CaKQRp9de2zOviQ5FnJTyD/6sRAixTEr5WGCapC3MZjMcDgc8Hg+czvZzh7nd3lFYp06RyTs+nlxEXTFvB6K4cyTi64Z1uejmfuwY7dv+/Sl7et++kZMbKNKwWin/1MGDJLrVFEztdFLh86++onPp17+m8ynQHD9uQWFhLjIyoptd2w0NZK3PyKBYqfj4wG9biyQnk9ssM5M8GSdP0v3V7aZ9lJionnOICRwOB3lioqOBKVPoWjx6tGcKOKBrKUY6YmqA1qM5cnJyYDQaYbfbW+VckpKeEIuK6CRzOKiDutS4tmAUd45EDAavm0nJUr91Kx2PxER6ih84kJbhp/bw4nJRx5ufT8dHbYlPi4spdcjZs+Squeee4LiplGvb5bLDYMjFkiUbkJFhQlkZPZwo+a2Y1uh03uD8SZO8Xo7z5+k7JZ0Ro20cDupPo6Io32BGhnrCLMKJisczagOTyYQvvtiAlSvfxB13LEB2tgk1NRSvdfgwCYioKHpq7O7IsWAXd45EdLrmCUdtNrL47NlDYu+yy8j12qcPC7pQogRff/YZxb317q2ukZUeD/Dvf1PJpIQE4MknyeoTLLzXtgculwMFBWbEx5swZgxtlzsr/zAaKYxi2DBKwaTch8+fp+s9JUVd5xnTOU4nWd6ioijJ8PDhfD34okkRJ4SYBeAlAHoA/5RSPhvO9lx9tQkNDXYkJprw+edk2lfKygQypidUxZ0jGSW5MEBWoMJCemo3GoERI+jmH8kjgMONy+WtLakEIKsp7g2gOJt//pNcuxMn0uhT5SEgWHivbTv0eiPGjMnB97/fMwYuBIuEBKqNecUVdEwVj4jdTtbUlBS+ztWMr3hTSnv1lMEaXSFQIi5kkQdCCD2AVwDcCOAMgB1CiE+llAdD1YaWVFfTtGcP3Tj69w9OLEY4ijtHMgaDV0i4XBRDd+gQ3Sguv5zyVfFI18DgcFAsaEEBxXgplmm1xXft3g28+Sa1d9484PrrQxNXNWKECT/+8QaUlLyJqVMX4Ac/MKnSMpyfb9FcxRIh6DpOTaWULKWl5MIvKuL4OTWiiDeDgcWbP3RJxAkh4qSUDW189VKA2uMPkwEcl1KebGzTagA/ABA2EWe3UzxPKDLIh7O4cyRjMJAFDqCbyKFDVDYsLo4SYw4ZQlm/+UbfNWw26jD37qX92qsXCbjjxy349tvWQfzhwuEA8vKAb74h8b54cegshEqm/pEjTZg1y45x48K/P9oiP9+COXO8Mbn+jMRXG3q9twrBVVd54+eUkl/JyRw/Fy6UdFs6HcU2jhjB4s0f/C27dQ2AfwJIADBUCDEewCNSykcBQEq5MmgtbM0gUPJhhTMApoRw+0yEExVFMXIA3Vj27SMLUmIiCbpBg8gVw4KufRoayHW1fz/Fl6WleWOR2griD6eQO30aePVVyj92003Af/xH6OKmlFHrGRkUrH38eGi2eylYLBS353a7AbQ/El8rREdTPGx6euv4OSVPZ3w8W+KDictFrm6nk+6vkyZRnDKnhfIffy1x/wvgJgCfAoCUco8Q4rqgtaqbCCEWAVgEAP369YPZbA7q9pxOIDq6DkZjcLejJfT6OiQmmsPdjG6TmOi10Hk8dMOpqKAneqXkmb9uL5utDgcOmIPWVjXgdpP1zWajjvCyy1qL3cLCXLhcdng8HrhcdhQWvokJE+xtrzCISAl8/vkgvPPOcMTHO/Gb3xxGVlZlSLfvcgFjxpD15/hxdZ8j/fsnwWAwQEoJg8GA/v2Tgt7WUO+Pyy/3Drqx2731qHU6mtTw4BYJ91a321uBp08fuo8aDHRNnDjR9fWF8jzxeOiBTy347U6VUp4Wzc/gDiq0BZWzAIb4fB7cOK8JKeXrAF4HgOzsbJmTkxPUBp0+DWzZYobdHtztaInERHOHpWG0js1GcZAeD7lchw8nl2taWvtP7gcOmJGZmRPSdoaKigoaCHDihDfFS3viNiMjGgZDbqMlLhoZGQtQWxtai05NDQ1e2LeP4qQWLIhGUtJ41NaGbvs2G8XcDfG5m6n5HMnMzEFGxsSQxsSFe3/YbJRkuaiIBkG53d4qEeEaIanFe6uUZO2sq6P747Bh5C7t2zcwNb9DeZ7U1ZFHRi34XTu10aUqhRBRAJ4AcCh4zeqQHQAuF0JkgMTbvQDuD1NbNAMnCQ4svqNcHQ4SMPv3041dEXS9ewfmBqVGpCQRe+ECDQhRkm/6k+dNGaBTWPgmMjIWhPx83LePBJzVCjzwAOV/C6WF5eJFEv633EIxgloiO9ukaRdqV4mJAQYPpunqq+lh5cwZSixcXu5NYRQbqw4rndqwWuk+ocSMT5xI8YjsLg0c/nYxi0GDFwaBhNNXAH4SrEZ1hJTSJYT4KYB/g1KMvCmlPBCOtmgFThIcXIxGeqIEyP1y/DiJOoOBYp2GDfPG2GkZKcmdfO4cxbvV1no7sYEDu7auESNMmDDBHlILnG/lhUGDgCVLqHMOFR4Pid4hQ6h0Fndk2sJgoOu8b19gwgQKwi8poeu9pISWEYIEXVxc5D7AdYTHQ/GwDQ1ktUxKAiZPpuvMt/41Ezj8Os2klGUA5ga5LX4jpfwcwOfhbodW4CTBoSMqqnnaktOnyVKl15P74OhR7whNLSSs9HjI4nD6NLlKrVbvKD615XfriLNngddeo//xve9R5YVQJn2128ktl5VF7lsOltc2QtB13KsX5aKz2chFrjzklJTQQ4OUdKzj40ncRdJxl5L+d0ND89jBvn3J3ThwII/oDwX+jk59C8ATUsqqxs+9APyPlHJBMBvHBAZOEhweDAavu8zjIVG3fTvd/AAa4TpwILkgk5PVk6vK5SLBUVxMbiOHg27OKSnUTi3h8QBff02VF2JigJ/9jERUKKmtpY5u+nQaCclEHkp4Rd++lFxYSqC+nlyJZWU04rWsjM5HgB7g4uIoBCEU13wgwmmcTjqPrVbvvJQUb2xbcjJZ3iJJqGoBfw2+4xQBBwBSykohxIQgtYkJMJwkOPzodGSN69ePPitPsceOkesVIMtQv37kegiFtc7joRuyMtXVecWby0XbTkrSbpmi0lKKfTt2jCxg8+eHXoSWlVFHfcstXH2hJyEEuQ8TEsh1P348uRdramgqKSFhV1pK9wJFyCmWO6ORrjtl1GZ3hJ6/4TQuV9uTQkwMWd/796f7U1KSNrwJkY6/Ik4nhOglpawEACFEahd+y6iAQCcJdrvpAnc6aVIueOV9fHxiU54lNViX1IYSO+ObWNTlosDps2ebW+tSU71P+jExdHM3GGjyfa989u0Q7PbmQq2y0lthpL6+dbuio+kGreV4Ho8H2LgRWLOGhPPDDwPXXBPa89DjoY560CCKf+MEsoxe73XBDhtG8zweukaVtDx2Oz1M1dTQq2LF9RV6MTEUW+mLEN7vlffK5927veE0TqcD+flmJCSYmt0nAK91UEl4HBdHIjQ+nubFxfn3P7VY1UPL+Hur/h8AFiHE+6ASW3cB+K+gtYoJG04nWWIKCykG6vRpurn4ijQl1qNjJgGgC18JBm45JSez6d0Xg4GebpOS6LNirTt/nkSzMim0FCW+N2OjkUSb203LKd9FRXnz20WiwC4ro7JZhw5R/rWHHgq9BczpJAvL2LEUAK/G8lmMOtDpWj/MtUQRelYrvV64AFx7LV3TyuTx0LUupTcHmzLf5crB+vUUThMVZcQdd+Rg4kTvfcJopHtCIO7FkVDVQ2v4O7BhlRBiJ4DpjbPuDGetUiYwKKPlTp4k0XbyJIk2RSgkJ1MMT3x8a+tPe5+VV6NxP4qLx+D/t3fmYXKU953//qanj7k10kijA12DkISEQEhjocGAB8QlIAGHYDt2bCdW1k7WOTYH+8RrZzdPHGeT2Mk6m2Q3cRw7eIONgh1jEh8KCCs2IBsJzCUhISEECBCS0Dmaq493//hWpXs0M1LPqLurqvv7eZ56urt6jrferq73W7/z0CFaJPbvB7Zvz8eEALx4TJ8+UtjNm8eMTi18Y1vrisG/iDc2VkYkh6F8jXPAD34A3HcfX//CLwDXXFN5kTo4SGvq1VczVkiIsZiItepMoXf0KItoF8vq1T24+OLNFbGOVVtXjygwkWK/O8zsMIAUAJjZPOfcq2UbmSg5x47lxdrLL3Pzg1RTKQq2m26iiOrqOr/MopaWI7j44pH7fHfhoUP57a23+Pj88/kMp4YG4OKLgeXLgWXLGCdWbRajchKLVU4Eh6F8zbFjwJe/zPpvS5cCGzbku2xUkr4+ur5uuIFuVCHGIghrVaXq+/X0MIkO4LH19PSW/X/WOsVmp/406FKdDeAQgPlgsd/l5RuaKAXHjgGPPQb88IcUSwAX+LlzWbyyq4uibdas8lttCussnUkux7pLL70E7NjB7amn+N60aRR0vqhTvaHwEGT5GueArVuBe+/lDcIHPsDCvUG46I953bpuuUUJDOLsVLO1qru7Bxs3VsbqJ0ixlrhPA1gL4GHn3OVmdi2Any/fsMT5kMkAzz5L99Kzz3KxW7oUWLeO3QTmzQtfxmFdHRe/qVOBd7yDYz50KC/otm3j8ZgxKHjZMsY8LVoUvmOpJYIqX3PiBHDPPcBPfsJzYMMGZs0FweHDvLG47jqWiRHibJyPtWr79q144IF7cccdydAKpFrr6hE0xYq4tHPubTOrM7M659z3zezzZR2ZmDAHD9Li9uijzG6aMoWWgauvzpe2iApmHHNnJxfHbJbuX1/UbdoEfOc7jKtbsoTtXLq7ZaWrNEGUr9m2DfjKVxh/9p73MAQgCOubcwwHmDmT8XfqwCCKYbLWKt8NOzw8hK997V4lDQgAxYu442bWDOAHAO41s0MAxihQICrN0BDw5JO0Uu3ezcXs0kvZWHvFiupJEPA7HixaBNx+O2P5du+moHvuOVpl/vEfecxr17IuWDIZ9Khrg1KXrxmPw4eZuPDUUwwB2LAhuNizXI5Zw4sXA1dcEe2SLKLyTMZa5bthc7kc0unqcsOKyXPWS4+ZJZ1zQwBuBzAI4DfB9lttAP6g/MMT47F/P61uW7dS0MyYAfzsz7Im1ZQpQY+u/DQ0UKitXEmLyCuvAD/6ETsiPP00BdyqVRR0y5ZpkY0yg4PAt78NfO97vEm5805g/frgblD8EiJqoSUqie+GdW6oIkkDqvcWDc61tG0FsArA3zjnPujtu6e8QxLjkctRqGzaxFpu8TiwejWtbosX1+5iYsbM2gUL6F578UWK2+3b+djSwji7tWtpyav2TNcwlPwoBbkcP7+vf51JLz09wF135VuZBcHgIHvJvvOd/M4JUSl8N+wDD3wJd9zxkbIKK9V7iw7nEnEJM3s/gCvN7GfOfNM598/lGZYoxDmW4Lj/ftZxmzMH+Pmfpyhpagp6dOGiro5JHEuXco6ee47C94c/BB55hKUnrriCc3fBBUGPtvSEoeRHKXjpJeCrX2U5nK4u4Fd/lUk5QdLXxy4XN96oEiIiGLq7e9DQMITly8v7na7mDNpq41wi7pdB9+kUAD91xnsOgERcmdm/n+Jt504Wxv3Yx4A1a2rX6jYR4nG6VFetosv5Jz+hoPvud+memzuXFpW1a6PX2H08giz5UQqOHaPl7fHHGRbwS79EC1zQ5/vx47QM3nILS94IUc2o3lt0OKuIc849amaPAzjgnFObrQpy+DDwjW8wxqu5Gfi5nwOuvVblNCZLQwP7Z155JTN3t22jULjvPvbYXLGCrWwuuyzacxxUyY/zZXiYYQL/+q8US7fdxqbxYcj4PHKE58/116uEiKgNVO8tOpwz3Ns5lzMz9UqtEKdOAf/yL3T9xWJczNavL775sDg3ra2smbduHfDGGyyG/PjjwDPP0D29di0F3fz50YufC6Lkx/ngHLOrN26kWFq9mnGNYxWEDmJsb73FMjfXXKMm9qK2UL23aFBszt5mM7sTwD87d+7W52LiDA0BDz3E2meDg6ztdscdwQZx1wKzZzNY/s47Wa7ksceAf/93YPNmxj1ddRXdeVFyt1aq5Mf5sn9/E774RWDXLsYn/tf/ilGt2oIim2XdxcWLGb4QZeusEKJ6KVbEfQzAbwHImNkgAAPgnHOtZRtZjZDNsjjvAw8w7mblSpYKUeB0Zamro0t1xQr2v/zxjynoNm5kTGK1uFuDxjlgzx7esDz1VDcaG4EPfYiWrrDUNBweZjjDqlWsuRg1a6wQonYoSsQ55xQJUmKco/vu/vvp0rvwQuBXfkVlC8JAYyPjD6+9lgVdH310pLt1zRputVzWZaKk0xTGDz3E8jhNTcDtt7+KdevmhyrDur+fLb16e1lQWAghwkxRIs7Mrhlrv3PuB6UdTm1w7Bg7DDzzDONtPv5xxgLpjj98zJqVd7fu3ElB99hjwPe/TxdrdzcF3aJFEnRjceIE4zu3bGFCyezZwIc/TBd1R8fLOHVqftBD/A+OH6dl/JZbmAkuhBBhp1h36t0Fz1MA1gB4EsB1JR9RFeMc22Nt3MjF4r3vZcabugmEn7o64JJLuA0NUYBv28bPc/NmlsPwBd2FF0rQvfwyrW5PPMFs00svBW64gd0zwnizcugQM0/VxF4IESWKdaeOqBFnZnMBfL4sI6pSDh8G/uEfaM1ZsgT4xV+MXlN6QZLJvEt1cJCC7oknaG16+GEmo7zjHdy6umpH0GUyzDR96CEW602l6JK+/vrwnuu5HDNQ/SSWMJQ0EUKIYpmsDegAgJDkkYWbXI7upK9/nRaID36Q8Ta1srBXO6kUO0BccQULCj/9NC10jzwC/Nu/AVOn5gXdZZcFPdry0NdHAfvIIwwVmDGDdQ2vvjrcZTkyGQq4ZctoRQ1LYoUQIpwMDITPc1ZsTNxfgh0aAKAOwEoAT5VrUNXCwYPAl77EbLxLLmEsUEdH0KMS5aKhgbFePT35DhHbttHdumkT0NLyTlx0ERMiLrqIdeiiKBzSaVraXniB2759DA9YtoyZppdeGv6blKEhWsfXrmVZkzC6eIUQ4SCXY8u9xkZg5sygRzOSYjXl9oLnGQBfc849VobxVAXZLK0w3/wmy1Fs2MD2TlooaofCDhH9/bTQ7d17BDt2zMJT3u1PKsX4ucWLuXV1AYlEsOMei1yO7d980bZnD8twmDGD86abKFyj0ou2sAdqVMYshAiGoSFe72bMYOxz2NbxYmPi7vGfm1k7gLllG1HEee01Wt/27wcuv5yWiSlTgh5V+cjlaJnJZvmYyXBLpRgsnkzyeSJRvpN/796toe5Q0NhIMXfTTbtx6tQsHDsGvPgitz17WCPQOVrlFi7Mi7qLLgqmU4dzLHuzcydF265dtCwCjB275hpa3ZYsiV4nkaNHaSW89Va6uoUQYiyc481eIgEsWMC1LIwU607dAuCnvZ9/EsAhM3vcOfebZRxbpMhkDN/6FltmNTQAv/zLDHwPm2qfKH19tCTlcvljcY7PfeHR2Mht6lTW/2psZMD/qlVspfT224yV8n/HjMKuFOJu796t+Oxn1yGdZqPmu+/eHEohV0h7ez6ODuCFYs+evLDbtImdO8womjo6eCPQ1sbHwuetrRN3yeZy/Ez7+kZu+/fvx4svnsbRoxfh9GmaBKdP53l88cXA0qXR6lxRiHN0n06ZwmSLMNWmE0KEi3SaN64dHVzXwhweUqw7tc05d9LMfgnAV5xz/8PMni3nwKLEs88Cn/rUarz6Khfm97+fi2tUSacpurJZnsSLF1OYJZN0DycS+W28IM8dO2it8clkKFb6+tgf1hd3hw5xgQUoWpJJlngoVpjs2rUF6fQwnMsikxnGrl1bQi/izqSpiZ06Vq7k66Ehxpy9+CLjzY4cAfbu5dydiRnnq60tL+za2iiQ/fk+czt9Oj/nI1kA4E3U1X0Dt976DvT2LqqKGE4/gaGri27fMLqshRDh4PRprj/z54c7McunWBFXb2azALwHwCfLOJ7I8dprwLvfDTQ1xfFrv0brUxTJ5ViYdXCQAmDFCp7EpXIF19fnhUYh2Sy/NKdPU9wdPMiK/tksBV1r69mzgZYu7UU8nkAmM4z6+gSWLu0tzYADJJmkAC4UwQDFyMmTLEp74sTYj6+/zp/JZim4m5u5tbQAc+fysbmZwtF/3twMPPHE/8H3vvdJAMcBxJBKfRodHZ8I4vBLin/TsHYt3b9hvqMWQgRHJkMPRXs7jRdRSTorVsT9AYBNAB51zm0zsy4Ae8o3rOgwdy7wZ38GpFLbEI9fFfRwJowvnszo91+0iAGclTqBYzEKtdZWdkdYvJiWwMOHGVf48sv8ciUSFIBnCrpFi3pw992bQx0TVyrq62naP1csVy6Xn7NiyWYvxyOPDCGTiVWFGHaOFsxUivFv1WBRFEKUh/5+Ps6dG71Qi2ITG+4HcH/B630A7vRfm9knnHP/s/TDiwZ33AE8+mgGQ0NBj6Q4Ct2l06Yxc3bOnPAUOo3H2Z5p9mzGYxUKunR6tKBbtKinqsXbRKmrm7jLsJrEcDpNN31XF8MbwnJeCyHChe8JamlhQfKw1YArhlIN+S4ANSviooBzdLf57tJLLqG7tL096JGdnfp6WuhmzWLB3CNH6G7du5eLdTxOQRePBz3S6FMNYvjUKbpQr7ySVt2oJxYJIcrDwABF3OzZFHFRvVaUSsRF9PBrg4EBllaYP59ZhpV0l5aS+noWWpw5E1i9Oi/oXnqJyQCpFAWd4p5qDz/7tLERuO02WpiFEOJMCgv3dnZGP9GpVCJuzFw3ESy5HBe2RILNx6upsGksxi9gZyeTSQ4dAnbvBl55hSJuypTw1vURpcV3ny5aRGut3KdCiLEYHGS8cGcnb/ijan0rRJa4KqWvjxmLy5axZ2c1L2yxWN7levo04+eef57Wx8ZGJk1Uw5dVjObkSX7mV11FEafPWQhxJr71raGBxoyoW98KKZWIu//cPyIqQTZL61tzM3DLLbzjqCWamoDly1mY9uBBdhs4cIBCr71dsXPVgu8+bWoCfuqn1H1BCDE2vvUtrG2zzpdiOzakAGwAsBzAf9h0nHMf8R7/qCyjExPixAmmSl92GYVMLQuWWIwZt3PmMNh93z62kBoczNdNq7Yvc63gu08XLwa6u+U2Lxe5HG8KMxk+dy7fcWUs6ur4nhmf+6/9RyEqSS5Hj1Q1Wt8KKdYS9/8A7AJwE1gz7gMAXijXoMTESKcZ5D91KrBuXXStEtu3b8XWrVvQ09OL7u7SZUm2tOSF7cGD7CZx8KCsc1HDOXb5yGbZv7WrS+LgfHAu3+vYF2n+fjMmEsXj/P4kk3wdi3FzLv87uVy+NmE2mxd+hX+78O/HYvy79fX6/ER5GBzkutjZWZ3Wt0KKFXGLnHN3mdntzrl7zOyrAH5YzoGJ4jh6FBgeZkD3kiXRzDoFKODe+958D9SNGzeXVMgBXDQuuIDbiRPMat29m/Pndy8Q4cSP8ezqYiJLS0vQI4omuRwzuTMZWshSKZ73fks9X6TFYqXN8naO4m54mAtsYfs3XzAmEsosF+eHH/uWStELUwtW+mJFXNp7PG5mlwA4CGBGeYYkimF4mDFBc+awpVCUe7UCwNat7IGazWYBDGPr1i0lF3GFtLVRDKxYAbzxBq1zb77JhWzKlGgWfaxGhodpZW5vr80Yz1LgC7dsloKprS1vXauUaPKFWn09k42mTuW40ml+xqdPMxQk7a00vrVOiGIZGuK5VK2xb+NR7FL1BTNrB/B7AB4E0Azgv5dtVOKsHDnCiJaPegAAGg9JREFUx2uuARYurI67154e9kAFaInr6emtyP+Nx1k/b/58WjX9xvOZDBe7xsaKDEOcQS7H89yMhXsvvDC6VuYgGE+4pVLhWdzq6igkk8m8ZTWTGWmty+UY01pfX1nRKaKD7673W0fWgvWtkGLbbn3Re/rvALrKNxxxNnI54K23WGH6yiuj1+PtbHR392Djxs1liYkrFr8v6aWXMqO10DrX3i4RUSn8BJ2LL6altKEh6BFFgzNdpWEUbufiTGvdvn3AvHl5d3o2y+9hKiVBJ/JdF+rreSMelfO8lBSbndoJ4I8AzHbOrTezZQB6nHN/X9bRif8gk2Ew/vLlzMirRkHR3d0TiHg7k2SSlp+uLlrn9uzh5hwXRomK8uB3Fpk9O9oJOpXGD+KOxRhWETXhdi4aGrh1dFCk9vWxPqB/zMlkdV4Pxfik07xetLbyvHjjjeo53ydKse7UfwDwZQCf9F6/CGAjAIm4CjA4yKy8nh7WP6vVk7XSmLF907RpwMqVtM4995xi50pNJsPzO5kErruOlhed4+fGjwFqamKh62oSbmNhxmNMpfidHBpiLN2JE1zQzfJZtKI68RMX4nFeJxTuUryI63DO/ZOZfQIAnHMZM8uWcVzC4+RJirgbb2QSgwiGVIodAS68kAklr7xC61w6rZ6t54OfXX3ZZXSfVmstp1IyPEwB09DAhawWLcOFgm7qVM5Jfz9w7BgFnR9vJ0FXHTjHzxdgclNrq663PsWe4qfNbBq8HqlmthbAibKNSgBgYHciAdx6K2OygqBctduiihmzn2bMAC6/nEVn9+1jq69slneGLS26wJyNTIaLbSoFzJzJLOGoZ1dXgnSaN3TJJDB3LsVbNVveisW3wCWTvE762a6+oFNSRLQZGuLW3k4LrIT5SIqdjt8Cs1IvNLPHAEwH8LNlG1WN4xwTGGbMYAZqUCbjStRuizL19Yzfmj2bdfoOHWJm64EDXFiam+nq0kLLc7qvj1siwbCAeJxJJOLsZDIUI/E4axw2NuqcOhuJBLcpU7j4nzoFHD9OV1w8TkGn+Qs/2Wy+3+mCBdXd//t8KDY79SkzexeAJWCz+93OufQ5fk1MAj+BYckSYM2aYGslVbp2W5TxrSNz53LBfeMNFhI+eHBkiYdaI52mRSSbpdWtu5uiNx5n9q8Yn2yWLqR4nHPW3CzxMRHOjKEbGGD83KlTeeudatGFD991asbzXi0Sz06x2al3Afiec26HmX0KwCoz+0Pn3FPlHV5tMTREF+qaNcxCDfrEDap2W9RpaGDs3IUXcsF4/XVg1y4KmoMH+X5TU/W6BZxjLKdfOf2SS1jPsK0t6JFFg1yOi1gsRuEr9/z5U1fH71xTE2+U+/sZj3nqVL5zhTJcg2dggJ/PtGkq61QsxS4jv+ecu9/MrgKwDsDnAPxfAFeUbWQ1Rl8fLyjXXcd6N2FgsrXbtm/figceuBd33JGsectdSwtdh0uWMLP16qvpbj1wgKLOLC/qon7BGh6m1S2Xy9cy7OysXrFaavw+pAMDwPTpSpYpF/X1jMFsbeU529dHd2t/P9+r9izfsOFcvkxOSwvPfSU4FU+xl1c/E/VWAH/nnPu2mf1hOQZkZr8P4D8BOOzt+m/Oue94730CwAZvPL/unNtUjjFUmmPH+HjbbbwDCRMTrd3mx9ENDw/ha1+7V3F0HmYUaQsXcsvl6Np5+23gtdfofs163zLfYhD2BdwXHP39vHtOpZjsMX9+bbqOzwe/S0EsxvqEURf0USGRYHZrezvn34+fc07xc+XGv37kcrxhaW+vvW4LpaBYEfe6mf0tgBsA/ImZJQGUc4n5X865zxXu8AoMvw/AcgCzATxsZoudc5EtdeIcg+Hb24He3upowO7H0eVyOaTTiqMbj7o6fu7t7Sxdks1S1B05Arz6Kt2ufiuZhoZ8/E5Qws538fkV0v1jmD6dorSzk4k4Eh8Twzm6nRMJBm8fPKg5DAL/e+YXFfbj5/r6+BklEhIYpSKbzdf1a2+ngFNs4uQpVsS9B8DNAD7nnDtuZrMA3F2+YY3J7QDuc84NAXjZzPYCWANga4XHURJyOV6wFy5kEd9qMR/7cXTODSmObgLEYvm2X4sX07J1/DhF/ltvcUE5cSIv7Jzj7yWT+Wy8Ul0ICwVbJpO3InZ00ErU0UFXVHNz+K2FYcav9zZtGj93zWU4KIyf8wXH8eMUdED++yYmhl8ip76eN3wtLbphKQXm/NUgJHju1F8AcBLAdgC/7Zw7ZmZ/BeBHzrl/9H7u7wF81zn39TH+xkcBfBQAOjs7V993331lHXM6DZw61Qfnijel+UViq6n/qc/OnTvw5JNPYPXqNVi2bHnQwwkNg4N9SKXOz9yay+Vjp7LZkVsuV5pxmvFCG4/zsa6uPBfbUsxHVPHFeDw+0l3X19eH5mowyZeIsM1HLscbG3/ZNKu8uzVq3xv/elVXl7+elJqwnSfny7XXXvukc667mJ8NJOTYzB4GMHOMtz4JJkx8Giws/GkAfwbgIxP5+865LwD4AgB0d3e73t7e8xnuOXntNeDRR7dgaOjc/yebpQVuxQoWOa3Gu+/ly3uxbNlyLF/eG/RQQsWOHVvKOif+na5vRXOO55dZ/rHweV3dyP2+xa1S8Xjlno8wci7r25YtW1Du61WUCOt8pNP5DNfh4cqWLInK98a35Psu6nIWpw7reVIJAhFxzrnri/k5M/s7AP/qvXwdwNyCty/w9kWGTIauscsvZ5shBcyKUhKPc1NSQfjwY9/icSZ+qHBptInHGcvV1kZR3t/PpIhTp/h+LXaJcI6CNu1VkG1pYcybzvXyErrkfzOb5Zx703v5bgDPe88fBPBVM/tzMLHhIgBPBDDESZFOM77JrwEnhKgNFPtW3RS2/MpmaQ0/fZq1ErNZft7xeHXG0WWzPL/92NnmZiY7pVIqLVQpwjjNf2pmK0F36n4AHwMAr9DwPwHYCSAD4ONRyUz1Bdw738mgdSFEZQmiB7Csb7WHH47Q1EQx45eOOXkynxgRi1HQRTWoP53mcfltzNraeLy1ZnkMC6ETcc65D57lvc8A+EwFh3PeDA2xFti73sXMPiFEZQmiB7C/eHd0yPpWq/hxcskkhU42m3e9njyZby3lJw7V14dT2DnHcWcyfJ1MMrs0laIYVVhQsIROxFUTfmr6unXsqSmEqDyV7gF8+jQXYzXtFoXEYkBjI7eOjnz82PAw14rBwXxCki/u6uvLl9F5Jn7mrZ/tDuTH0tzMskJyk4YPfRxlor+fdb2uvx6YMyfo0QhRu1SqB3AuR5dZW5sKH4tz49eba2piPB3AcyidppgaGqKoGxrKt+jzf2ZggM8LM8vP9do5irNMJl+iyN/vlxVKJjkev+5kLBZskXFxbiTiykBfH0Xc+vW8mAshgmOyPYAngp+8MGsWRZwQk6GuLu+CLawh6ouvdJolrVpbKb78Gmx+/cjCGpJ+LTv/fT/BorExXw7Fd+HW18stGlUk4krMyZO8oK9fH74+qELUKhPtATwR+vu5QC5YoNZMojzEYtySST5Onx70iERYkIgrIX5bpPXr8+ZxIUR1IvepECJoJOJKhN9G5+ab5U4RotqR+1QIEQYk4kpAfT0DQW++mVk8QojqRe5TIURYkIgrAVOnssWIBJwQ1Yvcp0KIsCERVyKU2SPKQRCdBsRo5D4VQoQRiTghQkoQnQbEaOQ+FUKEFZXwEyKkFHYaSKfZaUBUjlyOJYOamtj7VAJOCBE2ZIkTIqRUqtOAGE2h+7S1VeESQohwIhEnREipRKcBMZqBAYq2+fPV+1QIEW4k4oQIMeXsNCBG4meftrQAM2cq+1QIEX4k4oQQNU86DQwOAp2dwJQpcp8KIaKBRJwQoqYZGODjvHlAQ0OwYxFCiIkgESeEqEmco/u0qYnu03pdDYUQEUOXLSFEzZHJsP7bjBlAe7vcp0KIaCIRJ4SoKQYGaIWbNw9obAx6NEIIMXkk4oQQNYHvPm1spPs0Hg96REIIcX5IxAkhqh7ffdrRAUybJvepEKI6kIgTQlQ1/f18nDuXSQxCCFEtSMSJs5JO04rhXH4rBr9wKsDficW41deriKqoDLkccPo00NzM+m/KPhVCVBu6rIkRZDLsG5nN0uWUSrF3ZF1dXoiZjd7q6ka+fuMNoKuLf8//m0NDLKjqB5b7Lq1YjL8vgSdKxeAgb0BmzWIHBrlPhRDViERcjZPNUmBlMlzoEgmWXGhsBJJJiqvJUl8/tvXDuby4y2bz4m5oKO/6qq/nWCTqxETwkxdSKeCCC3gOCSFEtSIRV2PkchRL2Sxf19fT0tbUxAWvEi4nM2YG+tmBzc2jx3f6NHDyZL6afiLBn5dFRYyHb+2dPl2134QQtYFEXA3gHMVQNkvLVksLRVsyGb4yC3V1bH3U0MBMwuFhWulOnqSwc05WOjES52jBjcWA+fNphRNCiFpAIq6K8Rc359jUu62N4idKFopEgltrK610g4MUc6dOjbTSyW1Wm/ilQ9rbKfol7IUQtYREXBWSy+WTB6ZOpXgLm8VtMtTVMVavsZEus+FhHuepU9zMwmldFOVBpUOEELWORFwVkc1S1JixoGlra3WXVfAtcG1teYvM0aMUdHV1dKvJMlN9+OVrWlpUOkQIUdvo8lcFZDIUb7EYLVStrbUnXvwEjdZWWuj6+oDjxyns6usp6KLkRhZj099PETd7tkqHCCGERFyESacZIxaPsx5Wc/P5lQSpFhIJupHb2zk/p04BJ07QvRyP0+UqosXQEMV5WxutzHKZCyGERFwk8UspxOO0SDQ3yyIxFmYjM137+2md8+PnUim54sJONstEloYGZZ4KIcSZaAmLEM5xQYvHWci0sVHirVjq6ih2m5tpwezvB95+m25ouVvDh59ZbSbXqRBCjIdEXETw3UnTp7NciNymkycep1uurS1fg853t6pcSfAMDDDOc9o0usRrLb5TCCGKRSIu5PhNvJNJYMECxXOVmlSKm+9u9bNbYzHOtQRE5UinKeBaWnizIjEthBBnRyIuxPhNvGfMoPVN7qTyUehu9bNbjx2jRchPhtD8lwf/RiUeB+bNY5iAEEKIcyMRF0L8Ra2hQU28g6Awu3VgQMkQ5cTvKDJzJi1wChMQQoji0XIUMvx4oM5OxmzJ+hMcZvkOEZkMhfXRo/mafKmURMdk8Nun5XKs6zd1qoSxEEJMBl06Q4Jfhb6piW2EZH0LF/X1I5Mh+vqYDJHN5rtDSNCdHb80TixG4fbaawwVEEIIMTkk4kLAwADFwKxZtEzI+hZu/GSIadMoSvzuENksxV4yKUHn4xxFbybDOZszh5bNujqd50IIcb5IxAVINsuYoOZmWiRUhT5a+DFyvqAr7A6Ry9V2/Tm/jy+Qt2CqUK8QQpQWibiA8Ou+zZqlQqbVQGF3iOnTxxZ0tYB/XsfjjOtsaqqdYxdCiEqjy2sA9PczLkh136qTM9t9+QWFczkKu7o6CptEojrEezZL8ZbNUrR1dvLYq+HYhBAizEjEVRA/eaGtje5TFZKtfurq8hmuu3ZRuPuJEX55DYCWq6iIukyG9QszGb6Ox1nHsLVVCTlCCFFJJOIqRDrNxbuzU4V7a5lkkltbGwWcn7HZ18cSJrkczw3fUheGBIl0mls2y9f++BsbOUa5S4UQIhh0+a0AAwNcsOfNo5tJCIBizRd1ra08R9Lp0aLOJxbj78Ri3Moh8PwxpNP5/93QwJIgDQ0UbbIgCyFEOJCIKyPO5TsvzJoli4U4O2YUSYkEk118QZXJ0Ao2PMwtneaj7870XbJmFHb+5r/ni7FcLv+zZxtDQwO7VSST4bEGCiGEGI1kRZnIZingpk3jpoVQTJRCUTcWzvE887dcLi/whofzYq6+fuSjX6NtvEe5+oUQIhoEJuLM7C4Avw/gYgBrnHPbC977BIANALIAft05t8nbfzOAvwAQA/BF59wfV3rcxeCXWZgzhxYVIcqBHztXSxberVu3YsuWLejt7UVPT0/QwxFCiEAJ8vL/PICfAfC3hTvNbBmA9wFYDmA2gIfNbLH39l8DuAHAAQDbzOxB59zOyg353PT306Kh8iFClJatW7di3bp1GB4eRiKRwObNmyXkhBA1TWBOPufcC8653WO8dTuA+5xzQ865lwHsBbDG2/Y65/Y554YB3Of9bGg4eZJ1subPl4ATotRs2bIFw8PDyGazGB4expYtW4IekhBCBEoYI7XmAHit4PUBb994+0NDZycwc6ay94QoB729vUgkEojFYkgkEujt7Q16SEIIEShldaea2cMAZo7x1iedc98q4//9KICPAkBnZ2dF7tiHh/vwzDPl/z9Roa+vT5aSM9CcjGQy8/HZz34WTz/9NFauXImhoaGqm0+dIyPRfIxGczKaWp6Tsoo459z1k/i11wHMLXh9gbcPZ9l/5v/9AoAvAEB3d7erxB27H2wtiOZjNJqTkUxmPqp9/nSOjETzMRrNyWhqeU7C6E59EMD7zCxpZgsBXATgCQDbAFxkZgvNLAEmPzwY4DiFEEIIIQIjyBIj7wbwlwCmA/i2mT3tnLvJObfDzP4JwE4AGQAfd85lvd/5VQCbwBIjX3LO7Qho+EIIIYQQgRKYiHPOfRPAN8d57zMAPjPG/u8A+E6ZhyaEEEIIEXrC6E4VQgghhBDnQCJOCCGEECKCSMQJIYQQQkQQiTghhBBCiAgiESeEEEIIEUEk4oQQQgghIohEnBBCCCFEBDHnXNBjKCtmdhjAKxX4Vx0AjlTg/0QFzcdoNCcj0XyMRnMyEs3HaDQno6m2OZnvnJtezA9WvYirFGa23TnXHfQ4woLmYzSak5FoPkajORmJ5mM0mpPR1PKcyJ0qhBBCCBFBJOKEEEIIISKIRFzp+ELQAwgZmo/RaE5GovkYjeZkJJqP0WhORlOzc6KYOCGEEEKICCJLnBBCCCFEBJGIGwcz+5KZHTKz5wv2XWZmW83sOTP7FzNr9fbHzeweb/8LZvaJgt+52cx2m9leM/vdII6lVJRwTvZ7+582s+1BHEspmOB8JMzsy97+Z8yst+B3Vnv795rZ/zYzC+BwSkIJ52SL97152ttmBHA4542ZzTWz75vZTjPbYWa/4e2famYPmdke77Hd22/eObDXzJ41s1UFf+vD3s/vMbMPB3VM50uJ5yRbcI48GNQxnQ+TmI+l3vdpyMx+54y/VRXrTYnnpCrWm3FxzmkbYwNwDYBVAJ4v2LcNwLu85x8B8Gnv+fsB3Oc9bwSwH8ACADEALwHoApAA8AyAZUEfW5Bz4r3eD6Aj6OOp8Hx8HMCXveczADwJoM57/QSAtQAMwHcBrA/62EIwJ1sAdAd9PCWYj1kAVnnPWwC8CGAZgD8F8Lve/t8F8Cfe81u8c8C8c+LH3v6pAPZ5j+3e8/agjy/IOfHe6wv6eAKYjxkA3gHgMwB+p+DvVM16U6o58d7bjypYb8bbZIkbB+fcDwAcPWP3YgA/8J4/BOBO/8cBNJlZPYAGAMMATgJYA2Cvc26fc24YwH0Abi/32MtFieakapjgfCwD8Ij3e4cAHAfQbWazALQ6537keMX5CoA7yj32clGKOanAMCuGc+5N59xT3vNTAF4AMAe8Dtzj/dg9yH/mtwP4iiM/AjDFO0duAvCQc+6oc+4YOI83V/BQSkYJ56QqmOh8OOcOOee2AUif8aeqZr0p4ZxUPRJxE2MH8l+KuwDM9Z5/HcBpAG8CeBXA55xzR8GT7rWC3z/g7asmJjonAAXev5nZk2b20UoOtgKMNx/PAPhpM6s3s4UAVnvvzQHPC59aOkfGmxOfL3sukN+LsovZx8wWALgcwI8BdDrn3vTeOgig03s+3jWjKq8l5zknAJAys+1m9iMzi+zNj0+R8zEetXyOnI1qXm8k4ibIRwD8ZzN7EjTxDnv71wDIApgNYCGA3zazrmCGWHEmMydXOedWAVgP4ONmdk2Fx1xOxpuPL4EX1e0APg/gcXB+aoHJzMkHnHMrAFztbR+s6IhLjJk1A/gGgP/inBthkfYssDVXJqBEczLfsVL/+wF83swuLP1IK4POkdGUaE6qeb2RiJsIzrldzrkbnXOrAXwNjD8AeAH5nnMu7bmFHgPdQq9jpGXhAm9f1TCJOYFz7nXv8RCAb4KCryoYbz6ccxnn3G8651Y6524HMAWM83gdPC98auYcOcucFJ4jpwB8FRE+R8wsDi5E9zrn/tnb/ZbvEvQeD3n7x7tmVNW1pERzUnie7APjKC8v++DLwATnYzxq+RwZl2pebwCJuAlhXoacmdUB+BSAv/HeehXAdd57TWDw7S4woPsiM1toZgkA7wMQyQyq8ZjonJhZk5m1FOy/EcDzZ/7dqDLefJhZo3e8MLMbAGScczs918BJM1vruQw/BOBbwYy+PEx0Tjz3aoe3Pw7gNkT0HPE+078H8IJz7s8L3noQgJ9h+mHkP/MHAXzIy8hcC+CEd45sAnCjmbV7GXk3evsiR6nmxJuLpPc3OwC8E8DOihxECZnEfIxH1aw3pZqTal9vACg7dbwNtBi8CQZKHgCwAcBvgJaCFwH8MfLFkpsB3A/G/uwEcHfB37nF+/mXAHwy6OMKek7AzKlnvG1HlOdkgvOxAMBuMED3YdAN5P+dbvDC8hKAv/J/J4pbKeYEQBOYqfqsd478BYBY0Mc2yfm4CnT5PAvgaW+7BcA0AJsB7PGOfar38wbgr71z4TkUZOiCbum93vaLQR9b0HMC4Erv9TPe44agj61C8zHT+26dBJOBDoDJUUCVrDelmhNU0Xoz3qaODUIIIYQQEUTuVCGEEEKICCIRJ4QQQggRQSTihBBCCCEiiEScEEIIIUQEkYgTQgghhIggEnFCCCGEEBFEIk4IIYQQIoJIxAkhhBBCRJD/D3XQFtwBIQitAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 5))\n",
    "ax.plot(annual_mean_df.year, annual_mean_df.height, 'k.')\n",
    "ax.set_ylabel(\"{} [{}]\".format(ssh.standard_name, 'mm'))\n",
    "ax.set_title('Annual mean')\n",
    "conf_int = prediction.conf_int()\n",
    "pred_int = prediction.conf_int(obs=True)\n",
    "ax.fill_between(annual_mean_df.year, pred_int[:, 0], pred_int[:, 1], alpha=0.1, color='blue', label='prediction interval')\n",
    "ax.fill_between(annual_mean_df.year, conf_int[:, 0], conf_int[:, 1], alpha=0.3, color='blue', label='confidence interval')\n",
    "ax.plot(annual_mean_df.year, fit.predict(), color='blue', label='fitted trend')\n",
    "ax.legend()\n",
    "ax.grid(True)\n",
    "fig.savefig('nodal.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnEAAAE/CAYAAAApN5W5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsvXl8W9WZ//8+3tfEzuYsTkggQIKz4wYMhBpStk6H0n6BAKE0ZEqgtKXTTlPodDpDCzOlhWkplBboAIFpALMV+EFpKSmmLC7FSQNDFkggJomdPXFieZEs6fz+OD6+V7Icy7YkS/Lzfr3uy5Lu1b2Pj2Xpo2dVWmsEQRAEQRCE1CJjqA0QBEEQBEEQ+o+IOEEQBEEQhBRERJwgCIIgCEIKIiJOEARBEAQhBRERJwiCIAiCkIKIiBMEQRAEQUhBRMQJgiD0E6VUtVJq51DbIQjC8EZEnCAISY9SqlYpdUgplTvUtgiCICQLIuIEQUhqlFJTgUWABi4cUmMEQRCSCBFxgiAkO1cBfwVWAV9271BKrVJK3aOUelEp1aKUelspdZxrv1ZKXaeU2qKUau46VnXtu1kp9VvXsVO7js/qun+1UmpT13k/VkpdG63BXee5vuu6LUqpW5RSxyml3lJKHVFKPaGUynEd/zml1PouG99SSs1x7btJKfVR13k2KqW+4Nq3TCn1hlLqji5P5Tal1AX9WVxBEFIXEXGCICQ7VwGru7bzlFJlYfsvA34IlAJbgf8M2/854FPAHOBS4Lwor7u367kjgKuBnyulFvTD7vOAk4FTge8C9wNXApOBWcDlAEqp+cCDwLXAaOA+4HlX6PgjjCdyZNfv+Vul1ATXdU4BPgDGAD8FHrBCVRCE9EZEnCAISYtS6gzgGOAJrfVajKC5Iuyw32mt/6a19mOE3ryw/bdprZu11tuBVyPsj4jW+kWt9Ufa8BrwMkZMRctPtdZHtNYbgPeBl7XWH2utDwMvAfO7jlsB3Ke1fltrHdBaPwx4MeIPrfWTWusmrXVQa10DbAEWuq7zidb6N1rrAPAwMAEIF7qCIKQhIuIEQUhmvowRP/u77j9KWEgV2O263QYU9XN/RJRSFyil/qqUOqiUagY+i/F2Rcse1+32CPetHccA/9IVSm3uutZkYGKXHVe5Qq3NGC+e247u309r3dZ1M6rfURCE1CZrqA0QBEGIhFIqHxP+zFRKWaGSC5QopeZqrd8d5CVagQLX/fGua+cCT2NCuc9prTuVUs8C8QhT7gD+U2sdHgZGKXUM8BtgMVCntQ4opdbHyQ5BEFIM8cQJgpCsXAQEgJMwIdB5wEzgdYy4GizrgTOVUlOUUiOB77n25WAE4z7A31UscG4MrhmJ3wDXKaVOUYZCpdQ/KKWKgUJMVe4+MMUWGE+cIAiCiDhBEJKWLwMPaa23a6132w34JbDUVpEOFK31n4Aa4D1gLfCCa18LcAPwBHAIk4f3/GCudxQ76oFrML/XIUxxxrKufRuB/wbqMOHY2cCb8bBDEITUQ2mth9oGQRAEQRAEoZ+IJ04QBEEQBCEFEREnCIIgCIKQgoiIEwRBEARBSEFExAmCIAiCIKQgIuIEQRAEQRBSkLRv9jtmzBg9derUuF+ntbWVwsLCuF8nVZD16ImsSSiyHj2RNQlF1qMnsiY9Sbc1Wbt27X6t9dhojk17ETd16lTq6+vjfp3a2lqqq6vjfp1UQdajJ7Imoch69ETWJBRZj57ImvQk3dZEKfVJtMdKOFUQBEEQBCEFEREnCIIgCIKQgoiIEwRBEARBSEHSPicuEp2dnezcuZOOjo6YnXPkyJFs2rQpZudLdVJpPfLy8igvLyc7O3uoTREEQRCEqBmWIm7nzp0UFxczdepUlFIxOWdLSwvFxcUxOVc6kCrrobXmwIED7Ny5k2nTpg21OYIgCIIQNcMynNrR0cHo0aNjJuCE1EUpxejRo2PqlRUEQRCERDAsRRwgAk7oRl4LgiAIQioybEXcUNPQ0MCsWbOiPv7ee+/lkUceOeoxq1at4utf/3rEff/1X//VL/sEQRAEQUhuRMSlCNdddx1XXXXVgJ8vIk4QBEEQ0gsRcUNIIBDgmmuuoaKignPPPZf29nY++ugjzj//fE4++WQWLVrE5s2bAbj55pu54447AHjnnXeYM2cO8+bNY+XKlSEevaamJs4//3yOP/54vvvd7wJw00030d7ezrx581i6dGnif1FBEARBEGKOiLgoqaur48c//jF1dXUxO+eWLVv42te+xoYNGygpKeHpp59mxYoV3H333axdu5Y77riD66+/vsfzrr76au677z7Wr19PZmZmyL7169dTU1PD//3f/1FTU8OOHTu47bbbyM/PZ/369axevTpm9guCIAhCuhMMgtcLHg8cPjzU1oQyLFuM9Je6ujoWL16Mz+cjJyeHNWvWUFVVNejzTps2jXnz5gFw8skn09DQwFtvvcUll1zSfYzX6w15TnNzMy0tLd3Xv+KKK3jhhRe69y9evJiRI0cCcNJJJ/HJJ58wefLkQdsqCIIgCOlOIACdneDzQXs7tLWZ+wBag1LQ9RGbFIiIi4La2lp8Ph+BQACfz0dtbW1MRFxubm737czMTPbs2UNJSQnr16+P2Tn9fv+gbBQEQRCEdMTvdwRbW5sRbfYjUynIzITsbHB9rOLxDI2tvSHh1Ciorq4mJyeHzMxMcnJyqK6ujst1RowYwbRp03jyyScB04j23XffDTmmpKSE4uJi3n77bQAef/zxqM6dnZ1Np/06IQiCIAjDBK0doXboEDQ2wtat8PHHsGMH7N0LHR2QkwNFRWYrLIS8PCPkkpkhF3FKqQeVUnuVUu+7HhullPqTUmpL18/SrseVUuoupdRWpdR7SqkFibCxqqqKNWvWcMstt8QslNobq1ev5oEHHmDu3LlUVFTw3HPP9TjmgQce4JprrmHevHm0trZ2h0+PxooVK5gzZ44UNgiCIAhpizt/bf9+2L7dCLaGBti5Ew4cMN63/PxQwZabCxlDroj6TzKEU1cBvwTcTdBuAtZorW9TSt3Udf9G4ALg+K7tFODXXT/jTlVVVUzF29SpU3n//W7dyne+853u23/4wx96HH/zzTd3366oqOC9994D4LbbbqOyshKAZcuWsWzZsu7j3LlyP/nJT/jJT34SK/MFQRAEYUix+WvBoPGmtbUZj5slMxOysqCgwIRH05EhF3Fa678opaaGPfx5oLrr9sNALUbEfR54RGutgb8qpUqUUhO01rsSY21y8OKLL/LjH/8Yv9/PMcccw6pVq4baJEEQBEGIG5Hy12yGUGcntLSY/LWioqG1M9EMuYjrhTKXMNsNlHXdngTscB23s+uxYSXilixZwpIlS4baDEEQBEGIKVobUdbZacKibW0mXy0YNPuVMt617GyTswYmDJqfP3Q2DyXJKuK60VprpZTuz3OUUiuAFQBlZWXU1taG7B85ciQtLS0xsxFM495YnzOVSbX16Ojo6PE6iTUejyfu10glZD16ImsSiqxHT9JtTbQ2WzBoNh32aa9U36HQjg4PGzbUxs1GN8EgNDUl5FJRkawibo8NkyqlJgB7ux5vBNxNz8q7HgtBa30/cD9AZWWlDq8m3bRpE8XFxTE1uKWlJebnTGVSbT3y8vKYP39+XK9RW1sbt8rmVETWoyeyJqHIevQkVdfEHQ5tbw8Nh4KTv5aV1f/8tQ0baqmoqI6pvb3h8cAJJyTkUlGRrCLueeDLwG1dP59zPf51pdTjmIKGw8MtH04QBEEQkhWtHcFmw6Ht7aYIAULDoe7+a8LAGHIRp5R6DFPEMEYptRP4D4x4e0Ip9U/AJ8ClXYf/HvgssBVoA65OuMGCIAiCIBAMOvlr1rvm9Toh0YwMI9jy8lKzfUcqMOTLqrW+XGs9QWudrbUu11o/oLU+oLVerLU+Xmv9Ga31wa5jtdb6a1rr47TWs7XW9UNt/0C56667mDlzJkuXLuX555/ntttuA+DZZ59l48aN3cetWrWKpn4G4BsaGpg1a1bExx999NHBGT6A6wqCIAipjd9vRNrhw7B7N2zbZvqvffKJaZ575IgRbwUFTv+1ggLTQFcEXPwYck/ccOVXv/oVr7zyCuXl5QBceOGFgBFxn/vc5zjppJMAI+JmzZrFxIkTB31NK+KuuOKKHvv8fj9ZWfJyEARBGM64q0Pd7TwkHJqcyKf2EHDdddfx8ccfc8EFF7B8+XJKS0upr6/niiuu4Pnnn+e1117j1ltv5fLLL6e+vp6lS5eSn59PXV0dGzdu5Nvf/jYej4cxY8awatUqJkyYwNq1a1m+fDkA5557bsTr3nTTTWzatIl58+bx5S9/mdLSUp555hk8Hg+BQIDXXnuN22+/nSeeeAKv18sXvvAFfvjDH9LQ0MAFF1zAGWecwVtvvcWkSZN47rnnyM/Pj+q6giAIQvJhm+V2dpo2Hm1tTjhUKQmHpgLyZxkC7r33XiZOnMirr77Kt771re7HTzvtNC688EJuv/121q9fz4033khlZSWrV69m/fr1ZGVl8Y1vfIOnnnqqWzx9//vfB+Dqq6/m7rvv7jFr1c1tt93GokWLWL9+ffd1161bx1NPPcVrr73Gyy+/zJYtW/jb3/7G+vXrWbt2LX/5y18A2LJlC1/72tfYsGEDJSUlPP3001FfVxAEQRg6rHetrQ2am2HXLjM31IZDd+0y4VClTBi0uFjCoanCsPfE/fM/w/r1gz9PIJDfPSh33jy4887BnzOcDz74gPfff59zzjmn65oBJkyYQHNzM83NzZx55pkAfOlLX+Kll16K6pznnHMOo0aNAuDll1/m5Zdf7m614fF42LJlC1OmTGHatGnMmzcPgJNPPpmGhoZBXVcQBEGIPe5ig44OEwq1zXJtz7WsLCPObLNcIXUZ9iIuldBaU1FRQV1dXcjjzc3NAz5nYWFhyPm/973vce2114Yc09DQQK4r+SEzM5P29vYBX1MQBEEYPL31XrPh0OEwO3S4M+xFXKw8Zi0t7TFpbltcXBwy6cB9/8QTT2Tfvn3U1dVRVVVFZ2cnH374IRUVFZSUlPDGG29wxhlnsHr16qjOHc55553HD37wA5YuXUpRURGNjY1kZ2f3enxJSUlU1xUEQRAGTjDoCLZAAHbudLxrVrBJscHwZNiLuGTjsssu45prruGuu+7iqaeeYtmyZVx33XXdhQ1PPfUUN9xwA4cPH8bv9/PP//zPVFRU8NBDD7F8+XKUUr0WGMyZM4fMzEzmzp3LsmXLKC0tDdl/7rnnsmnTJqqqqgAoKirit7/9LZk2ThyBaK4rCIIgRIcVa7b3Wltb6GQDv99s+fniXRNA6fBBZWlGZWWlrq8PbSe3adMmZs6cGdPrpNqYqXiTausRj9dEOKk6LideyHr0RNYklHRej/BB7zYcGgg4uWvuUVSWRI6YShXSbeyWUmqt1roymmPFEycIgiAIcSTcu9bebvLYwIg5K9aklYfQX0TECYIgCEIMcOeu9eZdy8gwuWtFRUNtrZAOiIgTBEEQhH5ytNw18a4JiUJEnCAIgiD0gu275veHTjUIrwzNypLK0HSnvr6O116r5eKLq7sLAIcaEXGCIAjCsEdrp/LTzgzt6OjZdy0zUypDhyP19XUsWbIYn8/Hr3+dw5o1a5JCyImIEwRBEIYVgYATDnVPNXA3a5C+a4KburpaOjt9BIMBfD4ftbW1SSHiJFKfJhR1Zck2NTVx8cUXH/XYO++8k7a2tu77n/3sZwc19cFy77338sgjjxz1mPXr1/P73/9+0Nfqi4aGBmbNmhX36wiCkLxo7XjVDh2CpqbQmaFNTXD4sAmNFhSYYgO75eXBUVpkCsOMqqpqsrNzyMzMJCcnJ2la34gnDti+3XwLGwwej+quNsrLgylTBm9XIBA4aqPdSEycOJGnnnrqqMfceeedXHnllRQUFADETFRdd911fR6zfv166uvr+exnPxv1ef1+P1lZ8lIVBKF3bCjU3cbD63X2y8xQYTBUVlZRU7Mm6XLixBOHEXDub2CD3foShA0NDcyYMYOlS5cyc+ZMLr744m7P2NSpU7nxxhtZsGABTz75JB999BHnn38+J598MosWLWLz5s0AbNu2jaqqKmbPns2//du/hZzbeqACgQDf+c53mDVrFnPmzOHuu+/mrrvuoqmpibPOOouzzjqr+5r79+8H4Gc/+xmzZs1i1qxZ3Nk1k6yhoYGZM2dyzTXXUFFRwbnnnhtxdurNN9/MHXfcARjv3o033sjChQs54YQTeP311/H5fPz7v/87NTU1zJs3j5qaGlpbW1m+fDkLFy5k/vz5PPfccwCsWrWKCy+8kLPPPpvFixdz2WWX8eKLL3Zfa9myZTz11FM0NDSwaNEiFixYwIIFC3jrrbf6/fcXBCF1CAaNOPN44MABM4Jq61b46CPzhXz3bmhpMaLN/b5cWGhCo1IpKgyUysoqrr32e0kj4EA8cUPGBx98wAMPPMDpp5/O8uXL+dWvfsV3vvMdAEaPHs26desAWLx4Mffeey/HH388b7/9Ntdffz1//vOf+eY3v8lXv/pVrrrqKu65556I17j//vtpaGhg/fr1ZGVlcfDgQUaNGsXPfvYzXn31VcaMGRNy/Nq1a3nooYd4++230Vpzyimn8OlPf5rS0lK2bNnCY489xm9+8xsuvfRSnn76aa688sqj/o5+v5+//e1v/P73v+eHP/whr7zyCj/60Y+or6/nl7/8JQD/+q//ytlnn82DDz5Ic3MzCxcu5DOf+QwA69at47333mPUqFH87ne/44knnuAf/uEf8Pl8rFmzhl//+tdorfnTn/5EXl4eW7Zs4fLLLyd8QocgCKmH1mazrTva2ox3ze93Cg0yMox3TQoNhOGKfCcZIiZPnszpp58OwJVXXskbb7zRvW/JkiUAeDwe3nrrLS655BLmzZvHtddey65duwB48803ufzyywH40pe+FPEar7zyCtdee213KHLUqFFHtemNN97gC1/4AoWFhRQVFfHFL36R119/HYBp06Yxb948AE4++WQaGhr6/B2/+MUv9nn8yy+/zG233ca8efOorq6mo6OD7du3A3DOOed023zBBRfw6quv4vV6eemllzjzzDPJz8+ns7OTa665htmzZ3PJJZewcePGPu0SBCG5sO07Wlpg716Tr7Z1q8ln27HDPNbR4TTJLS42PwsKTHhUBJwwXBFP3BChwt513PcLCwsBCAaDlJSUsH79+qjOEU9yXSVamZmZEcOpvT0nMzMTv98f8RitNU8//TQnnnhiyONvv/129zoA5OXlUV1dzR//+Edqamq47LLLAPj5z39OWVkZ7777LsFgkDxJdhGEpMX2XHNPNOjocCYagCkmyM42Ai0jwwg2QRAiI564IWL79u3U1dUB8Oijj3LGGWf0OGbEiBFMmzaNJ598EjCC59133wXg9NNP5/HHHwdg9erVEa9xzjnncN9993ULqIMHDwJQXFxMS0tLj+MXLVrEs88+S1tbG62trfzud79j0aJFg/xNQwm/9nnnncfdd9+N7qrt//vf/97rc5csWcJDDz3E66+/zvnnnw/A4cOHmTBhAhkZGfzv//4vgUAgpvYKgtB/wqtCd+1yqkIbGkxVaHOz8cDl5zuetaIicz8rS7xrghANIuKGiBNPPJF77rmHmTNncujQIb761a9GPG716tU88MADzJ07l4qKiu7E/1/84hfcc889zJ49m8bGxojP/cpXvsKUKVOYM2cOc+fO5dFHHwVgxYoVnH/++d2FDZYFCxawbNkyFi5cyCmnnMJXvvIV5s+fH8PfGs466yw2btzYXdjwgx/8gM7OTubMmUNFRQU/+MEPen3uueeey2uvvcZnPvMZcnJyALj++ut5+OGHmTt3Lps3bw7x3gmCEF9sg9z2dtOqY88eI9K2bDE/d+6E/fsjh0Lz8yUUKgiDRWl3d8M0pLKyUocnum/atImZM2d2349NixFPd6+2vlqMNDQ08LnPfY73339/cBdNYlpaWihOoThI+GsiHtTW1iZNb6FkQNajJ8m8JoGAM37K3cLDFiC4x0/FqiPQhg21VFRUx+ZkaYKsSU8SuSYeD5xwQnyvoZRaq7WujOZYyYkjNj3dWlq05G4IgpDyhM8KtRMN3LNC7XD3ggLxpAnCUCIibgiYOnVqWnvhBEFIfoJBpzmuz+eItUgtPPLypL+aICQjIuIEQRDSGJu3ZqtCbSjUPdjdhkKzs2WagSCkEiLiBEEQ0gCtnbw192B3n88Z7G5DoTLYXRDSAxFxgiAIKYZ7TqgVa+4ig4wMs9l+a5K3JgjpiYg4QRCEJMVdEeoWa+HNcaXIQBCGJ5KqmuKsXLmSiooKVq5cyb333ssjjzzS45iGhgZmzZo1BNbBaaed1ucxd955J21tbXG3ZdmyZTz11FNxv44g9JdAwIizYNAMdW9sDB3q3tRk+rAFApGb42Zni4AThOGIeOKAZ54xTSoHg9eb1Z1jUlYGXWND487999/PwYMHyczMTMwF+8lbb73V5zF33nknV155JQUFBVGfNxAIJO3vLAi94W7fYcdOeb3mPph9hw5JRaggCNEhbxEYAVdePrht4kTdfTsaQfjII490T1KwA+wbGho4++yzmTNnDosXL+4eBL9s2TJuuOEGTjvtNI499thub9KFF16Ix+Ph5JNPpqamhptvvpk77rgDgLVr1zJ37lzmzp3LPffc033dQCDAypUr+dSnPsWcOXO47777AKfJ6MUXX8yMGTNYunRp9yisd955h9NOO425c+eycOFCWlpaej1POLYBcm/nv+uuu2hqauKss87qniDx8ssvU1VVxYIFC7jkkkvweDyAac1y4403smDBAm6//XYWLlzYfZ2GhgZmz54NwI9+9CM+9alPMWvWLFasWEG6N7QWko9g0Iiz1lYjypqaQsdONTbCwYOm6CAnx/GqZWQ4Q91FwAmC0BfyNjEEbNiwgVtvvZU///nPvPvuu/ziF78A4Bvf+AZf/vKXee+991i6dCk33HBD93N27drFG2+8wQsvvMBNN90EwPPPP09+fj7r169nyZIlIde4+uqrufvuu7tnrVoeeOABRo4cyTvvvMM777zDb37zG7Zt2waYuaV33nknGzdu5OOPP+bNN9/E5/OxZMkSfvGLX/Duu+/yyiuvkJ+ff9Tz9Eak899www1MnDiRV199lVdffZX9+/dz66238sorr7Bu3ToqKyv52c9+1n2O0aNHs27dOm666SZ8Pl/3NWtqarrX4Otf/zrvvPMO77//Pu3t7bzwwgsD+TMJQp8Eg0aIWbHmnhH6ySdm7NSBA0bQWbFmQ6EFBaZCVMSaIAgDRcKpQ8Cf//xnLrnkEsaMGQPAqFGjAKirq+OZZ54B4Etf+hLf/e53u59z0UUXkZGRwUknncSePlx9zc3NNDc3c+aZZ3af66WXXgKMl+u9997r9uYdPnyYLVu2kJOTw8KFCykvLwdg3rx5NDQ0MHLkSCZMmMCnPvUpAEaMGHHU80ybNq1XuyKd/4wzzgg55q9//SsbN27k9NNPB8Dn81FVVdW93y1WL730UmpqarjpppuoqamhpqYGgFdffZWf/vSntLW1cfDgQSoqKvjHf/zHo66ZIByN3hrjdnaaXDStTYGBbd8hvdYEQUgEIuJShFxXU6fBhAe11tx9992cd955IY/X1taGXCMzMxO/TdTpx3mORjTn11pzzjnn8Nhjj0U8h3vA/ZIlS7jkkkv44he/iFKK448/no6ODq6//nrq6+uZPHkyN998Mx2DHYwrDBv6EmsgjXEFQUgexJE/BJx99tk8+eSTHDhwAICDBw8CppLz8ccfB2D16tUsWrRoQOcvKSmhpKSEN954o/tclvPOO49f//rXdHZ2AvDhhx/S2tra67lOPPFEdu3axTvvvAOYwfZ+v7/f5zkaxcXFtLS0AHDqqafy5ptvsnXrVgBaW1v58MMPIz7vuOOOIzMzk1tuuaXbQ2cF25gxY/B4PFKNKkSktzDoli1Oztr+/UbEZWeHVoMWFpowqNTVCIIw1IgnbgioqKjg+9//Pp/+9KfJzMxk/vz5rFq1irvvvpurr76a22+/nbFjx/LQQw8N+BoPPfQQy5cvRynFueee2/34V77yFRoaGliwYAFaa8aOHcuzzz7b63lycnKoqanhG9/4Bu3t7eTn5/PKK6/0+zxHY8WKFZx//vnduXGrVq3i8ssvx+v1AnDrrbdywgknRHzukiVLWLlyZXduXElJCddccw2zZs1i/Pjx3WFgYXjirgYNnw8KEgYVBCG1UeleuVdZWanr6+tDHtu0aRMzZ87svh+bFiPt5ObmA4ltMZKstLS0UFxcPNRmRE34ayIe2ApdwRDL9QgEnCkGPp8Rah0dkcVaVlbyetE2bKiloqJ6qM1IGmQ9eiJr0pNEronHA734FGKGUmqt1roymmPFE0dsBFdLi58U0iyCkJJYsWZz1tx91uww94wMI9JycsSzJghCeiMiThCEpMM9bsrrNV4161mzBQYZGcarJmJNEIThiog4QRCGDPcgd/cEg2DQOcaGQXNzzYgpQRCEwVBfX0ddXS1VVdVUVlb1/YQkZtiKOK01SoYNCgyuZYvQN1o7Yi18kPtHH/Uc5J6fL3NABUGID/X1dSxZspjOTh/Z2TnU1KxJaSE3LEVcXl4eBw4cYPTo0SLkhjlaaw4cOECexOMGjdZOCLSz0wmB+nyOZy0jwwmDZmQgeaSCICSUurpaOjt9BAIBwEddXa2IuFSjvLycnTt3sm/fvpids6OjQ4SAi1Raj7y8vO5JEkLfuBvidnY6laBdLQO7CwyyssxWUCCeNUEQkoOqqmqys3MA44mrqqoeapMGxbAUcdnZ2UcdDzUQamtrmT9/fkzPmcrIeqQ+R6sEtSjl9FhzDeQQBEFISiorq6ipWSM5cYIgpD5ah/ZYc+eruYsLpG2HIAjpQmVlVb/Fm9cLzc1w8GD8+8T1BxFxgjAMsMUFNgTqzldz13XYStC8PCPcBGG4kE4Vi8Lg8XrNSL59+2DHDjhwwJkAs3jxUFvnICJOENKISGOm3PlqFlsJKvlqgpB+FYu6NmglAAAgAElEQVRC/+noMKJt714j2rpGmpORYWYmjxtn7m/fPnQ2RkJEnCCkIO6WHbYZbqR8NVtcIPlqgtA76VaxKPRNR4cRanv3GmF26JAz8aW42IzPDP+Cm4zdqI4q4pRSo6I4R1Br3RwjewRB6MLdX80dArX5au7+apKvJggDJ90qFoWetLcb0bZnj/G0NTc7hVlFRTB+fGpGJfryxDV1bUf71TKBKTGzyIVSqgFoAQKAX2td2SUsa4CpQANwqdb6UDyuLwiJIHx4u4RABSGxpFvFogBtbca71tRkRFtLixOdKCyECROG2sLY0JeI26S1PmqfCKXU32NoTyTO0lrvd92/CVijtb5NKXVT1/0b42yDIAwK64a38z+tUOvoMCIOQoe3SwhUEBLLQCoWheTB4zGibdcuI9paW837bk6O8bSli2gLpy8RF80rOtGv+s8D1V23HwZqEREnJAm2EW6kqQVer5MUK1WggiAIA0Nr41mznrZAANauNfuys41oGy7TYI4q4rTWHfa2UqoUmOx+jtZ6nfuYOKCBl5VSGrhPa30/UKa13tW1fzdQFsfrC0JEoi0scIdAbZWTIAiCED3BIBw5YnLampqgsdG834KJWEyebHLahiMqmuHfSqlbgGXARxhhBaC11mfHzzRQSk3SWjcqpcYBfwK+ATyvtS5xHXNIa10a9rwVwAqAsrKykx9//PF4mgmAx+OhSD6hu0mX9dA6dAsGI1coKdV3nlpHh4e8vNRfk1gh69ETWZNQZD16MlzWxOYK+3wmqmHfd+38Zff7bWamh0AgMWvi88U/NHvWWWet1VpXRnNstC1GLgWO01r7Bm5W/9FaN3b93KuU+h2wENijlJqgtd6llJoA7I3wvPuB+wEqKyt1dXV13G2tra0lEddJFVJpPXrzqrkLC9xetaysgRUWbNhQS0VFdczsTnVkPXoiaxKKrEdP0nFN/H44fNg01N2xw+S12Ykx+fmmECGrF7WydWsd27atZtq05UyfHt/sLq1NSszll8f1Mv0iWhH3PlBCBMEUL5RShUCG1rql6/a5wI+A54EvA7d1/XwuUTYJqcvRctXc46Vsrlp2trTrEARBiAc+n2nxsX8/7NxperXZ9+GCAhg92rwP98XWrXXcfvti/H4vWVmrWblyTdyFXLIRrYj7MfB3pdT7gNc+qLW+MC5WGcqA3ynj8sgCHtVa/0Ep9Q7whFLqn4BPMF5CQegxBzQar5q06xAEQYgvdhrCvn1GtO3v6jehlPGyjRkzsAKvzZtNk2atg/j9PjZvrhUR1wsPAz8B/g8I9nFsTNBafwzMjfD4ASCJJpcJieZoXrVIc0ClCa4gCELisD3awhvr2uKuSNMQBsKMGaZJs/HE5TBjRvXgT5piRCvi2rTWd8XVEkFw4Z5WEO5Vs33VQLxqgiAIQ4lt99Hc7FSOejyhjXXjNQ1h+vQqVq5cw7ZtDyYkJy4ZiVbEva6U+jEmH80dTl0XF6uEYcPRphVYr5oVapmZppxc+qoJgiAMDcGgEWy2R1tTk9PuIycn8dMQpk+vYv58Ly0tw0/AQfQizk5tONX1mAbi2mJESA/Cw5/WqxZeVCDTCgRBEJKLzk5TObp/v/GyhVeOFhdDaenRzyHEj6hEnNb6rHgbIqQ24eHPQMBpyGgb4GodGv7Mz5fwpyAIQjJhixBs5eiBA0a0KWXes6OtHBUSQ1QiTilVAlyFGTrvnthwQ3zMEpIVt0g7WvjTet2kqEAQBCE50drMGG1uNh62nTvNZARbhDCYylEhMUQbTv098FcSWJ0qDB02Ty0QiNxTzXrPjhb+zMgwAk4QBEFIDoJBExq1+Wzu8VXZ2YnPZ0sV2tqMwN250xRtJBPRirg8rfW342qJkFBsnpoVauHVn1aoSfWnIAhCamKb6h44YESbO58tL8+0+0iWfLatW+vYvLmWGTOqh6zKtLPTrNPOnUbgWuF26JBzzDHHDIlpvRKtiPtfpdQ1wAuEVqcejItVQkwIz1Oz4U9381sryqxQy8sT17kgCEIqYkOje/YYEWLFR7Lns9nJC52dPrKzc+I+eSEYNFMirEizgm3PHictKCsLJk6EGTOgvBwmTTI/W1riZtaAiFbE+YDbge9jqlLp+nlsPIwSoid8SoHP53jVemvTIXlqgiAIqY07NLprl9na2sw+258tVk11440zeSEQ08kLWpv1cXvVGhuNt83tyBg3zoi0hQuNUCsvN4+FC16tUzec+i/AdK31/ngaI/SODX1G8qhZoaa1I9SkTcfRqa+vo66ulqqqaiorh2d/IUEQBk+i3ku8XiPaOjrglVdg924nNJqba0KjI0bE7fJxxZm84Bvw5IXW1lCxZgWbFbYAJSVGrJ19tiPWJkxI7c/KaEXcVqCtz6OEQeH2qIX3U3OPk8rIMJvkqQ2M+vo6lixxXPc1NWtEyAmC0G/i9V5ipyAcPmzEWlOTM7pq6lTzeLKGRgeCnbwQTU6cz+fkrTmetSoOupK78vONQDvlFCcMOmmSEbrpRrQirhVYr5R6ldCcOGkx0k8iVX56vaGNb7UOrfwUoRZb6uqM6z4QCAA+6upqRcQJMSPRXl7xKg8dsXovsQ11Dx40omT3bvNYpNFV2dmmwW66MX16VYh4CwQi563t3dszb23WrEOUlY3vFmujRg2fz8xoRdyzXZsQBYGAI9asR81ukSo/MzOl8W0iqaoyrnsw356rqqqH2iQhTUi0l1e8ykPLQN5LbG+2w4eNIGlsJMSLVFAAI0cagTIc0NoRr26x1tTkNIq3eWvWuxaet1ZcvJmWlvFD+4sMEdFObHg43oakMlaw7dhhhFowrJOeTChILiorq6ipWSPeCyHmJNrLmype5XT1FkbzXuLzmQa6tjfbrl3mMXAiLalSgDBYPJ6eYm3nTmhvd44pKTECbeZMR6xNnCh9R3vjqCJOKXW/1nrFYI9Jd2weWzAoLTpShcrKqrT6MBGSg0R7eVPBq5zu3kL3e4k7l23PHieXzZKfb4oP0t3L5vU6EyDc4dDwtSgvh1NPdXLW0jVvLZ709VK6SCnVcZT9CpC5qji5CoIgDF8S7eVNBa9yqngLB0JHh/GyHTxoBNuePaG5bOnuZQsEnJ50brEWKW/tpJNC+62VlqbvuiSSvkTcyijO8XosDBEEQUgHEu3lTXavcip4C6PB7zeCzTbT3bXLeN2Ucprppmsum81bCy8y2LUrNG+trAwmT46ctybEh6O+3CQXLn6ka46IIAiCm1TwFoZjm7oePgz79hmxcuCA413KzTVetnScM+rx9BRrjY2heWulpUagzZrlhEElb21oSMPvDMlPuueICIIguEl2b2FbmxMW3bXLeNqsh8mGRceOTa98Z6+355zQ8Ly1wkIj1qqqQvutFRYOnd1CKCLihoB0zhERIlNfX8ezz67mooty5W+dQojHPP2weWzNzY5ga28PDYuWlKRPWNTmrYVXhO7b53gWs7ONOKuoCM1bKymRvLVkJ6qXqVLqEq31k309JkRHuuSICNFhPa8+n5fHHlstntcUQTzmqY/XGzr5YPduEy60gi0vz3iVSkqG2tLBE563ZkVbeN7a+PEwZQqcdpoj1saNSy8v43Ai2u8a3wPCBVukx4QoSMUcEWHgWM9rMBiks1M8r6mCeMxTC9uPzZ3H1tLi7M/LM2HRdJh20NISeU5oh6uXxKhRRqTNmhU6J1S6KKQXffWJuwD4LDBJKXWXa9cIwB9Pw9KdZM8REWKH9bxq7RXPawohHvPkpaPDCJkjR0yo0O1h09opPEj19h5eb6RJBlUR89ZOO80Ra5Mmmd9fSH/68sQ1AfXAhcBa1+MtwLfiZZQgpBPW8/rssw9y0UXLRbynCOIxTw7a20Mb6O7ZY8ZWhQu28Sk8dcnv75m31tgYOW9t3jwzJ1Ty1gTou8XIu8C7SqlHtdadCbJJENKOysoq8vO9VFSIEEglxGOeOOxM0ZYWU3Swd6/pL/bEE45Iyc83W6qGRLU2rUrCxZo7by0jwwjSY45xvGuTJjl5a8N5TqjQk2hz4hYqpW4Gjul6jgK01vrYeBkmCIIgpCd+vxFrHg/s3288Tvv2mUpKcKpEJ05M3V5sR470DIWG562NHm0E2uzZTpGB5K0J/SFaEfcAJny6FgjEzxxBEAQhnWhvN2LtyBHjXdu3z3jarHctM9MIttGje1ZIpkKYsKOjZ7+1nTvN72uxeWunnx7ab03y1pKLYNCZhd7Z6Wzu12Wy9ciLVsQd1lq/FFdLBEEQhJSls9OINY/HtLqwHjav1xFjublGsEVTcLB1ax3btq1m2rRcpk8f+pC2328KKMLnhO7b5xyTk2O8h3PmhBYZjByZGoI0XdG6pzBzizObd6iUEWmFhTBmDBQVmdB9Xp752+bmOmHvZKGv6tQFXTdfVUrdDjwDeO1+rfW6ONomCIIgJBnBoJlwYKtDbSjUzhHV2vGujRgxsKa5W7fWcfvti/H7vWRlrWblyjUJE3LBoMlbC2/hsXu3E+61eWtTp8IZZzhiLd2mOqQCfr/5mx054ogzcESz1k54vqjIEWdFReYxK85yc83tvsS2xxPf36e/9PXv9d9h9ytdtzVwdmzNEQRBEJIBrY1Ys961AweMd+3gQbPPfjjm5Zktlu08Nm82Pfq0DuL3+9i8uTYuIu7IkZ5zQpuaeuatlZfD3LmOWJO8tfgTDPb0nAWDoa8xW51cUmIEtPWc5ec7wsyKs3QV131Vp56VKEMEQRCExGPFmq0MPXjQeNasWLPH5OQYsTZmTPw/EGfMMD36jCcuhxkzqgd1vo6OyEUG7ry1oiIj0E4/3RFrkrcWH/x+05zZLdDsFA1wvLmFhSYUbcVZQYEjzPLyzM/MTNiwwYwMG45EO3br2xEePgys1Vqvj61JgiAIQqwJBo1Qs5vNWzt0yAkTJlqs9cb06VWsXLmGbdseZNq05VF74WzeWqR+a5acHCPO3Hlr5eUm9Ct5a4PDFga4BVq49wyM+CoqMlMliovN2ltR5vaeCX0TbbZCZdf2/3Xd/xzwHnCdUupJrfVP42GcIAiC0D9sgUFbm/E0HThgBNvhw2a/9a5ZsTZqVHKGmqZPr2L+fC8tLT0FnM1bC68I7Stvrbx8aMVpKhMe3vT5nJC63Z+Z6Yw2s96zoqLI3jMhNkQr4sqBBVprD4BS6j+AF4EzMW1HRMQJgiAkiGDQhAhbW41YO3TICLWDB819CE3ozs1N3aR7d96aFW1Hy1vLzt6Cx/MXKisrmDHj1KEzPIXoTaCBeQ1ZgVZUZMKbI0aY2wUFjjCzP8WbmViiFXHjcFWlAp1Amda6XSnl7eU5giAIwiDo6DCizFaDWq/akSNOcQGYD9i8PKciNBVpb3f6rVmx1th4Wo+8tfLy0IrQSZPM7w1OVWtnp4/XX89JaFVrsmLba/h8TpjTvm7s/owMJ//MCrTCQqdoJdrKTSHxRCviVgNvK6We67r/j8CjSqlCYGNcLBMEQRgG+HyOUPN64W9/M0Lt0KGe7RJsv6pUDgn6/WbMVHihwf79zjG5uUacVVYeoKxsQki/taPhVLUG4lrVmkwEAo448/mcPmZuwVVYaITZ+PFmDQsLHe+ZeNBSm6hEnNb6FqXUS8DpXQ9dp7Wu77q9NC6WCUISUl9fJwPRhX5jPWrt7SYEeuiQ2ZqbQ5vhTp1qBE5urvmwHUiPtWjYurWOzZtrmTGjOm4iJxg0wixcrLnz1jIzjbA47jg480zHu2ZFanHxB7S0RD93y6lq9cWkqjUZsOLMbnl5sGeP2ae1eY0UF5tw+YgRZsvPdzyzubmpK/iFvumr2e8IrfURpdQo4OOuze4bpbU+GG8DBSFZqK+vY8kSE6rJzs6hpmaNCDkBcHLU2tud3mpWpDU3G++IzS3KyHAajIY3w83OhtLS+NrqDjlmZ8cm5Hj4cM9JBo2NRqBaxowxIm3evNB+a7EUqraqNd4CNVZo3VOkudtsgBFiI0c6XjSvF847z/GiSRXn8Kavf59HMZWoazHNfVXYz2Pjap0gJBF1dSZUEwgEAB91dbUi4oYRPp8j0trbjXCxIs3dxd1++NqKvNLS5KrGG0zIMTxvzQq2lhbnmOJiI9IWLXIqQidOdPLW4s306VVJI95sPprX64Q8wyksNGtWUmJEvc1Fs9608NfOhg0wblxi7BeSn76a/X6u6+e0xJgjCMlLVZUJ1YDxYFRVVQ+1SUIM8fuNSLEeNY/HEWlHjpgPYTdZWY5HLZUqP6MJOdq8tfBpBgcOOMfYvLX5852h7rbf2nDCFg1YoWaxYr6gwHjQSkqcfDS3SEuV142QnETb7Fdhct+mdeXHTQHGa63/FlfrBCGJqKysoqZmjeTEpSiBQGSRduSI2drbndmf4IQ9c3IGPgM0GXGHHE88sZoRI6pYty5UrO3Z0zNvbfp0qK52BNvo0cNDgNj2G1ak2dC42+M6cqQZO1ZS4szk7M2TJgixJNq3pV8BQcys1FuAFuBp4FNxsksQkpLKyioRb0lKZ6cj0Do6jEg7fNgZ1G77p4GTd2RFmvWWpCtah88JraKxsYoXXgjNWxs71oi0BQucvLXx49NHwPaGrfD0ekPXw1JUZERrSYnZrEgrKJAZqsLQEu2/5ila6wVKqb8DaK0PKaUknVIQhISgtflw7ehwhJoVZ3YLTwrPyDAfsDk55gO3uHh4tFFob+85dmrnztC8vREjjEBzV4ROmmQ8R+lKIOCINLdQ09q8TkaONCK2tNRpZCshTyHZiVbEdSqlMjHFDCilxmI8c4IgCIPGetHsZkdG2a211YS1wAllZWamZ7gzWjo7nbw1Ww0anreWl9fTs5bOeWvBoCP2rVCzr5esLONFKyszQq242PGmSZ80IVWJ9m3vLuB3wDil1H8CFwP/FjerUoy//rWOxx9fjdebK6E2QQgjEAgVaB0djufsxReNR83tRQNnEHtOjvGSpHJz28ESDJoB7uHetd27HWGbmWnaddi8NVtkkKxzUQeD1k7os6PDWQO3UJsyxfzuNj/NjocShHQj2ma/q5VSa4HFmPYiF2mtN8XVshShrq6OCy5YjNfr5fHHV0vvMGFYYUNU4QLN4zE/W1pCm9naZPCsLPNB6/cPTy9aJLSGQ4dy2Lw5tOdaU1No1eNwyVvz+51B9+GTK4qKnNBnSYmp+LRCTTxqwnCiP//2W4Aj9jlKqSla6+1xsSqFqK2txefzEQwG6eyU3mFC+mD7W7kFms1D83gcgRZOVpbxnmVnm5BVb81rMzMT1zss2Whr69kcd/v2TtrbT+s+ZsQII9LcFaETJ6afR6mz0+Txeb3mS4Ftimx77B1/vPlZWOiINan4FARDtC1GvgH8B7AHCOA0+50TP9NSg+rqanJycvB6vQnpHSZjn4TB4g5H2a29PVSceTzmcRuKc3vQohFogsGdt+YOhx50zbrJy4PRo1vwep9EqXfJzPyAr3/9FubOTZ/if/uas18G3APYCwpM6HP06FCxlpdnGttWVAyd3YKQ7ETrifsmcKLW+kCfRw4zqqqqeOmlNaxa9SBf+MLyuAorGfsk9IU7vGkFWltbqAfNFgm4e10pFSrQJMTZP4JB2Lu3Z0Xonj0989ZOOCG0Oe7o0fDii7/kmWd+gNYBgsFMduz4dEqKuN7EmlImBFpWZvIbR4404q2wUMZGCcJgiPZtegdwOJ6GpDKnnlpFR4eX2bPjK6hk7NPwxbbYCPeeWa+Z3dz5Z/Z5ttWG3YZzkcBg0do0CA4Xa5Hy1srL4eSTHbFWVta7MHamKHhTZnC7DYO6iwugp1grKjJiTb4UCELsOeq/lVLq2103PwZqlVIvAt1ZMFrrn8XRtqPZdT7wCyAT+B+t9W1DYUeikbFP6Ul4aNPrNeKstdURZ62tRkCECzS3OJPwZmyJlLe2c6f5W1hGjjReNXfe2qRJJp+rP9gpCtu2Pci0acuTZvYnhE66cE8rKCgwXsSxY53igqIiaX4rCImkr+9GxV0/t3dtOV3bkNHVr+4e4BxgJ/COUup5rfXGobQrEcjYp9TCjuqxm88H778f6j1zhzbBCT9lZjrizE4UEO9ZfOjsNJ40t2DbuRMOHXKOyc834qyy0vGsTZpkhHOsmD69ivnzvbS0DM3/tQ2F2iID+1hWlvGqTZsGTU11bNhQy5lnVnPqqfL+IwhDzVFFnNb6h9GcRCl1t9b6G7ExqU8WAlu11h93Xftx4PNA2os4kLFPyUC4OPN6jRhze8/czWnBfBhOnQrbt0toE2Dr1jq2bVvNtGm5CfM62by1cM/anj2hhRsTJsCMGaF5a6NGpVfrCr/fmXwRDDq/v62IHTvW3LYNcZUyObk33GBycu+7T3JyBSEZiFWWwukxOk80TMLk6Fl2Aqck8PpCGmK9EOGhzfCQZiRxBj3zziINB8/OhnHjEvc7JStbt9Zx++2Lu/K/VrNy5ZqYCjmbt+YWa42NoXlrSjl5awsXOp61srL0a1/R2WlCw+5Cg6ws81qcPt0I1MJCI9iOlrcmObmCkHykZaqpUmoFsAKgrKyM2trauF7PJJ172LAhvtdJJTo6kmc9tHa8DfZ2MGhyfexttzfCjVJms5V0g/HGZGZ6KC6uHfgJ0oRt21bj93sJBoP4/V62bXuQ+fMjNJyLgtbWLLZvLwzZduwoxONxErNKS71MntzK7NmtTJnSyuTJZsvNTb7JgYN9jdjXsfuLRkGB8fhmZRmBmpkZ+gXj8GGz9cX48SPIyspCa01WVhbjx4+I+/94Mr2PJAuyJj1J5JoEg+YLYbKQiiKuEZjsul/e9Vg3Wuv7gfsBKisrdXV1dVwN8npN09+KivheJ5XYsCG+6xE+zNrm8liPmcdjvA9tbeZ4m4xtG4laj5ltq5GVFf+wZnFxLS0t1fG9SAowbVouWVmruzxxuUybtrzPPDCbt+YOg0bKWysvd/LWbDi0qCgXyAVGdR9rva7JRrSvkWDQCYd2djpfLkaONB62sWOdytBYNQeuqKhm2rQFCc3Jjff7SCoia9KTRK6Jx2PaBCULsRJxicwWeQc4Xik1DSPeLgOuSOD1U5JUaBKstZNv5g5rtrX1FGeRZm1CqDjLzzchonTKZUoHjlaJGW3e2sSJJm/NLdbSLW/NYgVbW5vJZQPzhWPUqJ75a/Fu4yE5uYKQXPTrX14pVaC1bouw6xcxsqdPtNZ+pdTXgT9iWow8qLXekKjrpyJD3SQ4EHAagFoPSEdHaJ5Za6v5kLKVmtZjBqH5ZllZ0og2HTjuuCqmTNF88EEVL73kCLamptA5mePGGZFm89bKy81j6Za35sZ6kMMF2+TJRrAVF5v/gXReA0EQoiPasVunAf8DFAFTlFJzgWu11tcDaK1Xxc3CCGitfw/8PpHXTGXikZAcqRDA53O8ZlrD8887XjPo6SWRNhrDg9bWnv3WGhuhtdWZE1pSYsTa2Wc7Ym3ChP73W0sltHY8bFa45ueb/4UpUxzBVlwsgk0QhMhE68/4OXAe8DyA1vpdpdSZcbNKiCn9aRLs90f2mrnbZ7S1mVwcd/PZ8CrNSZPMffGaDR98vp55a42NkfPWFi6EY4/9kLFjT2DSJJO7le7Y1ACvq4bDDngfN878rzQ2wqJFQ2ejkJ6kQjqNMDCi/njVWu9Qoa6UQOzNEeLBggVVPPLIGt56q5Z586oZP76Kjz6KHM60Hdkt7lwzWwSQn28+dI/mNcvMjF1CtZBcBAKR89b27u2ZtzZzZugkA3feWnFxEy0tSZQhHENsW4/2dqeoprgYjjnGCLaSEnM/fLpBY2Pk8wnCQBnqdBohvkQ9O7UrpKqVUtnAN4FN8TNL6As7S9Md0nSHM93CzOeDYLCKsrIqdu+G3budCs2sLEeclZZK2EZw0BoOHow8J9Tma9m8tfJyOOWU4ZO35saOpbKj0bQ27WjGjzdbSYnxsqVzaFhIXhLd30+8foklWhF3HaZ4YRKmIvRl4GvxMmq4Yqszw/PM7AeEu22Gu3UGOAUBNpxpxVkiKtaE1Mfj6SnWdu40rz1Laanxpp10kuNdmzgRtm+vY/PmWmbMqE6qmZ/xIFIeW2amaRI8Y4ZZo5IS460WhN5IpNBJ5Mxt8folnqg+3rXW+4GlcbYlbYkkzNzVmW5hFgg4YUobmlLK8ZbZcKa0zhAGgtcLu3aFzghtbDQTDiw2b+3UUx3P2sSJkfPW7PQF+6Yd6+kLQ01np/kf7egI7cV2/PFGuNnWHlKQI0RLooVOImduy1SPxBNtderDwDe11s1d90uB/9ZaL4+ncamCLQbYsiVUmFnPmXvAOUTOM8vNlepMIXYEAqa3WvhQ9337euatnXRS6FD30tLovyBs3mzetLUO4Pf72Ly5NmVFXDDoeL1te5v8fFMlO3GiExYNz2MThP4wFEInUf39Eun1EwzRBtrmWAEHoLU+pJSaHyebUo6DB0315oYNocJM2mYI8cbmrUWaE+rOWysrM20rqqocsRaLvLUZM8ybtt/vIysrhxkzqgf9OyUKn88INlstmpFhxlMdd5z5WVJi/n8FIZaks9BJpNdPMEQr4jKUUqVa60MASqlR/XjusEApk8QsCPHC4+lZEdrY2DNvrbzcVIW6Q6E5OfGxyU5fSPacuGDQSVmw80ULC00DXXfxgeSPCvFmMEKnvr6OZ59dzUUX5SatQJKpHokl2res/wbqlFJPYkZsXQz8Z9ysEoRhjNcbud+aO2+toMAItKqq0BYehYWJt3f69KqkE2/uXDYwXraxY83Mw9GjpfhAGFoGInRsLp3P5+Wxx1ZL0YAARF/Y8IhSai1wVtdDX9Rab4yfWYKQ/ti8tfCKUHfemm2cXFHhiLXyciNCpLDFoLXxRrpHVeXlmfUaP954J0eOHD4tT4T0xObSBYNBOjulaEAw9KfZ7wal1D4gD0ApNUVrvT1ulglCmmDz1j78cBRbtjiibSUaObAAAB0/SURBVNeu0Ly18eNN3tpppzmCbdw4yakMx+/v2WantNQUaIwd6+SyicgV0gmbS6e1NyG5dNLvLTWItjr1QkxIdSKwFzgG0+y3In6mCULq4c5bcxcamLy1OYCZWjBpEsyaFTonNFZVj1u3plffNq/XhEbtDN7sbCN458wxazlypFSMCumPzaV79tkHueii5XFv2Cv93lKDaD1xtwCnAq9orecrpc4CroyfWYKQ3Hi9jkBzh0MPH3aOKSx08tbKy+GEE9ZRWrogrhWPqd63LRh0eijaAoTiYjj2WKcAQXokCsOVysoq8vO9VFTE939a+r2lDtGKuE6t9QGlVIZSKkNr/apS6s64WiYISYDfH9pvzf50563l5JgK0NmzQ4sMwvPWiouP0NISX3tTrW+bOzRq12rMGCPabGhUChAEIbGkcxuUdCNaEdeslCoC/gKsVkrtBVrjZ5YgJBat4cCB0AKD8Ly1jAyn39ppp4X2W0uWvLVk79vm7s2Wl2c8lxIaFYTkQvq9pQ5HFXFKqVyttRf4PNABfAszfmsk8KP4mycIsaelpeckg8ZGpx0FGEFRXh6/vLV4kUx92+yc0dZWRwgXFsIxx5i1PHwYzjgjeQSwIAgO0u8tNejLE1cHLADu1Vp/qeuxh+NrkiDEBpu3Ft5vLVLe2umnh4ZCU7lT/1D1bXM31A0ETHi0pMQMhi8rMxWk7nXdsEEEnCAIwmDoS8TlKKWuAE5TSn0xfKfW+pn4mCUI0WPz1sLFWm95a+45oSNHSpL8QAlv9ZGRYULLM2Y4DXVzc4fWRkEQhHSmLxF3HSZ8WgL8Y9g+DYiIExJGMNhzTujOnSZvLRAwx2RkmByrY44x3jUr1saOFa/PYOnsNFWjdtZoVpYRxpMmOflsMrZKEAQhcRz1LVdr/YZS6i1gp9ZaxmwJCePIkZ6etfC8tdGjjUibM8cRa6mQt5YqeL1GtHV2mvt5eaGVtyNGiDAWBEEYSvr83qy1DiqlZFaqEBc6OiLnrR054hxTVOTkrVkRkep5a8mGLULweIxX0/Znmz7dGV1VWCihZ0EQhGQi2uDHGqXU/wOe0dpmGQlC9Pj9sHt35Lw1i81bs541u40YIeIh1gSDZopEa6u5DUaoVVSYvLbSUunPJgiCkOxEK+KuBb4N+JVSHYACtNZ6RNwsE1KSYND0W9u8eTRbtzr5a7t398xbmzrVtJiwYm3MGAnPxYtg0Ai29nZzOyPDhKOPP95pqitFCIIgCKlFVCJOa10cb0OE1MPmrbkLDZqabN7abMDJW5s71xFr48dL3lq8CZ+EoJRZ94oKp3JU/gaCIAipTVQiTil1ZqTHtdZ/ia05QjLizltzh0PdI6Rs3toZZ5h8teOPX8eoUQsSEpJLt4HvA8HvdzxtYKpEJ0ww4nn0aFM5mpk5tDYKgiAIsSXacOpK1+08YCGwFjg75hYJQ0Y0eWu5uUakzZvnFBlEyltLxJxQSP2B7wMlvN1HTo75W0yc6LT7kNC0IAjCwNHafC76/SYNJRhMvghGtOHUkB5xSqnJwJ1xsUiIO8Eg7N/fsyrUnbeWmWnCb9OmwaJFjmBLtry1VBv4PlB8PiPafD5zPy8PJk8Obfcx2OKP+vo6mZUoCMKwIxAwQi0QcAq9wHzW5eaa99e8PBPhSEkRF4GdwMxYGiLEHq1777dmPThghFl5eah3bcKE1GjcmuwD3weK7dHm8xlxVlgIxx7rtPsoKoptxW59fR1LljgezZqaNSLkBEFIG9xeNeussGRnG5GWl2eiGllZZkuFFJRoc+LuxkxoAMgA5gHr4mWU0H/a2x2B5i428HicY2zemvWslZeb8Fsqt5JIpoHvg8H2aPP7Q3u0TZjg9GiLJ3V1xqMZCAQAH3V1tSLiBEFIOYJBR6hZsaa18arl5Jj31txcI9yys41YS+UWVtH6Wupdt/3AY1rrN+Ngj9AHfr8ZMxVeZHDggHOMzVubP7/nnNB0ZKgGvg8Ud2Ndv988NnIknHii42lLdCPjqirj0QTjiauqqk6sAYIgCP3A7VULBh2HRVaW+Qy0Ys3tVUtlsdYb0ebEPWxvK6VKgclxs0gAnLy1cLG2Z0/PvLXp06G62ikyGD06ufLWhjtaG09pfr7JOwSTx3bSSVBWlhyNdSsrq6ipWSM5cYIgJA19hUDz8822Y4fpO5qVNfw++6INp9YCF3YdvxbYq5R6S2v9rTjaNixw562F567ZJHYwDVknTYIFCxyxNn58auStDTesaPN4nCTZUaPMm8055xjRlpc3tDZGorKySsSbIAgJJzwEaudCZWaaEOiIEU4I1HrW3F41GyodjkQrAUZqrY8opb4CPKK1/g+l1HvxNCwdsXlrbrEWnrc2YoQRaJ/+dOic0GT80BcMVrS1tpo3IKWMN3TOHGeEVW4ubNhgctwEQRCGI7YK1LbssEKstxCoOCn6JtolylJKTQAuBb4fR3vSgs5Ok7cWLtjceWt5eaaoYMGCnv3WhOQm0tzRMWNMWFtGWAmCMJzROlSsae2ItVSuAk1WohVxPwL+CLyhtX5HKXUssCV+ZqUOHR3whz/Ayy8fwyefGMG2e7fz4Z6ZabwvNm/NFhqMGjX8YvepihVtHo/zhjRmjJk7Om6cEW3D1ZUvCMLwJDxfzb43KmXeDwsLjVhzh0DlMy/2RFvY8CTwpOv+x8D/s/eVUt/TWv849uYlP8EgXHstaD0tJG/NirWyMnEJpxrBoJk52toaKtpOPNH8LC1NvoaPgiAI8cCdr2ar6ZWK3Ag3Ur6aEF9iJS8uAYaliCsogBdfhN27X0epRUNtjjAAwkVbRoYJi86YIaJNEIThQfjUAltckJVlRJrkqyUnsfozDGvdPWsWNDcHQqYgCP0nUYPs3aINzLdGK9psTpuINkEQ0pFAwORtu6tAtXby1fLzJV8tlYiViNN9HyIIvRPPQfYi2gRBGG6EN8N1FxfY/mpusSb5aqmJeOKEpCCWg+yPFh4V0SYIQrpgK0GtZ80+1ltxQXa25KulG7EScU/2fYgg9M5gBtmLaBMEIZ1xV4IGg9DSEloJWlTUs22HiLXhQbQTG/KAfwIqgO62s1rr5V0//ysu1gnDhv4Mso9UPTpunIg2QRBSm97GTCnlVIJmZ8OUKeZnus4DFaInWk/c/wKbgfMwPeOWApviZZQwPOltkH1vfdqkelQQhFQkkljT2ogyd9sOGwZ1i7WMjKGftSwkD9GKuOla60uUUp/XWj+slHoUeD2ehgnDF60dT5tNyB0zBk44wXjaRLQJgpAKaO3kq0Uj1qRth9Bfon3JdHb9bFZKzQJ2A+PiY5Iw3AifPQpGrMlEBEEQUoFwseYe4J6b6+SsSY81IdZE+1K6XylVCvwAeB4oAv49blYJSUM8erdZ0ebxOJ620aOd2aOlpSLaBEFIPnobNWWnF4hYExJNtGO3/qfr5mvAsfEzR0gmYtW7TWszY7alxZkpO3o0zJ3riDYZGC8IQrLQV87ayJEi1oTkINrq1DLgv4CJWusLlFInAVVa6wfiap0wpAy0d5vWRqzt2+d8Wx01CmbPNrNkRbQJgpAM9LfAQMSakGxE+5JcBTwEfL/r/odADSAiLo3pT+82m9Pm95s3wenToaLCEW15eb0+VRAEIa6E91mzaRzu1h0i1oRUJNqX6hit9RNKqe8BaK39SqlAX08SUpuj9W7r6DA5bVa0lZSYlh/jx8O2bXW89NJqLroolwkT4jcDVRAEwU2kCQbgNMUtLjbtOez0AumzJqQ60Yq4VqXUaLpmpCqlTgUOx8MgpdTNwDXAvq6H/lVr/fuufd/DNB0OADdorf8YDxsEB9u7zeuFAwfMm6PW5pvriSca0TZqlNO3qL6+jquuWozP5+Wxx1ZTU7OGykoRcoIgxJbwAgOLLTDIzw/1rIlYE9KRaEXctzFVqccppd4ExgIXx80q+LnW+g73A115eJdhpkZMBF5RSp2gtRaPYBzweo2nzeczb35FRSZEOnGiCY8WFER+Xl2dyaMLBoN0dvqoq6sVEScIwoAJBBzBprVTEZqdHTob1Ao2EWvCcCLa6tR1SqlPAydiht1/oLXu7ONpsebzwONaay+wTSm1FVgI1CXYjrSks9NUj/p85n5hIRx3nBFto0b1LtrCqaoyeXRae8nOzqGqqjpuNguCkD4Eg6G91pRyHsvPN5udDSqD3AXBEG116iXAH7TWG5RS/wYsUErdqrVeFye7vq6UugqoB/5Fa30ImAT81XXMzq7HhAHQ2Wk8bV6vuV9QAFOnwqRJxtNWVDSw81ZWVlFTs4Znn32Qiy5aLl44QRBCCAadMKg7by0ry4RCCwqMWMvOhqYmOFaaWglCryjtTibo7SCl3tNaz1FKnQHcAtwB/LvW+pQBXVSpV4DxEXZ9HyPU9mPy724BJmitlyulfgn8VWv9265zPAC8pLV+KsL5VwArAMrKyk5+/PHHB2Jm1BgvlgetB6h8EoBt+2H7tGVkmDdJ+802IyO21+vo8JCXl7zrMRTImoQi69GTdFsTG/4MJyPDbLZCtDevmsfjoWig3yjTFFmTnqTbmpx11llrtdaV0RwbbU6c/b70D8BvtNYvKqVuHZB1gNb6M9Ecp5T6DfBC191GYLJrd3nXY5HOfz9wP0BlZaWurq4eqKlRsWMHvPFGLV5vfK/TH/x+0/Kjvd28iebmwuTJMGWKCY8WF8c3HLFhQy0VFdXxu0AKImsSiqxHT1J1TWzOms1bc7fvsKFQd0VotNTW1hLv9+9UQ9akJ8N5TaIVcY1KqfuAc4CfKKVygRj7bgxKqQla611dd78AvN91+3ngUaXUzzCFDccDf4uHDalIIGDCox0d5k00O9uERsvLzXSEESMkh0QQhMHhLjKwvdbAePQLC52KUGnfIQiJIVoRdylwPnCH1rpZKTUBWBknm36qlJqHCac2ANcCdOXjPQFsBPzA14ZzZWowaERbW5u5n5lpRNvkycbTNnJk7EOkgiAMD2zemhVrlqwsI9QKCpxqUCkyEIShI9rq1DbgGdf9XcCu3p8xcLTW/397dx9jWXkXcPz7m50Zdtju7Oyyy0ApBdZlY6BWu6zIai0bsUCJEV/SlGBaYpvwh5hUY01oqn8ZghrT1JcmSiKkTazEN1KM1QropIkOLYsBykIpCyJli6JBhd2mO8vM4x/POd57587dzMzee889534/yc2997lv5/z2zDm/fV4/fIbX7gLuGsTvjrqVldw8evJkfj4xARdemNcfPe+8PNmuSZu0tiNHFllcXODQocMOtmlzpgXdZ2ZyDf4557SSNc8x0mhxcZERtbKSa9m+853WifWCC+Bd74Ldu3PStpG+JdK4OnJkkQ996DpOn15iamp6bCegbu+3Vmpfdmqz/dYkVcckbkSklBO2kydzAjcxAXv25KWszj8/J22u5ydtXDkB9fLyMtD8Cajb51trbwptnxy3HJXu5LhSvZkWVCSlPHL0xInWxJZ79sD+/fl+58580pV0dsoJqGGpURNQr24KLcumprrnWxvENEKSqmcSN0TtSVtKuS/bu98N8/M5aZuernoLpeYpJ6Cuc5+4Xk2hW7fmQUxbt9oUKo0jk7gBOnUqL2VVzp00NwdXXtlK2rZurXoLpfFw8OCh2iRvp051Lj0F3U2hTuEhCUzi+qp90XjInYX378+jSHftyh2HJanXqNCVlZyk2RQqaT1M4vqg/N/w0hLs29dK2ta7aLyk5lprgtyInKTNzubatXKC3G9/O8/3KEnrYRLXB+WUH9deW/WWSKrKWgMNoHuC3LJ2zaZQSWfLJK4PyskxpX5zktrRtLycp/FYa6DB3FweHVombJ4bJA2KSZw0opyktnoptZK1cs61lFr91mZmHGggqTomcdKIGrdJaqvWPo1HSrlsy5Zcu7Z9e+c0HtauSRoFJnHSiGrqJLVVW6t2DVrTeMzMdPZdk6RR5SlKGlFNmKS2au0DDcqEbXXtmtN4SKorkzhphNVpktoqWbsmaRx5OpNUK+0jQ8tJcicmrF2TNH5M4iSNpF6rGkxNdY8MtXZN0jjy1Cepcu2rGqSUb2XtWrnAu/OuSVInkzhJQ5NSqzl0ZSWvNQyuaiBJm2ESJ2kgVlY61wwtnXNOXjN0agre+c7WRLmSpI0xiZN01lYvQ5VSTsxmZnJzaNl3bWqqVbtWNpdKkjbHJE7Sup1psIFTeUjScHmalbSmlZXOudciWou8O9hAkqpnEiepq3YNWisbzM3lfmwONpCk0WISJ42RtZpDIfdZK5tDy/5rDjaQpNFmEic1VHtzaJmsReRaNZtDJan+TOKkBli9FBW0RofaHCpJzWQSJ9XMW2/lhG15uVU2OelSVJI0bjzNSyMqpVbtWpmwlc2hs7Od03nYHCpJ48ckThoBq6fzgNZkuNu3d/ZfszlUkgQmcdLQrV7dAOy/JknaOJM4aYDa+6+VC77bf02S1A9eOqQ+WD3/Wqm9/9rx47Bvn/3XJEn9YRInbVD7gIOy/1q5HNXcXKv/2uRkZ8IWYQInSeofkzjpDNaaMHdiItes7dhh/zVJUnVM4qTCWgMOJidzwnbuufZfkySNFi9HGkvrmTB3etr1QyVJo8skTo3XnrC1L/i+fXvnhLkmbJKkOjGJU2OsZ4SoKxxIkprCJE61tHpJqojWklQ7drQSttUjRCVJagqTOI281VN6pOSSVJIkmcRppKxeQ7ScW80pPSRJ6mQSp8qcaQ421xCVJOnMTOI0FKtr2KA1B9uuXc7BJknSRnnJVN+V/dZOnGiVOWmuJEn95WVUZ6Vc5aB9Dray+XN+3oRNkqRB8dKqdVudsKWUE7RylYOyD9uWLfDyy3kggiRJGgyTOK2pPWEr+7CVCVt7k6irHEiSVA2TOHUlbBGtdURN2CRJGk0mcWOmV8K2bVtr4XcTNkmSRp9JXIP1ahI1YZMkqf4qW1UyIj4YEUcjYiUiDq567ZMRcSwinouIG9rKbyzKjkXEncPf6tG1vAzf/S6cPAlvvplvp0/n5tD5ebjkEti3D/buzc9nZ/NyVSZwqpPFxUXuvvtuFhcXq94USapclTVxTwM/A/xRe2FEXAHcAlwJvB14OCL2Fy9/Fng/8ArwWEQ8mFJ6ZnibPBrONErUPmxqqsXFRa677jqWlpaYnp7mkUce4dChQ1VvliRVprIkLqX0LEB0r6d0M3B/SukU8K8RcQy4unjtWErpxeJz9xfvbXQSt7ICS0t5pYNSOehg9bQeUpMtLCywtLTE8vIyS0tLLCwsmMRJGmuj2CfuIuDRtuevFGUA31pV/kNrfUFE3A7cDjA/P8/CwkL/t7JNSnDq1AmOHj373ylr1koRrUXgJybqs4boiRMnBh73ujEmnTYaj9nZWSYnJ0kpMTk5yezsbOPi6THSyXh0MybdxjkmA03iIuJh4II1XvpUSumLg/rdlNI9wD0ABw8eTIcPHx7UTwFw6lSuJbjyyvX/Tq+1RFc3idZ1pYOFhQUGHfe6MSadNhqPw4cPc+DAgf//XBNr4TxGOhmPbsak2zjHZKApQkrpxzfxsePAxW3P31GUcYbykZZSq0m0rGWbmMjJmou/S+t36NChRiZvkrQZo5g2PAh8ISI+TR7YcDnwNSCAyyPiMnLydgtwa2Vb2UNKrRq25eVWU+jMDMzN5RGhJmySJOlsVZZKRMRPA78P7AH+JiKeSCndkFI6GhF/Rh6w8BZwR0ppufjMLwJfBrYA96aUjla0+V1SytN6lDVsO3a0Bh1MTVW9dZIkqWmqHJ36APBAj9fuAu5ao/xLwJcGvGkbNjWVm0T37s01bHUZfCBJkurLRr0+KEeOWuMmSZKGpbIVGyRJkrR5JnGSJEk1ZBInSZJUQyZxkiRJNWQSJ0mSVEMmcZIkSTVkEidJklRDJnGSJEk1ZBInSZJUQyZxkiRJNRQppaq3YaAi4j+BfxvCT+0G/msIv1MXxqObMelkPLoZk07Go5sx6da0mFySUtqznjc2Pokblog4klI6WPV2jArj0c2YdDIe3YxJJ+PRzZh0G+eY2JwqSZJUQyZxkiRJNWQS1z/3VL0BI8Z4dDMmnYxHN2PSyXh0MybdxjYm9omTJEmqIWviJEmSasgkroeIuDciXouIp9vKvj8iFiPi6xHx1xExW5RPRcTnivJnI+KTbZ+5MSKei4hjEXFnFfvSL32MyUtF+RMRcaSKfemHDcZjOiLuK8qfjIjDbZ+5qig/FhG/FxFRwe70RR9jslD83TxR3M6vYHfOWkRcHBH/GBHPRMTRiPh4Ub4rIh6KiOeL+51FeRTHwLGIeCoiDrR9123F+5+PiNuq2qez1eeYLLcdIw9WtU9nYxPx+N7i7+lURHxi1Xc14nrT55g04nrTU0rJ2xo34H3AAeDptrLHgGuLxx8FfqN4fCtwf/H4XOAl4FJgC/ACsBeYBp4Erqh636qMSfH8JWB31fsz5HjcAdxXPD4feByYKJ5/DbgGCOBvgQ9UvW8jEJMF4GDV+9OHeFwIHCgebwe+CVwB/DZwZ1F+J/BbxeObimMgimPiq0X5LuDF4n5n8Xhn1ftXZUyK105UvT8VxON84AeBu4BPtH1PY643/YpJ8dpLNOB60+tmTVwPKaWvAK+vKt4PfKV4/BDws+XbgW0RMQnMAEvAG8DVwLGU0osppSXgfuDmQW/7oPQpJo2xwXhcAfxD8bnXgP8BDkbEhcBsSunRlM84nwd+atDbPij9iMkQNnNoUkqvppT+pXj8JvAscBH5PPC54m2fo/VvfjPw+ZQ9CswVx8gNwEMppddTSv9NjuONQ9yVvuljTBpho/FIKb2WUnoMOL3qqxpzveljTBrPJG5jjtL6o/ggcHHx+C+Ak8CrwMvA76SUXicfdN9q+/wrRVmTbDQmkBO8v4+IxyPi9mFu7BD0iseTwE9GxGREXAZcVbx2Efm4KI3TMdIrJqX7iiaQX69zE3MpIi4F3gN8FZhPKb1avPTvwHzxuNc5o5HnkrOMCcDWiDgSEY9GRG3/81NaZzx6Gedj5EyafL0xidugjwK/EBGPk6t4l4ryq4Fl4O3AZcCvRMTeajZx6DYTk/emlA4AHwDuiIj3DXmbB6lXPO4ln1SPAJ8B/pkcn3GwmZj8XErp+4AfLW4fHuoW91lEvA34S+CXUkodNdJFDezYTRPQp5hckvJM/bcCn4mI7+n/lg6Hx0i3PsWkydcbk7iNSCl9I6V0fUrpKuBPyf0PIJ9A/i6ldLpoFvoncrPQcTprFt5RlDXGJmJCSul4cf8a8AA54WuEXvFIKb2VUvrllNIPpJRuBubI/TyOk4+L0tgcI2eISfsx8ibwBWp8jETEFPlC9Ccppb8qiv+jbBIs7l8rynudMxp1LulTTNqPkxfJ/SjfM/CNH4ANxqOXcT5Gemry9QZM4jYkihFyETEB/Brwh8VLLwM/Vry2jdz59hvkDt2XR8RlETEN3ALUcgRVLxuNSURsi4jtbeXXA0+v/t666hWPiDi32F8i4v3AWymlZ4qmgTci4pqiyfAjwBer2frB2GhMiubV3UX5FPAT1PQYKf5N/xh4NqX06baXHgTKEaa30fo3fxD4SDEi8xrgf4tj5MvA9RGxsxiRd31RVjv9ikkRi3OK79wN/AjwzFB2oo82EY9eGnO96VdMmn69ARyd2utGrjF4ldxR8hXgY8DHyTUF3wR+k9ZkyW8D/pzc9+cZ4Ffbvuem4v0vAJ+qer+qjgl55NSTxe1onWOywXhcCjxH7qD7MLkZqPyeg+QTywvAH5SfqeOtHzEBtpFHqj5VHCO/C2ypet82GY/3kpt8ngKeKG43AecBjwDPF/u+q3h/AJ8tjoWv0zZCl9wsfay4/XzV+1Z1TIAfLp4/Wdx/rOp9G1I8Lij+tt4gDwZ6hTw4ChpyvelXTGjQ9abXzRUbJEmSasjmVEmSpBoyiZMkSaohkzhJkqQaMomTJEmqIZM4SZKkGjKJkyRJqiGTOEmSpBoyiZMkSaqh/wOgfc75VhBvSwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "exog = fit.model.exog.copy()\n",
    "exog[:, 2:4] = 0 \n",
    "prediction = fit.get_prediction(exog=exog)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 5))\n",
    "ax.plot(annual_mean_df.year, annual_mean_df.height, 'k.')\n",
    "ax.set_ylabel(\"{} [{}]\".format(ssh.standard_name, 'mm'))\n",
    "ax.set_title('Annual mean')\n",
    "conf_int = prediction.conf_int()\n",
    "pred_int = prediction.conf_int(obs=True)\n",
    "ax.fill_between(annual_mean_df.year, pred_int[:, 0], pred_int[:, 1], alpha=0.1, color='blue', label='prediction interval')\n",
    "ax.fill_between(annual_mean_df.year, conf_int[:, 0], conf_int[:, 1], alpha=0.3, color='blue', label='confidence interval')\n",
    "ax.plot(annual_mean_df.year, fit.predict(exog=exog), color='blue', label='fitted trend')\n",
    "ax.legend()\n",
    "ax.grid(True)\n",
    "fig.savefig('linear.pdf')"
   ]
  }
 ],
 "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}