Pan_Arctic_discharge_v7_regional_error_calculations_v2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "8574dad9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "from scipy.spatial.distance import pdist\n",
    "import skgstat\n",
    "import xdem\n",
    "from typing import Callable\n",
    "\n",
    "# date_string='Fall2021_v11'\n",
    "date_string='Summer2022_v13'\n",
    "pd.set_option('display.max_columns', 250)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f5c3102d",
   "metadata": {},
   "outputs": [],
   "source": [
    "Arctic_discharge_fluxes_with_area = gpd.read_file('Arctic_discharge_fluxes_with_area_epsg3995_'+date_string+'_V4.gpkg', driver=\"GPKG\")#change file name based on distance\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "596e3efc",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_early_D = pd.read_pickle('df_all_reg_early_'+date_string+'.pkl')\n",
    "df_late_D = pd.read_pickle('df_all_reg_late_'+date_string+'.pkl')\n",
    "\n",
    "df_early_T = pd.read_pickle('df_all_reg_early_TC_'+date_string+'.pkl')\n",
    "df_late_T = pd.read_pickle('df_all_reg_late_TC_'+date_string+'.pkl')\n",
    "\n",
    "df_late_T['Tg_errH_corr'].fillna(0,inplace=True)#fix a problem where iceland is nan in this column (no area gained)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "df8d5541",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Alaska\n",
      "BBI\n",
      "FJL\n",
      "Greenland\n",
      "Iceland\n",
      "NZ\n",
      "QEI\n",
      "SZ\n",
      "Svalbard and Jan Mayen\n",
      "Alaska\n",
      "BBI\n",
      "FJL\n",
      "Greenland\n",
      "Iceland\n",
      "NZ\n",
      "QEI\n",
      "SZ\n",
      "Svalbard and Jan Mayen\n"
     ]
    }
   ],
   "source": [
    "regions = df_early_T['region'].unique()\n",
    "# regions = regions[regions!='Greenland'] #uncomment to exclude periphery Greenland from the total numbers\n",
    "for period in ['early','late']:\n",
    "    list_df_reg = []\n",
    "    list_df_uncertain = []\n",
    "    for region in regions:\n",
    "        print(region)\n",
    "        if period == 'early':\n",
    "            df_T = df_early_T[df_early_T['region']==region]\n",
    "            df_D = df_early_D[df_early_D['region']==region]\n",
    "            time='2000_to_2010'\n",
    "            time2='2000_2010'\n",
    "            time3 = 'early_to_middle_years_elapsed'\n",
    "        if period == 'late':\n",
    "            df_T = df_late_T[df_late_T['region']==region]\n",
    "            df_D = df_late_D[df_late_D['region']==region]\n",
    "            time='2010_to_2020'\n",
    "            time2='2010_2020'\n",
    "            time3 = 'middle_to_late_years_elapsed'\n",
    "            \n",
    "        df_FA = Arctic_discharge_fluxes_with_area[Arctic_discharge_fluxes_with_area['region']==region]\n",
    "        \n",
    "        Tg_errA = np.sum(df_FA['perimeter_error_gained_'+period]/df_FA[time3]*df_FA['thickness_for_area_gained_'+period]*0.9/1000)\n",
    "        Tl_errA = np.sum(df_FA['perimeter_error_lost_'+period]/df_FA[time3]*df_FA['thickness_for_area_lost_'+period]*0.9/1000)\n",
    "\n",
    "        sig_ID_raw_reg = df_D['V_tot_err_corr']\n",
    "        sig_TC_raw_reg = np.sqrt(df_T['Tg_errH_corr']**2 + df_T['Tl_errH_corr']**2 + Tg_errA**2 + Tl_errA**2)\n",
    "        sig_ID_clim_reg = np.sqrt(sum(df_FA['melt_correct_gt_'+period+'_err_discharge']**2))\n",
    "        sig_TC_clim_reg = np.sqrt(sum(df_FA['melt_correct_gt_'+period+'_err_terminus_change']**2))\n",
    "\n",
    "        sig_FA_reg = np.sqrt(sig_ID_raw_reg**2 + sig_ID_clim_reg**2 + sig_TC_raw_reg**2 + sig_TC_clim_reg**2) *2 #convert to 2 sigma error\n",
    "        sig_FA_reg=sig_FA_reg[0]\n",
    "#         print('--',sig_FA_reg)\n",
    "        \n",
    "        sig_D = np.sqrt(sig_ID_raw_reg**2 + sig_ID_clim_reg**2) *2\n",
    "        sig_TC = np.sqrt(+ sig_TC_raw_reg**2 + sig_TC_clim_reg**2) *2\n",
    "        \n",
    "        df_reg = pd.DataFrame()\n",
    "        df_reg = df_reg.assign(region=[region],FA=[sum(df_FA['Frontal_ablation_'+time+'_gt_per_yr_mean'])],FA_err=[sig_FA_reg],\n",
    "                               Discharge=[sum(df_FA[time2+'_mean_flux_gt'])],D_err=[sig_D[0]],\n",
    "                               Terminus_change=[sum(df_FA['terminus_gt_change_per_year_'+period])],TC_err=[sig_TC[0]]\n",
    "                              )\n",
    "        list_df_reg.append(df_reg)\n",
    "        \n",
    "        df_uncertain = pd.DataFrame()\n",
    "        df_uncertain = df_uncertain.assign(region=[region],\n",
    "                                           D_raw=[sig_ID_raw_reg[0]**2],\n",
    "                                           TC_raw=[sig_TC_raw_reg[0]**2],\n",
    "                                           D_clim=[sig_ID_clim_reg**2],\n",
    "                                           TC_clim=[sig_TC_clim_reg**2]\n",
    "                                          )\n",
    "        list_df_uncertain.append(df_uncertain)\n",
    "        \n",
    "    if period=='early':\n",
    "        df_reg_early = pd.concat(list_df_reg)\n",
    "        df_uncertain_early = pd.concat(list_df_uncertain)\n",
    "    if period=='late':\n",
    "        df_reg_late = pd.concat(list_df_reg)\n",
    "        df_uncertain_late = pd.concat(list_df_uncertain)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3e9b7eed",
   "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>region</th>\n",
       "      <th>D_raw</th>\n",
       "      <th>TC_raw</th>\n",
       "      <th>D_clim</th>\n",
       "      <th>TC_clim</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Alaska</td>\n",
       "      <td>0.019243</td>\n",
       "      <td>0.006858</td>\n",
       "      <td>0.012096</td>\n",
       "      <td>0.000409</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>BBI</td>\n",
       "      <td>0.000287</td>\n",
       "      <td>0.001110</td>\n",
       "      <td>0.000017</td>\n",
       "      <td>0.000002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>FJL</td>\n",
       "      <td>5.603590</td>\n",
       "      <td>0.772033</td>\n",
       "      <td>0.000418</td>\n",
       "      <td>0.000015</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Greenland</td>\n",
       "      <td>0.147894</td>\n",
       "      <td>0.462761</td>\n",
       "      <td>0.002228</td>\n",
       "      <td>0.000131</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Iceland</td>\n",
       "      <td>0.000846</td>\n",
       "      <td>0.000131</td>\n",
       "      <td>0.001534</td>\n",
       "      <td>0.000004</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>NZ</td>\n",
       "      <td>0.176438</td>\n",
       "      <td>0.081965</td>\n",
       "      <td>0.000313</td>\n",
       "      <td>0.000013</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>QEI</td>\n",
       "      <td>0.101518</td>\n",
       "      <td>0.206036</td>\n",
       "      <td>0.002656</td>\n",
       "      <td>0.000069</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>SZ</td>\n",
       "      <td>0.278130</td>\n",
       "      <td>0.046234</td>\n",
       "      <td>0.002624</td>\n",
       "      <td>0.000013</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Svalbard and Jan Mayen</td>\n",
       "      <td>0.973204</td>\n",
       "      <td>0.783041</td>\n",
       "      <td>0.005015</td>\n",
       "      <td>0.000064</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   region     D_raw    TC_raw    D_clim   TC_clim\n",
       "0                  Alaska  0.019243  0.006858  0.012096  0.000409\n",
       "0                     BBI  0.000287  0.001110  0.000017  0.000002\n",
       "0                     FJL  5.603590  0.772033  0.000418  0.000015\n",
       "0               Greenland  0.147894  0.462761  0.002228  0.000131\n",
       "0                 Iceland  0.000846  0.000131  0.001534  0.000004\n",
       "0                      NZ  0.176438  0.081965  0.000313  0.000013\n",
       "0                     QEI  0.101518  0.206036  0.002656  0.000069\n",
       "0                      SZ  0.278130  0.046234  0.002624  0.000013\n",
       "0  Svalbard and Jan Mayen  0.973204  0.783041  0.005015  0.000064"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_uncertain_early"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6d0dc89e",
   "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>region</th>\n",
       "      <th>FA_early</th>\n",
       "      <th>FA_early_err</th>\n",
       "      <th>FA_late</th>\n",
       "      <th>FA_late_err</th>\n",
       "      <th>D_early</th>\n",
       "      <th>D_early_err</th>\n",
       "      <th>D_late</th>\n",
       "      <th>D_late_err</th>\n",
       "      <th>TC_early</th>\n",
       "      <th>TC_early_err</th>\n",
       "      <th>TC_late</th>\n",
       "      <th>TC_late_err</th>\n",
       "      <th>gate_length_km</th>\n",
       "      <th>num_glaciers</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>total</td>\n",
       "      <td>44.468614</td>\n",
       "      <td>6.225412</td>\n",
       "      <td>51.980958</td>\n",
       "      <td>4.620898</td>\n",
       "      <td>35.398313</td>\n",
       "      <td>5.414075</td>\n",
       "      <td>40.415414</td>\n",
       "      <td>2.61204</td>\n",
       "      <td>-9.070301</td>\n",
       "      <td>3.073036</td>\n",
       "      <td>-11.565544</td>\n",
       "      <td>3.811816</td>\n",
       "      <td>3802.047382</td>\n",
       "      <td>959</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  region   FA_early  FA_early_err    FA_late  FA_late_err    D_early  \\\n",
       "0  total  44.468614      6.225412  51.980958     4.620898  35.398313   \n",
       "\n",
       "   D_early_err     D_late  D_late_err  TC_early  TC_early_err    TC_late  \\\n",
       "0     5.414075  40.415414     2.61204 -9.070301      3.073036 -11.565544   \n",
       "\n",
       "   TC_late_err  gate_length_km  num_glaciers  \n",
       "0     3.811816     3802.047382           959  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#global total\n",
    "df_global = pd.DataFrame()\n",
    "df_global['region'] = 'total'\n",
    "df_global = df_global.assign(FA_early=[sum(df_reg_early['FA'])],\n",
    "                             FA_early_err=[np.sqrt(sum(df_reg_early['FA_err']**2))],\n",
    "                             FA_late=[sum(df_reg_late['FA'])],\n",
    "                             FA_late_err=[np.sqrt(sum(df_reg_late['FA_err']**2))],\n",
    "                             \n",
    "                             D_early=[sum(df_reg_early['Discharge'])],\n",
    "                             D_early_err=[np.sqrt(sum(df_reg_early['D_err']**2))],\n",
    "                             D_late=[sum(df_reg_late['Discharge'])],\n",
    "                             D_late_err=[np.sqrt(sum(df_reg_late['D_err']**2))],\n",
    "                             \n",
    "                             TC_early=[sum(df_reg_early['Terminus_change'])],\n",
    "                             TC_early_err=[np.sqrt(sum(df_reg_early['TC_err']**2))],\n",
    "                             TC_late=[sum(df_reg_late['Terminus_change'])],\n",
    "                             TC_late_err=[np.sqrt(sum(df_reg_late['TC_err']**2))],\n",
    "                             \n",
    "                             gate_length_km = [sum(Arctic_discharge_fluxes_with_area['gate_length_km'])],\n",
    "                             num_glaciers = [Arctic_discharge_fluxes_with_area['RGIId'][Arctic_discharge_fluxes_with_area['region']!='Greenland'].count()]\n",
    "                             )\n",
    "df_global['region'] = 'total'\n",
    "df_global"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3a87abcb",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "48.22  ±  5.42\n"
     ]
    }
   ],
   "source": [
    "print(str(round((df_global['FA_early']+df_global['FA_late'])/2,2)[0]),' ± ',str(round((df_global['FA_early_err']+df_global['FA_late_err'])/2,2)[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "bf69abd0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#build table 1\n",
    "df_gate = Arctic_discharge_fluxes_with_area.groupby('region').sum()[['gate_length_km']].round().reset_index()\n",
    "df_num = Arctic_discharge_fluxes_with_area.groupby('region').count()[['RGIId']].reset_index()\n",
    "df_num.rename(columns={'RGIId':'num_glaciers'},inplace=True)\n",
    "\n",
    "df_final = df_reg_early.merge(df_gate,left_on='region',right_on='region').merge(df_num,left_on='region',right_on='region')\n",
    "df_final.rename(columns={'FA':'FA_early','FA_err':'FA_early_err',\n",
    "                        'Discharge':'D_early','D_err':'D_early_err',\n",
    "                        'Terminus_change':'TC_early','TC_err':'TC_early_err'},inplace=True)\n",
    "\n",
    "df_final = df_final.merge(df_reg_late)\n",
    "df_final.rename(columns={'FA':'FA_late','FA_err':'FA_late_err',\n",
    "                        'Discharge':'D_late','D_err':'D_late_err',\n",
    "                        'Terminus_change':'TC_late','TC_err':'TC_late_err'},inplace=True)\n",
    "\n",
    "df_final = df_final[['region', 'FA_early', 'FA_early_err','FA_late','FA_late_err',\n",
    "'D_early', 'D_early_err','D_late', 'D_late_err',\n",
    "'TC_early', 'TC_early_err','TC_late', 'TC_late_err','gate_length_km', 'num_glaciers',]]\n",
    "\n",
    "# df_final = df_final.append(df_global).round(2)\n",
    "df_final['gate_length_km']=df_final['gate_length_km'].astype('int')\n",
    "df_final.reset_index(drop=True,inplace=True)\n",
    "df_final.to_pickle('FA_regional_results'+date_string+'.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "37d862e2",
   "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>region</th>\n",
       "      <th>FA_early</th>\n",
       "      <th>FA_early_err</th>\n",
       "      <th>FA_late</th>\n",
       "      <th>FA_late_err</th>\n",
       "      <th>D_early</th>\n",
       "      <th>D_early_err</th>\n",
       "      <th>D_late</th>\n",
       "      <th>D_late_err</th>\n",
       "      <th>TC_early</th>\n",
       "      <th>TC_early_err</th>\n",
       "      <th>TC_late</th>\n",
       "      <th>TC_late_err</th>\n",
       "      <th>gate_length_km</th>\n",
       "      <th>num_glaciers</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Alaska</td>\n",
       "      <td>11.590140</td>\n",
       "      <td>0.392966</td>\n",
       "      <td>10.680662</td>\n",
       "      <td>0.332951</td>\n",
       "      <td>11.494634</td>\n",
       "      <td>0.354058</td>\n",
       "      <td>9.791141</td>\n",
       "      <td>0.180708</td>\n",
       "      <td>-0.095507</td>\n",
       "      <td>0.170487</td>\n",
       "      <td>-0.889521</td>\n",
       "      <td>0.279645</td>\n",
       "      <td>80</td>\n",
       "      <td>42</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>BBI</td>\n",
       "      <td>0.088878</td>\n",
       "      <td>0.075245</td>\n",
       "      <td>0.087500</td>\n",
       "      <td>0.075406</td>\n",
       "      <td>0.043441</td>\n",
       "      <td>0.034889</td>\n",
       "      <td>0.028674</td>\n",
       "      <td>0.024221</td>\n",
       "      <td>-0.045437</td>\n",
       "      <td>0.066668</td>\n",
       "      <td>-0.058826</td>\n",
       "      <td>0.071411</td>\n",
       "      <td>40</td>\n",
       "      <td>86</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>FJL</td>\n",
       "      <td>10.458657</td>\n",
       "      <td>5.050171</td>\n",
       "      <td>7.438550</td>\n",
       "      <td>3.293486</td>\n",
       "      <td>8.782888</td>\n",
       "      <td>4.734557</td>\n",
       "      <td>4.680333</td>\n",
       "      <td>2.212169</td>\n",
       "      <td>-1.675769</td>\n",
       "      <td>1.757326</td>\n",
       "      <td>-2.758217</td>\n",
       "      <td>2.439950</td>\n",
       "      <td>1377</td>\n",
       "      <td>328</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Greenland</td>\n",
       "      <td>4.305632</td>\n",
       "      <td>1.565904</td>\n",
       "      <td>3.176132</td>\n",
       "      <td>1.090514</td>\n",
       "      <td>2.251698</td>\n",
       "      <td>0.774911</td>\n",
       "      <td>1.881270</td>\n",
       "      <td>0.430582</td>\n",
       "      <td>-2.053934</td>\n",
       "      <td>1.360723</td>\n",
       "      <td>-1.294862</td>\n",
       "      <td>1.001908</td>\n",
       "      <td>697</td>\n",
       "      <td>537</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Iceland</td>\n",
       "      <td>0.101617</td>\n",
       "      <td>0.100303</td>\n",
       "      <td>0.026796</td>\n",
       "      <td>0.088645</td>\n",
       "      <td>0.098955</td>\n",
       "      <td>0.097580</td>\n",
       "      <td>0.000601</td>\n",
       "      <td>0.073061</td>\n",
       "      <td>-0.002662</td>\n",
       "      <td>0.023212</td>\n",
       "      <td>-0.026196</td>\n",
       "      <td>0.050200</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>NZ</td>\n",
       "      <td>2.665301</td>\n",
       "      <td>1.017308</td>\n",
       "      <td>4.148618</td>\n",
       "      <td>0.933137</td>\n",
       "      <td>2.016913</td>\n",
       "      <td>0.840834</td>\n",
       "      <td>3.144788</td>\n",
       "      <td>0.371454</td>\n",
       "      <td>-0.648388</td>\n",
       "      <td>0.572637</td>\n",
       "      <td>-1.003829</td>\n",
       "      <td>0.856018</td>\n",
       "      <td>182</td>\n",
       "      <td>39</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>QEI</td>\n",
       "      <td>4.136396</td>\n",
       "      <td>1.114052</td>\n",
       "      <td>4.275447</td>\n",
       "      <td>1.181872</td>\n",
       "      <td>2.681971</td>\n",
       "      <td>0.645519</td>\n",
       "      <td>2.244291</td>\n",
       "      <td>0.328192</td>\n",
       "      <td>-1.454426</td>\n",
       "      <td>0.907974</td>\n",
       "      <td>-2.031156</td>\n",
       "      <td>1.135391</td>\n",
       "      <td>481</td>\n",
       "      <td>252</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>SZ</td>\n",
       "      <td>3.504501</td>\n",
       "      <td>1.143682</td>\n",
       "      <td>5.330368</td>\n",
       "      <td>0.875906</td>\n",
       "      <td>3.145425</td>\n",
       "      <td>1.059725</td>\n",
       "      <td>4.229354</td>\n",
       "      <td>0.594152</td>\n",
       "      <td>-0.359076</td>\n",
       "      <td>0.430105</td>\n",
       "      <td>-1.101014</td>\n",
       "      <td>0.643579</td>\n",
       "      <td>283</td>\n",
       "      <td>45</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>Svalbard and Jan Mayen</td>\n",
       "      <td>7.617491</td>\n",
       "      <td>2.654298</td>\n",
       "      <td>16.816885</td>\n",
       "      <td>2.481380</td>\n",
       "      <td>4.882389</td>\n",
       "      <td>1.978100</td>\n",
       "      <td>14.414962</td>\n",
       "      <td>1.051844</td>\n",
       "      <td>-2.735102</td>\n",
       "      <td>1.769864</td>\n",
       "      <td>-2.401923</td>\n",
       "      <td>2.247415</td>\n",
       "      <td>657</td>\n",
       "      <td>166</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   region   FA_early  FA_early_err    FA_late  FA_late_err  \\\n",
       "0                  Alaska  11.590140      0.392966  10.680662     0.332951   \n",
       "1                     BBI   0.088878      0.075245   0.087500     0.075406   \n",
       "2                     FJL  10.458657      5.050171   7.438550     3.293486   \n",
       "3               Greenland   4.305632      1.565904   3.176132     1.090514   \n",
       "4                 Iceland   0.101617      0.100303   0.026796     0.088645   \n",
       "5                      NZ   2.665301      1.017308   4.148618     0.933137   \n",
       "6                     QEI   4.136396      1.114052   4.275447     1.181872   \n",
       "7                      SZ   3.504501      1.143682   5.330368     0.875906   \n",
       "8  Svalbard and Jan Mayen   7.617491      2.654298  16.816885     2.481380   \n",
       "\n",
       "     D_early  D_early_err     D_late  D_late_err  TC_early  TC_early_err  \\\n",
       "0  11.494634     0.354058   9.791141    0.180708 -0.095507      0.170487   \n",
       "1   0.043441     0.034889   0.028674    0.024221 -0.045437      0.066668   \n",
       "2   8.782888     4.734557   4.680333    2.212169 -1.675769      1.757326   \n",
       "3   2.251698     0.774911   1.881270    0.430582 -2.053934      1.360723   \n",
       "4   0.098955     0.097580   0.000601    0.073061 -0.002662      0.023212   \n",
       "5   2.016913     0.840834   3.144788    0.371454 -0.648388      0.572637   \n",
       "6   2.681971     0.645519   2.244291    0.328192 -1.454426      0.907974   \n",
       "7   3.145425     1.059725   4.229354    0.594152 -0.359076      0.430105   \n",
       "8   4.882389     1.978100  14.414962    1.051844 -2.735102      1.769864   \n",
       "\n",
       "    TC_late  TC_late_err  gate_length_km  num_glaciers  \n",
       "0 -0.889521     0.279645              80            42  \n",
       "1 -0.058826     0.071411              40            86  \n",
       "2 -2.758217     2.439950            1377           328  \n",
       "3 -1.294862     1.001908             697           537  \n",
       "4 -0.026196     0.050200               5             1  \n",
       "5 -1.003829     0.856018             182            39  \n",
       "6 -2.031156     1.135391             481           252  \n",
       "7 -1.101014     0.643579             283            45  \n",
       "8 -2.401923     2.247415             657           166  "
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_final"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "ad21841a",
   "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>region</th>\n",
       "      <th>FA_early</th>\n",
       "      <th>FA_late</th>\n",
       "      <th>D_early</th>\n",
       "      <th>D_late</th>\n",
       "      <th>TC_early</th>\n",
       "      <th>TC_late</th>\n",
       "      <th>gate_length_km</th>\n",
       "      <th>num_glaciers</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Alaska</td>\n",
       "      <td>11.59 ± 0.39</td>\n",
       "      <td>10.68 ± 0.33</td>\n",
       "      <td>11.49 ± 0.35</td>\n",
       "      <td>9.79 ± 0.18</td>\n",
       "      <td>-0.1 ± 0.17</td>\n",
       "      <td>-0.89 ± 0.28</td>\n",
       "      <td>80</td>\n",
       "      <td>42</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>BBI</td>\n",
       "      <td>0.09 ± 0.08</td>\n",
       "      <td>0.09 ± 0.08</td>\n",
       "      <td>0.04 ± 0.03</td>\n",
       "      <td>0.03 ± 0.02</td>\n",
       "      <td>-0.05 ± 0.07</td>\n",
       "      <td>-0.06 ± 0.07</td>\n",
       "      <td>40</td>\n",
       "      <td>86</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>FJL</td>\n",
       "      <td>10.46 ± 5.05</td>\n",
       "      <td>7.44 ± 3.29</td>\n",
       "      <td>8.78 ± 4.73</td>\n",
       "      <td>4.68 ± 2.21</td>\n",
       "      <td>-1.68 ± 1.76</td>\n",
       "      <td>-2.76 ± 2.44</td>\n",
       "      <td>1377</td>\n",
       "      <td>328</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Greenland</td>\n",
       "      <td>4.31 ± 1.57</td>\n",
       "      <td>3.18 ± 1.09</td>\n",
       "      <td>2.25 ± 0.77</td>\n",
       "      <td>1.88 ± 0.43</td>\n",
       "      <td>-2.05 ± 1.36</td>\n",
       "      <td>-1.29 ± 1.0</td>\n",
       "      <td>697</td>\n",
       "      <td>537</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Iceland</td>\n",
       "      <td>0.1 ± 0.1</td>\n",
       "      <td>0.03 ± 0.09</td>\n",
       "      <td>0.1 ± 0.1</td>\n",
       "      <td>0.0 ± 0.07</td>\n",
       "      <td>-0.0 ± 0.02</td>\n",
       "      <td>-0.03 ± 0.05</td>\n",
       "      <td>5</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>NZ</td>\n",
       "      <td>2.67 ± 1.02</td>\n",
       "      <td>4.15 ± 0.93</td>\n",
       "      <td>2.02 ± 0.84</td>\n",
       "      <td>3.14 ± 0.37</td>\n",
       "      <td>-0.65 ± 0.57</td>\n",
       "      <td>-1.0 ± 0.86</td>\n",
       "      <td>182</td>\n",
       "      <td>39</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>QEI</td>\n",
       "      <td>4.14 ± 1.11</td>\n",
       "      <td>4.28 ± 1.18</td>\n",
       "      <td>2.68 ± 0.65</td>\n",
       "      <td>2.24 ± 0.33</td>\n",
       "      <td>-1.45 ± 0.91</td>\n",
       "      <td>-2.03 ± 1.14</td>\n",
       "      <td>481</td>\n",
       "      <td>252</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>SZ</td>\n",
       "      <td>3.5 ± 1.14</td>\n",
       "      <td>5.33 ± 0.88</td>\n",
       "      <td>3.15 ± 1.06</td>\n",
       "      <td>4.23 ± 0.59</td>\n",
       "      <td>-0.36 ± 0.43</td>\n",
       "      <td>-1.1 ± 0.64</td>\n",
       "      <td>283</td>\n",
       "      <td>45</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>Svalbard and Jan Mayen</td>\n",
       "      <td>7.62 ± 2.65</td>\n",
       "      <td>16.82 ± 2.48</td>\n",
       "      <td>4.88 ± 1.98</td>\n",
       "      <td>14.41 ± 1.05</td>\n",
       "      <td>-2.74 ± 1.77</td>\n",
       "      <td>-2.4 ± 2.25</td>\n",
       "      <td>657</td>\n",
       "      <td>166</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   region      FA_early       FA_late       D_early  \\\n",
       "0                  Alaska  11.59 ± 0.39  10.68 ± 0.33  11.49 ± 0.35   \n",
       "1                     BBI   0.09 ± 0.08   0.09 ± 0.08   0.04 ± 0.03   \n",
       "2                     FJL  10.46 ± 5.05   7.44 ± 3.29   8.78 ± 4.73   \n",
       "3               Greenland   4.31 ± 1.57   3.18 ± 1.09   2.25 ± 0.77   \n",
       "4                 Iceland     0.1 ± 0.1   0.03 ± 0.09     0.1 ± 0.1   \n",
       "5                      NZ   2.67 ± 1.02   4.15 ± 0.93   2.02 ± 0.84   \n",
       "6                     QEI   4.14 ± 1.11   4.28 ± 1.18   2.68 ± 0.65   \n",
       "7                      SZ    3.5 ± 1.14   5.33 ± 0.88   3.15 ± 1.06   \n",
       "8  Svalbard and Jan Mayen   7.62 ± 2.65  16.82 ± 2.48   4.88 ± 1.98   \n",
       "\n",
       "         D_late      TC_early       TC_late  gate_length_km  num_glaciers  \n",
       "0   9.79 ± 0.18   -0.1 ± 0.17  -0.89 ± 0.28              80            42  \n",
       "1   0.03 ± 0.02  -0.05 ± 0.07  -0.06 ± 0.07              40            86  \n",
       "2   4.68 ± 2.21  -1.68 ± 1.76  -2.76 ± 2.44            1377           328  \n",
       "3   1.88 ± 0.43  -2.05 ± 1.36   -1.29 ± 1.0             697           537  \n",
       "4    0.0 ± 0.07   -0.0 ± 0.02  -0.03 ± 0.05               5             1  \n",
       "5   3.14 ± 0.37  -0.65 ± 0.57   -1.0 ± 0.86             182            39  \n",
       "6   2.24 ± 0.33  -1.45 ± 0.91  -2.03 ± 1.14             481           252  \n",
       "7   4.23 ± 0.59  -0.36 ± 0.43   -1.1 ± 0.64             283            45  \n",
       "8  14.41 ± 1.05  -2.74 ± 1.77   -2.4 ± 2.25             657           166  "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_show = pd.DataFrame()\n",
    "for i in range(0,len(df_final)):\n",
    "    df_show.loc[i,'region']=df_final['region'][i]\n",
    "    df_show.loc[i,'FA_early']=str(round(df_final['FA_early'][i],2))+' ± '+str(round(df_final['FA_early_err'][i],2))\n",
    "    df_show.loc[i,'FA_late']=str(round(df_final['FA_late'][i],2))+' ± '+str(round(df_final['FA_late_err'][i],2))\n",
    "    \n",
    "    df_show.loc[i,'D_early']=str(round(df_final['D_early'][i],2))+' ± '+str(round(df_final['D_early_err'][i],2))\n",
    "    df_show.loc[i,'D_late']=str(round(df_final['D_late'][i],2))+' ± '+str(round(df_final['D_late_err'][i],2))\n",
    "    \n",
    "    df_show.loc[i,'TC_early']=str(round(df_final['TC_early'][i],2))+' ± '+str(round(df_final['TC_early_err'][i],2))\n",
    "    df_show.loc[i,'TC_late']=str(round(df_final['TC_late'][i],2))+' ± '+str(round(df_final['TC_late_err'][i],2))\n",
    "    \n",
    "    df_show.loc[:,'gate_length_km']=df_final['gate_length_km']\n",
    "    df_show.loc[:,'num_glaciers']=df_final['num_glaciers']\n",
    "df_show"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "9c9e4625",
   "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>region</th>\n",
       "      <th>FA_total</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Alaska</td>\n",
       "      <td>11.14 ± 0.36</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>BBI</td>\n",
       "      <td>0.09 ± 0.08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>FJL</td>\n",
       "      <td>8.95 ± 4.17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Greenland</td>\n",
       "      <td>3.74 ± 1.33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Iceland</td>\n",
       "      <td>0.06 ± 0.09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>NZ</td>\n",
       "      <td>3.41 ± 0.98</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>QEI</td>\n",
       "      <td>4.21 ± 1.15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>SZ</td>\n",
       "      <td>4.42 ± 1.01</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>Svalbard and Jan Mayen</td>\n",
       "      <td>12.22 ± 2.57</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   region      FA_total\n",
       "0                  Alaska  11.14 ± 0.36\n",
       "1                     BBI   0.09 ± 0.08\n",
       "2                     FJL   8.95 ± 4.17\n",
       "3               Greenland   3.74 ± 1.33\n",
       "4                 Iceland   0.06 ± 0.09\n",
       "5                      NZ   3.41 ± 0.98\n",
       "6                     QEI   4.21 ± 1.15\n",
       "7                      SZ   4.42 ± 1.01\n",
       "8  Svalbard and Jan Mayen  12.22 ± 2.57"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_show = pd.DataFrame()\n",
    "for i in range(0,len(df_final)):\n",
    "    df_show.loc[i,'region']=df_final['region'][i]\n",
    "    df_show.loc[i,'FA_total']=str(round((df_final['FA_early'][i]+df_final['FA_late'][i])/2,2))+' ± '+str(round((df_final['FA_early_err'][i]+df_final['FA_late_err'][i])/2,2))\n",
    "\n",
    "df_show"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "3a7288c4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.1059368556443\n",
      "0.22294763124865832\n"
     ]
    }
   ],
   "source": [
    "#sea level rise of discharge:\n",
    "print((df_global['D_early'][0]+df_global['D_late'][0])/2*20/360)\n",
    "print((df_global['D_early_err'][0]+df_global['D_late_err'][0])/2*20/360)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d32ffb04",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10c59bd8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "7080d97d",
   "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>region</th>\n",
       "      <th>D_raw</th>\n",
       "      <th>TC_raw</th>\n",
       "      <th>D_clim</th>\n",
       "      <th>TC_clim</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Alaska</td>\n",
       "      <td>0.024324</td>\n",
       "      <td>0.006858</td>\n",
       "      <td>0.014131</td>\n",
       "      <td>0.000409</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>BBI</td>\n",
       "      <td>0.000253</td>\n",
       "      <td>0.001110</td>\n",
       "      <td>0.000018</td>\n",
       "      <td>0.000002</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>FJL</td>\n",
       "      <td>11.128534</td>\n",
       "      <td>0.772036</td>\n",
       "      <td>0.000389</td>\n",
       "      <td>0.000015</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Greenland</td>\n",
       "      <td>0.176968</td>\n",
       "      <td>0.462840</td>\n",
       "      <td>0.002108</td>\n",
       "      <td>0.000131</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Iceland</td>\n",
       "      <td>0.001920</td>\n",
       "      <td>0.000131</td>\n",
       "      <td>0.001702</td>\n",
       "      <td>0.000004</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>NZ</td>\n",
       "      <td>0.259470</td>\n",
       "      <td>0.081965</td>\n",
       "      <td>0.000302</td>\n",
       "      <td>0.000013</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>QEI</td>\n",
       "      <td>0.138093</td>\n",
       "      <td>0.206036</td>\n",
       "      <td>0.002484</td>\n",
       "      <td>0.000069</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>SZ</td>\n",
       "      <td>0.497362</td>\n",
       "      <td>0.046234</td>\n",
       "      <td>0.002447</td>\n",
       "      <td>0.000013</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Svalbard and Jan Mayen</td>\n",
       "      <td>2.040199</td>\n",
       "      <td>0.783037</td>\n",
       "      <td>0.004710</td>\n",
       "      <td>0.000064</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                   region      D_raw    TC_raw    D_clim   TC_clim\n",
       "0                  Alaska   0.024324  0.006858  0.014131  0.000409\n",
       "0                     BBI   0.000253  0.001110  0.000018  0.000002\n",
       "0                     FJL  11.128534  0.772036  0.000389  0.000015\n",
       "0               Greenland   0.176968  0.462840  0.002108  0.000131\n",
       "0                 Iceland   0.001920  0.000131  0.001702  0.000004\n",
       "0                      NZ   0.259470  0.081965  0.000302  0.000013\n",
       "0                     QEI   0.138093  0.206036  0.002484  0.000069\n",
       "0                      SZ   0.497362  0.046234  0.002447  0.000013\n",
       "0  Svalbard and Jan Mayen   2.040199  0.783037  0.004710  0.000064"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_uncertain_early"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "f6c5cfb0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<BarContainer object of 9 artists>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAD4CAYAAADCb7BPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAT8klEQVR4nO3de5RlZX3m8e8j7YWLqJGKy4Ch1CG4jPHaE+9KVBId4ngJk4URgcSEiSbeVhiHjCYhE82QOCYaMUwQFHUQFES8MKMQtIOJSKyG5u4lYmMITCgvATGogL/5Y78lp4uqru73nK4+lf5+1qpV++yzz96/Onvv87zve/Y5lapCkqTtdY+dXYAkaW0yQCRJXQwQSVIXA0SS1MUAkSR1WbeaG9tnn31qdnZ2NTcpSWvexo0bv1FVMzu7jsVWNUBmZ2eZm5tbzU1K0pqX5LqdXcNSHMKSJHUxQCRJXQwQSVIXA0SS1MUAkSR1MUAkSV0MEElSFwNEktTFAJEkdVnVT6JrsmaPPXdVt7f5+ENWdXuSpps9EElSFwNEktTFAJEkdTFAJEldDBBJUhcDRJLUxQCRJHUxQCRJXQwQSVIXA0SS1MUAkSR1MUAkSV0MEElSFwNEktTFAJEkdTFAJEldDBBJUpcVAyTJu5PclOTKkXk/luT8JF9pvx+wY8uUJE2bbemBnAo8d9G8Y4ELquoA4IJ2W5K0C1kxQKrqQuBbi2a/AHhvm34v8MLJliVJmna974E8qKpuBGi/f3y5BZMcnWQuydz8/Hzn5iRJ02aHv4leVSdV1fqqWj8zM7OjNydJWiW9AfLPSR4M0H7fNLmSJElrQW+AfAw4sk0fCXx0MuVIktaKbbmM93TgIuDAJNcneTlwPHBwkq8AB7fbkqRdyLqVFqiqlyxz17MnXIskaQ3xk+iSpC4GiCSpiwEiSepigEiSuhggkqQuBogkqYsBIknqYoBIkroYIJKkLgaIJKmLASJJ6mKASJK6GCCSpC4GiCSpiwEiSepigEiSuhggkqQuBogkqYsBIknqYoBIkroYIJKkLgaIJKmLASJJ6mKASJK6GCCSpC4GiCSpiwEiSepigEiSuowVIElel+SqJFcmOT3JfSZVmCRpunUHSJJ9gVcD66vqUcBuwGGTKkySNN3GHcJaB+yeZB2wB3DD+CVJktaC7gCpqn8C/ifwdeBG4OaqOm/xckmOTjKXZG5+fr6/UknSVBlnCOsBwAuAhwI/AeyZ5PDFy1XVSVW1vqrWz8zM9FcqSZoq4wxhPQf4WlXNV9XtwNnAUyZTliRp2o0TIF8HnpRkjyQBng1cM5myJEnTbpz3QC4GzgIuAa5o6zppQnVJkqbcunEeXFV/APzBhGqRJK0hfhJdktTFAJEkdTFAJEldDBBJUhcDRJLUxQCRJHUxQCRJXQwQSVIXA0SS1MUAkSR1MUAkSV0MEElSFwNEktTFAJEkdTFAJEldDBBJUhcDRJLUxQCRJHUxQCRJXQwQSVIXA0SS1MUAkSR1MUAkSV0MEElSFwNEktTFAJEkdTFAJEldDBBJUpexAiTJ/ZOcleSLSa5J8uRJFSZJmm7rxnz824FPVtWhSe4F7DGBmiRJa0B3gCTZG3gGcBRAVf0A+MFkypIkTbtxhrAeBswD70lyaZKTk+y5eKEkRyeZSzI3Pz8/xuYkSdNknABZBzweOLGqHgd8Fzh28UJVdVJVra+q9TMzM2NsTpI0TcYJkOuB66vq4nb7LIZAkSTtAroDpKr+H/CPSQ5ss54NXD2RqiRJU2/cq7BeBZzWrsC6FvjV8UuSJK0FYwVIVW0C1k+mFEnSWuIn0SVJXQwQSVIXA0SS1MUAkSR1MUAkSV0MEElSFwNEktTFAJEkdTFAJEldDBBJUhcDRJLUxQCRJHUxQCRJXQwQSVIXA0SS1MUAkSR1MUAkSV0MEElSFwNEktTFAJEkdTFAJEldDBBJUhcDRJLUxQCRJHUxQCRJXQwQSVIXA0SS1MUAkSR1GTtAkuyW5NIkn5hEQZKktWESPZDXANdMYD2SpDVkrABJsh9wCHDyZMqRJK0V4/ZA3ga8HvjhcgskOTrJXJK5+fn5MTcnSZoW3QGS5BeBm6pq49aWq6qTqmp9Va2fmZnp3ZwkacqM0wN5KvAfk2wGzgCeleR/T6QqSdLU6w6QqvrdqtqvqmaBw4BPV9XhE6tMkjTV/ByIJKnLukmspKo2ABsmsS5J0tpgD0SS1MUAkSR1MUAkSV0MEElSFwNEktTFAJEkdTFAJEldDBBJUhcDRJLUxQCRJHUxQCRJXQwQSVIXA0SS1MUAkSR1MUAkSV0MEElSFwNEktTFAJEkdTFAJEldDBBJUhcDRJLUxQCRJHUxQCRJXQwQSVIXA0SS1MUAkSR1MUAkSV0MEElSl+4ASfKQJJ9Jck2Sq5K8ZpKFSZKm27oxHnsH8DtVdUmS+wIbk5xfVVdPqDZJ0hTr7oFU1Y1VdUmb/g5wDbDvpAqTJE23cXogP5JkFngccPEk1idJO9rsseeu6vY2H3/Iqm5vNYz9JnqSvYAPA6+tqluWuP/oJHNJ5ubn58fdnCRpSowVIEnuyRAep1XV2UstU1UnVdX6qlo/MzMzzuYkSVNknKuwApwCXFNVfza5kiRJa8E4PZCnAi8DnpVkU/v5DxOqS5I05brfRK+qvwUywVokSWuIn0SXJHUxQCRJXQwQSVIXA0SS1MUAkSR1MUAkSV0MEElSFwNEktTFAJEkdTFAJEldDBBJUhcDRJLUxQCRJHUxQCRJXQwQSVIXA0SS1MUAkSR1MUAkSV0MEElSFwNEktTFAJEkdTFAJEldDBBJUpd1O7sASbuO2WPPXdXtbT7+kFXd3q7GHogkqYsBIknqYoBIkrr4Hoi0gzjer3/r7IFIkrqM1QNJ8lzg7cBuwMlVdfxEqtKaY2t7uq3m/nHf7Dq6eyBJdgPeCTwPeCTwkiSPnFRhkqTpNk4P5GeBf6iqawGSnAG8ALh6EoUtZgtXkqZLqqrvgcmhwHOr6tfb7ZcBT6yq31603NHA0e3mgcCX+svtsg/wjVXe5nKmpZZpqQOmp5ZpqQOsZSnTUgfsnFr2r6qZVd7misbpgWSJeXdLo6o6CThpjO2MJclcVa3fWdsfNS21TEsdMD21TEsdYC3TXAdMVy072zhXYV0PPGTk9n7ADeOVI0laK8YJkC8AByR5aJJ7AYcBH5tMWZKkadc9hFVVdyT5beBTDJfxvruqrppYZZOz04bPljAttUxLHTA9tUxLHWAtS5mWOmC6atmput9ElyTt2vwkuiSpiwEiSeqyZgIkyYuSVJJHtNuzSa7sXNetk63uR+u9M8mmJJcluSTJU9r82SS3jdz3uSQHtvsOSvKJHVzPws/s6PaSHJXkhK08/kFJPpDk2iQbk1yU5EU7otaRbW5Isl2XSC63P5NsTrLPVh53XJJjtrfGZdZ1avts1PY8ppK8deT2MUmOa9OfWrTvbkhy8SRqXaKO/ZJ8NMlX2r4+Icm927Fy86I6ntMec+vI49+Q5Kokl7dlnthZx4Yk69v67hxZ3zMX1bApyS1J/mQ71r2jzvnljr2JbK8do5Xk343Me12bt9MvJV4zAQK8BPhbhqu9ptVtVfXYqnoM8LvA/xi576sj970X+G+rWM/Cz+ZtfWCSAOcAF1bVw6rqCQzP/X6LlvMbnft9H3jxUiFXVb+wsN+ApwK3AG+cdAFtP58NnFNVBwAHALsDf9oW+eyiY+ivFz3+ycAvAo+vqkcDzwH+cYySfqat77aR9X1ltAbgZcDNwNuW+Zv+rR2TV7Dl696h7KBv/NheayJAkuzFcBK9nCUCpLWsP9ta/aMt/wcnubC1WK5M8vRFj9untaoPWW4dY9gb+HbHfdPiWcAPqup/Lcyoquuq6h2t53Jmko8D5yXZM8m7k3whyaVJXgDD96UleUubf3mS/9zmH9Ram2cl+WKS09oL2RaSnJhkrrVu/3Bk/uYkf9j20xW0D7UmeWCSL2fo7X0TeECb//Akn2y9qM+m9WIXbes3Wp2XJflwkj3a/FOT/EWGXuO1ab2MDE5IcnWSc4Ef73iO72C4oud1Kyz3duD/VNX5HdtYybOA71XVewCq6s5WzxHAXtvw+AcD36iq77fHf6OqbkjyvCQfWlio7fOPt+kl92vzo095Z+idnQe8P8lMm/cK4GKG5+6ERfvpz5J8BviTDB8vuKjt0z9arvgk57Tj4qoM35qxMP/WJG9ux8Pnkzyozd+m9Y6sZ68kFywcqyPnxmySa5K8q237vCS7L7Oacxi+JookD2MIz/mRbdzt+Uzy7CQfGVnm4CRnt+mfb3/DJe083qvN3+K8Wuo8uZuqmvof4HDglDb9OeDxwCxwZZu3B3CfNn0AMNemfwd4Q5veDbhvm74VeBDDgXjw1taxnXXeCWwCvth28hPa/FngtnbfV4EbgZ9s9x0EfGIHPW8L9WwCPrJ4e8BRwAnLPPbVwJ8vc99RDB8k/bF2+4+Bw9v0/YEvA3syfIXNG9v8ewNzwENbDTcz9GbuAVwEPK0ttwFY36YX1r9bm//odnsz8Ko2/Urg9jb9UeDrbV8ewvDNCPsAFwAHtGWeCHy6TR8HHNOmHzjy971pZP2nAme2Oh/J8P1vAC8Gzm+1/QTwL8Ch27l/bmVoTGwG7gccAxy3aJkXtf137x10jCy5n4FLgde2/bRp5OfhC7W333u1+V8G/hJ4Zpu/ru2LPdvtE0eOkeX26wbg6W191Z7fZwK/TztOGQL3xGX20yeA3drtjwFHtOnfWqh3ib9zoZbdgSsXjoO2/ee36T/lruN4W9d768jzsHeb3gf4B4YGzyxDCD623fehhedn0XqOa8fF2cCjgDcAR7LCedK28UVgpt33AeD5rYYLR/bLfwV+f5nz6uSVjp810QNhGL46o02f0W6PuifwrtYaPZPhRIfhw46/mmFc+Weq6jsjy18AvL7uatUtt47tsTBk9AjgucD7RlrWC0NYD2c4MVfjWvLRIayx3rtI8s7WGvtCm3V+VX2rTf88cGySTQwH8H2An2zzj2jzLwYeyBDOAH9fVddX1Q8ZXjBml9jsLye5hOHF7KfZcp+c3X5v5K6v1Xkyw4vLv1bVuQy9vD2BpwBntjr+iqHVvNijWu/kCuClbXsLzqmqH1bV1QwND4BnAKdX1Z1VdQPw6SXWuaKqugV4H8ML+RaS7Av8BfAr1Vr4O0BY4iuIuOs5XTyE9dXRharqVuAJDI2FeeCDSY6qqjuATwLPzzCkdAhDwMPW9+ttbX0/ZGhYfLDNe1qS5zEcU49eZj+dWUMPCoYRi9Pb9Pu38ve/OsllwOcZvllj4fj8AUMgwXCMzW7nehcE+OMklwN/DezLXcfQ16pq0xLbWMoZDKMvLwQ+sui+uz2fNaTA+4HDk9yf4dz4v8CTGJ7vv2vnw5HA/iPrGj2vtlYPsAb+I2GSBzJ0sx+VpBhSthhaOwteB/wz8BiGluL3AKrqwiTPYDh435/kLVX1Pobk3wj8AvA3W1tHr6q6KMPY9lJfgPYx4D3jrH8VXAX80sKNqvqt9vfMtVnfHVk2wC9V1RZflNnC81VV9alF8w9iGP9fcCeLjsUkD2Voef37qvp2klMZgmnB90ceOzr8tfjFMMC/1DB2vjWnAi+sqsuSHMXQS1q8rYX1LbetXm8DLmHkmGjP3XuB41tw7Shb7Oe27b0ZXuS+xPAexFa1F+0NwIb2wn4kw/P5QYZW+reAL1TVd7Zhv1JVd7Z2138HLgdewXDe/1X7/cpl9tPoMQkr7J92HD4HeHJV/WuSDSO13N5ehOHux+f27PeXMrwGPKGqbk+yeWQbi8+B5YawAD4OvIVhZOSWhXbpCs/ne9rjvscQrne04+r8qlrcCF8wel6tmA9roQdyKPC+qtq/qmar6iHA19jyzdz7ATe21uzLGA4ykuwP3FRV7wJOYRj6guEA+DXgEUmO3do6erXxw92Aby5x99MYhrKm2aeB+7Qx5wV7LLPsp4BXLfS2kjxuZP4rktyzzf+pJHtu4/b3ZnhBuLmNPz9vGx7zOeCVSfZordUHMAwTfS3Jf2o1JMljlnjsfYEbW60v3YZtXQgcluF9ngcDP7cNj1lS68l9iOE9vgXHMLw38c7e9W6jC4A9khwBP/o/P28FTmBo+W9VkgOTHDAy67HAdW16A8M59xsMYQIr79f92/ruwXDuP5ahF7kX8A6GodBt2U9/x13vly633P2Ab7fweARD63wl27Lexdu4qYXHz7Fla3+bVdVtDMNNb15017LPZ+sZ38Bw8cWpbfbngaemXdXVzpWf6qkJ1kAPhGG4avF/OvwwW17F9JfAh9uLxGe4qyVyEPBfktzO8EJyxMIDWivnMODjSW7Zyjq2x+6tWwhDS/XIkdbUw9t9Yege/3rH+idhHVu2fI5K8sKR209qQ0vV5v95ktczDE98l+EgXtxS+iOGVvTlLUQ2M1xJczJDN/iSNn+eoQu+otbCvJShhXwtw4m7kpczDHt8E/gOw5VLMJzoJyZ5I8NQ5RnAZYse+3sMw2zXMVz1ct8VtvURhp7xFQzj/3+z9cVX9FZg9F8hvAm4fuR4guHFrjuoltL284uAdyb5PYbW8ger6s2thf70RTW8qarOGrm9F/CONkxyB8MY/9Ft3XdmuGT8KIZeybbs190Zel7FcPHAvdpyswzvhd7GcJXX7QzP+z8t86e9BvhAktcwvF4s5ZPAb7bhpS8xvLiuZMX1tiG7hXPsNIbXmDnuen+0S1WdscS8lZ7P0xjeB7m6LT/fem6nJ7l3W+aNDM/ldvOrTHYx7cDft6pev7Nr0fTJcPXh6cCLq2rjzq5nLWo93HdV1c9OQS0nAJdW1Sk7ZP0GyK4jySkMV3L8clVdt9LykrZPkt9kuCDitVV13k6uZSPDqMHBO+oiDANEktRlLbyJLkmaQgaIJKmLASJJ6mKASJK6GCCSpC7/HwVB12L1hxoMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.bar(df_uncertain_early['region'],df_uncertain_early['D_raw'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "id": "94767d4d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:xlabel='region'>"
      ]
     },
     "execution_count": 175,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAF1CAYAAAD4PxH2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtjUlEQVR4nO3deXRUZbb+8e8mCJFBCZN6iXbAhCEkhBmcQEQFhwZFLqA4IcilG2fbJX1bbFu7f9Li1RYnGhtEkQsqgtNFxYmpG4UgsyKDoERQEBSJCJiwf39UUp2EDFUhUOTk+ayVBXXOqTq7IDx5ees9+5i7IyIilV+1WBcgIiIVQ4EuIhIQCnQRkYBQoIuIBIQCXUQkIBToIiIBUT1WJ27YsKEnJSXF6vQiIpXS0qVLv3P3RsXti1mgJyUlkZmZGavTi4hUSmb2ZUn7NOUiIhIQCnQRkYAoM9DNbJKZbTez1SXsH2xmK/O+/mVmGRVfpoiIlCWSOfTJwBPA8yXs3wR0d/fvzewiYALQpWLKEznUL7/8QlZWFvv27Yt1KSJHTHx8PImJiRx33HERP6fMQHf3+WaWVMr+fxV4+BGQGPHZRcohKyuLunXrkpSUhJnFuhyRCufu7Ny5k6ysLJo2bRrx8yp6Dn0o8FYFv6ZIIfv27aNBgwYKcwksM6NBgwZR/y+0wpYtmlkPQoF+dinHDAeGA5x22mkVdWqpghTmEnTl+R6vkBG6mbUB/gH0dfedJR3n7hPcvaO7d2zUqNh18SIiUk6HPUI3s9OAmcA17r7u8EuS8kga9X+xLuGI2jzmkliXcMwwM+644w7+53/+B4CHH36Y7Oxs7rvvviN2zqSkJOrWrQtAbm4u/fr1Y/To0dSsWfOInVOiF8myxWnAIqCFmWWZ2VAzG2FmI/IOuRdoADxlZsvNTJd/ihxBNWvWZObMmXz33XdH9bwffvghq1atYvHixXzxxRcMHz48oufl5uYe4cokX5mB7u5Xuvsp7n6cuye6+0R3H+/u4/P2D3P3BHdvm/fV8ciXLVJ1Va9eneHDh/Poo48esu/LL7+kZ8+etGnThp49e/LVV18BcP3113PLLbdw5pln0qxZM2bMmBF+ztixY+nUqRNt2rThj3/8Y5nnr1OnDuPHj+fVV19l165dxR4zd+5cevTowVVXXUV6ejoAl112GR06dKB169ZMmDABgJdeeok77rgDgMcee4xmzZoBsHHjRs4+u8SP46QEulJUpBIaOXIkU6dOZffu3YW233TTTVx77bWsXLmSwYMHc8stt4T3bdu2jYULF/Lmm28yatQoAObMmcP69etZvHgxy5cvZ+nSpcyfP7/M859wwgk0bdqU9evXl3jM4sWL+ctf/sKnn34KwKRJk1i6dCmZmZmMGzeOnTt30q1bNxYsWADAggULaNCgAV9//TULFy7knHPOifrPpapToItUQieccALXXnst48aNK7R90aJFXHXVVQBcc801LFy4MLzvsssuo1q1aqSmpvLtt98CoUCfM2cO7dq1o3379qxdu7bUkC6orBvMd+7cudAa6nHjxpGRkUHXrl3ZsmUL69ev5+STTyY7O5s9e/awZcsWrrrqKubPn8+CBQsU6OUQs26LInJ4brvtNtq3b8+QIUNKPKbg0reCH2Dmh7G78/vf/57/+q//iurce/bsYfPmzTRv3rzEY2rXrh3+/dy5c3nvvfdYtGgRtWrV4txzzw2vsT7jjDN49tlnadGiBeeccw6TJk1i0aJF4Q99JXIaoYtUUvXr12fAgAFMnDgxvO3MM89k+vTpAEydOrXMeehevXoxadIksrOzAfj666/Zvn07AD179uTrr78+5DnZ2dn89re/5bLLLiMhISGiWnfv3k1CQgK1atVi7dq1fPTRR+F93bp14+GHH6Zbt260a9eODz/8kJo1a3LiiSdG9Nrybxqhi1Rid955J0888UT48bhx47jhhhsYO3YsjRo14tlnny31+RdeeCGfffYZZ5xxBhD6wPOFF16gYcOGbNiwgfr164eP7dGjB+7OwYMHufzyyxk9enTEdfbu3Zvx48fTpk0bWrRoQdeuXcP7zjnnHLZs2UK3bt2Ii4vj1FNPpWXLlhG/tvyblTUPdqR07NjRdYOLilOV1qF/9tlntGrVKobVBN/q1auZNGkSjzzySKxLqdKK+143s6UlrSbUlIuIHCItLU1hXglpykVEym3VqlVcc801hbbVrFmTjz/+OEYVVW0KdBEpt/T0dJYvXx7rMiSPplxERAJCgS4iEhAKdBGRgFCgi4gEhAJdKr2kUf9XoV+R+Pnnn+nevTsrVqygbdu2tG3blvr169O0aVPatm3L+eefzy+//MKoUaNISUkhLS2Nzp0789Zbsb1D44MPPkhycjItWrTgnXfeKfaYl19+mdatW1OtWjUKXiuyePHi8HvNyMhg1qxZQKgNQP72tm3b0rBhQ2677TYA9u/fz8CBA0lOTqZLly5s3rw5/HrPPfccKSkppKSk8Nxzz4W3P/HEEyQnJ2NmhVoEjx07NnyOtLQ04uLiwt0ef/jhB/r370/Lli1p1aoVixYtKvSeHn744UNeb+XKlZxxxhm0bt2a9PT0Q2731qdPH9LS0gpte+mll0hNTaV169bhnjkffvhhofcfHx/Pq6++CoRaK/zhD3+gefPmtGrVKtx7580334yos2W0tMpFpBwmTZpEv379yMjICK/yuP7667n00kvp378/AKNGjWLbtm2sXr2amjVr8u233zJv3ryIXj83N5e4uLgKrfnTTz9l+vTprFmzhq1bt3L++eezbt26Q86TlpbGzJkzD+nvkpaWRmZmJtWrV2fbtm1kZGTw61//mrp16xZa6dKhQwf69esHwMSJE0lISGDDhg1Mnz6du+++mxdffJFdu3bxpz/9iczMTMyMDh060KdPHxISEjjrrLO49NJLOffccwud/6677uKuu+4C4I033uDRRx8NX8l666230rt3b2bMmMGBAwfYu3dv+Hlbtmzh3XffLXTby5ycHK6++mqmTJlCRkYGO3fu5LjjjgvvnzlzJnXq1Cl0/vXr1/Pggw/yz3/+k4SEhHCLhB49eoTf/65du0hOTubCCy8EYPLkyWzZsoW1a9dSrVq18HMuueQSRo8ezd13302tWrXK/suLkEboIuUwdepU+vbtW+L+vXv38swzz/D444+Hm2KddNJJDBgwoMTn1KlTh3vvvZcuXbqwaNEi7r//fjp16kRaWhrDhw/H3dm+fTsdOnQAYMWKFZhZuOf56aefXijIinrttdcYNGgQNWvWpGnTpiQnJ7N48eJDjmvVqhUtWrQ4ZHutWrWoXj00Bty3b1+x97xcv34927dvD3dKfO2117juuusA6N+/P++//z7uzjvvvMMFF1xA/fr1SUhI4IILLuDtt98GoF27diQlJZX4PgCmTZvGlVdeCcCPP/7I/PnzGTp0KAA1atSgXr164WNvv/12HnrooUL1zpkzhzZt2pCRkQFAgwYNwj/YsrOzeeSRR7jnnnsKnfOZZ55h5MiR4f41jRs3PqSuGTNmcNFFF4VD+umnn+bee++lWrVqhZ5jZpx77rm8+eabpb7PaCnQRaJ04MABvvjii1JDZ8OGDZx22mmccMIJEb/uTz/9RFpaGh9//DFnn302N910E0uWLGH16tX8/PPPvPnmmzRu3Jh9+/bx448/smDBAjp27MiCBQv48ssvady4camjva+//ppTTz01/DgxMbHY5lul+fjjj8NTFOPHjw8HfL5p06YxcODAcHgWPGf16tU58cQT2blz52HVsnfvXt5++22uuOIKAL744gsaNWrEkCFDaNeuHcOGDeOnn34C4PXXX6dJkybh4M63bt06zIxevXrRvn17HnroofC+0aNHc+eddx7yZ7lu3TrWrVvHWWedRdeuXcM/gAqaPn16+AcNhG7U8eKLL9KxY0cuuuiiQq2J8//uKpICXSRK3333XaERYEWJi4sLhxSE5ma7dOlCeno6H3zwAWvWrAFCHRX/+c9/Mn/+fP77v/874v7hxfVtivbO8l26dGHNmjUsWbKEBx988JB556KBVtI5D6eWN954g7POOis83ZKTk8Mnn3zCb37zG5YtW0bt2rUZM2YMe/fu5S9/+Qv333//Ia+Rk5PDwoULmTp1KgsXLmTWrFm8//77LF++nA0bNnD55ZcX+5z169czd+5cpk2bxrBhw/jhhx/C+7dt28aqVavo1atXeNv+/fuJj48nMzOTG2+8kRtuuCG8r3HjxmzdujWi9xwpBbpIlI4//vhDgqyo5ORkvvrqK/bs2RPx68bHx4f/279v3z5++9vfMmPGDFatWsWNN94YPuc555wTHpX37duXFStWsHDhQrp161bq6ycmJrJly5bw46ysLP7jP/4j4voKatWqFbVr12b16tXhbStWrCAnJyc8JVT0nDk5OezevZv69esfVi1Ff2gkJiaSmJhIly5dgNDUzieffMLGjRvZtGkTGRkZJCUlkZWVRfv27fnmm29ITEyke/fuNGzYkFq1anHxxRfzySefsGjRIpYuXUpSUhJnn30269atC8/lJyYm0rdvX4477jiaNm1KixYtCo24X3rpJS6//PJCc/GJiYnhH9KXX345K1euDO/bt28fxx9/fETvOVIKdJEoJSQkkJubW2qo16pVi6FDh3LLLbdw4MABIDSCe+GFFyI6R/5rN2zYkOzs7EL3AO3WrRsvvPACKSkpVKtWjfr16zN79mzOOussAGbNmsXvf//7Q16zT58+TJ8+nf3797Np0ybWr19P586dI37fmzZtIicnBwjdu/Tzzz8vNO1UcF674DnzV7DMmDGD8847LzzVMWfOHL7//nu+//575syZU2hkW5Ldu3czb968Qp9fnHzyyZx66ql8/vnnALz//vukpqaSnp7O9u3b2bx5M5s3byYxMZFPPvmEk08+mV69erFy5Ur27t1LTk4O8+bNIzU1ld/85jds3bqVzZs3s3DhQpo3b87cuXOB0B2fPvzwQyD0v7R169aF74Fa0vu/7LLL+OCDDwCYN29eoRuCrFu37pBVNIdLq1yk0ivYWvdoufDCC1m4cCHnn39+icf8+c9/5p577iE1NZX4+Hhq165d7H//i1OvXj1uvPFG0tPTSUpKolOnTuF9+SGaPyI/++yzycrKCn9Yt3HjxmLn7lu3bs2AAQNITU2levXqPPnkk+H/EQwbNowRI0bQsWNHZs2axc0338yOHTu45JJLaNu2Le+88w4LFy5kzJgxHHfccVSrVo2nnnqKhg0bhl//pZdeYvbs2YXOOXToUK655hqSk5OpX79++OYb9evXZ/To0eH3de+994anUMaNG8dDDz3EN998Q5s2bbj44ov5xz/+AYR+WF144YWF7oYE8PjjjzN48GAOHDhAs2bNyuwDn5CQwB133EGnTp0wMy6++GIuuaT076P8H0KpqanExcUxduxYGjRoAMDmzZvZsmUL3bt3L/ScUaNGMXjwYB599FHq1KkTfh8QmlJ78MEHSz1ntNQPPSDUD/3oWrZsGY888ghTpkyJaR3Fufrqq3n00Udp1KhRrEuREnz77bdcddVVvP/++6UeF20/dI3QRcqhXbt29OjR44isFz9ckU7rSOx89dVXR+SeqQp0kXIquGIhGl26dGH//v2Ftk2ZMoX09PSKKEsqgYJTaBVJgS5ylOnmD3KkaJWLiEhAlBnoZjbJzLab2eoS9puZjTOzDWa20szaV3yZIiJSlkhG6JOB3qXsvwhIyfsaDjx9+GWJiEi0ygx0d58P7CrlkL7A8x7yEVDPzE6pqAJFRCQyFTGH3gTYUuBxVt62Q5jZcDPLNLPMHTt2VMCpRYD7TqzYrwjk90PPzc1l8+bNHH/88bRr145WrVrRuXPnQv29o5GUlBTu2X3mmWeW6zWKoz7owe6Dnq8iVrkU11Gn2KuV3H0CMAFCFxZVwLlFYiK/H3r+GvTTTz+dZcuWAaHuf/369ePgwYMMGTKk3Of417/+VSG1qg968Pug56uIEXoWcGqBx4lAxbYQEznGlNYPvVmzZjzyyCPhUVlxsrOzGTJkCOnp6bRp04ZXXnnlkGPyg2Xu3Ll0796dAQMG0Lx5c0aNGsXUqVPp3Lkz6enpbNy4sdRa1Qc9+H3Q81VEoL8OXJu32qUrsNvdt1XA64ockyLph96+fXvWrl1b4v4HHniAE088kVWrVrFy5UrOO++8Us+5YsUKHnvsMVatWsWUKVNYt24dixcvZtiwYTz++OOlPld90P8tqH3Q80WybHEasAhoYWZZZjbUzEaY2Yi8Q2YDXwAbgGeA3x6RSkWOEZH0Qy+rR9J7773HyJEjw4/zR38l6dSpE6eccgo1a9bk9NNPD//XPj09vdD8dKS1qA96sPqg54tklcuV7n6Kux/n7onuPtHdx7v7+Lz97u4j3f10d093d3XckkCLpB/6smXLSm0g5u5RhWr+bewAqlWrFn5crVq1cEvbkqgPevD7oOfTlaIiUSqrH/rmzZv53e9+x80331zia1x44YU88cQT4cfff//9YdelPuhVtw96PvVykcrvvt1H/ZRF+6Fv3LiRdu3asW/fPurWrcvNN99c6gqXe+65h5EjR4aX4f3xj38MrxApL/VBr7p90POpH3pAqB/60XUs9kNXH/RjX6R90POpH7rIUXAs9kNXH/Rj35Hqg55PgS5STpH0Q3/22Wd57LHHCm0766yzePLJJ49UWXIMO1J90PMp0EWOoCFDhhzW1aIi0dAqFxGRgFCgi4gEhAJdRCQgNIculV76cxV7c+VV162q0NcTOVo0Qhcph/x+6CtWrAj3w65fvz5Nmzalbdu2nH/++fzyyy+MGjWKlJQU0tLS6Ny5M2+99VZU55k8eTI33XQTAOPHj+f555+vkPo3bdpEly5dSElJYeDAgRw4cKDY43r37k29evW49NJLC20fOnQoGRkZtGnThv79+5OdnQ2U3rv87bffpkWLFiQnJzNmzJjwa+3atYsLLriAlJQULrjggvBVszt37qRHjx7UqVMn/GcApfdhh+J7l+f78ccfadKkSaHXK6l3eb4lS5YQFxfHjBkzwttK6sM+cODAcF1JSUm0bds2/JyS+rCff/75FXKlMCjQRcolvx96RkYGy5cvZ/ny5fTp04exY8eyfPly3nvvPUaPHs22bdtYvXo1q1ev5o033mDPnj3lPueIESO49tprK6T+u+++m9tvv53169eTkJDAxIkTiz3urrvuKvbiqUcffZQVK1awcuVKTjvttHAbg7vuuiv85/Hggw/SvXt36tevT25uLiNHjuStt97i008/Zdq0aXz66acAjBkzhp49e7J+/Xp69uwZDvv4+HgeeOABHn744ULnzu/Dnv/1q1/9KnyVbcHe5WvWrOFvf/tboeeOHj36kKs6C/Yu/+yzzxg0aFB4X25uLnffffch7Qny+7CvXbuWFStWhC/+efHFF8N1XXHFFeG68vuwjx8/njVr1jB37txw75drrrmGp556qoS/qego0EXKobR+6BBq9frMM8/w+OOPhxtpnXTSSQwYMKDE57z99tu0b9+ejIwMevbsecj+++67Lxxu5557LrfffjvdunWjVatWLFmyhH79+pGSknJIP++i3J0PPviA/v37A3DdddeF77RTVM+ePalbt+4h2/NbDLg7P//8c7GNxgr2OFm8eDHJyck0a9aMGjVqMGjQIF577TWgcO/0grXUrl2bs88+m/j4+BLfS9E+7KX1Ll+6dCnffvttuFNlvpJ6l0OotcAVV1xRaFtZfdjz/1xeeuml8PsvrQ97nz59mDZtWonvMRoKdJEoRdIPfcOGDZx22mnF9lYpzo4dO7jxxht55ZVXWLFiBS+//HKZz6lRowbz589nxIgR9O3blyeffJLVq1czefJkdu7cWeLzdu7cSb169cI9zcvTHx1Ca+xPPvlk1q5de0gjsqK9y0vrg/7tt99yyimh2xCfcsop4bv7RKJoH/aSepcfPHiQO++8k7Fjxx7yGiX1Lv/666+ZNWsWI0aMKHR8aX3Y8y1YsICTTjqJlJSUcF0l9WFPSEhg//79pf6dRUqBLhKlSPqhR+ujjz6iW7duNG3aFCDcsKo0ffr0AUI90Vu3bh3ul96sWbNCLWqLqoj+6BC6Cnbr1q20atWKF198sdC+or3LK+qcRRVtqVtS7/KnnnqKiy++uNAPlXwl9S6/7bbb+Otf/3pIa4eS+rAXVLQDY0l92PNVVI90rXIRiVIk/dCTk5P56quv2LNnT7FTFkVF2x8dKNQTvWi/9NJ6pDds2JAffviBnJwcqlevflj90ePi4hg4cCBjx44tdEVscb3LS+qDftJJJ7Ft2zZOOeUUtm3bVuwt3opTUh/2rl27HtK7fNGiRSxYsICnnnqK7OxsDhw4QJ06dRgzZswhvcvz30dmZmZ4Pv27775j9uzZVK9ena5dux7Sh71goOfk5DBz5kyWLl1aqK78PuxAuA97/tRaRfVIV6BLpXe0lxkW7Ide0vxurVq1GDp0KLfccgt///vfqVGjBtu2beP999/n6quvPuT4M844g5EjR7Jp0yaaNm3Krl27Ihqll6Vnz548//zzNGnSJLzNzOjRowczZsxg0KBBPPfcc6V+HlCUu7Nx40aSk5Nxd9544w1atmwZ3p/fu7xgs7BOnTqxfv16Nm3aRJMmTZg+fTr/+7//C/y7d/qoUaOiqqWkPuTTpk3j+uuvL9S7fOrUqeFjJk+eTGZmZjiE83uX33DDDYV6l2/atCn8nOuvv55LL72Uyy67DCDch71FixbhPuz53nvvPVq2bEliYmJ4W69evXjooYfYu3cvNWrUYN68edx+++3hP89vvvmmzPupRkKBLlIORfuhF+fPf/4z99xzD6mpqcTHx1O7du1ib4sG0KhRIyZMmEC/fv04ePAgjRs35t133z2sGg8ePMiGDRuK/cHw17/+lUGDBnHPPffQrl278Ad8mZmZjB8/Pty/+5xzzmHt2rVkZ2eTmJjIxIkTueCCC7juuuv48ccfcXcyMjJ4+umnw69dXO/y6tWr88QTT9CrVy9yc3O54YYbaN26NRDqHT5gwAAmTpzIaaedVujzg6SkJH788UcOHDjAq6++Gu5LDsX3YS+td3lJSutdXpLS+rAX/d8JlN6HfenSpXTt2vWQ+7SWh/qhB4T6oR9dx2I/9KJWr17NpEmTeOSRR2JdipTi1ltvpU+fPsWubIq2H7o+FBUph4L90I9VaWlpCvNKIC0trdgwLw9NuYiUUyT90IvTpUsX9u/fX2jblClTSE+v2BYGUjnceOONFfZaCnSplMqzKuRY8fHHH8e6BKkEyjMdrikXqXTi4+PZuXNnub7hRSoDd2fnzp2lXiVbHI3QpdJJTEwkKyuLHTt2xLoUkSMmPj6+0NLHSCjQpdLJv2hERArTlIuISEBEFOhm1tvMPjezDWY2qpj9J5rZG2a2wszWmJnuiisicpSVGehmFgc8CVwEpAJXmllqkcNGAp+6ewZwLvA/ZlajgmsVEZFSRDJC7wxscPcv3P0AMB0o2mzBgboWWkdWB9gFlNwdSEREKlwkgd4EKNiLMytvW0FPAK2ArcAq4FZ3P1j0hcxsuJllmlmmViiIiFSsSAK9uKs3ii4A7gUsB/4DaAs8YWaHdPZ39wnu3tHdOzZq1CjKUkVEpDSRBHoWULArfCKhkXhBQ4CZHrIB2AS0REREjppIAn0JkGJmTfM+6BwEvF7kmK+AngBmdhLQAviiIgsVEZHSlXlhkbvnmNlNwDtAHDDJ3deY2Yi8/eOBB4DJZraK0BTN3e7+3RGsW0REiojoSlF3nw3MLrJtfIHfbwUuLPo8ERE5enSlqIhIQCjQRUQCQoEuIhIQCnQRkYBQoIuIBIQCXUQkIBToIiIBoUAXEQkIBbqISEAo0EVEAkKBLiISEAp0EZGAUKCLiASEAl1EJCAU6CIiAaFAFxEJCAW6iEhAKNBFRAJCgS4iEhAKdBGRgFCgi4gEhAJdRCQgFOgiIgGhQBcRCQgFuohIQCjQRUQCIqJAN7PeZva5mW0ws1ElHHOumS03szVmNq9iyxQRkbJUL+sAM4sDngQuALKAJWb2urt/WuCYesBTQG93/8rMGh+hekVEpASRjNA7Axvc/Qt3PwBMB/oWOeYqYKa7fwXg7tsrtkwRESlLJIHeBNhS4HFW3raCmgMJZjbXzJaa2bXFvZCZDTezTDPL3LFjR/kqFhGRYkUS6FbMNi/yuDrQAbgE6AWMNrPmhzzJfYK7d3T3jo0aNYq6WBERKVmZc+iERuSnFnicCGwt5pjv3P0n4Cczmw9kAOsqpEoRESlTJCP0JUCKmTU1sxrAIOD1Ise8BpxjZtXNrBbQBfisYksVEZHSlDlCd/ccM7sJeAeIAya5+xozG5G3f7y7f2ZmbwMrgYPAP9x99ZEsXERECotkygV3nw3MLrJtfJHHY4GxFVeaiIhEQ1eKiogEhAJdRCQgFOgiIgGhQBcRCQgFuohIQCjQRUQCQoEuIhIQCnQRkYBQoIuIBIQCXUQkIBToIiIBoUAXEQkIBbqISEAo0EVEAkKBLiISEAp0EZGAUKCLiASEAl1EJCAU6CIiAaFAFxEJCAW6iEhAKNBFRAJCgS4iEhAKdBGRgFCgi4gEhAJdRCQgIgp0M+ttZp+b2QYzG1XKcZ3MLNfM+ldciSIiEokyA93M4oAngYuAVOBKM0st4bi/Au9UdJEiIlK2SEbonYEN7v6Fux8ApgN9iznuZuAVYHsF1iciIhGKJNCbAFsKPM7K2xZmZk2Ay4Hxpb2QmQ03s0wzy9yxY0e0tYqISCkiCXQrZpsXefw34G53zy3thdx9grt3dPeOjRo1irBEERGJRPUIjskCTi3wOBHYWuSYjsB0MwNoCFxsZjnu/mpFFCkiImWLJNCXAClm1hT4GhgEXFXwAHdvmv97M5sMvKkwFxE5usoMdHfPMbObCK1eiQMmufsaMxuRt7/UeXMRETk6Ihmh4+6zgdlFthUb5O5+/eGXJSIi0dKVoiIiAaFAFxEJCAW6iEhAKNBFRAJCgS4iEhAKdBGRgFCgi4gEhAJdRCQgFOgiIgGhQBcRCQgFuohIQCjQRUQCQoEuIhIQCnQRkYBQoIuIBIQCXUQkIBToIiIBoUAXEQkIBbqISEAo0EVEAkKBLiISEAp0EZGAUKCLiASEAl1EJCAU6CIiAaFAFxEJiIgC3cx6m9nnZrbBzEYVs3+wma3M+/qXmWVUfKkiIlKaMgPdzOKAJ4GLgFTgSjNLLXLYJqC7u7cBHgAmVHShIiJSukhG6J2BDe7+hbsfAKYDfQse4O7/cvfv8x5+BCRWbJkiIlKWSAK9CbClwOOsvG0lGQq8VdwOMxtuZplmlrljx47IqxQRkTJFEuhWzDYv9kCzHoQC/e7i9rv7BHfv6O4dGzVqFHmVIiJSpuoRHJMFnFrgcSKwtehBZtYG+AdwkbvvrJjyREQkUpGM0JcAKWbW1MxqAIOA1wseYGanATOBa9x9XcWXKSIiZSlzhO7uOWZ2E/AOEAdMcvc1ZjYib/944F6gAfCUmQHkuHvHI1e2iIgUFcmUC+4+G5hdZNv4Ar8fBgyr2NJERCQaulJURCQgFOgiIgGhQBcRCQgFuohIQCjQRUQCQoEuIhIQCnQRkYBQoIuIBIQCXUQkIBToIiIBoUAXEQkIBbqISEAo0EVEAkKBLiISEAp0EZGAUKCLiASEAl1EJCAU6CIiAaFAFxEJCAW6iEhAKNBFRAJCgS4iEhAKdBGRgFCgi4gEhAJdRCQgFOgiIgERUaCbWW8z+9zMNpjZqGL2m5mNy9u/0szaV3ypIiJSmjID3czigCeBi4BU4EozSy1y2EVASt7XcODpCq5TRETKUD2CYzoDG9z9CwAzmw70BT4tcExf4Hl3d+AjM6tnZqe4+7YKr1iqpvtOjHUFR9Z9u2NdgQRAJIHeBNhS4HEW0CWCY5oAhQLdzIYTGsEDZJvZ51FVK1WWHf1TNgS+O2pn+1MM3qFUVr8qaUckgV7cd5qX4xjcfQIwIYJzisSUmWW6e8dY1yESjUg+FM0CTi3wOBHYWo5jRETkCIok0JcAKWbW1MxqAIOA14sc8zpwbd5ql67Abs2fi4gcXWVOubh7jpndBLwDxAGT3H2NmY3I2z8emA1cDGwA9gJDjlzJIkeFpgal0rHQwhQREansdKWoiEhAKNBFRAJCgS4iEhCRrEMXCTwzaw7cReiijfC/C3c/L2ZFiURJH4qKAGa2AhgPLAVy87e7+9KYFSUSJQW6CGBmS929Q6zrEDkcCnQRwMzuA7YDs4D9+dvdfVesahKJlgJdBDCzTcVsdndvdtSLESknBbqISEBo2aIIYGa1zOweM5uQ9zjFzC6NdV0i0VCgi4Q8CxwAzsx7nAX8OXbliERPgS4Scrq7PwT8AuDuPxOT+2qIlJ8CXSTkgJkdT96NWczsdAqsdhGpDHSlqEjIfcDbwKlmNhU4C7g+lgWJREurXETymFkDoCuhqZaP3P3o3VNUpAJohC4CmNkMYBLwlrsfjHU9IuWhOXSRkPHAYGC9mY0xs5axLkgkWppyESnAzE4ErgT+AGwBngFecPdfYlqYSAQ0QhfJkzeHfj0wDFgGPAa0B96NYVkiEdMIXQQws5lAS2AKMNndtxXYl+nuHWNWnEiEFOgigJmd5+4fxLoOkcOhQBfJY2ZpQCoQn7/N3Z+PXUUi0VGgiwBm9kfgXEKBPhu4CFjo7v1jWZdINPShqEhIf6An8I27DwEygJqxLUkkOgp0kZCf8y4oyjGzEwjdvUg3t5BKRVeKioRkmlk9QuvOlwLZwOKYViQSJc2hixRhZknACe6+Mta1iERDgS5Vmpm1L22/u39ytGoROVwKdKnSzOwgsAbYkb+pwG539/OOflUi5aM5dKnq7gSuAH4GpgOz3D07tiWJlI9G6CKAmTUl1JSrL/Al8P/cfXlMixKJkpYtigDuvgl4DZgDdAaax7YikehphC5Vmpk1AwYRGplvITTt8qa774tpYSLloECXKi3vQ9GVhEbnP5J3k+h87v5ILOoSKQ99KCpV3f38O8TrxLIQkcOlEbqISEDoQ1ERkYBQoIuIBIQCXUQkIPShqAhgZjUJXTGaRIF/F+5+f6xqEomWAl0k5DVgN6HWuftjXItIuWiViwhgZqvdPS3WdYgcDs2hi4T8y8zSY12EyOHQCF0EMLNPgWRgE6EpFyPUPrdNTAsTiYICXQQws18Vt93dvzzatYiUlz4UFeHfwW1mjYH4GJcjUi6aQxcBzKyPma0nNOUyD9gMvBXTokSipEAXCXkA6Aqsc/emQE/gn7EtSSQ6CnSRkF/cfSdQzcyqufuHQNsY1yQSFc2hi4T8YGZ1gPnAVDPbDuTEuCaRqGiViwhgZrWBfYSWKw4GTgSm5o3aRSoFBbqISEBoykWqNDPbQ5HbzhWwH9gI/MHd3z96VYmUj0boIiUwszggjdDUi/q8yDFPq1xESuDuue6+Ang81rWIREIjdBGRgNAIXUQkIBToIiIBoVUuUqWVscoFdz/hKJYjclgU6FKluXtdADO7H/gGmMK/Ly6qG8PSRKKmD0VFADP72N27lLVN5FimOXSRkFwzG2xmcWZWzcwGA7mxLkokGgp0kZCrgAHAt3lf/5m3TaTS0JSLiEhA6ENREcDMGgE3AkkU+Hfh7jfEqiaRaCnQRUJeAxYA76G5c6mkNOUiApjZcndvG+s6RA6HPhQVCXnTzC6OdREih0MjdBHCV4zWJtQD/RdCFxe5rhSVykSBLiISEPpQVCSPmSUAKUB8/jZ3nx+7ikSio0AXAcxsGHArkAgsB7oCi4DzYliWSFT0oahIyK1AJ+BLd+8BtAN2xLYkkego0EVC9rn7PgAzq+nua4EWMa5JJCqachEJyTKzesCrwLtm9j2wNaYViURJq1xEijCz7sCJwNvufiDW9YhESoEuIhIQmkMXEQkIBbqISEAo0EVEAkKrXKRKy+vhUuIHSerlIpWJAl2qNHevC2Bm9wPfAFMINeYaDNSNYWkiUdMqFxHAzD529y5lbRM5lmkOXSQk18wGm1mcmVUzs8HozkVSySjQRUKuAgYA3+Z9/WfeNpFKQ3PoUuWZWRww0t37xroWkcOhEbpUee6eC3SIdR0ih0sjdJGQZWb2OvAy8FP+RnefGbuSRKKjQBcJqQ/spPANLRxQoEuloWWLIiIBoRG6CGBm8cBQoDWF7yl6Q8yKEomSPhQVCZkCnAz0AuYRurfonphWJBIlTbmIAGa2zN3bmdlKd29jZscB77i7bhItlYZG6CIhv+T9+oOZpRG6Y1FS7MoRiZ7m0EVCJphZAjAaeB2ok/d7kUpDUy4iIgGhKRcRwMwamNnjZvaJmS01s7+ZWYNY1yUSDQW6SMh0YDtwBdAf+A54MaYViURJUy4igJktdfcORbZlunvHWNUkEi2N0EVCPjSzQXm90KuZ2QDg/2JdlEg0NEKXKq3APUUNqM2/b2oRB2TrnqJSmSjQRUQCQuvQpUozs/al7Xf3T45WLSKHSyN0qdLM7MNSdrsu/ZfKRIEuIhIQmnIRyZPXwyWVwu1zn49dRSLR0QhdBDCzPwLnEgr02cBFwEJ37x/LukSioXXoIiH9gZ7AN+4+BMgAasa2JJHoKNBFQn5294NAjpmdQKgNQLMY1yQSFc2hi4Rkmlk94BlgKZANLI5pRSJR0hy6SBFmlgSc4O4rY12LSDQU6CKAmb1GqLvia+7+U6zrESkPzaGLhDwCnA18amYvm1l/M4sv60kixxKN0EUKMLM44DzgRqC3mnNJZaIPRUXymNnxwK+BgUB74LnYViQSHY3QRQAzexHoArwNvATMzVvGKFJpKNBFADPrDbzr7rllHixyjFKgS5VmZp2ALe7+Td7jawndV/RL4D533xXL+kSioVUuUtX9HTgAYGbdgDHA88BuYEIM6xKJmj4UlaoursAofCAwwd1fAV4xs+WxK0skehqhS1UXZ2b5A5uewAcF9mnAI5WKvmGlqpsGzDOz74CfgQUAZpZMaNpFpNLQh6JS5ZlZV+AUYE7+Zf9m1hyoo3uKSmWiQBcRCQjNoYuIBIQCXUQkIBToIiUws/vN7PxY1yESKc2hS5VgZkbo+139WSSwNEKXwDKzJDP7zMyeAj4BRpvZEjNbaWZ/KnDcaDNba2bvmtk0M/td3vbJZtY/7/c9zWyZma0ys0lmVjNv+2Yz+5OZfZK3r2Us3qsIKNAl+FoQupT/bqAJ0BloC3Qws25m1pFQ75Z2QD+gY9EXyLvRxWRgoLunE7p+4zcFDvnO3dsDTwO/O2LvRKQMCnQJui/d/SPgwryvZYRG6y2BFEJ3KXrN3X929z3AG8W8Rgtgk7uvy3v8HNCtwP6Zeb8uBZIq/B2IREhXikrQ5d8f1IAH3f3vBXea2e0RvIaVsX9/3q+56N+UxJBG6FJVvAPcYGZ1AMysiZk1BhYCvzaz+Lx9lxTz3LVAUl47AIBrgHlHo2iRaGg0IVWCu88xs1bAotCCF7KBq919iZm9Dqwg1AM9kyI9XNx9n5kNAV7Oa+S1BBh/VN+ASAS0bFGqPDOr4+7ZZlYLmA8MVw8XqYw0QheBCWaWCsQDzynMpbLSCF1EJCD0oaiISEAo0EVEAkKBLiISEAp0EZGAUKCLiASEAl1EJCD+P7JymEF+x95mAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_uncertain.pivot('region','D_raw').plot(kind='bar')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "id": "fcb4180d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:title={'center':'2000-2010'}, xlabel='region', ylabel='squared unceratinties (Gt2 a-2)'>"
      ]
     },
     "execution_count": 196,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAGFCAYAAAD6j5gmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxU0lEQVR4nO3de5xVdb3/8deb4TIEjhcUNUkHDMUUQYQ0U1Ixb4kd07JQk8xIT5qmntQux0snS/P0qzRDungpNI8ZamXeL0Ba3CSvaApqaApiwqCggJ/fH2sNbsaZPXtm9tprZvb7+Xjsx+y1vnuv74cBPvu7v+t7UURgZmbVo0feAZiZWWU58ZuZVRknfjOzKuPEb2ZWZZz4zcyqjBO/mVmVceI3M6syTvzW5UnqI+mXkp6X1CDpYUmHFJSPk7RA0puS7pO0XUGZJF0saVn6uESSCsrr0/e8mV7jgCJxDJR0vaSXJC2X9BdJezR5zYQ0zjck3Sxps4Kyz0h6MK3r/mauP1LS3LR8rqSR7f+tWTVz4rfuoCfwT+BjwMbAt4H/S5P25sDv03ObAXOAGwreOwn4D2AEsCtwGPDlgvLrgYeBAcA3gd9J2qKFOPoDs4Hd07quAf4kqT+ApJ2BK4HjgC2BN4ErCt7/GvAj4PtNLyypN3AL8Btg0/Tat6TnzdpEnrlr3ZGkR4ALSBL2xIjYKz3fD3gV2C0iFkh6ELg6Iqak5V8EvhQRe0raAXgU2DwiGtLyGcDUiJhcYhwrgP0iYq6ki4D6iJiQlm0PPAkMaLx+ev5E4NiI2Lfg3IHAVcCgSP/TSnoBmBQRt7fz12RVyi1+63YkbQnsADwO7Az8vbEsIt4Ank3P07Q8fV5YtrAwKTcpby2OkUBv4Jnm6oqIZ4G301hbszPwSGzYUnuk1FjMCjnxW7ciqRcwFbgmIhaQdL8sb/Ky5cBG6fOm5cuB/mk/f2vvLRZHHfBr4IKIaLxGu6/XwfeabcCJ37oNST1Iku3bwCnp6ZVAXZOX1gENLZTXASvTlnXR90p6XNLK9LFPQRx9gT8Af42I7xW8t7VYiunIe8024MRv3ULaQv8lyU3TIyNiTVr0OMmN28bX9QO2T8+/pzx9Xlg2RNJGzZVHxM4R0T99zEiv3we4GXiRDW8SNxfLEKAP8HQJf8THgV0LRxyR3Ix+vIXXm7XIid+6i58BOwHjI2JVwflpwC6SjpRUC/w3SV/5grT8WuAMSdtIej9wJnA1QEQ8DcwHzpNUK+kIkmR7U3MBpN1MvwNWAZ+PiHeavGQqMF7SPukH0IXA7wtuHNekMfYEeqR19krfez+wDvhqOny18RvNvW37NZkBEeGHH136AWwHBLCapEuk8XFMWn4AsIAkId9PMrKm8b0CLiEZSvla+lwF5fXpe1YBTwEHFInjY2kcbzaJY5+C10wAXgDeIBmeuVlB2cT0/YWPqwvKdwPmprHMIxmZlPvv34+u9/BwTjOzKuOuHjOzKpNZ4pf0K0lLJD1WcG4zSXdJ+kf6c9Os6jczs+Zl2eK/Gji4yblzgHsiYihwT3psZmYVlGkfv6R64I8RsUt6/BSwb0T8S9LWwP0RsWNmAZiZ2XtUuo9/y4j4F0D6c2CF6zczq3o98w6gJZImkaycSL9+/XYfNmxYzhGZmXUtc+fOfTUi3rOabKUT/yuSti7o6lnS0gsjWS1xCsDo0aNjzpw5lYrRzKxbkPR8c+cr3dVzK3B8+vx4kgksZmZWQVkO57weeAjYUdLidJ3z7wMfl/QP4OM0s+GEmZllK7Ounoj4XAtF47Kq08zMWtdpb+6amQGsWbOGxYsXs3r16rxD6bRqa2sZNGgQvXr1av3FOPGbWSe3ePFiNtpoI+rr69lwVWqDZKHNZcuWsXjxYgYPHlzSe7xWj5l1aqtXr2bAgAFO+i2QxIABA9r0jciJ38w6PSf94tr6+3HiNzOrMu7jr0bnb1ym6zTd+9sse/Xn/Kms13vu+59o9TU1NTUMHz6cNWvW0LNnT44//nhOP/10evTomm1nJ34zs1b07duX+fPnA7BkyRImTJjA8uXLueCCC1p979q1a+nZs3Ol2q75cWVmlpOBAwcyZcoULr/8clpa3fjqq6/m05/+NOPHj+fAAw9k5cqVjBs3jlGjRjF8+HBuuSVZtOCSSy7hJz/5CQBf+9rX2H///QG45557OPbYYzP7Mzjxm5m10ZAhQ3jnnXdYsqTF5cZ46KGHuOaaa7j33nupra1l2rRpzJs3j/vuu48zzzyTiGDs2LHMmDEDgDlz5rBy5UrWrFnDzJkz2WeffTKLv3N9/zAz6yJa28vk4x//OJttttn6137jG99g+vTp9OjRgxdffJFXXnmF3Xffnblz59LQ0ECfPn0YNWoUc+bMYcaMGeu/CWTBid/MrI0WLlxITU0NAwe2vKVIv3791j+fOnUqS5cuZe7cufTq1Yv6+npWr169/vlVV13FXnvtxa677sp9993Hs88+y0477ZRZ/O7qMTNrg6VLl3LSSSdxyimnlDx+fvny5QwcOJBevXpx33338fzz766WPHbsWC699FLGjh3LPvvsw+TJkxk5cmSmcxfc4jezLqWU4ZfltmrVKkaOHLl+OOdxxx3HGWecUfL7jznmGMaPH8/o0aMZOXIkhRtL7bPPPnz3u9/lIx/5CP369aO2tjbT/n1w4jcza9W6deva9PqJEycyceLE9cebb745Dz30ULOvHTduHGvWrFl//PTTT7crxrZwV4+ZWZVxi9/MrJ3uuOMOzj777A3ODR48mGnTpuUUUWmc+M3M2umggw7ioIMOyjuMNnNXj5lZlXHiNzOrMk78ZmZVxonfzKzK+OZuFapffV1ZrvNcWa5i1kbl2k9i/fWK7yuxbNkyxo0bB8DLL79MTU0NW2yxBQATJkzgV7/6FT179qSmpoYzzzyTz3/+8+WNLwNO/GZmRQwYMGD9Wvznn38+/fv356yzzmLy5MlMmzaNWbNmUVdXx/Lly7n55ptLuua6deuoqanJLuhWuKvHzKwdLrroIq644grq6uoA2HjjjTn++ONbfH19fT0XXnghe++9NzfeeCM///nPGTNmDCNGjODII4/kzTffZN26dQwZMoSI4PXXX6dHjx5Mnz4dSJZ2eOaZZ8oSuxO/mVkbNTQ00NDQwPbbb9+m99XW1jJz5kw++9nP8qlPfYrZs2fz97//nZ122olf/vKX1NTUsMMOO/DEE08wc+ZMdt99d2bMmMFbb73F4sWL+eAHP1iW+J34zczaKCLatXrm0Ucfvf75Y489xj777MPw4cOZOnUqjz/+OJC07KdPn8706dM599xzmTlzJrNnz2bMmDFli9+J38ysjerq6ujXrx8LFy5s0/sK1+ifOHEil19+OY8++ijnnXceq1evBpLEP2PGDGbNmsWhhx7K66+/zv3338/YsWPLFn/RxC9pkKSzJN0iabak6ZKukPQJSf7QMLOqde655/KVr3yFFStWALBixQqmTJlS8vsbGhrYeuutWbNmDVOnTl1/fo899uDBBx+kR48e1NbWMnLkSK688sqyLtXc4qgeSVcB2wB/BC4GlgC1wA7AwcA3JZ0TEdPLFo2ZWWtaGX5ZKSeffDIrV65kzJgx9OrVi169enHmmWeW/P7vfOc77LHHHmy33XYMHz6choYGAPr06cMHPvAB9txzTyD5BnD99dczfPjwssWulvaNlLRLRDzW4hul3sC2EVGe28xFjB49OubMmZN1NVWj/pw/leU6eWyIYdXnySefzHQbwu6iud+TpLkRMbrpa1ts8RdL+mn520DmSd/MzMqrWFdPHfANku6eP0fEdQVlV0TEf1YgPjOzLuWII45g0aJFG5y7+OKLO9XyzcVm7l4F/AO4CThB0pHAhIh4C9izEsGZmXU1nX0TFig+qmf7iDgnIm6OiMOBecC9kgZUKDYzM8tAsRZ/H0k9IuIdgIj4rqTFwHSgf0WiMzOzsivW4v8DsH/hiYi4BjgTeLsjlUr6mqTHJT0m6XpJtR25npmZla7FxB8RX4+Iu5s5f3tEDG1vhZK2Ab4KjI6IXYAa4LPtvZ6ZmbVNm5ZllvTHiDisTPX2lbQGeB/wUhmuaWZVYPg15ZvIBPDo8Y+2+pqamhqGDx/OmjVr6NmzJ8cffzynn346PXq0bQGDfffdl0svvZTRo0dz6KGHct1117HJJpu0M/L2a+t6/Nt0tMKIeFHSpcALwCrgzoi4s+nrJE0CJgFsu+22Ha3WzKzd+vbtu35N/iVLljBhwgSWL1/OBRdc0O5r3nbbbWWKru3aut7Owx2tUNKmwCeBwcD7gX6Sjm36uoiYEhGjI2J04243ZmZ5GzhwIFOmTOHyyy+npZUP1q1bx1lnncXw4cPZddddueyyy97zmvr6el599VWee+45hg0bxoknnsguu+zCMcccw913381HP/pRhg4dyqxZs8r+Z2hT4o+IE8pQ5wHAoohYGhFrgN8De5XhumZmFTFkyBDeeecdlixZ0mz5lClTWLRoEQ8//DCPPPIIxxxzTNHrPfPMM5x22mk88sgjLFiwgOuuu46ZM2dy6aWXctFFF5U9/la7eiQNBb4HfIhkkTYAImJIO+t8AdhT0vtIunrGAV6Ix8y6lJZa+wB33303J510Ej17Jil2s802K3qtwYMHr1+Ebeedd2bcuHFIYvjw4Tz33HNli7lRKS3+q4CfAWuB/YBrgV+3t8KI+BvwO5IJYY+mMZS+lqmZWc4WLlxITU0NAwcObLa8rRu19OnTZ/3zHj16rD/u0aMHa9eu7ViwzSgl8feNiHtIVvJ8PiLOp8n4/raKiPMiYlhE7BIRx6XLQJiZdXpLly7lpJNO4pRTTmkxuR944IFMnjx5fdJ+7bXXKhliq0oZ1bM63XTlH5JOAV4Emv+YMzPLWCnDL8tt1apVjBw5cv1wzuOOO44zzjijxdefeOKJPP300+y666706tWLL33pS5xyyikVjLi4FtfjX/8CaQzwJLAJ8B2gDvhBRPw18+hSXo+/vLwev3UlXo+/NGVZj79RRMxOn64EvlCWCM3MLDdtncBlZmapO+64g7PPPnuDc4MHD+70SzM78ZuZtdNBBx3UqTZYKVVbZ+6amVkXVzTxSzpI0hcl1Tc5X44ZvGZmloMWE7+ki4BvAsOBeySdWlDcecYlmZlZmxTr4x8P7BYRayWdD1wnaUhEfA0ofUqamVkZPTmsvEM7d1rwZFmv1xUU6+rpGRFrASLidZIPgjpJNwK9KxCbmVnuli1bxsiRIxk5ciRbbbUV22yzzfrjSy65hGHDhrHLLrswYsQIrr322jZd+/777+eww5ItTm699Va+//3vZ/FHeI9iLf5nJX0sIh4AiIh1wBcl/Q9wZEWiMzPL2YABA9avxX/++efTv39/zjrrLCZPnsy0adOYNWsWdXV1LF++nJtvvrnd9Rx++OEcfvjh5Qm6FcVa/J8G3rMQdER8C/hAZhGZmXUBF110EVdccQV1dXUAbLzxxhx//PEtvn727NnstddejBgxgg9/+MM0NDRsUH711VevX9Zh4sSJnHzyyey3334MGTKEBx54gBNOOIGddtqJiRMndjj2Flv8EbEKQNI9ETGuSfG1JMspm5lVnYaGBhoaGth+++1Lev3bb7/N0UcfzQ033MCYMWNYsWIFffv2Lfqef//739x7773ceuutjB8/nr/85S/84he/YMyYMcyfP5+RI0e2O/5io3pqJW0GbC5pU0mbpY96kp2zzMyqUluXXX7qqafYeuutGTNmDAB1dXXr1+pvyfjx49evyb/lllsyfPhwevTowc4779zhNfqL1fxl4HSSJD+Xd0fyrAB+2qFazcy6sLq6Ovr168fChQsZMqT1Pana+kEBbLAmf9P1+ju6Rn+xrp4fAz+WdGpEvHfDSDOzHHSW4ZfnnnsuX/nKV7jhhhuoq6tjxYoV/Pa3v2XSpEnvee2wYcN46aWXmD17NmPGjKGhoaHVrp4stZj4JdUBWzYmfUmfBhojvSMiXqlAfGZmndLJJ5/MypUrGTNmDL169aJXr16ceeaZzb62d+/e3HDDDZx66qmsWrWKvn37cvfdd1c44ne1uB6/pCnAgxFxdXr8DPBnkuS/NiJOqlSQXo+/vLwev3UlXo+/NOVaj38MST9/o4aIODW92MxyBGpmZpVXLPH3jA2/DhxX8HyTbMIxM+vajjjiCBYtWrTBuYsvvrhTLd9cLPG/I2mriHgZICIeA5C0DfBOJYIzM4P2jYrJSx6bsLS2hW5TxWbu/gD4g6SxkjZKHx8Dbk7LzMwyV1tby7Jly9qc3KpFRLBs2TJqa2tLfk+x4Zy/kfQq8D/AzkAAjwP/HRF/7miwZmalGDRoEIsXL2bp0qV5h9Jp1dbWMmjQoJJfX3TqWETcDtze0aDMzNqrV69eDB48OO8wupViSzZ8K12yoaXy/SUdlk1YZmaWlWIt/kdJ+vhXA/OApUAtMBQYCdwNXJR1gGZmVl7F+vhvAW6RNBT4KLA1yTo9vwEmNa7eaWZmXUvx5eGAiPgH8I8KxGJmZhVQbDinmZl1Q078ZmZVxonfzKzKtJr4JV0iqU5SL0n3SHpV0rGVCM7MzMqvlBb/gRGxAjgMWAzsAPxXplGZmVlmSkn8vdKfhwLXR8RrGcZjZmYZKyXx/0HSAmA0cI+kLYDVHalU0iaSfidpgaQnJX2kI9czM7PStZr4I+Ic4CPA6IhYA7wJfLKD9f4YuD0ihgEjgM6xiaaZWRUo5ebu+4CvAD9LT72fpPXfLulevmOBXwJExNsR8Xp7r2dmZm1TSlfPVcDbwF7p8WKSpZrbawjJuj9XSXpY0i8k9evA9czMrA1KSfzbR8QlwBqAdI2ejmyF0xMYBfwsInYD3gDOafoiSZMkzZE0x+twm5mVTymJ/21JfUk2YkHS9sBbHahzMbA4Iv6WHv+O5INgAxExJSJGR8ToLbbYogPVmZlZoVIS/3kkm7F8QNJU4B7g6+2tMN3D95+SdkxPjQOeaO/1zMysbUpZnfMuSfOAPUm6eE6LiFc7WO+pwFRJvYGFwBc6eD0zMytRi4lf0rCIWCCpsRvmX+nPbSVtGxHz2ltpRMynAyODzMys/Yq1+M8AJgH/20xZAPtnEpGZmWWq2A5ck9Knh0TEBjN1JdVmGpWZmWWmlJu7D5Z4zszMuoBiffxbAdsAfSXtxrtj9+uA91UgNjMzy0CxPv6DgInAIOCHBecbgG9kGJOZmWWoWB//NcA1ko6MiJsqGJOZmWWo1XH8wB8lTQDqC18fERdmFZSZmWWnlMR/C7AcmEvHlmowM7NOoJTEPygiDs48EjMzq4iShnNKGp55JGZmVhGltPj3BiZKWkTS1SMgImLXTCMzM7NMlJL4D8k8CjMzq5hiE7jqImIFybh9MzPrJoq1+K8DDiMZzRNsuOtWkGyhaGZmXUyxCVyHpT8HVy4cMzPLWqujeiTdU8o5MzPrGor18deSLMa2uaRN2XCRtvdXIDYzM8tAsT7+LwOnkyT5ubyb+FcAP802LDMzy0qxPv4fAz+WdGpEXFbBmMzMLEOlbLZ+maRdgA8BtQXnr80yMDMzy0ariV/SecC+JIn/NpIJXTMBJ34zsy6olLV6jgLGAS9HxBeAEUCfTKMyM7PMlJL4V0XEO8BaSXXAEjx5y8ysyyplrZ45kjYBfk4yumclMCvLoMzMLDtFE78kAd+LiNeByZJuB+oi4pFKBGdmZuVXtKsnIgK4ueD4OSd9M7OurZQ+/r9KGpN5JGZmVhGl9PHvB5wk6TngDbwRi5lZl+aNWMzMqkyrXT0R8TzwAWD/9PmbpbzPzMw6p1KWZT4POBs4Nz3VC/hNlkGZmVl2Smm5HwEcTtK/T0S8BGyUZVBmZpadUhL/2+mwzgCQ1C/bkMzMLEulJP7/k3QlsImkLwF3k8ziNTOzLqiUZZkvlfRxkg1YdgT+OyLuyjwyMzPLRCnLMg8GZjQme0l9JdVHxHMdqVhSDTAHeLFxY3czM8teKV09NwLvFByvS8911GnAk2W4jpmZtUEpib9nRLzdeJA+792RSiUNAj4B/KIj1zEzs7YrJfEvlXR444GkTwKvdrDeHwFfZ8NvEhuQNEnSHElzli5d2sHqzMysUSmJ/yTgG5JekPRPkslcX25vhZIOA5ZExNxir4uIKRExOiJGb7HFFu2tzszMmihlVM+zwJ6S+gOKiIYO1vlR4HBJh5Js3l4n6TcRcWwHr2tmZiUoZVRPH+BIoB7omezNAhFxYXsqjIhzSZd/kLQvcJaTvplZ5ZSyOuctwHKSbRffyjYcMzPLWimJf1BEHJxF5RFxP3B/Ftc2M7PmlXJz90FJwzOPxMzMKqKUFv/ewERJi0i6erwDl5lZF+YduMzMqkwpiT8yj8LMzCqmlMT/J5LkL5Jx94OBp4CdM4zLzMwyUsoErg1u7EoaRQdm7pqZWb7avGl6RMwDxmQQi5mZVUApM3fPKDjsAYwCvGqamVkXVUoff+HG6mtJ+vxvyiYcMzPLWil9/BdUIhAzM6uMNvfxm5lZ1+bEb2ZWZZz4zcyqTIt9/JIuo8is3Yj4aiYRmZlZpoq1+OeQrMFfSzKE8x/pYySwLvPIzMwsEy22+CPiGgBJE4H9ImJNejwZuLMi0ZmZWdmV0sf/fjYcy98/PWdmZl1QKRO4vg88LOm+9PhjwPmZRWRmZpkqZQLXVZL+DOyRnjonIl7ONiwzM8tKq109kgQcAIyIiFuA3pI+nHlkZmaWiVL6+K8APgJ8Lj1uAH6aWURmZpapUvr494iIUZIeBoiIf0vqnXFcZmaWkVJa/Gsk1ZBO5pK0BfBOplGZmVlmSmnx/wSYBgyU9F3gKOBbmUZlZtYdnb9xma6zvENvL5r4JfUAFgFfB8aR7Lv7HxHxZIdqNTOz3BRN/BHxjqT/jYiPAAsqFJOZmWWolD7+OyUdmQ7rNDOzLq6UPv4zgH7AWkmrSbp7IiLqMo3MzMwyUcrM3Y1ae42ZmXUdpbT4kbQpMJRkiWYAImJ6VkGZmVl2Wk38kk4ETgMGAfOBPYGHgP0zjczMzDJRys3d04AxwPMRsR+wG7A006jMzCwzpST+1RGxGkBSn4hYAOyYbVhmZpaVUvr4F0vaBLgZuEvSv4GXsgzKzMyyU8qoniPSp+enm7FsDNze3golfQC4FtiKZM2fKRHx4/Zez8zM2qaUm7vbFhwuSn9uBbzQzjrXAmdGxDxJGwFzJd0VEU+083pmZtYGpXT1/IlkZU6RDOccDDwF7NyeCiPiX8C/0ucNkp4EtgGc+M3MKqCUrp7hhceSRgFfLkflkupJRgn9rRzXMzOz1pUyqmcDETGPZHhnh0jqD9wEnB4RK5opnyRpjqQ5S5d69KiZWbmU0sd/RsFhD2AUHRzHL6kXSdKfGhG/b+41ETEFmAIwevTo6Eh9ZmadQf3q68pynec6+P5S+vgL1+pZS9Lnf1N7K0xX+fwl8GRE/LC91zEzs/YppY//gjLX+VHgOOBRSfPTc9+IiNvKXI+ZmTWjlK6eW4uVR8ThbakwImaSjBAyM7MclNLVs4hk3P5v0uPPkXQx3ZFRTGZmlqFSEv9uETG24PgPkqZHxDeyCsrMzLJTynDOLSQNaTyQNBjYIruQzMwsS6W0+L8G3C9pYXpcT5kmcJmZWeWVMqrndklDgWHpqQUR8Va2YZmZWVZa7eqR9Gmgd0T8HRgPXJ8u22BmZl1QKX38304XU9sbOAi4BvhZtmGZmVlWSkn869KfnwB+FhG3AL2zC8nMzLJUSuJ/UdKVwGeA2yT1KfF9ZmbWCZWSwD9DMlnr4Ih4HdgM+K8sgzIzs+yUMqrnTeD3BcfrN1IxM7Oux102ZmZVxonfzKzKOPGbmVUZJ34zsyrjxG9mVmWc+M3MqowTv5lZlXHiNzOrMk78ZmZVxonfzKzKOPGbmVUZJ34zsyrjxG9mVmWc+M3MqowTv5lZlXHiNzOrMq1uxGJWEedvXKbrLC/Pdcy6Mbf4zcyqjBO/mVmVceI3M6sy7uM3s+7J941a5MRvZt1S/errynKd58pylc7Fid/MOqZcLWvolq3rzsiJ36wrcZK1Msjl5q6kgyU9JekZSefkEYOZWbWqeItfUg3wU+DjwGJgtqRbI+KJSsdi1tWUq98aumfftZUmj66eDwPPRMRCAEm/BT4JOPFXMd+I67r8YdT1KCIqW6F0FHBwRJyYHh8H7BERpzR53SRgUnq4I/BUmULYHHi1TNcqF8dUGsdUus4Yl2MqTTlj2i4itmh6Mo8Wv5o5955Pn4iYAkwpe+XSnIgYXe7rdoRjKo1jKl1njMsxlaYSMeVxc3cx8IGC40HASznEYWZWlfJI/LOBoZIGS+oNfBa4NYc4zMyqUsW7eiJiraRTgDuAGuBXEfF4BUMoe/dRGTim0jim0nXGuBxTaTKPqeI3d83MLF9endPMrMo48ZuZVRknfjOzKuNF2qwoSe+PCA+3NcuIpB2A/wK2oyAnR8T+mdXpm7uVJekPNDNhrVFEHF7BcFol6YWI2DaHejcrVh4Rr1UqlkaSPlWsPCJ+X6lYCkl6CPhmRNzbTNk9ETEuh7A6HUnfA/5fRCxppuziiDg7h7CQ9HdgMjAXWNd4PiLmZlZnd0/8kvYELgN2AnqTDCF9IyLqcornY8XKI+KBSsVSCkn/jIgPtP7Kste7iOQDUsC2wL/T55sAL0TE4Bxiuip9OhDYC2hMtPsB90dE0Q+GDON6EVgO3AacGxFrCsoejojdKhxPA803bgREjv/3lgOvA/8ZEX9qUjYvIkblFNfciNi9knVWQ1fP5SSTxG4ERgOfBz6YVzCdLbGXIJeWQWNilzQZuDUibkuPDwEOyCmmL6Qx/BH4UET8Kz3emmTF2by8AuwN/AT4m6TPRUTj2lYV//uLiI0qXWeJFgLHAVPTf0dnRcTqtKy5pWQq5Q+S/hOYBrzVeDLLb7XVkPiJiGck1UTEOuAqSQ/mFYukocA3SFqwPwR+DuwDPAt8MSLm5BDTZbTcQtukstG8x5iIOKnxICL+LOk7eQYE1Dcm/dQrwA55BQMQEW8CJ0o6ArhL0kURMZkcEpqk/Ru7nSQNjohFBWWfyqtLjOTbxmOS9gAuAeZImhARj+QUT6Pj05//VXAugCFZVVgNif/NdGmI+ZIuAf4F9MsxnquAa4E64G/A6cARJMn/p8AeOcRU7MOm4h9ETbwq6VvAb0j+MxwLLMs3JO6XdAdwPUlMnwXuyzekRERMkzQLuFrSoUD/HMK4FGjsNrmp4DnAt4C8Ej8AaSv/q2mr/w+SfpJzPJXvtqyCPv7tSFpkvYGvARsDP42IZ3OKZ35EjEyfPxMRH2yuzBLpTd7zgLHpqenABXnc3C2U3ujdJz2cHhHTcozl9og4uJnz/wV8JyJqKxzP+vsKTe8x5HHPoVjdkjYHfgV8IiJqcorrfcAZwLYRMSntFdgxIv6YVZ3V0OL/UEQ8D6wGLgCQdBJJ10oe3il4vqJIWcW0MtLoLZLf1U8j4p+ViyqRJvjTKl1va9Luilxbro2aS/rp+R8AP6hwOLDhv6Wm/67ybGmeJmmriHgZQNLngSNJ9n8ZW+yNGbuKZETPXunxYpJ7kpkl/mpo8T8IfKugz/HrwH4RcUhO8bwJPEPS97p9+pz0eEhEVLwbqpWRRj2BnYHPRcRHKhTSeukY57OAeio0xrmEmD4FXEwyukfkP1rlv4sUR0RU9J6IpNdJvpmJ5FvR9MYiYO+I2LSS8RTENQ84ICJekzQW+C1wKjAS2CkijsoprjkRMbrJN6W/R8SIrOqshhb/4cAf06+9BwPD0nN52SnHuluyKCJeKFJ+j6RdKxbNhm4kGeP8CwrGOOfsEmB8RDyZdyCpN5o59z7gRGAAUOmb4Z9Mf/YD7iZp5T9F8q370grHUqimoIvwaGBKRNwE3CRpfn5h8bakvqTfhiRtT8HonkxERLd/kLTMHiH5SqW842kmvs3zjAuYV/D8prx/H01im5t3DM3E9Je8YygS20YkN1AXkX4rySGG3sCPgNeAecDDwIvAOWn5bjn9bh4DeqbPFwBjC8ty/Ds7EHgAWApMJel62jfLOrtti79gEonSn71JhkcdJSkiv6/lewLfJ/lP8R3g1ySJv4ekz0fE7XmEVfA8syFk7VTxMc4lmCPpBuBmNowptz7/9Cb4GcAxwDXAqIj4d07hXEryjWO7iGhI46sDLpX0M5Jv3hUfyUIyCusBSa8Cq4AZaWwfJJkAl4uIuFPSXGBPkv+Lp0VEpvsAd/s+/s5G0hyScfwbk2y4cEhE/FXSMOD6yGHEQ+GsxTxnMDYnncHbVEREbh9QBTN4C0VEnFDxYABJPwA+RfLv6acRsTKPOArieQYYGk2Si6Qakk3ED4mIv+YU257A1sCdEfFGem4HoH9EzMsppt+RjCy6PSIqMsCj2yd+SR8F5kfEG5KOJRlT/KMo3qedZTzz493hnE9GxE4FZQ/nlPjXkfQTC+gLvNlYRI43La00kt4h+eaxlg1HzeTy9yfp6YhodkJbsbJqJekA4AskLf4bgasjYkGWdXbbrp4CPwNGSBoBfB34JUn3StE1czJU+Im+qklZXssj5DJ+uVSSdgE+BKwfjx4R1+YYTy3wRZLRToUx5dLij4jOtrz6E2m35QZ/R2nDq7PcEO80IuJu4G5JGwOfI5l5/U+SWf2/iYK1l8qlGlr88yJiVDrk7cWI+GXOCzIVa13XRkSvPOLqrCSdB+xLkvhvAw4BZkZOQ+/SmG4kuTk4AbiQpF/9yYjodPMN8iBpG5I5DqtIxqcHMIbk3/sREfFijuF1SpIGkMxKPw54ieQm797A8IjYt+z1VUHifwC4neSr1FiSO+fzI2J4roFZSSQ9CowAHo6IEZK2BH4REeNzjOnhiNhN0iMRsaukXsAdkePcgs5I0v4k34oEPB4R9+QcUqck6fckw8x/TdLN86+CsjkRMbrcdVZDV8/RJC2zL0bEy5K2JZ/ZjNY+qyLiHUlr05EhS8h/5FHjV+/X026ol0kmmFmBSCZNvmePAHuPy6OZvRQAskj6UAWJP5Lp2T8sOH6BZJE06xrmSNqEpL9zLrASmJVrRDBF0qbAt4FbSRZCKzZ71qxFEXFvpe9jVUNXT3MbsayMiI1zDczaTFI9UBf5L6NrVjZ53Mfq9i1+mt+IZWiuEVmrJLV4813SqDzGXEs6o1h5RPywWLlZC47i3ftYX2i8j5VlhdWQ+IlOtBGLlex/i5QFkMeN1M66s5R1bRW/j1UNib+zbcRiJYiI/fKOoamIuCDvGKxbqvh9rGro49+O5BO0F+9uxHJFRDxT9I3WaUjai/cuy5znBK4dSCYGbhkRu6Qrlx4eEf+TV0zWPVTqPla3T/zWtUn6Ncm+BfN5d1nmiIiv5hjTAyT7o14Z766f/lhE7JJXTNb1FLuPBZDlfaxu29WTTvxp6VMtIsNNDqysRpPsotaZWijvi4hZ0gb7mK/NKxjrsuYAj5NMKoUNV8nN9D5Wt038wGHNnBMwiGR1TOsaHgO2Irk301m8mm6W0bhxxlF0rvisaziTZOvHVSS7gU2r1MqqVdHVI2kkyezdz5BsUHFTRFyea1BWEkn3kWyNN4sN177PbRc1SUNIlkDeC/g3yb+pYyPiubxisq5L0mCSxdk+CTwPXBQR87Oss9u2+NMbcJ8l+YUuA24g+aDrdKNFrKjz8w6gqYhYCBwgqR/Qo3GzEbP2iIhFkm4hWcTuOGAHkntamem2Lf50jfIZJGv0PJOeW5jnBh7WPunIrKERcbek95HsnZpbspV0EXBJRLyeHm8KnBkR38orJut60m+OnyVp6f+TpLvnjxGxOvO6u3HiP4Lkl7oXyeqcvyVZ1TGPLd+snSR9CZgEbBYR20saCkyOiHE5xvSeDXM6285l1vmljdNHgFuAFTQZjJLlTPBu29UTEdOAaenX8f8gGcO/Zbrn57SIuDPP+KxkXwE+DPwNICL+IWlgviFRI6lPRLwFIKkv0CfnmKzruZB3k33/SlbcbRN/o3RfzanA1HRD6k8D5wBO/F3DWxHxduPQSUk9yWmnsgK/Ae5J994N4ASSDc7NShYR5+dVd7ft6rHuIV1m43WSxfVOBf4TeCIivplzXIcA40iGCN8ZEXfkGY9ZWzjxW6cmqQfJ/rYHkiTZO0ju1fgfrlk7OfFbp5f2oW8bEU/lHEcDzXcziWQ2eF2FQzJrFyd+69QkHU6yVWbviBicTsa7MM8JXGblJKkPyQzeejZciPDCrOrskdWFzcrkPJJRPa8DpDMa6/MLx6zsbiEZy78WeKPgkZluP6rHury1EbG8yYJoZt3JoIg4uJIVusVvnd1jkiaQjJ0fKukywDuoWXfyoKThlazQffzWqaVLNHyTZFQPJKN6/qcS09rNKkHSE8AHSRb7e4t3BwvsmlmdTvzWWUmqAe6IiAPyjsUsK+laVO8REc9nVae7eqzTioh1JHsmb5x3LGZZiYjn0yS/imS4cOMjM765a53dauBRSXdRMNIhz60XzcopHbL8v8D7SfYH3w54Etg5qzqd+K2z+1P6gHdbQR7iY93Jd4A9gbsjYjdJ+5HsI5IZJ37rlCR9kmSY20/T41nAFiTJ/+w8YzMrszURsUxSD0k9IuI+SRdnWaETv3VWXyfZT6FRb2B3kuVrrwJuzCMoswy8Lqk/MJ1kFeElJJO5MuObu9ZZ9Y6IfxYcz4yI1yLiBaBfXkGZZeCTJDd2v0ayadSzwPgsK/RwTuuUJD0TER9soezZiNi+0jGZdRdu8Vtn9bd028UNSPoyMCuHeMzKSlKDpBUtPJZK+qukTLYYdYvfOqV0e8WbSWYyzktP706yxeF/RMQrOYVmlrl08uIuwNSI2KXs13fit85M0v68O5758Yi4N894zCpJ0pcj4sqyX9eJ38ysuriP38ysyjjxm5lVGU/gMjPLQZE9nAHIcg9nJ34zsxxExEYAki4EXgZ+TbIO1THARlnW7Zu7ZmY5kvS3iNijtXPl5D5+M7N8rZN0jKSadKG2Y4B1WVboxG9mlq8JwGeAV9LHp9NzmXFXj5lZlfHNXTOzHEnaAvgSUE9BTo6IE7Kq04nfzCxftwAzgLvJuG+/kbt6zMxyJGl+RIysZJ2+uWtmlq8/Sjq0khW6xW9mlqN0Bm8/kiXI15BM4oosZ+468ZuZVRnf3DUzy5mkTYGhQG3juYiYnlV9TvxmZjmSdCJwGjAImA/sCTwE7J9Vnb65a2aWr9OAMcDzEbEfsBuwNMsKnfjNzPK1OiJWA0jqExELgB2zrNBdPWZm+VosaRPgZuAuSf8GXsqyQo/qMTPrJCR9DNgYuD0i3s6sHid+M7Pq4j5+M7Mq48RvZlZlnPjNzKqMR/WYmeUgXaOnxZusWa7V48RvZpaDiNgIQNKFwMvAr0kWaDsG2CjLuj2qx8wsR5L+FhF7tHaunNzHb2aWr3WSjpFUI6mHpGPIeCcuJ34zs3xNAD4DvJI+Pp2ey4z7+M3MciKpBvhKRHyykvW6xW9mlpOIWAfsXul63eI3M8vXw5JuBW4E3mg8GRG/z6pCJ34zs3xtBixjw41XAsgs8Xs4p5lZlXGL38wsR5JqgS8CO7PhnrsnZFWnb+6ameXr18BWwEHAAyR77zZkWaG7eszMciTp4YjYTdIjEbGrpF7AHRHhzdbNzLqpNenP1yXtQrIDV32WFbqP38wsX1MkbQp8G7gV6J8+z4y7eszMqoy7eszMciRpgKTLJM2TNFfSjyQNyLJOJ34zs3z9FlgCHAkcBbwK3JBlhe7qMTPLkaS5EbF7k3NzImJ0VnW6xW9mlq/7JH02XYu/h6TPAH/KskK3+M3MclCw566Afry7+UoNsDLLPXed+M3MqozH8ZuZ5UDSqGLlETEvs7rd4jczqzxJ9xUpjiyXbHDiNzOrMu7qMTPLWbpGz4fYcFnmazOrzy1+M7P8SDoP2Jck8d8GHALMjIijsqrT4/jNzPJ1FDAOeDkivgCMAPpkWaETv5lZvlZFxDvAWkl1JMs3DMmyQvfxm5nla46kTYCfA3OBlcCsLCt0H7+ZWSchqR6oi4hHMq3Hid/MLD+SbiFZjfOWiHijEnW6j9/MLF8/BPYGnpB0o6SjJNW29qaOcIvfzKwTkFQD7A98CTg4y0XafHPXzCxnkvoC44GjgVHANZnW5xa/mVl+JN0A7AHcDvwfcH86vDO7Op34zczyI+lg4K6IWNfqi8tVpxO/mVnlSRoD/DMiXk6PP0+y7+7zwPkR8VpWdXtUj5lZPq4E3gaQNBb4PnAtsByYkmXFvrlrZpaPmoJW/dHAlIi4CbhJ0vwsK3aL38wsHzWSGhvf44B7C8oybZS7xW9mlo/rgQckvQqsAmYASPogSXdPZnxz18wsJ5L2BLYG7mxcrkHSDkB/77lrZmZl4z5+M7Mq48RvZlZlnPjNOkjShZIOyDsOs1K5j9+sgCSR/L/IdK0Uszy5xW9VT1K9pCclXQHMA74tabakRyRdUPC6b0taIOkuSddLOis9f7Wko9Ln4yQ9LOlRSb+S1Cc9/5ykCyTNS8uG5fFnNQMnfrNGO5JMlz8b2Ab4MDAS2F3SWEmjSdZR2Q34FDC66QXSzTOuBo6OiOEk82ROLnjJqxExCvgZcFZmfxKzVjjxmyWej4i/Agemj4dJWv/DgKEkOyTdEhGrIqIB+EMz19gRWBQRT6fH1wBjC8p/n/6cC9SX/U9gViLP3DVLNO51KuB7EXFlYaGkr5VwDbVS/lb6cx3+v2c5covfbEN3ACdI6g8gaRtJA4GZwHhJtWnZJ5p57wKgPp1yD3Ac8EAlgjZrC7c6zApExJ2SdgIeSgb4sBI4NiJmS7oV+DvJeulzaLKeSkSslvQF4MZ08a3ZwOSK/gHMSuDhnGYlktQ/IlZKeh8wHZiU5XoqZllxi9+sdFMkfQioBa5x0reuyi1+M7Mq45u7ZmZVxonfzKzKOPGbmVUZJ34zsyrjxG9mVmWc+M3Mqsz/B9tdI/YMrkb3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_uncertain_early.plot.bar(x='region',stacked=True,ylim=(0,10),ylabel='squared unceratinties (Gt2 a-2)', title='2000-2010')\n",
    "# plt.set_xlim(0,10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "id": "466aadca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:title={'center':'2010-2020'}, xlabel='region', ylabel='squared unceratinties (Gt2 a-2)'>"
      ]
     },
     "execution_count": 197,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAGFCAYAAADw//QxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAy0UlEQVR4nO3deZyVdd3/8debARkCR3NBTVTAME1RRHAryS23RG/T0jCXzEhv7dbUbtMWl+6806zfXZoilVupeXsbSmaaO1AZm7hrobjgipgwKCjg5/fHdQ0expkz18C5znVmzvv5eJzHnGv9fmbE8znf67spIjAzs/rVo+gAzMysWE4EZmZ1zonAzKzOORGYmdU5JwIzszrnRGBmVuecCMzM6pwTgXUrknpL+rWk5yU1S3pI0v4lx/eS9JSkdyTdJ2mzkmN7pPsWSHqujXsPTI+/k95j7zJx9Jd0g6SX0/v9RdJOrc4Zk8b5tqRbJK1TcuxiSf9Mf4enJB3d6tphkmakscyQNGzV/mJmTgTW/fQEXgQ+A6wFfA/43/RDfD3g9+m+dYDpwI0l174NXAl8q5173wA8BKwLfAf4P0nrt3NuP2AasENa1jXAHyX1A5C0NXAFcBSwAfAOcFmrWEanv8MxwM8k7ZpeuwZwK/Bb4KPpvW9N95t1mjyy2Lo7SY8A55F8gB8bES0fqH2BN4DtI+KpkvP3Bn4VEQNL9m0BPAqsFxHN6b7JwHURMS5jHAuBPSJihqQLgIERMSY9tjnwJLBuy/1bXTsReCAifiJpH+AqYECk/wNLegEYGxF3dOZvYwauEVg3J2kDYAvgcWBr4OGWYxHxNvBMur8jWwPPtvqQfjjjtaSPbtYAZpfcrzSWZ4D30lhbX9sHGJn+Di3XPhIrf4t7JGssZq05EVi3JakXcB1wTfqNvx+woNVpC4A1M9xula+V1AT8BjgvIlru0Zn7jSNJGneubixmbXEisG5JUg+SD9/3gJPT3YuAplanNgEfehTThrLXSnpc0qL0tVtJHH2APwAPRsR/Z71fyfU/BrYBvlhSA1id38PsQ5wIrNuRJODXJI2wh0bE0vTQ48B2Jef1BTbng0cu5TwODJZU+q17u5ZrI2LriOiXvian9+8N3AK8BHy9jfuVxjIY6A38o2TfecD+wD4RsbDVtdumv2eLbTP+HmYf4kRg3dHlwFbA6IhYXLJ/ArCNpEMlNQLfJ3nW/hQktYh0f69kU40tPXEi4h/ALOCcdP8hJB++N7cVQPpY6v+AxcDREfF+q1OuA0ZL2i1NSOcDvy9piD4LGAN8NiLmt7r2fmA58B9pd9mWGs+9nfgbmX0gIvzyq9u8gM2AAJaQPEJpeR2ZHt8beIrkA/p+kp47Ldfunl5b+rq/5PjA9JrFwNPA3mXi+Ex6/Tut4tit5JwxwAskXUVvBdYpORbAu62uPbvk+PbAjDSWmSQ9nwr/+/vVNV/uPmpmVuf8aMjMrM7llggkbZIOx38y7VFxShvn7J4Ov5+Vvr6fVzxmZta2njneexlwekTMTHtazJB0V0Q80eq8yRFxYI5xmJlZGbnVCCLilYiYmb5vJhk+v3Fe5ZmZ2aqpShuBpIEkvRz+3sbhXSQ9LOlP6URcZmZWRXk+GgIgnW3xZuDUWHlQDCTd3jaLiEWSDiAZfDOkjXuMBcYC9O3bd4ctt9wy36DNzLqZGTNmvBERbc6Wm2v30XRQzW3AnRHx0wznPweMiIg32jtnxIgRMX369MoFaWZWByTNiIgRbR3Ls9dQyzD/J9tLApI2bBkmL2nHNJ7WoyjNzCxHeT4a+hTJohuPSpqV7jsb2BQgkjncDwNOlLSMZITkEeERbmZmVZVbIoiIKYA6OOdS4NK8YjAzs47l3lhsZlYpS5cuZe7cuSxZsqToUGpWY2MjAwYMoFevXpmvcSIwsy5j7ty5rLnmmgwcOJCVZ+E2SCYRnT9/PnPnzmXQoEGZr/NcQ2bWZSxZsoR1113XSaAdklh33XU7XWNyIjCzLsVJoLxV+fs4EZiZ1Tm3EZhZlzXw23+s6P2e+9HnOjynoaGBoUOHsnTpUnr27MkxxxzDqaeeSo8eXfd7tROBmVkn9OnTh1mzZgHw+uuvM2bMGBYsWMB5553X4bXLli2jZ8/a+9jtuinMzKxg/fv3Z/z48Vx66aW0Nxb26quv5gtf+AKjR49mn332YdGiRey1114MHz6coUOHcuuttwJw0UUX8fOf/xyAb37zm+y5554A3HPPPXz5y1/O9fdwIjAzWw2DBw/m/fff5/XXX2/3nL/97W9cc8013HvvvTQ2NjJhwgRmzpzJfffdx+mnn05EMGrUKCZPngzA9OnTWbRoEUuXLmXKlCnstttuuf4OtVdHMTPrYjqaGeezn/0s66yzzopzzz77bCZNmkSPHj146aWXeO2119hhhx2YMWMGzc3N9O7dm+HDhzN9+nQmT568oqaQFycCM7PV8Oyzz9LQ0ED//v3bPadv374r3l933XXMmzePGTNm0KtXLwYOHMiSJUtWvL/qqqvYdddd2Xbbbbnvvvt45pln2GqrrXL9HfxoyMxsFc2bN48TTjiBk08+OXP//QULFtC/f3969erFfffdx/PPP7/i2KhRo7j44osZNWoUu+22G+PGjWPYsGG5j51wjcDMuqws3T0rbfHixQwbNmxF99GjjjqK0047LfP1Rx55JKNHj2bEiBEMGzaM0oW2dtttN374wx+yyy670LdvXxobG3NvH4CcF6bJgxemMatfTz75ZO6PSbqDtv5OhSxMY2ZmXYMfDZmZVcCdd97JmWeeudK+QYMGMWHChIIiys6JwMysAvbdd1/23XffosNYJX40ZGZW55wIzMzqnBOBmVmdcyIwM6tzbiw2s67r3LUqfL8FZQ/Pnz+fvfbaC4BXX32VhoYG1l9/fQDGjBnDlVdeSc+ePWloaOD000/n6KOPrmx8OXEiMDPLaN11112xFsG5555Lv379OOOMMxg3bhwTJkxg6tSpNDU1sWDBAm655ZZM91y+fDkNDQ35BZ2BHw2Zma2mCy64gMsuu4ympiYA1lprLY455ph2zx84cCDnn38+n/70p7npppv45S9/yciRI9luu+049NBDeeedd1i+fDmDBw8mInjrrbfo0aMHkyZNApKpKGbPnl2x+J0IzMxWQ3NzM83NzWy++eaduq6xsZEpU6ZwxBFH8PnPf55p06bx8MMPs9VWW/HrX/+ahoYGtthiC5544gmmTJnCDjvswOTJk3n33XeZO3cuH//4xyv2OzgRmJmthohYpdlBDz/88BXvH3vsMXbbbTeGDh3Kddddx+OPPw4k3/wnTZrEpEmTOOuss5gyZQrTpk1j5MiRFYsfnAjMzFZLU1MTffv25dlnn+3UdaVrFBx77LFceumlPProo5xzzjksWbIESBLB5MmTmTp1KgcccABvvfUW999/P6NGjaro71A2EUgaIOkMSbdKmiZpkqTLJH1OkpOImRlw1llncdJJJ7Fw4UIAFi5cyPjx4zNf39zczEYbbcTSpUu57rrrVuzfaaed+Otf/0qPHj1obGxk2LBhXHHFFRWfmrrdXkOSrgI2Bm4DLgReBxqBLYD9gO9I+nZETKpoRFZdlex+10HXO7OKq5F/cyeeeCKLFi1i5MiR9OrVi169enH66adnvv4HP/gBO+20E5ttthlDhw6lubkZgN69e7PJJpuw8847A0kN4YYbbmDo0KEVjb/d9QgkbRMRj7V7obQGsGlEVK7pOgOvR1BhTgTWhXg9gmw6ux5BuzWCckkgPf4eUNUkYGZmlVfu0VATcDbJ46E/RcT1Jccui4h/r0J8ZmZd1iGHHMKcOXNW2nfhhRfW3HTV5UYWXwX8E7gZOE7SocCYiHgX2LkawZmZdWVdYVEaKN9raPOI+HZE3BIRBwEzgXslrVul2MzMrArK1Qh6S+oREe8DRMQPJc0FJgH9qhKdmZnlrlyN4A/AnqU7IuIa4HTgvY5uLGkTSfdJelLS45JOaeMcSfq5pNmSHpE0vLO/gJmZrZ5yvYb+s539dwBDMtx7GXB6RMyUtCYwQ9JdEfFEyTn7p/caAuwEXJ7+NDOzKunUNNSSbouIA7OcGxGvAK+k75slPUnSA6k0ERwMXBvJYIYHJa0taaP0WjOzsoZeU9mBVY8e82iH5zQ0NDB06FCWLl1Kz549OeaYYzj11FPp0aNzky3svvvuXHzxxYwYMYIDDjiA66+/nrXXXnsVI189nV2PYONVKUTSQGB74O9t3O/Fku256b6VEoGkscBYgE033XRVQjAzq4g+ffqsWJPg9ddfZ8yYMSxYsIDzzjtvle95++23Vyi6VdPZ+YIe6mwBkvqRdEE9NSIWtj7cxiUfGuocEeMjYkREjGhZDcjMrGj9+/dn/PjxXHrppbQ3S8Py5cs544wzGDp0KNtuuy2XXHLJh84ZOHAgb7zxBs899xxbbrklxx9/PNtssw1HHnkkd999N5/61KcYMmQIU6dOzeX36FSNICKO68z5knqRJIHrIuL3bZwyF9ikZHsA8HJnyjAzK9LgwYN5//33ef3119lggw0+dHz8+PHMmTOHhx56iJ49e/Lmm2+Wvd/s2bO56aabGD9+PCNHjuT6669nypQpTJw4kQsuuCDzymed0WGNQNIQSf8n6QlJz7a8Mlwn4NfAkxHx03ZOmwgcnfYe2hlY4PYBM+tq2qsNANx9992ccMIJ9OyZfO9eZ511yt5r0KBBDB06lB49erD11luz1157IYmhQ4fy3HPPVTLsFbLUCK4CzgH+H7AH8BXafqTT2qeAo4BHJc1K950NbAoQEeOA24EDSOYseie9t5lZl/Hss8/S0NBA//792zze2YVrevfuveJ9jx49Vmz36NGDZcuWrV6w7ciSCPpExD2SFBHPA+dKmkySHNoVEVPoIGGkvYVOyhytmVkNmTdvHieccAInn3xyux/2++yzD+PGjWP33Xdf8Wioo1pBtWVJBEvSRWj+Kelk4CWg7dRnZlZFWbp7VtrixYsZNmzYiu6jRx11FKeddlq75x9//PH84x//YNttt6VXr1587Wtf4+STT65ixB1rdz2CFSdII4EngbWBHwBNwI8j4sHco2uD1yOoMK9HYF2I1yPIpmLrEbSIiGnp20X4Gb6ZWbfT2QFlZmbWhjvvvJMzzzxzpX2DBg3qElNROxGYmVXAvvvuW3MLzmTV2ZHFZmbWzZRNBJL2lfTVdK6g0v2dGmFsZma1q91EIOkC4DvAUOAeSd8oOVxbfZ/MzGyVlWsjGA1sHxHLJJ0LXC9pcER8k2wji83McvXklpXtSrrVU09W9H5dRblHQz0jYhlARLxFkhiaJN0ErFGF2MzMasr8+fMZNmwYw4YNY8MNN2TjjTdesX3RRRex5ZZbss0227Dddttx7bXXdure999/PwcemCz3MnHiRH70ox/l8Su0qVyN4BlJn4mIBwAiYjnwVUn/BRxalejMzGrIuuuuu2ItgnPPPZd+/fpxxhlnMG7cOCZMmMDUqVNpampiwYIFqzVL6EEHHcRBBx1UmaAzKFcj+ALwocmvI+K7rDx1tJlZXbvgggu47LLLaGpqAmCttdbimGOOaff8adOmseuuu7Lddtux44470tzcvNLxq6++esU0FMceeywnnngie+yxB4MHD+aBBx7guOOOY6uttuLYY4+tSPzl1ixeDCDpnojYq9Xha4HW+8zM6k5zczPNzc1svvnmmc5/7733OPzww7nxxhsZOXIkCxcupE+fPmWv+de//sW9997LxIkTGT16NH/5y1/41a9+xciRI5k1axbDhg1brd+hXK+hRknrAOtJ+qikddLXQOBjq1WqmVk30dlppp9++mk22mgjRo4cCUBTU9OKtQraM3r06BVrEmywwQYrrVdQiTUKypX+deBUkg/9GXzQU2gh8IvVLtnMrBtoamqib9++PPvsswwePLjD8zubOICV1iRovV5BJdYoKPdo6GfAzyR9IyI+vMimmVnBaqW751lnncVJJ53EjTfeSFNTEwsXLuR3v/sdY8eO/dC5W265JS+//DLTpk1j5MiRNDc3d/hoKG/tJgJJTcAGLUlA0heAlmjvjIjXqhCfmVnNO/HEE1m0aBEjR46kV69e9OrVi9NPP73Nc9dYYw1uvPFGvvGNb7B48WL69OnD3XffXeWIV9buegSSxgN/jYir0+3ZwJ9IksGyiDihWkGW8noEFeb1CKwL8XoE2VRyPYKRJO0ELZoj4hvpDaesbqBmZlYbyiWCnrFydeGokvdr5xOOmVn3ccghhzBnzpyV9l144YU1N111uUTwvqQNI+JVgIh4DEDSxsD71QjOzKy1Vel1U5QiFqXpaPnhtpQbWfxj4A+SRklaM319BrglPWZmVlWNjY3Mnz9/lT7s6kFEMH/+fBobGzt1Xbnuo7+V9AbwX8DWQACPA9+PiD+tTrBmZqtiwIABzJ07l3nz5hUdSs1qbGxkwIABnbqm7HC2iLgDuGN1gjIzq5RevXoxaNCgosPodspNMfHddIqJ9o7vKenAfMIyM7NqKVcjeJSkjWAJMBOYBzQCQ4BhwN3ABXkHaGZm+SrXRnArcKukIcCngI1I5hn6LTC2ZXZSMzPr2spPeQdExD+Bf1YhFjMzK0C57qNmZlYHnAjMzOqcE4GZWZ3rMBFIukhSk6Reku6R9IakL1cjODMzy1+WGsE+EbEQOBCYC2wBfCvXqMzMrGqyJIJe6c8DgBsi4s0c4zEzsyrLkgj+IOkpYARwj6T1gSUdXSTpSkmvS3qsneO7S1ogaVb6+n7nQjczs0roMBFExLeBXYAREbEUeAc4OMO9rwb26+CcyRExLH2dn+GeZmZWYVkaiz8CnARcnu76GEntoKyImAT4MZKZWY3L8mjoKuA9YNd0ey7J1NSVsIukhyX9SdLWFbqnmZl1QpZEsHlEXAQsBUjnGKrE8kAzgc0iYjvgEpIFb9okaayk6ZKmex5yM7PKypII3pPUh2RhGiRtDry7ugVHxMKIWJS+vx3oJWm9ds4dHxEjImLE+uuvv7pFm5lZiQ4nnQPOIVmcZhNJ15HMRHrs6hYsaUPgtYgISTuSJKX5q3tfMzPrnCyzj94laSawM8kjoVMi4o2OrpN0A7A7sJ6kuSQJpVd6z3HAYcCJkpYBi4EjwguRmplVXbuJQNKWEfGUpOHprlfSn5tK2jQiZpa7cUR8qYPjlwKXdipaMzOruHI1gtOAscBP2jgWwJ65RGRmZlVVboWysenb/SNipZHEkhpzjcrMzKomS6+hv2bcZ2ZmXVC5NoINgY2BPpK254OxA03AR6oQm5mZVUG5NoJ9SbqJDgB+WrK/GTg7x5jMzKyKyrURXANcI+nQiLi5ijGZmVkVZRlQdpukMcDA0vM9W6iZWfeQJRHcCiwAZlCBqSXMzKy2ZEkEAyKio3UFzMysi8rUfVTS0NwjMTOzQmSpEXwaOFbSHJJHQwIiIrbNNTIzs+7s3LUqeK8Fq3V5lkSw/2qVYGZmNa3cgLKmiFhIMm7AzMy6qXI1guuBA0l6CwUrr0oWwOAc4zIzsyopN6DswPTnoOqFY2Zm1dZhryFJ92TZZ2ZmXVO5NoJGksnl1pP0UVaedO5jVYjNzMyqoFwbwdeBU0k+9GfwQSJYCPwi37DMzKxayrUR/Az4maRvRMQlVYzJzMyqKMvi9ZdI2gb4JNBYsv/aPAMzM7Pq6DARSDoH2J0kEdxOMsBsCuBEYGbWDWSZa+gwYC/g1Yj4CrAd0DvXqMzMrGqyJILFEfE+sExSE/A6HkxmZtZtZJlraLqktYFfkvQeWgRMzTMoMzOrnrKJQJKA/46It4Bxku4AmiLikWoEZ2Zm+Sv7aCgiArilZPs5JwEzs+4lSxvBg5JG5h6JmZkVIksbwR7ACZKeA97GC9OYmXUrXpjGzKzOdfhoKCKeBzYB9kzfv5PlOjMz6xqyTEN9DnAmcFa6qxfw2zyDMjOz6snyzf4Q4CCS9gEi4mVgzTyDMjOz6smSCN5Lu5EGgKS++YZkZmbVlCUR/K+kK4C1JX0NuJtklLGZmXUDWaahvljSZ0kWpPkE8P2IuCv3yMzMrCqyTEM9CJjc8uEvqY+kgRHxXAfXXQkcCLweEdu0cVzAz4ADSHoiHRsRMzv/K5iZ2erI8mjoJuD9ku3l6b6OXA3sV+b4/sCQ9DUWuDzDPc3MrMKyJIKeEfFey0b6fo2OLoqIScCbZU45GLg2Eg+StEFslCEeMzOroCyJYJ6kg1o2JB0MvFGBsjcGXizZnpvu+xBJYyVNlzR93rx5FSjazMxaZEkEJwBnS3pB0oskg8u+XoGy1ca+aOvEiBgfESMiYsT6669fgaLNzKxFll5DzwA7S+oHKCKaK1T2XJKpK1oMAF6u0L3NzCyjLL2GegOHAgOBnklnH4iI81ez7InAyZJ+B+wELIiIV1bznmZm1klZZh+9FVhAskzlu1lvLOkGYHdgPUlzgXNI5ikiIsYBt5N0HZ1N0n30K50J3MzMKiNLIhgQEeW6gbYpIr7UwfEATursfc3MrLKyNBb/VdLQ3CMxM7NCZKkRfBo4VtIckkdDXqHMzKwb8QplZmZ1LksiaLNvv5mZdQ9ZEsEfSZKBgEZgEPA0sHWOcZmZWZVkGVC2UkOxpOFUZmSxmZnVgE4vQp9OFT0yh1jMzKwAWUYWn1ay2QMYDnjmNzOzbiJLG0HpQvXLSNoMbs4nHDMzq7YsbQTnVSMQMzMrRqfbCMzMrHtxIjAzq3NOBGZmda7dNgJJl1BmVHFE/EcuEZmZWVWVqxFMJ1mDoJGky+g/09cwYHnukZmZWVW0WyOIiGsAJB0L7BERS9PtccCfqxKdmZnlLksbwcdYeSxBv3SfmZl1A1kGlP0IeEjSfen2Z4Bzc4vIzMyqKsuAsqsk/YlkgXmAb0fEq/mGZWZm1ZJlriEBewODI+J8SZtK2jEipuYfnuVt4JLrK3av5yp2JzOrpixtBJcBuwAti9E3A7/ILSIzM6uqLG0EO0XEcEkPAUTEvyStkXNcZmZWJVlqBEslNZAOLpO0PvB+rlGZmVnVZEkEPwcmAP0l/RCYAlyQa1RmZlY1ZR8NSeoBzAH+E9iLZN3if4uIJ6sQm5mZVUHZRBAR70v6SUTsAjxVpZjMzKyKsjwa+rOkQ9NupGZm1s1k6TV0GtAXWCZpCcnjoYiIplwjMzOzqsgysnjNjs4xM7OuK0uNAEkfBYaQTEkNQERMyisoMzOrnixTTBwPnAIMAGYBOwN/A/bMNTIzM6uKLI3FpwAjgecjYg9ge2BerlGZmVnVZEkESyJiCYCk3hHxFPCJfMMyM7NqydJGMFfS2sAtwF2S/gW8nGdQZmZWPR3WCCLikIh4KyLOBb4H/Br4tyw3l7SfpKclzZb07TaO7y5pgaRZ6ev7nYzfzMxWU5bG4k1LNuekPzcEXujgugaS6ao/C8wFpkmaGBFPtDp1ckQcmD1kMzOrpCyPhv5IMvOoSLqPDgKeBrbu4LodgdkR8SyApN8BBwOtE4GZmRUoy4CyoaXbkoYDX89w742BF0u25/LBcpeldpH0MEm7wxkR8XiGe5uZWYVkGlBWKiJmShqZ4dS25iaKVtszgc0iYpGkA0gapId86EbSWGAswKabbtr6sJmZrYYsbQSnlWz2AIaTbRzBXGCTku0BtOptFBELS97fLukySetFxButzhsPjAcYMWJE62RiZmarIcs4gjVLXr1J2gwOznDdNGCIpEHp0pZHABNLT5C0YcusppJ2TOOZnz18MzNbXVnaCM5blRtHxDJJJwN3Ag3AlRHxuKQT0uPjgMOAEyUtAxYDR0SEv/GbmVVRlkdDE8sdj4iDyhy7Hbi91b5xJe8vBS7tOEwzM8tLlsbiOSTjBn6bbn8JeI7km76ZmXVxWRLB9hExqmT7D5ImRcTZeQVlZmbVk6WxeH1Jg1s2JA0C1s8vJDMzq6YsNYJvAvdLejbdHki2AWVmZtYFZOk1dIekIcCW6a6nIuLdfMMyM7Nq6fDRkKQvAGtExMPAaOCGdJoJMzPrBrK0EXwvIpolfRrYF7gGuDzfsMzMrFqyJILl6c/PAZdHxK3AGvmFZGZm1ZQlEbwk6Qrgi8DtknpnvM7MzLqALB/oXyQZPLZfRLwFrAN8K8+gzMyserL0GnoH+H3J9ivAK3kGZWZm1dPp9QjMzMo6d60K3WdBZe5jHfKzfjOzOudEYGZW55wIzMzqnBOBmVmdcyIwM6tz7jVktadSvU7APU/MMnCNwMyszjkRmJnVOScCM7M650RgZlbn3Fhs1pV5OgerANcIzMzqnGsEZtbtDb1maMXu9egxj1bsXrXCicDMur1H57xQdAg1zYnAzKwAA5dcX7F7Pbea17uNwMyszjkRmJnVOScCM7M650RgZlbnnAjMzOqcE4GZWZ1z91Ez6/ZqqatmLXKNwMyszuWaCCTtJ+lpSbMlfbuN45L08/T4I5KG5xmPmZl9WG6PhiQ1AL8APgvMBaZJmhgRT5Sctj8wJH3tBFye/lx9npXRzCyTPNsIdgRmR8SzAJJ+BxwMlCaCg4FrIyKAByWtLWmjiHglx7isxvl5bnaV+ls9V5G7JGoxJitPyWdwDjeWDgP2i4jj0+2jgJ0i4uSSc24DfhQRU9Lte4AzI2J6q3uNBcamm58Anq5QmOsBb1ToXpXimLKpxZigNuNyTNl095g2i4j12zqQZ41AbexrnXWynENEjAfGVyKolQqXpkfEiErfd3U4pmxqMSaozbgcUzb1HFOejcVzgU1KtgcAL6/COWZmlqM8E8E0YIikQZLWAI4AJrY6ZyJwdNp7aGdggdsHzMyqK7dHQxGxTNLJwJ1AA3BlRDwu6YT0+DjgduAAYDbwDvCVvOJpR8UfN1WAY8qmFmOC2ozLMWVTtzHl1lhsZmZdg0cWm5nVOScCM7M650RgZlbnPPuoZSbpYxHh7r1mOZG0BfAtYDNKPp8jYs9cy3VjcXEk/YE2BtC1iIiDqhhOhyS9EBGbFlDuOuWOR8Sb1YqllKTPlzseEb+vViwtJP0N+E5E3NvGsXsiYq9qx1SLJP038P8i4vU2jl0YEWcWEBaSHgbGATOA5S37I2JGruXWUyJIxypcAmwFrEHSrfXtiGgqKJ7PlDseEQ9UK5YsJL0YEZt0fGbFy51DkjAFbAr8K32/NvBCRAyqdkxpXFelb/sDuwItH757APdHRNlEkVNMLwELSLpmnxURS0uOPRQR21c5nmba/rIjIAr8f28B8Bbw7xHxx1bHZkZEITMhS5oRETtUu9x6ezR0KcnAtpuAEcDRwMeLCqbWPugzKORbQ8sHvaRxwMSIuD3d3h/Yu4iY0ri+ksZxG/DJlsGQkjYimXm3CK8BnwZ+Dvxd0pciomVurqr/94uINatdZkbPAkcB16X/js6IiCXpsbamvqmWP0j6d2AC8G7LzrxrvfWWCIiI2ZIaImI5cJWkvxYVi6QhwNkk33B/CvwS2A14Bvhq68n3qhTTJbT/DW7t6kbzISMj4oSWjYj4k6QfFBlQamCrEfGvAVsUFUxEvAMcL+kQ4C5JF6QDOKv+ASdpz5bHVJIGRcSckmOfL+LxWSoi4jFJOwEXAdMljYmIRwqKp8Ux6c9vlewLYHCehdZbIngnne5ilqSLgFeAvgXGcxVwLdAE/B04FTiEJBn8gkqtzdA55ZJP1RNTK29I+i7wW5L/Ob4MzC82JADul3QncANJXEcA9xUbEkTEBElTgaslHQD0KyCMi4GWxyw3l7wH+C5QVCIAIK0F/EdaK/iDpJ8XHE8xjznrrI1gM5Jva2sA3wTWAn4REc8UFM+siBiWvp8dER9v65gl0kbjc4BR6a5JwHlFNRaXShuOd0s3J0XEhILiuCMi9mtj/7eAH0REY5XjWdEu0bqNoog2i3JlS1oPuBL4XEQ0FBTXR4DTgE0jYmz61OATEXFbnuXWW43gkxHxPLAEOA8gnfuokEQAvF/yfmGZY1XTQU+md0n+Vr+IiBerF1Ui/cA/pdrlZpE+4ij0220ax4eSQLr/x8CPqxwOrPxvqfW/qyK/hZ4iacOIeBVA0tHAoSTr4Ywqd2HOriLpMbRruj2XpE0z10RQbzWCvwLfLXlm+Z/AHhGxf0HxvEMy4Z6AzdP3pNuDI6Lqj6066MnUE9ga+FJE7FKlkFZI+1ifAQykin2sO5LWBi4k6T0kCuwRI+n7ZQ5HRFS1TUXSWyQ1N5HUmCa1HAI+HREfrWY8JXHNBPaOiDcljQJ+B3wDGAZsFRGHFRTX9IgY0aom9XBEbJdnufVWIzgIuC2tJu8HbJnuK8pWBZbdnjkR8UKZ4/dI2rZq0azsJpI+1r+ipI91DbgIGB0RTxYdCPB2G/s+AhwPrAtUu3H94PRnX+BuklrA0yS18ourHEuphpJHiocD4yPiZuBmSbOKC4v3JPUhrS1J2pyS3kO5iYi6epF8a3uEpAqmouNpI771iowLmFny/uai/x6tYptRdAztxPWXomNoJ641SRpk55DWWAqIYQ3gf4A3gZnAQ8BLwLfT49sX9Ld5DOiZvn8KGFV6rMD/ZvsADwDzgOtIHlXtnne5dVEjKBnUovTnGiTdsQ6TFFHcoJadgR+R/E/yA+A3JImgh6SjI+KOIsIqeZ9rl7VVUEgf6wymS7oRuIWV4yqkzSBtVD8NOBK4BhgeEf8qIhaSb/0fIVkvtzmNrwm4WNLlJDXzInrK3AA8IOkNYDEwOY3t4yQD8goREX+WNAPYmeT/xVMiIvd1lOuqjaDWSJpOMo5gLZIFKPaPiAclbQncEAX0qCgdVVnkCMu2pCOMW4uIKDRhlYwwLhURcVwBsfwY+DzJv6dfRMSiasfQKp7ZwJBo9UEjqYFkUfb9I+LBgmLbGdgI+HNEvJ3u2wLoFxEzC4rp/0h6Lt0REVXrMFJXiUDSp4BZEfG2pC+T9Gn+nyj/TDzPeGbFB91Hn4yIrUqOPVRQIlhO8pxZQB+SleOg4CkBLBtJ75PUSpaxcq+cQv77SfpHRLQ5uK7csXolaW+SlRp3JmkTuzoinsq73Lp4NFTicmA7SdsB/wn8muRxTNk5f3JUmvEXtzpW1HQOhfSfzkrSNsAngRX94SPi2uIiAkmNwFdJelSVxlX1GkFE1NrU8k+kjzlX+m+UfhGrhcb1mhIRdwN3S1oL+BLJyPAXSWYd+G2UzB1VSfVWI5gZEcPTLnYvRcSvC55gqty378aI6FVEXLVK0jnA7iSJ4HZgf2BKFNTVrySum0gaHMcA55M8m38yImpyzEM1SdqYZHzFYpL+8QGMJPn3fkhEvFRgeDVJ0roko+aPAl4maTT+NDA0InbPpcw6SwQPAHeQVL1GkbTMz4qIoYUGZplIehTYDngoIraTtAHwq4gYXXBcD0XE9pIeiYhtJfUC7oyCxzfUEkl7ktSYBDweEfcUHFJNkvR7km7tvyF5LPRKybHpETEij3Lr7dHQ4STf2r4aEa9K2pRiRlvaqlkcEe9LWpb2PHmd2ujZ1FJdfyt9dPUqyaA3S0UyiPNDayTYh1wabawlAZBXEoA6SwSRDCf/acn2CySTvlnXMF3S2iTPS2cAi4CphUaUGC/po8D3gIkkk7uVG+Fr1qaIuLeIdrB6ezTU1sI0iyJirUIDs06TNBBoiuKnDTarmKLaweqqRkDbC9MMKTQi65CkdhvzJQ0vsM/3aeWOR8RPyx03a8NhfNAO9pWWdrC8C623REDU0MI0ltlPyhwLoKhG2Vpdfcu6rkLaweotEdTawjSWQUTsUXQMbYmI84qOwbqdQtrB6q2NYDOSDNuLDxamuSwiZpe90GqGpF358DTURQ8o24JksOIGEbFNOjvrQRHxX0XGZV1bNdvB6ioRWNcm6Tck6zbM4oNpqCMi/qOwoFgxPuVbwBXxwRzyj0XENkXGZV1HuXYwgLzbweri0VA6EKm9jBeR86IPVjEjSFaZq7VvLx+JiKnSSmvDLysqGOuSpgOPkwxyhZVnAc69HawuEgFwYBv7BAwgmf3TuobHgA1J2nZqyRvpAiIti4kcRu3FaLXtdJKlMheTrJY2oZozx9bdoyFJw0hGF3+RZMGOmyPi0kKDskwk3UeylOBUVp73v8hV5pA0mGTa512Bf5H8u/pyRDxXZFzW9UgaRDLZ3MHA88AFETEr73LrokaQNuYdQfIHng/cSJIEa7I3irXr3KIDaEtEPAvsLakv0KNlARazzoqIOZJuJZmU7yhgC5I2sVzVRY0gnaN9MskcQ7PTfc8WvaCJdV7a82tIRNwt6SMka88W+sEr6QLgooh4K93+KHB6RHy3yLis60hrlUeQ1AReJHk8dFtELKlK+XWSCA4h+SPvSjL76O9IZq0sYok8W0WSvgaMBdaJiM0lDQHGRcReBcf1oUWEam11N6tt6ZfVR4BbgYW06tyS9yj1ung0FBETgAlp1f3fSMYQbJCumTohIv5cZHyW2UnAjsDfASLin5L6FxsSAA2SekfEuwCS+gC9C47Jupbz+eDDv1+1C6+LRNAiXZf0OuC6dIHvLwDfBpwIuoZ3I+K9lm6aknpS0EpurfwWuCdduziA40gWjTfLJCLOLbL8ung0ZN1DOi3IWySTBX4D+HfgiYj4TpFxAUjaH9iLpFvynyPizoJDMsvMicC6DEk9SNYG3ofkA/dOkrYe/yM2Ww1OBNalpM/fN42Ip2sglmbafjQlkhHrTVUOyWyVOBFYlyHpIJKlRdeIiEHp4MDzix5QZlYpknqTjDAeyMoTK56fZ7k98ry5WYWdQ9Jr6C2AdMTlwOLCMau4W0nGEiwD3i555aqueg1Zl7csIha0mtzNrDsZEBH7VbtQ1wisK3lM0hiSfvtDJF0CeIU5607+KmlotQt1G4F1GemUEt8h6TUESa+h/6rWMHyzvEl6Avg4ycSF7/JBx4Ntcy3XicC6AkkNwJ0RsXfRsZjlJZ1L60Mi4vk8y/WjIesSImI5yZrTaxUdi1leIuL59EN/MUnX5JZXrtxYbF3JEuBRSXdR0pOi6KUqzSol7SL9E+BjJOurbwY8CWydZ7lOBNaV/DF9wQffktyFyLqTHwA7A3dHxPaS9iBZRyVXTgRW8yQdTNKt7hfp9lRgfZJkcGaRsZlV2NKImC+ph6QeEXGfpAvzLtSJwLqC/yRZT6LFGsAOJNP1XgXcVERQZjl4S1I/YBLJLMmvkwwuy5Ubi60rWCMiXizZnhIRb0bEC0DfooIyy8HBJA3F3yRZROsZYHTehbr7qNU8SbMj4uPtHHsmIjavdkxm3YlrBNYV/D1dpnIlkr4OTC0gHrOKktQsaWE7r3mSHpSU25KsrhFYzUuXo7yFZKTlzHT3DiTLQf5bRLxWUGhmuUsHU24DXBcR2+RShhOBdRWS9uSD/tSPR8S9RcZjVk2Svh4RV+RybycCM7P65jYCM7M650RgZlbnPKDMzKxgZda/BiDv9a+dCMzMChYRawJIOh94FfgNyTxaRwJr5l2+G4vNzGqEpL9HxE4d7as0txGYmdWO5ZKOlNSQTjx3JLA870KdCMzMascY4IvAa+nrC+m+XPnRkJlZnXNjsZlZjZC0PvA1YCAln88RcVye5ToRmJnVjluBycDdVKFtoIUfDZmZ1QhJsyJiWLXLdWOxmVntuE3SAdUu1DUCM7MakY4w7ksy5fpSkkFlkffIYicCM7M658ZiM7MaIumjwBCgsWVfREzKs0wnAjOzGiHpeOAUYAAwC9gZ+BuwZ57lurHYzKx2nAKMBJ6PiD2A7YF5eRfqRGBmVjuWRMQSAEm9I+Ip4BN5F+pHQ2ZmtWOupLWBW4C7JP0LeDnvQt1ryMysBkn6DLAWcEdEvJdrWU4EZmb1zW0EZmZ1zonAzKzOORGYmdU59xoyMytYOsdQuw22ec815ERgZlawiFgTQNL5wKvAb0gmnDsSWDPv8t1ryMysRkj6e0Ts1NG+SnMbgZlZ7Vgu6UhJDZJ6SDqSKqxU5kRgZlY7xgBfBF5LX19I9+XKbQRmZjVAUgNwUkQcXO2yXSMwM6sBEbEc2KGIsl0jMDOrHQ9JmgjcBLzdsjMifp9noU4EZma1Yx1gPisvRBNAronA3UfNzOqcawRmZjVCUiPwVWBrVl6z+Lg8y3VjsZlZ7fgNsCGwL/AAydrFzXkX6kdDZmY1QtJDEbG9pEciYltJvYA7I8KL15uZ1Yml6c+3JG1DskLZwLwLdRuBmVntGC/po8D3gIlAv/R9rvxoyMyszvnRkJlZjZC0rqRLJM2UNEPS/0haN+9ynQjMzGrH74DXgUOBw4A3gBvzLtSPhszMaoSkGRGxQ6t90yNiRJ7lukZgZlY77pN0RLoWQQ9JXwT+mHehrhGYmRWsZM1iAX35YDGaBmBR3msWOxGYmdU5jyMwMyuYpOHljkfEzFzLd43AzKxYku4rczjynmLCicDMrM750ZCZWQ1J5xj6JCtPQ31trmW6RmBmVhsknQPsTpIIbgf2B6ZExGF5lutxBGZmteMwYC/g1Yj4CrAd0DvvQp0IzMxqx+KIeB9YJqmJZLqJwXkX6jYCM7PaMV3S2sAvgRnAImBq3oW6jcDMrAZJGgg0RcQjuZflRGBmVhsk3Uoy2+itEfF2tcp1G4GZWe34KfBp4AlJN0k6TFJjRxetLtcIzMxqjKQGYE/ga8B+eU8658ZiM7MaIqkPMBo4HBgOXJN7ma4RmJnVBkk3AjsBdwD/C9yfdifNt1wnAjOz2iBpP+CuiFje4cmVLNeJwMysWJJGAi9GxKvp9tEk6xY/D5wbEW/mWb57DZmZFe8K4D0ASaOAHwHXAguA8XkX7sZiM7PiNZR86z8cGB8RNwM3S5qVd+GuEZiZFa9BUssX872Ae0uO5f6F3TUCM7Pi3QA8IOkNYDEwGUDSx0keD+XKjcVmZjVA0s7ARsCfW6aXkLQF0M9rFpuZWa7cRmBmVuecCMzM6pwTgVmFSTpf0t5Fx2GWldsIzMqQJJL/T3Kf78WsKK4RmLUiaaCkJyVdBswEvidpmqRHJJ1Xct73JD0l6S5JN0g6I91/taTD0vd7SXpI0qOSrpTUO93/nKTzJM1Mj21ZxO9qBk4EZu35BMkQ/zOBjYEdgWHADpJGSRpBMhfM9sDngRGtb5AuKHI1cHhEDCUZt3NiySlvRMRw4HLgjNx+E7MOOBGYte35iHgQ2Cd9PURSO9gSGEKyitStEbE4IpqBP7Rxj08AcyLiH+n2NcCokuO/T3/OAAZW/Dcwy8gji83a1rJerID/jogrSg9K+maGe6iD4++mP5fj/xetQK4RmJV3J3CcpH4AkjaW1B+YAoyW1Jge+1wb1z4FDEynCQA4CnigGkGbdYa/hZiVERF/lrQV8LekAxGLgC9HxDRJE4GHSeaMn06rOWEiYomkrwA3pROKTQPGVfUXMMvA3UfNVpGkfhGxSNJHgEnA2LznhDHLg2sEZqtuvKRPAo3ANU4C1lW5RmBmVufcWGxmVuecCMzM6pwTgZlZnXMiMDOrc04EZmZ1zonAzKzO/X8Up/fxYk/vOQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_uncertain_late.plot.bar(x='region',stacked=True,ylim=(0,2.5),ylabel='squared unceratinties (Gt2 a-2)', title='2010-2020')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de746410",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30a544ea",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "16160fcb",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/willkochtitzky/opt/miniconda3/envs/gis/lib/python3.7/site-packages/ipykernel_launcher.py:22: UserWarning: Boolean Series key will be reindexed to match DataFrame index.\n",
      "/Users/willkochtitzky/opt/miniconda3/envs/gis/lib/python3.7/site-packages/ipykernel_launcher.py:23: UserWarning: Boolean Series key will be reindexed to match DataFrame index.\n"
     ]
    }
   ],
   "source": [
    "df_all_early = pd.read_pickle('df_all_early_'+date_string+'.pkl').sort_values('rgiid').reset_index(drop=True)\n",
    "df_all_late = pd.read_pickle('df_all_late_'+date_string+'.pkl').sort_values('rgiid').reset_index(drop=True)\n",
    "\n",
    "Arctic_discharge_fluxes_with_area = Arctic_discharge_fluxes_with_area.sort_values('RGIId')\n",
    "\n",
    "# #matusevich changes\n",
    "# subset_index = (Arctic_discharge_fluxes_with_area['RGIId']=='RGI60-09.00515_2')|(Arctic_discharge_fluxes_with_area['RGIId']=='RGI60-09.00498_3')|(Arctic_discharge_fluxes_with_area['RGIId']=='RGI60-09.00498_2')|(Arctic_discharge_fluxes_with_area['RGIId']=='RGI60-09.00498')|(Arctic_discharge_fluxes_with_area['RGIId']=='RGI60-09.00515')|(Arctic_discharge_fluxes_with_area['RGIId']=='RGI60-09.00968_2')|(Arctic_discharge_fluxes_with_area['RGIId']=='RGI60-09.00968')|(Arctic_discharge_fluxes_with_area['RGIId']=='RGI60-09.00967') #retained by terminus\n",
    "\n",
    "#get a subset of df_all_early and df_all_late here\n",
    "# subset_index = Arctic_discharge_fluxes_with_area['terminus_gt_change_per_year_late']>0 #retained by terminus\n",
    "# subset_index = Arctic_discharge_fluxes_with_area['terminus_gt_change_per_year_late']<0 #lost by terminus\n",
    "\n",
    "#Just look at frontal ablation of advancing glaciers\n",
    "# subset_index = (Arctic_discharge_fluxes_with_area['net_area_change_per_year_early']>0)&\\\n",
    "#     (Arctic_discharge_fluxes_with_area['net_area_change_per_year_early']>Arctic_discharge_fluxes_with_area['net_area_change_per_year_early_error'])   \n",
    "# subset_index = (Arctic_discharge_fluxes_with_area['net_area_change_per_year_late']>0)&\\\n",
    "#     (Arctic_discharge_fluxes_with_area['net_area_change_per_year_late']>Arctic_discharge_fluxes_with_area['net_area_change_per_year_late_error'])\n",
    "\n",
    "# subset_index = Arctic_discharge_fluxes_with_area['FA_not_calculated_in_normal_way_early'] == 1\n",
    "subset_index = Arctic_discharge_fluxes_with_area['FA_not_calculated_in_normal_way_late'] == 1\n",
    "\n",
    "df_all_early=df_all_early[subset_index]\n",
    "df_all_late=df_all_late[subset_index]\n",
    "\n",
    "df_all_early['region']='same'\n",
    "df_all_late['region']='same'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "9b1e1190",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running for period: early\n",
      "Working on region: same\n",
      "Working on global total\n",
      "Running for period: late\n",
      "Working on region: same\n",
      "Working on global total\n"
     ]
    }
   ],
   "source": [
    "#ERROR ANALYSIS from ROMAIN - but only for subset\n",
    "for period in ['early','late']:\n",
    "    if period == 'early':\n",
    "        period_date = '2000_2010'\n",
    "    else:\n",
    "        period_date = '2010_2018'\n",
    "\n",
    "    print('Running for period: '+period)\n",
    "\n",
    "    # Parameters estimated for ghe variograms (nb_model x (correlation range, partial sill)) for velocity and ice thickness:\n",
    "    # ranges in meters, partial sills are in percentage of the variance\n",
    "    params_velocity =((691 ,63), (18000, 26),\n",
    "                      (1040000, 11))\n",
    "    params_icethickness = ((2100, 28), (47000, 50),\n",
    "                           (1354000, 22))\n",
    "\n",
    "    def sum_vgm(h, params, error1, error2):\n",
    "        \"\"\"Sum of variograms models\"\"\"\n",
    "        v = 0\n",
    "        for k in range(len(params)):\n",
    "            v += skgstat.models.spherical(h, params[k][0], params[k][1] * error1 * error2)\n",
    "\n",
    "        return v\n",
    "\n",
    "    def product_vgms(h, mean_vn_err, mean_Hn_err, mean_d, mean_vn, mean_Hn):\n",
    "        \"\"\"Quadratic product of independent velocity and ice thickness variograms\"\"\"\n",
    "\n",
    "        v = (sum_vgm(h, params_velocity, mean_Hn_err, mean_Hn_err) * mean_d ** 2 * mean_vn ** 2 + \\\n",
    "             sum_vgm(h, params_icethickness, mean_vn_err, mean_vn_err) * mean_d ** 2 * mean_Hn ** 2) /\\\n",
    "            ((mean_Hn_err**2 * mean_vn**2 * mean_d**2 + mean_vn_err**2 * mean_Hn**2 * mean_d**2) * 100)\n",
    "        return v\n",
    "\n",
    "    # Function to integrate the variogram from a spatial support to another: a double sum of covariances\n",
    "    def double_sum_covar(coords: np.ndarray, areas: np.ndarray, errors: np.ndarray,\n",
    "                         vgm_func: Callable):\n",
    "        \"\"\"\n",
    "        Double sum of covariance for euclidean coordinates\n",
    "        :param coords: Spatial support (typically, pixel) coordinates\n",
    "        :param areas:  Area of supports\n",
    "        :param errors: Standard errors of supports\n",
    "        :param vgm_func: Variogram function\n",
    "        :return:\n",
    "        \"\"\"\n",
    "\n",
    "        n = len(coords)\n",
    "        pds = pdist(coords)\n",
    "        var = 0\n",
    "        for i in range(n):\n",
    "            for j in range(n):\n",
    "\n",
    "                # For index calculation of the pairwise distance, see https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.pdist.html\n",
    "                if i == j:\n",
    "                    d = 0\n",
    "                elif i < j:\n",
    "                    ind = n * i + j - ((i + 2) * (i + 1)) // 2\n",
    "                    d = pds[ind]\n",
    "                else:\n",
    "                    ind = n * j + i - ((j + 2) * (j + 1)) // 2\n",
    "                    d = pds[ind]\n",
    "\n",
    "                var += errors[i] * errors[j] * (1 - vgm_func(d)) * areas[i] * areas[j]\n",
    "\n",
    "        total_area = sum(areas)\n",
    "        se_dsc = np.sqrt(var / total_area ** 2)\n",
    "\n",
    "        return se_dsc\n",
    "\n",
    "    # 2/ Propagation glacier -> region\n",
    "    list_df_reg = []\n",
    "    for region in ['same']:\n",
    "\n",
    "        print('Working on region: '+region)\n",
    "\n",
    "        # Isolate region of interest\n",
    "        if period=='early':\n",
    "            df_reg = df_all_early[df_all_early.region == region]\n",
    "        if period=='late':\n",
    "            df_reg = df_all_late[df_all_late.region == region]\n",
    "        # Sum volume discharges\n",
    "        V_tot_reg = np.sum(df_reg.V_tot.values)\n",
    "\n",
    "        # Uncertainty on the sum, considering everything independent\n",
    "        V_tot_err_uncorr_reg = np.sqrt(np.nansum(df_reg.V_tot_err_uncorr.values ** 2))\n",
    "\n",
    "        # Uncertainty on the sum, considering everything correlated\n",
    "        V_tot_err_fullycorr_reg = np.nansum(df_reg.V_tot_err_fullycorr)\n",
    "\n",
    "        # Get the regional weighted mean of velocities, thicknesses, flux gate width, and related uncertainties\n",
    "        mean_vn_reg = np.nansum(df_reg.mean_vn.values * df_reg.nb_flux_gates.values)/np.nansum(df_reg.nb_flux_gates.values)\n",
    "        mean_Hn_reg = np.nansum(df_reg.mean_Hn.values * df_reg.nb_flux_gates.values)/np.nansum(df_reg.nb_flux_gates.values)\n",
    "        mean_d_reg = np.nansum(df_reg.mean_d.values * df_reg.nb_flux_gates.values)/np.nansum(df_reg.nb_flux_gates.values)\n",
    "        mean_vn_err_reg = np.nansum(df_reg.mean_vn_err.values * df_reg.nb_flux_gates.values)/np.nansum(df_reg.nb_flux_gates.values)\n",
    "        mean_Hn_err_reg = np.nansum(df_reg.mean_Hn.values * df_reg.nb_flux_gates.values)/np.nansum(df_reg.nb_flux_gates.values)\n",
    "        mean_d_err_reg = np.nansum(df_reg.mean_d_err.values * df_reg.nb_flux_gates.values)/np.nansum(df_reg.nb_flux_gates.values)\n",
    "\n",
    "        # Define the average ice discharge correlation function in the region accordingly\n",
    "        def vgm_func_gla(h):\n",
    "            v = product_vgms(h, mean_vn_err_reg, mean_Hn_err_reg, mean_d_reg, mean_vn_reg, mean_Hn_reg)\n",
    "            return v\n",
    "\n",
    "        # Uncertainty on the sum, accounting for spatial correlations in velocities and ice thickness\n",
    "        coords = np.array((df_reg['mean_easting'].values, df_reg['mean_northing'].values)).T\n",
    "        V_tot_err_corr_reg = double_sum_covar(coords=coords, areas=np.ones(len(coords)), errors=df_reg.V_tot_err_corr.values,\n",
    "                                          vgm_func=vgm_func_gla) * len(coords)\n",
    "\n",
    "        # Store results in a dataframe\n",
    "        df_res_reg = pd.DataFrame()\n",
    "        df_res_reg = df_res_reg.assign(region=[region], V_tot=[V_tot_reg],\n",
    "                               V_tot_err_uncorr=[V_tot_err_uncorr_reg], V_tot_err_corr=[V_tot_err_corr_reg],\n",
    "                               V_tot_err_fullycorr=[V_tot_err_fullycorr_reg],\n",
    "                               nb_flux_gates=[np.sum(df_reg.nb_flux_gates.values)],\n",
    "                               mean_vn=[mean_vn_reg], mean_Hn=[mean_Hn_reg], mean_d=[mean_d_reg],\n",
    "                               mean_vn_err=[mean_vn_err_reg], mean_Hn_err=[mean_Hn_err_reg], mean_d_err=[mean_d_err_reg])\n",
    "        list_df_reg.append(df_res_reg)\n",
    "\n",
    "    df_all_reg = pd.concat(list_df_reg)\n",
    "    if period=='early':\n",
    "        df_all_reg.to_pickle('df_all_reg_early_subset_'+date_string+'.pkl')\n",
    "    if period=='late':\n",
    "        df_all_reg.to_pickle('df_all_reg_late_subset_'+date_string+'.pkl')\n",
    "\n",
    "    # 3/ Propagation region -> global\n",
    "\n",
    "    print('Working on global total')\n",
    "\n",
    "    # Sum volume discharges\n",
    "    V_tot_global = np.sum(df_all_reg.V_tot.values)\n",
    "\n",
    "    # Uncertainty on the sum, considering everything independent\n",
    "    V_tot_err_uncorr_global = np.sqrt(np.nansum(df_all_reg.V_tot_err_uncorr.values ** 2))\n",
    "\n",
    "    # Uncertainty on the sum, accounting for spatial correlations in velocities and ice thickness\n",
    "    # For the global step, no correlations after >1000 km, so we can simply sum quadratically\n",
    "    V_tot_err_corr_global = np.sqrt(np.nansum(df_all_reg.V_tot_err_corr.values ** 2))\n",
    "\n",
    "    # Uncertainty on the sum, considering everything correlated\n",
    "    V_tot_err_fullycorr_global = np.nansum(df_all_reg.V_tot_err_fullycorr)\n",
    "\n",
    "    # Get the global weighted mean of velocities, thicknesses, flux gate width, and related uncertainties\n",
    "    mean_vn_global = np.nansum(df_all_reg.mean_vn.values * df_all_reg.nb_flux_gates.values) / np.nansum(df_all_reg.nb_flux_gates.values)\n",
    "    mean_Hn_global = np.nansum(df_all_reg.mean_Hn.values * df_all_reg.nb_flux_gates.values) / np.nansum(df_all_reg.nb_flux_gates.values)\n",
    "    mean_d_global = np.nansum(df_all_reg.mean_d.values * df_all_reg.nb_flux_gates.values) / np.nansum(df_all_reg.nb_flux_gates.values)\n",
    "    mean_vn_err_global = np.nansum(df_all_reg.mean_vn_err.values * df_all_reg.nb_flux_gates.values) / np.nansum(\n",
    "        df_all_reg.nb_flux_gates.values)\n",
    "    mean_Hn_err_global = np.nansum(df_all_reg.mean_Hn.values * df_all_reg.nb_flux_gates.values) / np.nansum(\n",
    "        df_all_reg.nb_flux_gates.values)\n",
    "    mean_d_err_global = np.nansum(df_all_reg.mean_d_err.values * df_all_reg.nb_flux_gates.values) / np.nansum(\n",
    "        df_all_reg.nb_flux_gates.values)\n",
    "\n",
    "    df_global = pd.DataFrame()\n",
    "    df_global = df_global.assign(region=['Global total'], V_tot=[V_tot_global],\n",
    "                               V_tot_err_uncorr=[V_tot_err_uncorr_global], V_tot_err_corr=[V_tot_err_corr_global],\n",
    "                               V_tot_err_fullycorr=[V_tot_err_fullycorr_global],\n",
    "                               nb_flux_gates=[np.sum(df_all_reg.nb_flux_gates.values)],\n",
    "                               mean_vn=[mean_vn_global], mean_Hn=[mean_Hn_global], mean_d=[mean_d_global],\n",
    "                               mean_vn_err=[mean_vn_err_global], mean_Hn_err=[mean_Hn_err_global], mean_d_err=[mean_d_err_global])\n",
    "\n",
    "    if period=='early':\n",
    "        df_reg_global_early = pd.concat([df_all_reg, df_global])\n",
    "    if period=='late':\n",
    "        df_reg_global_late = pd.concat([df_all_reg, df_global])\n",
    "    #df_reg_global.to_csv('/home/atom/ongoing/work_frontal_ablation/df_region_global_final_'+period+'.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "6473bab7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Running for period: early\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/willkochtitzky/opt/miniconda3/envs/gis/lib/python3.7/site-packages/geopandas/geodataframe.py:1322: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  super(GeoDataFrame, self).__setitem__(key, value)\n",
      "/Users/willkochtitzky/opt/miniconda3/envs/gis/lib/python3.7/site-packages/xdem/spatialstats.py:884: RuntimeWarning: invalid value encountered in true_divide\n",
      "  std_err = np.sqrt(2*np.pi*full_int / area)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Working on region: same\n",
      "Running for period: late\n",
      "Working on region: same\n",
      "Running for period: total\n",
      "Working on region: same\n"
     ]
    }
   ],
   "source": [
    "#define subset from Arctic_discharge_fluxes_with_area\n",
    "df_Subset = Arctic_discharge_fluxes_with_area[subset_index]\n",
    "df_Subset['region']='same'\n",
    "\n",
    "#calculate the spatial correlation of the terminus change volume error from Romain\n",
    "# Define the NMAD as a robust measure of dispersion\n",
    "def nmad(data: np.ndarray, nfact: float = 1.4826) -> float:\n",
    "    \"\"\"\n",
    "    Calculate the normalized median absolute deviation (NMAD) of an array.\n",
    "\n",
    "    :param data: input data\n",
    "    :param nfact: normalization factor for the data; default is 1.4826\n",
    "\n",
    "    :returns nmad: (normalized) median absolute deviation of data.\n",
    "    \"\"\"\n",
    "    data_arr = np.asarray(data)\n",
    "    return nfact * np.nanmedian(np.abs(data_arr - np.nanmedian(data_arr)))\n",
    "\n",
    "for period in ['early','late','total']:\n",
    "    print('Running for period: '+period)\n",
    "    # 1/ Derive uncertainty for pixel -> glacier based on the area of averaging (lost or gained)\n",
    "    # Get variables of interest: ice thickness\n",
    "\n",
    "    if period == 'early':\n",
    "        dt = df_Subset['early_to_middle_years_elapsed'].values\n",
    "    else:\n",
    "        dt = df_Subset['middle_to_late_years_elapsed'].values\n",
    "\n",
    "    # Lost\n",
    "    meanHl = df_Subset['thickness_for_area_lost_' + period].values\n",
    "    meanHl_err_fullycorr = df_Subset['thickness_for_area_lost_' + period + '_err'].values\n",
    "    dSldt = df_Subset['area_km2_lost_' + period].values / dt\n",
    "    \n",
    "    # Remove aberrantly large uncertainty values for ice thickness uncertainty\n",
    "    meanHl_err_fullycorr[np.abs(meanHl_err_fullycorr - np.nanmedian(meanHl_err_fullycorr)) > 5 * nmad(meanHl_err_fullycorr)] = np.nan\n",
    "    # Replace no data that might exist in the uncertainty from the average uncertainty\n",
    "    meanHl_err_fullycorr[~np.isfinite(meanHl_err_fullycorr)] = np.nanmean(meanHl_err_fullycorr)\n",
    "\n",
    "    # Gained\n",
    "    meanHg = df_Subset['thickness_for_area_gained_' + period].values\n",
    "    meanHg_err_fullycorr = df_Subset['thickness_for_area_gained_' + period + '_err'].values\n",
    "    dSgdt = df_Subset['area_km2_gained_' + period].values / dt\n",
    "\n",
    "    # Remove aberrantly large uncertainty values for ice thickness uncertainty\n",
    "    meanHg_err_fullycorr[np.abs(meanHg_err_fullycorr - np.nanmedian(meanHg_err_fullycorr)) > 5 * nmad(meanHg_err_fullycorr)] = np.nan\n",
    "    # Replace no data that might exist in the uncertainty from the average uncertainty\n",
    "    meanHg_err_fullycorr[~np.isfinite(meanHg_err_fullycorr)] = np.nanmean(meanHg_err_fullycorr)\n",
    "\n",
    "    # Parameters estimated for the ice thickness variogram (correlation range, partial sill):\n",
    "    # range in meters, partial sill is in percentage of the variance\n",
    "    params_icethickness = ((2100, 28), (47000, 50),\n",
    "                           (1354000, 22))\n",
    "    \n",
    "    # Looping for every glacier\n",
    "    unique_rgiid = list(np.unique(df_Subset['RGIId'].values))\n",
    "\n",
    "    list_Hl_err_corr, list_Hg_err_corr = ([] for i in range(2))\n",
    "    for rgiid in unique_rgiid:\n",
    "\n",
    "        # Select only a given RGIId = glacier\n",
    "        index_glacier = (df_Subset['RGIId'] == rgiid)\n",
    "\n",
    "#         print('Working on glacier: '+rgiid + ' with '+ str(np.count_nonzero(index_glacier))+ ' flux gates...')\n",
    "\n",
    "        # Variograms parameters formatted for xdem\n",
    "        list_vgm=[(params_icethickness[0][0], 'Sph', params_icethickness[0][1]),\n",
    "                  (params_icethickness[1][0], 'Sph', params_icethickness[1][1]),\n",
    "                  (params_icethickness[2][0], 'Sph', params_icethickness[2][1])]\n",
    "\n",
    "        # a/ Terminus lost\n",
    "        # Number of effective samples in the area from the variogram\n",
    "        dSl = dSldt[index_glacier] * dt[index_glacier]\n",
    "        neff_l = xdem.spatialstats.neff_circ(area=dSl*10**6, list_vgm=list_vgm)\n",
    "\n",
    "        # Uncertainty value for the average ice thickness of the glacier considering spatial correlation\n",
    "        if dSl != 0:\n",
    "            meanHl_err_corr = meanHl_err_fullycorr[index_glacier][0] / np.sqrt(neff_l)[0]\n",
    "        else:\n",
    "            meanHl_err_corr = 0\n",
    "\n",
    "#         print('For terminus lost:')\n",
    "#         print('Replacing fully correlated averaged ice thickness error {:.1f} m into spatially correlated {:.1f} m'.format(meanHl_err_fullycorr[index_glacier][0], meanHl_err_corr))\n",
    "        list_Hl_err_corr.append(meanHl_err_corr)\n",
    "\n",
    "        # b/ Terminus gained\n",
    "        # Number of effective samples in the area from the variogram\n",
    "        dSg = dSgdt[index_glacier] * dt[index_glacier]\n",
    "        neff_g = xdem.spatialstats.neff_circ(area=dSg*10**6, list_vgm=list_vgm)\n",
    "\n",
    "        # Uncertainty value for the average ice thickness of the glacier considering spatial correlation\n",
    "        if dSg != 0:\n",
    "            meanHg_err_corr = meanHg_err_fullycorr[index_glacier][0] / np.sqrt(neff_g)[0]\n",
    "        else:\n",
    "            meanHg_err_corr = 0\n",
    "\n",
    "#         print('For terminus gained:')\n",
    "#         print('Replacing fully correlated averaged ice thickness error {:.1f} m into spatially correlated {:.1f} m'.format(meanHg_err_fullycorr[index_glacier][0],meanHg_err_corr))\n",
    "        list_Hg_err_corr.append(meanHg_err_corr)\n",
    "\n",
    "    # Load under the names previously used\n",
    "    meanHl_err = df_Subset['thickness_for_area_lost_' + period + '_err'].values\n",
    "    meanHg_err = df_Subset['thickness_for_area_gained_' + period + '_err'].values\n",
    "\n",
    "    # 2/ Propagation glacier -> region\n",
    "\n",
    "    list_regions = list(np.unique(df_Subset.region.values))\n",
    "\n",
    "    list_df_reg = []\n",
    "    for region in list_regions:\n",
    "\n",
    "        print('Working on region: '+region)\n",
    "\n",
    "        # Isolate region of interest\n",
    "        ind_reg = df_Subset.region.values == region\n",
    "\n",
    "        # Terminus ice lost\n",
    "        Tl_reg = np.sum(meanHl[ind_reg] * dSldt[ind_reg]) * 0.9 / 1000\n",
    "        mean_meanHl_reg = np.nanmean(meanHl[ind_reg])\n",
    "        mean_dSldt_reg = np.nanmean(dSldt[ind_reg])\n",
    "        mean_meanHl_err_reg = np.nanmean(meanHl_err[ind_reg])\n",
    "\n",
    "        # Terminus ice gained\n",
    "        Tg_reg = np.sum(meanHg_err[ind_reg] * dSgdt[ind_reg]) * 0.9 / 1000\n",
    "        mean_meanHg_reg = np.nanmean(meanHg[ind_reg])\n",
    "        mean_dSgdt_reg = np.nanmean(dSgdt[ind_reg])\n",
    "        mean_meanHg_err_reg = np.nanmean(meanHg_err[ind_reg])\n",
    "\n",
    "        # Separating two sources of terminus ice loss: one based on area uncertainty, one based on ice thickness uncertainty\n",
    "        # Here we only compute only the ice thickness source (second term of Equation 6 in the paper)\n",
    "\n",
    "        # a/ Terminus mass lost\n",
    "        # Considering all ice thicknesses uncertainties fully independent\n",
    "        Tl_errH_uncorr_reg = np.sqrt(np.nansum((meanHl_err[ind_reg] * dSldt[ind_reg]) ** 2)) * 0.9 / 1000\n",
    "\n",
    "        # Uncertainty on the sum, considering everything correlated\n",
    "        Tl_errH_fulluycorr_reg = np.nansum(meanHl_err[ind_reg] * dSldt[ind_reg]) * 0.9 / 1000\n",
    "\n",
    "        # Uncertainty on the sum, accounting for spatial correlations in velocities and ice thickness\n",
    "        coords = np.array((df_Subset['easting'].values[ind_reg], df_Subset['northing'].values[ind_reg])).T\n",
    "        Tl_errH_corr_reg = double_sum_covar(coords=coords, areas=dSldt[ind_reg], errors=meanHl_err[ind_reg],\n",
    "                                          vgm_func=vgm_func_gla) * sum(dSldt[ind_reg]) * 0.9 / 1000\n",
    "\n",
    "        # b/ Terminus mass gained\n",
    "        # Considering all ice thicknesses uncertainties fully independent\n",
    "        Tg_errH_uncorr_reg = np.sqrt(np.nansum((meanHg_err[ind_reg] * dSgdt[ind_reg]) ** 2)) * 0.9 / 1000\n",
    "\n",
    "        # Uncertainty on the sum, considering everything correlated\n",
    "        Tg_errH_fulluycorr_reg = np.nansum(meanHg_err[ind_reg] * dSgdt[ind_reg]) * 0.9 / 1000\n",
    "\n",
    "        # Uncertainty on the sum, accounting for spatial correlations in velocities and ice thickness\n",
    "        coords = np.array((df_Subset['easting'].values[ind_reg], df_Subset['northing'].values[ind_reg])).T\n",
    "        Tg_errH_corr_reg = double_sum_covar(coords=coords, areas=dSgdt[ind_reg], errors=meanHg_err[ind_reg],\n",
    "                                          vgm_func=vgm_func_gla) * sum(dSgdt[ind_reg]) * 0.9 / 1000\n",
    "\n",
    "        # Store results in a dataframe\n",
    "        df_res_reg = pd.DataFrame()\n",
    "        df_res_reg = df_res_reg.assign(region=[region], Tl=[Tl_reg], Tg=[Tg_reg],\n",
    "                               Tl_errH_uncorr=[Tl_errH_uncorr_reg], Tl_errH_fullycor=[Tl_errH_fulluycorr_reg],\n",
    "                               Tl_errH_corr=[Tl_errH_corr_reg],\n",
    "                               Tg_errH_uncorr=[Tg_errH_uncorr_reg], Tg_errH_fullycor=[Tg_errH_fulluycorr_reg],\n",
    "                               Tg_errH_corr=[Tg_errH_corr_reg],\n",
    "                               mean_Hl=[mean_meanHl_reg], mean_meanHl_err=[mean_meanHl_err_reg] ,mean_dSl=[mean_dSldt_reg],\n",
    "                               mean_Hg=[mean_meanHg_reg], mean_meanHg_err=[mean_meanHg_err_reg] ,mean_dSg=[mean_dSgdt_reg])\n",
    "        list_df_reg.append(df_res_reg)\n",
    "\n",
    "    if period=='early':\n",
    "        df_all_reg_early = pd.concat(list_df_reg)\n",
    "        df_all_reg_early.fillna(0).to_pickle('df_all_reg_early_TC_subset_'+date_string+'.pkl')\n",
    "    if period=='late':\n",
    "        df_all_reg_late = pd.concat(list_df_reg)\n",
    "        df_all_reg_late.fillna(0).to_pickle('df_all_reg_late_TC_subset_'+date_string+'.pkl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "047127dc",
   "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>region</th>\n",
       "      <th>FA_early</th>\n",
       "      <th>FA_early_err</th>\n",
       "      <th>FA_late</th>\n",
       "      <th>FA_late_err</th>\n",
       "      <th>D_early</th>\n",
       "      <th>D_early_err</th>\n",
       "      <th>D_late</th>\n",
       "      <th>D_late_err</th>\n",
       "      <th>TC_early</th>\n",
       "      <th>TC_early_err</th>\n",
       "      <th>TC_late</th>\n",
       "      <th>TC_late_err</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>total</td>\n",
       "      <td>4.840986</td>\n",
       "      <td>1.876154</td>\n",
       "      <td>4.294005</td>\n",
       "      <td>1.737846</td>\n",
       "      <td>1.485013</td>\n",
       "      <td>1.020344</td>\n",
       "      <td>0.663295</td>\n",
       "      <td>0.425918</td>\n",
       "      <td>-3.355973</td>\n",
       "      <td>1.574437</td>\n",
       "      <td>-3.63071</td>\n",
       "      <td>1.684845</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  region  FA_early  FA_early_err   FA_late  FA_late_err   D_early  \\\n",
       "0  total  4.840986      1.876154  4.294005     1.737846  1.485013   \n",
       "\n",
       "   D_early_err    D_late  D_late_err  TC_early  TC_early_err  TC_late  \\\n",
       "0     1.020344  0.663295    0.425918 -3.355973      1.574437 -3.63071   \n",
       "\n",
       "   TC_late_err  \n",
       "0     1.684845  "
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Read in the data you just wrote out above for each subset\n",
    "df_early_D = pd.read_pickle('df_all_reg_early_subset_'+date_string+'.pkl')\n",
    "df_late_D = pd.read_pickle('df_all_reg_late_subset_'+date_string+'.pkl')\n",
    "\n",
    "df_early_T = pd.read_pickle('df_all_reg_early_TC_subset_'+date_string+'.pkl')\n",
    "df_late_T = pd.read_pickle('df_all_reg_late_TC_subset_'+date_string+'.pkl')\n",
    "\n",
    "regions = df_early_T['region'].unique()\n",
    "for period in ['early','late']:\n",
    "    list_df_reg = []\n",
    "    list_df_uncertain = []\n",
    "    for region in regions:\n",
    "#         print(region)\n",
    "        if period == 'early':\n",
    "            df_T = df_early_T[df_early_T['region']==region]\n",
    "            df_D = df_early_D[df_early_D['region']==region]\n",
    "            time='2000_to_2010'\n",
    "            time2='2000_2010'\n",
    "            time3 = 'early_to_middle_years_elapsed'\n",
    "        if period == 'late':\n",
    "            df_T = df_late_T[df_late_T['region']==region]\n",
    "            df_D = df_late_D[df_late_D['region']==region]\n",
    "            time='2010_to_2020'\n",
    "            time2='2010_2020'\n",
    "            time3 = 'middle_to_late_years_elapsed'\n",
    "            \n",
    "        df_FA = Arctic_discharge_fluxes_with_area[subset_index]\n",
    "        \n",
    "        Tg_errA = np.sum(df_FA['perimeter_error_gained_'+period]/df_FA[time3]*df_FA['thickness_for_area_gained_'+period]*0.9/1000)\n",
    "        Tl_errA = np.sum(df_FA['perimeter_error_lost_'+period]/df_FA[time3]*df_FA['thickness_for_area_lost_'+period]*0.9/1000)\n",
    "\n",
    "        sig_ID_raw_reg = df_D['V_tot_err_corr']\n",
    "        sig_TC_raw_reg = np.sqrt(df_T['Tg_errH_corr']**2 + df_T['Tl_errH_corr']**2 + Tg_errA**2 + Tl_errA**2)\n",
    "        sig_ID_clim_reg = np.sqrt(sum(df_FA['melt_correct_gt_'+period+'_err_discharge']**2))\n",
    "        sig_TC_clim_reg = np.sqrt(sum(df_FA['melt_correct_gt_'+period+'_err_terminus_change']**2))\n",
    "\n",
    "        sig_FA_reg = np.sqrt(sig_ID_raw_reg**2 + sig_ID_clim_reg**2 + sig_TC_raw_reg**2 + sig_TC_clim_reg**2) *2 #convert to 2 sigma error\n",
    "        sig_FA_reg=sig_FA_reg[0]\n",
    "#         print('--',sig_FA_reg)\n",
    "        \n",
    "        sig_D = np.sqrt(sig_ID_raw_reg**2 + sig_ID_clim_reg**2) *2\n",
    "        sig_TC = np.sqrt(+ sig_TC_raw_reg**2 + sig_TC_clim_reg**2) *2\n",
    "        \n",
    "        df_reg = pd.DataFrame()\n",
    "        df_reg = df_reg.assign(region=[region],FA=[sum(df_FA['Frontal_ablation_'+time+'_gt_per_yr_mean'])],FA_err=[sig_FA_reg],\n",
    "                               Discharge=[sum(df_FA[time2+'_mean_flux_gt'])],D_err=[sig_D[0]],\n",
    "                               Terminus_change=[sum(df_FA['terminus_gt_change_per_year_'+period])],TC_err=[sig_TC[0]]\n",
    "                              )\n",
    "        list_df_reg.append(df_reg)\n",
    "        \n",
    "        df_uncertain = pd.DataFrame()\n",
    "        df_uncertain = df_uncertain.assign(region=[region],\n",
    "                                           D_raw=[sig_ID_raw_reg[0]**2],\n",
    "                                           TC_raw=[sig_TC_raw_reg[0]**2],\n",
    "                                           D_clim=[sig_ID_clim_reg**2],\n",
    "                                           TC_clim=[sig_TC_clim_reg**2]\n",
    "                                          )\n",
    "        list_df_uncertain.append(df_uncertain)\n",
    "        \n",
    "    if period=='early':\n",
    "        df_reg_early = pd.concat(list_df_reg)\n",
    "        df_uncertain_early = pd.concat(list_df_uncertain)\n",
    "    if period=='late':\n",
    "        df_reg_late = pd.concat(list_df_reg)\n",
    "        df_uncertain_late = pd.concat(list_df_uncertain)\n",
    "        \n",
    "df_global = pd.DataFrame()\n",
    "df_global['region'] = 'total'\n",
    "df_global = df_global.assign(FA_early=[sum(df_reg_early['FA'])],\n",
    "                             FA_early_err=[np.sqrt(sum(df_reg_early['FA_err']**2))],\n",
    "                             FA_late=[sum(df_reg_late['FA'])],\n",
    "                             FA_late_err=[np.sqrt(sum(df_reg_late['FA_err']**2))],\n",
    "                             \n",
    "                             D_early=[sum(df_reg_early['Discharge'])],\n",
    "                             D_early_err=[np.sqrt(sum(df_reg_early['D_err']**2))],\n",
    "                             D_late=[sum(df_reg_late['Discharge'])],\n",
    "                             D_late_err=[np.sqrt(sum(df_reg_late['D_err']**2))],\n",
    "                             \n",
    "                             TC_early=[sum(df_reg_early['Terminus_change'])],\n",
    "                             TC_early_err=[np.sqrt(sum(df_reg_early['TC_err']**2))],\n",
    "                             TC_late=[sum(df_reg_late['Terminus_change'])],\n",
    "                             TC_late_err=[np.sqrt(sum(df_reg_late['TC_err']**2))],\n",
    "                             )\n",
    "df_global['region'] = 'total'\n",
    "df_global"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "0eaafd3f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/Users/willkochtitzky/bin/Jupyter_notebooks/Arctic/Arctic_discharge'"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pwd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "043baa48",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3f4a65e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81e0f4e8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90e30903",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2213254f",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}