RGI_13.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inversion with a corrected shapefile for Nordic Glacier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-06-04 13:46:05: oggm.cfg: Using configuration file: /home/pelto/anaconda2/envs/oggm_pub/lib/python3.7/site-packages/oggm/params.cfg\n",
      "2020-06-04 13:46:05: oggm.cfg: Multiprocessing switched ON according to the parameter file.\n",
      "2020-06-04 13:46:05: oggm.cfg: Multiprocessing: using all available processors (N=8)\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import xarray as xr\n",
    "\n",
    "import oggm\n",
    "import os\n",
    "from oggm import cfg, utils, workflow, tasks, graphics\n",
    "from oggm.core import inversion\n",
    "cfg.initialize(logging_level='WORKFLOW')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's read a file from the standard RGI:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/home/pelto/OGGM/rgi/RGIV61'"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "utils.get_rgi_dir(version='61')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'/home/pelto/OGGM/rgi/RGIV61/13_rgi61_CentralAsia/13_rgi61_CentralAsia.shp'"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sh = utils.get_rgi_region_file('13', version='61')\n",
    "sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Shapefiles are best read an manipulated with [geopandas](http://geopandas.org/) in python (see also oggm [working_with_rgi.ipynb](working_with_rgi.ipynb) tutorial):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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>RGIId</th>\n",
       "      <th>GLIMSId</th>\n",
       "      <th>BgnDate</th>\n",
       "      <th>EndDate</th>\n",
       "      <th>CenLon</th>\n",
       "      <th>CenLat</th>\n",
       "      <th>O1Region</th>\n",
       "      <th>O2Region</th>\n",
       "      <th>Area</th>\n",
       "      <th>Zmin</th>\n",
       "      <th>...</th>\n",
       "      <th>Lmax</th>\n",
       "      <th>Status</th>\n",
       "      <th>Connect</th>\n",
       "      <th>Form</th>\n",
       "      <th>TermType</th>\n",
       "      <th>Surging</th>\n",
       "      <th>Linkages</th>\n",
       "      <th>Name</th>\n",
       "      <th>check_geom</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>RGI60-13.00001</td>\n",
       "      <td>G079494E34471N</td>\n",
       "      <td>20000810</td>\n",
       "      <td>-9999999</td>\n",
       "      <td>79.4939</td>\n",
       "      <td>34.4713</td>\n",
       "      <td>13</td>\n",
       "      <td>5</td>\n",
       "      <td>0.432</td>\n",
       "      <td>5678</td>\n",
       "      <td>...</td>\n",
       "      <td>683</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>9</td>\n",
       "      <td>9</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>POLYGON ((79.49945 34.47555, 79.49880 34.47498...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>RGI60-13.00002</td>\n",
       "      <td>G077951E35545N</td>\n",
       "      <td>20020802</td>\n",
       "      <td>-9999999</td>\n",
       "      <td>77.9513</td>\n",
       "      <td>35.5452</td>\n",
       "      <td>13</td>\n",
       "      <td>5</td>\n",
       "      <td>0.367</td>\n",
       "      <td>5889</td>\n",
       "      <td>...</td>\n",
       "      <td>1307</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>9</td>\n",
       "      <td>9</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>POLYGON ((77.95519 35.55007, 77.95517 35.54953...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>RGI60-13.00003</td>\n",
       "      <td>G077930E35519N</td>\n",
       "      <td>20020802</td>\n",
       "      <td>-9999999</td>\n",
       "      <td>77.9295</td>\n",
       "      <td>35.5188</td>\n",
       "      <td>13</td>\n",
       "      <td>5</td>\n",
       "      <td>0.070</td>\n",
       "      <td>5646</td>\n",
       "      <td>...</td>\n",
       "      <td>426</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>9</td>\n",
       "      <td>9</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>POLYGON ((77.92863 35.52057, 77.92864 35.52030...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>RGI60-13.00004</td>\n",
       "      <td>G077924E35525N</td>\n",
       "      <td>20020802</td>\n",
       "      <td>-9999999</td>\n",
       "      <td>77.9237</td>\n",
       "      <td>35.5252</td>\n",
       "      <td>13</td>\n",
       "      <td>5</td>\n",
       "      <td>0.255</td>\n",
       "      <td>5698</td>\n",
       "      <td>...</td>\n",
       "      <td>605</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>9</td>\n",
       "      <td>9</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>POLYGON ((77.91713 35.52720, 77.91764 35.52719...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>RGI60-13.00005</td>\n",
       "      <td>G077914E35531N</td>\n",
       "      <td>20020802</td>\n",
       "      <td>-9999999</td>\n",
       "      <td>77.9141</td>\n",
       "      <td>35.5309</td>\n",
       "      <td>13</td>\n",
       "      <td>5</td>\n",
       "      <td>0.261</td>\n",
       "      <td>5644</td>\n",
       "      <td>...</td>\n",
       "      <td>389</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>9</td>\n",
       "      <td>9</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>POLYGON ((77.91939 35.52980, 77.91938 35.52958...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 24 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "            RGIId         GLIMSId   BgnDate   EndDate   CenLon   CenLat  \\\n",
       "0  RGI60-13.00001  G079494E34471N  20000810  -9999999  79.4939  34.4713   \n",
       "1  RGI60-13.00002  G077951E35545N  20020802  -9999999  77.9513  35.5452   \n",
       "2  RGI60-13.00003  G077930E35519N  20020802  -9999999  77.9295  35.5188   \n",
       "3  RGI60-13.00004  G077924E35525N  20020802  -9999999  77.9237  35.5252   \n",
       "4  RGI60-13.00005  G077914E35531N  20020802  -9999999  77.9141  35.5309   \n",
       "\n",
       "  O1Region O2Region   Area  Zmin  ...  Lmax  Status  Connect  Form  TermType  \\\n",
       "0       13        5  0.432  5678  ...   683       0        0     0         0   \n",
       "1       13        5  0.367  5889  ...  1307       0        0     0         0   \n",
       "2       13        5  0.070  5646  ...   426       0        0     0         0   \n",
       "3       13        5  0.255  5698  ...   605       0        0     0         0   \n",
       "4       13        5  0.261  5644  ...   389       0        0     0         0   \n",
       "\n",
       "   Surging  Linkages  Name  check_geom  \\\n",
       "0        9         9  None        None   \n",
       "1        9         9  None        None   \n",
       "2        9         9  None        None   \n",
       "3        9         9  None        None   \n",
       "4        9         9  None        None   \n",
       "\n",
       "                                            geometry  \n",
       "0  POLYGON ((79.49945 34.47555, 79.49880 34.47498...  \n",
       "1  POLYGON ((77.95519 35.55007, 77.95517 35.54953...  \n",
       "2  POLYGON ((77.92863 35.52057, 77.92864 35.52030...  \n",
       "3  POLYGON ((77.91713 35.52720, 77.91764 35.52719...  \n",
       "4  POLYGON ((77.91939 35.52980, 77.91938 35.52958...  \n",
       "\n",
       "[5 rows x 24 columns]"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gdf = gpd.read_file(sh)\n",
    "gdf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get the RGI geometry and plot it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>RGIId</th>\n",
       "      <th>GLIMSId</th>\n",
       "      <th>BgnDate</th>\n",
       "      <th>EndDate</th>\n",
       "      <th>CenLon</th>\n",
       "      <th>CenLat</th>\n",
       "      <th>O1Region</th>\n",
       "      <th>O2Region</th>\n",
       "      <th>Area</th>\n",
       "      <th>Zmin</th>\n",
       "      <th>...</th>\n",
       "      <th>Aspect</th>\n",
       "      <th>Lmax</th>\n",
       "      <th>Status</th>\n",
       "      <th>Connect</th>\n",
       "      <th>Form</th>\n",
       "      <th>TermType</th>\n",
       "      <th>Surging</th>\n",
       "      <th>Linkages</th>\n",
       "      <th>Name</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>52591</th>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>G080036E34291N</td>\n",
       "      <td>20110805</td>\n",
       "      <td>-9999999</td>\n",
       "      <td>80.036</td>\n",
       "      <td>34.291</td>\n",
       "      <td>13</td>\n",
       "      <td>8</td>\n",
       "      <td>23.097</td>\n",
       "      <td>5435</td>\n",
       "      <td>...</td>\n",
       "      <td>134</td>\n",
       "      <td>11760</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>9</td>\n",
       "      <td>9</td>\n",
       "      <td>CN5Z421I0054 Zhanjin Glacier</td>\n",
       "      <td>POLYGON ((80.02679 34.32982, 80.02684 34.32980...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1 rows × 23 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                RGIId         GLIMSId   BgnDate   EndDate  CenLon  CenLat  \\\n",
       "52591  RGI60-13.52594  G080036E34291N  20110805  -9999999  80.036  34.291   \n",
       "\n",
       "      O1Region O2Region    Area  Zmin  ...  Aspect   Lmax  Status  Connect  \\\n",
       "52591       13        8  23.097  5435  ...     134  11760       0        0   \n",
       "\n",
       "       Form  TermType  Surging  Linkages                          Name  \\\n",
       "52591     1         0        9         9  CN5Z421I0054 Zhanjin Glacier   \n",
       "\n",
       "                                                geometry  \n",
       "52591  POLYGON ((80.02679 34.32982, 80.02684 34.32980...  \n",
       "\n",
       "[1 rows x 23 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gl='RGI60-13.52594'\n",
    "rgidf = utils.get_rgi_glacier_entities([gl])\n",
    "rgidf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPkAAAD4CAYAAAAuE0A1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3RU1drH8e8zJT2BAAESCAQFpCTUSBOkShOw6xVRUATrvag0EUS9dlFQFAtSFS8KgjQBEaWK1FBCb9IJnfQ6s98/Mtw3NyZkSGZyksn+rJVFJnP23r8Bnpw5Z87ZW5RSaJrmuUxGB9A0zb10kWuah9NFrmkeThe5pnk4XeSa5uEsRge4EZUqVVIRERFGx9C0Emfbtm0XlVIheT1Xqoo8IiKCrVu3Gh1D00ocETme33P67bqmeThd5Jrm4XSRa5qH00WuaR5OF7mmeThd5Jrm4XSRa5qH00WuaR5OF7lGamoqS5cuZfiIkdzRvSerV69GzzPgOXSRl1EZGRnMnz+f7r3uokJIZR4b8gpT/zzNposW7n54AA0aNWXBggW62D2AlKZ/xOjoaKUvay26efPmMXDw01gqVEfVvh2/Oi0x+5X77/NK2Uk9tImMTbNpUr82P3z3DVWqVDEwsVYQEdmmlIrO67lSde26VnTvfTCOt97/iMA7R+Eddkue24iY8KvbGt+bo4n983saN7+Vtb+vpG7dusWcVnMFXeRlyJtvvc24z76i/EPvYwmqVOD2YrYS2PZRkoMq065DJw7s3U358uWLIanmSvqYvAxQSjF0+Eg+/HwKQfe95VSB5+TfqBu2qg0YP+FjNyXU3EkXuYdLSUmh76P9mfrDQoLufxtLYMVC9ePd7F4+nfQ5drvdxQk1d9NF7sF27txJg0ZNWbH7DEH3vfk/J9dulLVSOHaTlUOHDrkwoVYcdJF7oPPnzzN0+Aja3N6RpHq9COz+IiZvvyL3q+x2fHx8XJBQK066yD3IyZMneeb5f1Grdl1mrtlP8CMTCIjs5JK+lbKTEn9Jf5RWCumz6x4gMTGRocNHMus//8EvsgvBj04s9LF3blkJF0k59Cdmv/IEliuv9+SlkC7yUm7z5s30ufd+sqo0oOLjX2L2DXRJv1mJl0hbP4PUv2Lo2bMnP/7wFbc0jHJJ31rx0kVeis2e/T2Dnn4Wv87PElC3tcv6Td23hqTVX/PcM88wdswiAgIC+H3wk/j6+rpsDK346CIvpb79dhbPvjCMoPvewKvyTS7tW/21kcFPDGD0qJEEBAQA0KmTa47tteKnT7yVQrGxsTzzrxcIvHusywscwBTZk28XrKBS5SoEBVeiToNGDBg4iCNHjrh8LM399J68FBr/yad4N+qBV0hNt/TvUyMKanyAv7JjT44nKfECiw9t5cfmt7Jx/VoiIyPdMq7mHnpPXspcvnyZOXPm4BvZxe1jiZgwBwTjHVqXwDZ98WrQhe/+M9vt42qupYu8lHnvg3H41m2NJfDGrj93Ba+IZsyZr+8xL210kZciiYmJTPr8C3xufdCQ8b2rN+BSShaff/GFIeNrhVNgkYuIj4hsFpGdIrJHRN7I9fwwEVEi8rddy/Xaisg4EdkvIrtE5CcR0fcwFmD58uX4VauHpVxlQ8YXkxn/bi8xYtQYduzYYUgG7cY5sydPBzoppRoDTYDuItIKQETCgTuAEzfaFvgViFRKNQIOAqMK/zLKhrXr13Pp0DbsGWmGZbBWrI5Pq4d5YdhIwzJoN6bAIlfZkhwPrY6vawdlE4AROR473VYptUIpleV4biNQvVCvoAx5/913ua3d7STvXmloDv/IzmzZupU1a9YYmkNzjlPH5CJiFpEdwHngV6XUJhHpA5xWSu280bZ5bPYEsCyf9oNFZKuIbL1w4YIzcT2Wn58f4959m9SN32NLiTcsh8nqTUD3l+h9972sWLHCsByac5wqcqWUTSnVhOy9bQsRaQSMBsYWou3/fMgqIqOBLOC7fNpPVkpFK6WiQ0LyXGO9TAkICCA18SpZV+MMzeEb0QTfri9wf78nuKNHL/bt22doHi1/N3R2XSl1FVgN3AXUAnaKyDGyCzhGRKo60bb7tZ+JSH+gF/CI0p/LOKVBgwZUCQsHk9noKPjWakbwY5+xPT2E6NZtqVW3PqNeGc2xY8eMjqbl4MzZ9ZBrZ75FxBfoAmxXSlVWSkUopSKAU0AzpVScE233Ox53B0YCfZRSKS58TR5t69atxCfEYw0OMzoKAGKxEtDifioNnk5qi4F89ft+Ihs35aPxH2Oz2YyOp+HcnjwUWCUiu4AtZB9XL8lvYxEJE5GlTrT9DAgEfhWRHSLyZaFfhRulpaVx6tQpEhISDL0IRCnFt9/OonPX7vh3GETqmikk7lyBPTPdsEw5icmMT/X6BHV4gnIPvc/bk6bTvGUbDh48aHS0Mk8vrpBLfHw8y5cv5/sf57N27VoSrlzGO6AcmalJ2LKyqNsgil7d76BH9260b98ek8n91xNdunSJxwc9xZqNMfh1HQK2TOJmDad9565s2bqNwF4j8Kne0O05boRSdpK3/4x9+0+sW/07DRuWrHye5nqLK5TpIldKcezYMbZu3cqadX+wet0fHD64j8CakdhrROMT0RRLUAjiOP61Z6aTceYAGSd3wokY/M12hjz/DH0ffpiwMPe8fV60aBEDnhyM6aY2+N3WD5PVG1taEhcmDyQlKYEVK1Zw/z/6EnTP63hVudktGYoiee9qUtdMYdhLL/LsM09TubIxF/J4Ol3kuUz8bBIzZ81m/55YTF4++IbWJrNSbbzD6uFVtS4mr4KnOFJKkX5qD1n7V5F88E+iGjVm8OOP8cgjj7hkcoWNGzfywrCR7Dv8F76dns2+MyyHc5/15fSJv6hQoQLTp8/gpbHvUO6RCUUe1x0yr8aRvvkHkg9uZNTIEQwfNlRPQOFiushzuemWBlwJbUlAozuKNE3xNSorg9QjW7EfWEXWucO89OK/+OdzzxEcHOx8H0px+PBhli1bxoxZ33PoyFG8bn2AgKguiPnvdwRf+OJRjh7cx+nTp3niqWc5fj6Bcg+PK/JrcaeshAukrZ5MxtkDTP7ic/7xj4eMjuQxdJHnkJGRQfWaN2HpMdItEy5kXDhORsxPJB/cyO0dOtK7R1cqVqxIUFAQNpuN9PR0MjIySEhI4PTpMxw/dYaDR46we8d2zN5+eNWIwnRTa3wjmiIWa77jxH/zHH+sXMqgZ//JXsIJbNYLk1fp2Dumxx0mfv5rxG7fxs03l7xDjNJIF3kOI15+hamLVhPQ+xVExEXJ/s6elkTKwQ3IxaOYMpJRGSmI2QwmC5gs2C0+ZHqXw+wfjCUoBK+qtTH7O7/nT1oxka71Q5g3/yeqPD3DJfOqF6ektTN4KLoakybqpZdcQRe5w4kTJ7ilQRQVB3yGJaCCC5MVv6yEC6Rt+wkJuRn/yM5Gx7lhGReOwa8fcubEX0ZH8Qh66WKHxYsXE1CnRakvcABLUAgBHQcbHaPQrJVqcCbuDMnJyfj7+xsdx6OVmUkj/vrrL157821MddsbHUUDsNlQSuHnV7oOM0ojjy/yzMxMlixZwm3tOyKN78a3VlOjI2mAsmViMht//X1Z4LFFnpGRwZdffkXl0Or0HzKKjMb349/0TqNjaQ72zDRA9CqpxcDjjsmPHDnClGnTmDJ1BvZyYfj0HIl32C1Gx9Ic0uMOkxmzgJSj23ikXz+qVs33xkXNRTyuyGvXro1f6M2U6zGyRF7mWVbZ0pJIXTcd+4ntvPLySAYPWkhQUJDRscoEjyvyES+PYurKWF3gxciWEk/y7t8wXTiE2G3YqzXCWqkGypaVPYONspO+6Qf+cV8fPlo1p8jFferUKTZt2kTlypWpUqUKVapUISgoyK3XPZRmHlfkg58cyOeT2xDQ+Wmjo5QJmZdOkrjwTXp26UjXfgM4euwYq9f9SdLhHXh7eVO5cgiXLl1i1LQv6d27t0vGnPjZZ4yfMJHg8NpkJl0hNeEyymbj5lsa0KNrZ7p26Uzbtm3/u45bWedxRV6rVi3SkxOwZ6SWmss8S7PkZR/x79EjeGHIkGIZLzU1lY0bN1Ouw+P4N+0JQHnAnpHGlbhDzNi0i28Xvkz8yYMEV6xEWLXq1IqoSd2ba9G4URRNmzalTp06mPM5s5+ZmYnVmv/lxKWRR17x1qFLN3abahHQzDV7Di1v9vQUTn/al4z0tHyLxpWysrLo0fsuth2/QkC3FzB5538RjcrKJCvxAlkJF7AlXMSWcA5r/CnS446QlniZiJvrcvPNN2Gz2cnISOfEiZOcPnGM9NRkxrw6ljf//Ua+fZdEZe6Kt88+/ohWbdvjFR6FV0iE0XE8lnj5EFChCjt37qRZs2ZuHUspxeCnn2XbkXME3jUGMV9/bysWK9bgsL9Nk+VH9i+nhIsn2BJ/DixmxMuMuVV3KnevSvLe1Rw5dtyNr6T4eeTn5JGRkXw16VMS5r9GxrmjRsfxWCImrBHNWLt2LQA///wzN9drSNeevTl58qRLx3r3/Q+Yv3wVAXeOLLDAC2Ly9sO7Wj38G7THv15b/Oq2xrtqbUw+AYjFC7vd7qLUJYNH7skBHnmkL0rZeW74GKyPTEAsXkZH8kgS1pDX33qHL6fM4HTcOfw6PsXGbT+xdOlSnnrqqTzbJCYmsn37dhISEqhYsSJ79uzBZDLRuHFjmjVr9j9nyZVSxMbG8uZbb1HhkQnuv9tOTNjtpecQ1hkeW+QA/fr1Y/bc+WxaN5OAjoOMjuOR/Ou1JaNiOFdTE6jQox5itpK1dQ716tXLc/uYmBiaN29OxVoNwGQFeyZSPgwRM2knX6VqpQrMmjGVFi1aADD42eeZMXUqwV0GYw0OLYZXJHpPXtp8N3MakY2bkbRzOf6NuunPUt3AK6Tmf7/PvBpH5tXztGqVveTdgQMHmD9/Pj4+PiQmJTH+k0+pdOeLed4e66vsXNm3jg6duxLdoiUto5uwdPFCApv3xq9Rt2J5LXLhEJEtbi2WsYqLxxd5+fLlWbFsCT16303ioXVYm92Ld7X6pW6ShdIi/dQe2rVvj7e3NwDTZ8xk4uyf8ascgTJb8b1zFN6hdfNsK2LCv0F7vMMbsv/MQXb9eQKp0w3/em2LLX/W6T3cc/dbxTZecfD4IofsVUeOHNjLzJnf8MnnX3Fg8XuUa9QFv455HzNqhWetUJ1NK77nzJkzhIWF4e3thXe1+gS0e9TpPiyBlbDcUgk/2rgx6d8ppUi6eIbw8PBiHdfdPPLsel4sFgsDBz7Brm2buHj+HJazu0g5lNfai1pReIfdgjTsTqOm0cTFxVGzRg2sqZeNjuUUEaH8zU345ZdfjI7iUmWmyHMKDAxk9rczSVn1JbbURKPjeBz/W++D8CZMmz4jez761CtGR3KaNOjK8JdHk5SUVPDGpUSZLHKA9u3b81jff5C8fDzKllVwA+2GmOu0Zfo3s2jSpAkJJw+g7KVjXTS/2i1Jq1ibV8a8ZnQUlymzRQ4w8ePxNK1ZkZQ/vjU6isfxqRHFhcQ01q9fT5WqoWReOGZ0JKd5Ne7NgsWLjY7hMmW6yC0WCzOnfU1S7K+OmUo0VxEx4dNhMAMHP43YbWQlXDA6ktOsITW5fOUqGzZsMDqKS5TpIgeoXr06kVGNSTu+0+goHsenekO8Wj/K8b+OQCm6PkFMZrxbPszAp541OopLlPkiB7inVw/Sdy7DllR6ThCVFv5RXQgbPBnfCPfewOJKKau+JHHVZGrVqmV0FJfQRQ48//xz9O3amksznydp848o5VmXNRrNGhx23SWfSprM0/v4ad6PLF30k9FRXEIXORAUFMQXkyayf88uqsXvIXHR22TFnzc6lmYQa+VaXL5cOj7bd4Yu8hxq1KjB5g3reO7B7lz57iUS132LLVm/hS9LlC2LtJN7qF+/vtFRXEYXeS7e3t68+cbr7I3dQe96gVya8RxJv0wg5fBmVFam0fE0N0vet5Zbat9EdHSek6yUSrrI81GzZk1mTP2aE38d4dUBval+ciUXJj9O0u9fkX72IKVp2izNOSk7l5G2bhoff/i+0VFcqsA53kTEB1gLeJN9Q8uPSqnXcjw/DBgHhCilLjrbVkQqAD8AEcAx4EGl1HXfG7ti6eKiOHHiBFOnTWfK9JkkpqThG1obZfUjKygUn5tbYK0Yrm9lLaWStv+MLWYef6xZRYMGDYyOc8OKtHSxZP+v9VdKJYmIFVgPDFFKbRSRcGAKUA9onkeRX6/tB8BlpdR7IvIyEKyUGnm9LEYX+TVKKXbv3s3BgweJj49nw8bNLFq8hKSkJHwrVMHsVx6VkYItPRlLcDXEZMJu8cW7xYNYylUxOr7moOw2Ug9vxrZ7Kb7pV1i1cgV16tQxOlahFGkiR5X9W+Da1fpWx9e13wwTgBHAwkK0vQvo4Ph+JrAauG6RlxQiQlRUFFFRUQA88cQTKKU4f/48p06d4vz585QrVw5fX18OHTqEiBC7ew8ffTwUn4Zd8G7QCa9KNQx+FWWXPS2J5J3Lydy9nFo1azDy9WHcf//9eHl55hRhTk3JLCJmYBtQG5iklBopIn2AzkqpISJyDIjOvSfPr63j51eVUuVzbHdFKRWcR/vBwGCAGjVqND9+vPTOpHn48GEmffkVM2d+gwRURNVqjW/tllgrVDM6WpmRcngTKb9/Sc/uXXl52Etun2W2uBTp7XqujsoDPwFDgK+Brkqp+OsVeR5t/6mU2u1skedUUt6uF1VWVhYrV65k9pwfWbx4Mcrqh7lOO3yb9MTsG2h0PI9kS00kZc0UfC4fYfasmdx+++1GR3Kp6xX5DZ1dV0pdJftt9V1ALWCno8CrAzEiku8SlTnadnf86JyIhDoChgJl5uoTi8VC9+7dmTltChfPnWXFwjncUUO4NO0pklZPLVU3c5QGqcd3cuXbf/FAm3oc3Lfb4wq8IAUWuYiEOPbCiIgv0AXYrpSqrJSKUEpFAKeAZkqpOCfa7nc8vQjo7/i+P/kc13s6k8nErbfeyuxvZ3Jw3x4ebhXBlVkvkPTz+6Qe3aY/qisCpRTJMUtIXzGBhXNn8+Xnn+Hvn/+qK57KmT15KLBKRHYBW4BflVJL8ttYRMJEZKkTbd8D7hCRQ8AdjsdlWvXq1fn04/HEnT7Jm8/2JXDXDyR8P4zUo9uMjlYqJW+aS+Cx1cRs2UTnzn+fHbas8Mi10DyF3W5n3rx5vDj8ZVIDwvBtPwhLYEWjY5UKybtWYN/6A7E7YrKnoPJwLjsm14qXyWTigQce4PD+PTzZ+3auzh5G2qm9Rscq0ewZaSSumYbP3iVsWLemTBR4QcrElMylnY+PD++/+zbt293Gg337obq9hG+tpkbHKhGULYuM83+RfnIXlri9JBzfQ6fOXfhm8SYqVapkdLwSQRd5KdKzZ0+WLV5Iz953kRV9H/6Nu2Oy+hgdq1ilnz1E5p5f4PJJslLiSbl6ibDwGvTu0oluzw2nffv2urhz0UVeyrRr144N69Yw4pUxrJ36FEG9X8a7mufcFnk99vRk4hf8mzfGjKJ169ZUqVKF0NBQAgICjI5WoukTb6XYokWL6Pv4IMo/9AGWIM/feyVuX0pr37MsWTDf6Cgljj7x5qH69OnD6BHDOD/taY6/34vEDbOxpycbHcttzHF7+Mf99xodo9TRRV7KjRo5nDOnTjJo0CC61jBx5Zt/kXLYQ5d/MpnJytILYdwo/Xbdw6xatYp+/Z8gLTgCnxYPYa3oOYv3Je5cQXM5yoqlnrPwgavot+tlSMeOHTm0fw+PdGrGhVnDSs3yRM7wr9eWdWvXsm/fPqOjlCq6yD2Qn58fEyd8hI+vD7aUeKPjuIzJ24+A9o/TrWdvrl69anScUkMXuQer3yCSjHNHjI7hUv5Rd5AS0pAHHu6nb95xki5yD9bq1mZkXii9k2zkx7/942zZc5hJn39udJRSQRe5B6sSUglJ97z118Vsxb/bC4x8ZYw+PneCLnIPFhUVhfnyX0bHcAtrxXD82vanQ5euHD582Og4JZoucg/WpUsXbJdPkXZyt9FR3MI/6g5sje6mZZu27N2r787Ljy5yD+br68s7b75ByqqvPOqjtJwCGndHbv0HHbp09Zj1xF1NF7mHGzRoEPUiqpG8drrHno0OiLoDWg2ga88+vPPue8THe87Hhq6gi9zDmc1mfvl5ESHJf5H0+1fYM9KMjuQWfnVbU+7Bd5jw/QpCq4VzzwP/YM2aNUbHKhF0kZcBwcHB/LHmd9rX8ufKN8977LXt1grV8O8xlIpPfMmaq+Xpde+DrFy50uhYhtP3k5cRFSpUYP6c7/n99995dMBAkg6tx7fdE5gDsqe6z0q6TMre1ZgTzyFZqdj8KuLTsEupXPjB7FeOoOa9iU9P4udly+jSpYvRkQylb1Apg1JSUnjl1bF8PWUqXg06Y0k4S/KJ3dx77720bd0Sf39/duyKZfLXU/BpPwj/+qVznvL0Mwewr5zAyWNHPHYJpGtctoKK0XSRu9bRo0f5aMInNGvSiIceeuhvM6zExsbSsk07Kg6YhNm/fD69lGyJS97lrtYN+OKzT/H29jY6jtvoItcKrXqt2tg7vYhXSITRUQrFlppIyoqPqeqVwdzZs0rlssTO0LeaaoWyfv16riYkYg0ufcfl15h9AwnoM5pzVVsT3aoNAwYO4sKFsrUMld6Ta3lSSlEvsjGXb+6Of4P2RsdxCVtqAsl/zCJcXSBm85+YzWajI7mM3pNrN+zYsWOcPnMWv1J60i0vZt8gAjs/zckkxf0PPUx6errRkYqFLnItT3/++Sd+4Q0QEaOjuJSIiYBeo1i77yzde91VJuaM00Wu5SkxMRHl5ZkrgJqs3gTcOZxdp+Lpfc99pKV55lWA1+gi1/Jks9nA5DnHrLmJyUzAnSPZcjKJm+vWZ/PmzUZHchtd5FqeLBYLeOida9eIxUpA95fIbN6XDnd0p36jZowc9QqHDh0yOppL6SLX8uTv74/JVjZOTPnd0oZKT07hStRDTFlzkKbRLdi+fbvRsVxGF7mWp1tuuYW0MwdQym50lGJhsnrjEx5J4O2P49W4F1OnzzQ6ksvoItfy1Lx5c6qGVCTtaIzRUYqdtVp9fl+91ugYLqOLXMuTiDB21AgyNv3HY2eVyY9P9YYcPrjPYyaf0Fe8aflSStHm9o7sz6yEpXItbBf+wpKegM0rAEuNJvjUbIyYPfNu5aTl43mkUxM+Gf+R0VGcUqQr3kTER0Q2i8hOEdkjIm/ken6YiCgR+dvauSISLiKrRGSfo+2QHM81EZGNIrJDRLaKSIvCvDjNfUSEubNnUVviaJK1n2G9m/HBv/oy5M5mhBxcxOVpT5G0dYFHTivl2/ZxpkydRmxsrNFRiqzAPblkX/Lkr5RKEhErsB4YopTaKCLhwBSgHtBcKXUxV9tQIFQpFSMigcA24G6l1F4RWQFMUEotE5GewAilVIfrZdF78pJlx44dPPbEII4eO453wy74Ne7530koPEHKvjVk/jGDdat/JzIy0ug411WkPbnKluR4aHV8XfvNMAEYkeNx7rZnlVIxju8TgX3AtVuaFBDk+L4ccKbgl6KVJE2aNGHH1k2sXbmcexoEcWnGsyTNG0PKAc+YNdWvfnsstw2gXYdO7N+/3+g4hebUMbmImMneC9cGJimlRopIH6CzUmqIiBwDonPvyXP1EQGsBSKVUgkiUh/4BRCyf9m0UUr9bU0fERkMDAaoUaNG8+PHPW/ZH0+RkpLCypUreWLw05jbP4NvrWZGR3KJpNhf8YpdwNaNGwgLCzM6Tp6KfBeaUsqmlGoCVAdaiEgjYDQw1skAAcA84AWlVILjx88ALyqlwoEXgan5jD1ZKRWtlIoOCQlxZjjNIHFxcUz8/CsS4+PBg47TA6LuIL1WO7r3uiv7ct9S5oY+QlNKXQVWA3cBtYCdjr14dSBGRKrmbuM4jp8HfKeUmp/jqf7AtcdzAX3irRS7cOEC7TvdQUxaJSoP+hrfm5obHcmlAlo9xKmELP791ttGR7lhzpxdDxGR8o7vfYEuwHalVGWlVIRSKgI4BTRTSsXlaitk76H3KaXG5+r6DHBtNoJOgGddMFyGpKWl0bVnb1KrRRPU6kHMfuWMjuRyIoJftxf5cPzHLFmyxOg4N8SZPXkosEpEdgFbgF+VUvm+ShEJE5Gljoe3AY8CnRwfle1wnEkHGAR8JCI7gXdwHHdrpYtSin79H+dEmjf+bfsZHcetLIEVCerzCn0f7c+BAweMjuM0fTGMViRjX3+DidO/J+j+tzBZfYyOUyySd/2C954lrFv9G7Vq1TI6DqCnf9Lc5KvJk5kw6SsCer9SZgocwC+qK2l1OtOiTdtS8dGaLnKtUBYsWMDQkaMJuud1LAEVjI5TrESEgOi7odmDNG/ZmpeGDiMmpuTeyKOLXLthv/32G/0GDCSwz+hSuYySq/hHdaH8wx8y488TdOjWi/seetjoSHnSRa457ciRIzz19DP0vud+AnoOwzu0jtGRDGctX5Wg2/sT3H8SS5cu5ezZs0ZH+hvPvIVIc7mzZ8/Ssk1bbKGRVOj3MZagv92PVKaJxQvfcpW4dOkSoaGhRsf5H3pPrhXIbrfT7c4+UP8OynV/QRd4Pix+5Th37pzRMf5GF7lWoKSkJA7s24N/yweNjlKiqcCQEjkJpC5yrUCZmZmYrV4et9CCy4U344uvp5W4++t1kWsFMplM2ErwSiP2tCSSVk3m8pQnSfrpddJP7zMkh1/d1pw4f5X58+cXvHEx0kWuFah8+fKYRLClJRW8cTFTtkzifxxDzwaV2LD6V94b/jQZKz4i6ffJ2DOLd0ppMZnxajuAgYOf4fTp08U69vXoItcKJCKEhdck60rJm9cjedOPNG9Qm5nTphAZGcmTTw7k8P69tK5mIXHe2GL/xeRbszHm2q15570PinXc69FFrjklunlTMs4dMTrGfym7jaSNc7DtXcG306f8z/mC4OBgFs6by8O9O5M4v/gL3a/FA8z45huOHCkZf1+6yDWntGvdEtPFo0bHACDj3FES5rxMXfsJdm3fRrVqf7/qzmQyMWnixzx6T08S55bjbtcAABLsSURBVL+GvRgL3ewfjG/jnjz17D+x241fnEIXueaUhg0bQryxb9ez4s+T+NuXJC98g3FjhvLHmt8JDw/Pd3sRYeLH43m4d1cSfnode3pysWX1b/kA2/YeZvHixcU2Zn50kWtOqVu3LqkXTxk2fvLe1Vz9z1D6d6jPof17efLJgU59pCcifP7ZJ9zX7XaSV3xSDEkd45qtWKMfYOTosYZ/pKaLXHNKaGgo6ckJqKzMYh877dReMv+YwZ/r1zB+3AdUrlz5htqLCJ9/OhG/lHOkHNninpB58K3binPxKYbPJKOLXHOKyWTC1z8Qe0ZKsY5rS0siefl4vp0xjaioqEL34+3tzcTx48j487ti+0UlYsJ664MMHzXG0L25LnLNadVr1CQj7nCxjpm6bjoP3tOH3r17F7mvu+++m9ZNGpC88XsXJHOOb53WnItPZdmyZcU2Zm66yDWnvfDcM2TFLi14QxfJuHAM+8kdTPhoXNH6ychg7ty5DBz8NBv/3IAkxBXcyEVEBFNkd8ZPnFRsY+ami1xz2qOP9sM/7SLJe9cUy3hZCReo3yCSoKCggjfO5fLlyyxfvpwRI1+marVwnhn9LguOKnzuHIVftxfdkDZ//vVvZ8Mf6zlzxphPJ/T95JrTfH19WTBvDu06dsGram3S//yOhENbCLylNQHdhiAms0vHs1aoxsb5q2nToTMVK1Qgqv4thIVWxWw2czYujmMnT2Mxm7FaraSnp3Py9FnOxsVxPu4siQlXCap+C7aQOvj0GYtXSIRLs90Ik5cvfvXa8fXUqbz26qvFPr6erVW7Ye998AGjRo4kutVtLJw3h0bNovG6cwxeITVdPlbmlTNkXT6NLS0J25UzWDISQSmyvAIx+QcDCmW3ISYzZv8KmPyDMQcEYwkKcfkvnaJIjzuM7ZdxnDl5DIvF9fvW683Wqvfk2g0bPnQoNapXp1evXgQFBVG5SiiX091zRZk1OAxrsHvXH0s7sYusv7bi366/234xeFetTaJfBRYuXMh9993nljHyo/fkWpGsX7+edu3aUan3cLISL2JKuYQJBUohweH4Nb3T6IjXlXnpFFd/HINkpeMT1Y2A2we4bazk/eupdnoVO7ZsdHnfet51zS2UUrRr1w6AWxJj6BsZyKsPd+T1x7rx2mPdSN7wHRnn/zI4Zf4yL58mfv5YPvnwPY4eOYz5+Ea3LrvsV7c1R44eY+/evW4bIy/67bpWaCJCSkoKvr6+eT6fkprOu59/ifXBdxEpWfuTzMuniZ/3Kh9/8A5PDhwIwMJ5c+nc/U68wxu6ZT03MZnxrtOGr6ZM5ZPxH7m8//yUrL95rdTJr8ABhg19kSBTBumninfPVZDMq3EkzB/7PwUO0KpVKx7v/xipG75129jejXswffpMEhMT3TZGbrrINbcxmUw89/Rg0jfPQdlKxvRRWQkXSJw/lndef/V/Cvyaf7/2KhlHtpB56aRbxrcGh+FVvSEzZ850S/950UWuudVLL75A85uqkLzqS6OjkJVwgYT5Y3ll2Av88/nn8tymQoUKjBg+lLS17puQ0VyvE599NcUtfedFF7nmVhaLhYXz5+J/5TDJMcbdW5155Szxc0fz6rAhvDxi+HW3HTl8GFV9skjZ/rNbsvjUasqp06fZtWuXW/rPTRe55nb+/v6s/X0latcikvesKvbxMy6eIOHHMbz/7zEMHza0wO2tVis/zf2e1M0/kHnV9de5i8mMd1QP3h1XPCffdJFrxaJmzZqs+e1X1ObvSNkyr9jGTTu5m4QfxzBpwjiefeYZp9vVrVuXEUNfIm2Ne95W+zXpwaJFi4tlHjhd5FqxadiwIbt3bse+exlJKyZiS77qtrFsKfEkrp1B2vIPmTt7Fo899ugN9zHq5ZF4JZ0l7eRul+cz+wbh27AzEz751OV956aLXCtWVatWZW/sTh5sXYeEuaNIiv3NZfOjZyVeJGnXryT//D4Xpz5Fz7qB7I3dSY8ePQrVn9VqZewrL5OxZa5bTsJZa7dm+a+/ubzf3PRlrZohlFIsW7aMd8eNJyYmBv+IRmRVrI13tXp4Va2NmK1O9WNLvkLKgQ3IkXWkXTxFh46duP/u3tx9992UL1++yDkzMzOpUz+S1OjH8I1oUuT+crJnpnP+i8c4e/okwcHBReqrSJe1ioiPiGwWkZ0iskdE3sj1/DARUSLyt6UuRSRcRFaJyD5H2yG5nv+niBxwPFdyZqPX3E5E6NmzJ+tWrWTvru2MH/4k99xsxn/rTOI+60vSso9IP3swz7b29GSS960ladFbXJr+LO2CLjJz4ntcuXieJQvmMWDAAJcUOGTvzV/457PYDrj+HnqT1ZvAGvVZv369y/vOqcA9uWRPiemvlEoSESuwHhiilNooIuHAFKAe0FwpdTFX21AgVCkVIyKBwDbgbqXUXhHpCIwG7lRKpYtIZaXU+etl0XvysuHy5ctMmz6D1//9Jt6NeuLXpAcmv/JkxB0iM2Y+SUe306JVGwY+9ggPPPAAfn5+bs1z5swZ6tRrSPm+H2ItX9WlfSf8+QOPNK3Ipx9PKFI/RdqTq2zX7iO0Or6u/WaYAIzI8Th327NKqRjH94nAPuDaTPjPAO8ppdIdz1+3wLWyo0KFCgwb+hKxO2K4s7YPl2Y8x5VpT2FfOYHRA+/nQtxZ1v62gv79+7u9wAHCwsJ4fewYUn6Z4PIr96yVb2LtevfdFANOHpOLiJnsvXBtYJJSaqSI9AE6K6WGiMgxIDr3njxXHxHAWiBSKZUgIjuAhUB3IA0YppT623y5IjIYGAxQo0aN5sePH7+xV6iVehcvXuTUqVM0btzYsOWT7XY7Hbp0Y29WZQLa9HVdvxmpxH3+KBlpqUXqp8i3miqlbEqpJkB1oIWINCL7rfZYJwMEAPOAF5RSCY4fW4BgoBUwHJgjefwLKqUmK6WilVLRISEhzgyneZhKlSrRpEkTQ9dHN5lMfD9rJqk7lpKVcMFl/YrVB1tmJunp7luB9YY+QlNKXQVWA3cBtYCdjr14dSBGRP52wOI4jp8HfKeUyrlw8ylgvuNwYDNgB/528k7TSoqwsDAGPfkEaVtddzGPiBAYEsbRo+5bZ86Zs+shIlLe8b0v0AXYrpSqrJSKUEpFkF2wzZRScbnaCjAV2KeUGp+r6wVAJ8d2dQEvIN+3+5pWEowe9TJpB9aTleC6U0jWqnVYvXq1y/rLzZk9eSiwSkR2AVuAX5VS+a77IiJhInJtcu7bgEeBTiKyw/HV0/HcNOAmEdkNfA/0V6XpQ3utTAoJCeHOO3uRcvBP13Vaoxnfz1vguv5y0RfDaNoNWr9+PXfe9zDlH/vU6Yt2rseensy5Lx8nKeEqVmvh+tNzvGmaC7Vt25ao+nVI2vWrS/ozefvjH1KNLVvcsxijLnJNK4TxH7xHxtZ5LvvcXKrWY82atS7pKzdd5JpWCC1atOCWurVJOeiaC1ksNZsz56dFLukrN13kmlZILz7/DOrgapf05R0eyf69uzl79qxL+stJF7mmFVLv3r1JPX3AJRfHmLx88KvXlmnTZxQ9WO6+Xd6jppURQUFBPPboY6TFuOZttrVeR76YPMXl967rIte0Injt1VfIOLDGJVM4e1erT3Km4o8//nBBsv+ni1zTiqBq1aq8OORfpLvgUlcRQWq3Zca3s1yQ7P/pIte0Ihr60oukHtlCVkLRr8r2a9iZ2f+ZTUpKiguSZdNFrmlFFBQURL9+/UjdWfR52i1BlfCtehMrV650QbJsusg1zQVGDnuJ1NgV2DPTityXvUY08xa47jNzXeSa5gI33XQTLVu2IiW26Je6+tRoxG+/u24RCl3kmuYiH773Nmlb52NPTy5SP9bKEcQnJhMbG+uSXLrINc1Fmjdvzm1tWhV5KSgRE94NOvLZF65ZJFIXuaa50NuvjyV9649FXjDCu35H5syZi91uL3ImXeSa5kLR0dFEN2tG8t7VRerHEhwGPkFs2FD0G2B0kWuai419ZSS2nYuLdHmqiCA1m7Nocb6TMDlNF7mmuVjHjh3xtwrpJ4t24syrZhOWLC/62Xpd5JrmYiLCe2+9Qda2ol3q6lW1DkcO7ictrWifvesi1zQ3uO+++0g8sRd7RuEL1OTlS0BorSLP5KqLXNPcwM/Pj9Zt25FcxItjVHgzlixdXqQ+dJFrmptM/GgcaVvmFulSV6/qDVm9vmi3nuoi1zQ3adSoES1uvZXk2MLfbGIOrETcmTNFyqGLXNPc6NMJH5KxZQ7pp/cVqr0lsCJXLpwr0sdxusg1zY0iIyOZOfVrUn79tHBv200WrD5+XLp0qdAZdJFrmpvdc889tG9zKykb59xwWxEhoEoN9u/fX+jxdZFrmpuJCJ9P/JjU2BVkJd747DEmi5WsrMIv4qCLXNOKQfXq1en/2GOkbb/xy1SV3Y7ZbC702LrINa2YvDjkeVL3/IY9/Ubnb7OTvQp44egi17RiUqdOHXr26EHK1p9uqJ09I52AgIBCj6uLXNOK0Yfvv0PqzqXYkq843SYrPUUXuaaVFjVr1qRv34dJiVns1Pb29BTSEi4THh5e6DF1kWtaMRs5bKhjZteCZ49JP3uQeg2j8Pb2LvR4usg1rZjVrl2bli1bkrxrRYHbZp7eR+cOtxdpPF3kmmaAiePHkbZ5ToErosqpGLp3vaNIY+ki1zQDREVFMXrUSJKWvIeyZea5TcaFY9gSLtKpU6cijVVgkYuIj4hsFpGdIrJHRN7I9fwwEVEiUimPtuEiskpE9jnaDsljm3zba5onGzVyBC0j65C0ZvrfbkBRyk7a2qmMHfMKFoulSOM4sydPBzoppRoDTYDuItIKsosYuAM4kU/bLGCoUqo+0Ap4TkQaXHvSifaa5rFEhP98O4Py8YdI+u1LMq+cxZZ8hfSzh0j6eRx1Kgfwr38+X+RxCixylS3J8dDq+Lr2a2cCMCLH49xtzyqlYhzfJwL7gGo5Nrlue03zdBUrVmTLn+u5q2k10heMJX7WC/hs+IpHu7Vg7e+/FnkvDuBUDyJiBrYBtYFJSqlNItIHOK2U2unMJXciEgE0BTY5HjvVXkQGA4MBatSo4UxcTStVKlSowLSvv4Kvv3JL/04VuVLKBjQRkfLATyLSCBgNdHWmvYgEAPOAF5RSCSLi52x7pdRkYDJAdHS03uNr2g26obPrSqmrwGrgLqAWsFNEjgHVgRgRqZq7jYhYyS7w75RS8x0/vtnZ9pqmFU2Be3IRCQEylVJXRcQX6AK8r5SqnGObY0C0UupirrYCTAX2KaXGX/u5UioWKLC9pmlF58yePBRYJSK7gC3Ar0qpfG+KFZEwEVnqeHgb8CjQSUR2OL56Fjm1pmlOK3BPrpTaRfYJs+ttE5Hj+zNAT8f364ECz8rlbK9pmmvpK940zcPpItc0D6eLXNM8nBRl0vbiJiIXgOMGDF0JKE1n/nVe9yqJeWsqpULyeqJUFblRRGSrUira6BzO0nndq7Tl1W/XNc3D6SLXNA+ni9w5k40OcIN0XvcqVXn1MbmmeTi9J9c0D6eLXNM8XJkrchF50THf3G4Rme2Yw66CiPwqIoccfwbn07a7iBwQkcMi8nKOnzvVvgTlHSci+0Vkl4j85JgnoMTmzfG8y+cDdFdeEfmn47k9IvKBq/IWilKqzHyRPfXUX4Cv4/EcYADwAfCy42cvk30rbe62ZuAIcBPgBewEGjieK7B9CcvbFbA4vn+/pOd1PB8O/EL2xVCVSnJeoCOwEvB2PK5s5P/7MrcnJ/vOO18RsQB+wBmyJ8GY6Xh+JnB3Hu1aAIeVUkeVUhnA9452ONm+xORVSq1QSl1b8Hoj2ZN2lNi8Du6aD9AdeZ8B3lNKpQMopc67OPMNKVNFrpQ6DXxI9uywZ4F4pdQKoIpS6qxjm7PkmNAih2rAyRyPT/H/k1I6074k5c3pCWBZSc6bcz5AV+R0d16gLtBORDaJyBoRudWVuW9UmSpyx7HVtamrwgB/EennbPM8fubWzx/dnVdERpM9bfZ3RcmZoz+X580xH+BYV2T8nwHd9/drAYLJnoZ8ODBHnJnt1E3KVJGTPXXVX0qpC0qpTGA+0AY4JyKhAI4/83p7dYrs48JrqpP91g4n25ekvIhIf6AX8IhyHDiW0LzunA/QXX+/p4D5KttmwE72TS2GKGtFfgJoJSJ+jt+sncmeC34R0N+xTX9gYR5ttwB1RKSWiHgB/3C0w8n2JSaviHQHRgJ9lFIpLsrqlrxKqVilVGWlVITKnkHoFNBMKRVXEvM6nlsAdAIQkbpkn5gz7q41I8/6GfEFvAHsB3YD3wLeQEXgN+CQ488Kjm3DgKU52vYEDpJ9VnV0jp/n2b4E5z1M9vHkDsfXlyU5b67+j+Gis+tu/Pv1AmY5+owhewUiw/7P68taNc3DlbW365pW5ugi1zQPp4tc0zycLnJN83C6yDXNw+ki1zQPp4tc0zzc/wEVHMuBisItTgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# OGGM set-up\n",
    "cfg.PATHS['working_dir'] = utils.gettempdir(dirname='RGI_13_trial', reset=True)\n",
    "\n",
    "# Get the geometry and plot it\n",
    "# gl = utils.get_rgi_glacier_entities(['RGI60-13.00'], version='61') ##unless I put the version here it doesn't call V61\n",
    "rgidf.plot(edgecolor='k');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2020-06-04 13:46:19: oggm.workflow: Execute entity task define_glacier_region on 1 glaciers\n",
      "2020-06-04 13:46:21: oggm.workflow: Execute entity task define_glacier_region on 1 glaciers\n",
      "2020-06-04 13:46:22: oggm.workflow: Execute entity task glacier_masks on 1 glaciers\n",
      "2020-06-04 13:46:23: oggm.workflow: Execute entity task compute_centerlines on 1 glaciers\n",
      "2020-06-04 13:46:23: oggm.workflow: Execute entity task initialize_flowlines on 1 glaciers\n",
      "2020-06-04 13:46:23: oggm.workflow: Execute entity task compute_downstream_line on 1 glaciers\n",
      "2020-06-04 13:46:23: oggm.workflow: Execute entity task compute_downstream_bedshape on 1 glaciers\n",
      "2020-06-04 13:46:23: oggm.workflow: Execute entity task catchment_area on 1 glaciers\n",
      "2020-06-04 13:46:23: oggm.workflow: Execute entity task catchment_intersections on 1 glaciers\n",
      "2020-06-04 13:46:23: oggm.workflow: Execute entity task catchment_width_geom on 1 glaciers\n",
      "2020-06-04 13:46:24: oggm.workflow: Execute entity task catchment_width_correction on 1 glaciers\n",
      "2020-06-04 13:46:24: oggm.workflow: Execute entity task process_cru_data on 1 glaciers\n",
      "2020-06-04 13:46:25: oggm.workflow: Execute entity task local_t_star on 1 glaciers\n",
      "2020-06-04 13:46:25: oggm.workflow: Execute entity task mu_star_calibration on 1 glaciers\n",
      "2020-06-04 13:46:25: oggm.workflow: Execute entity task prepare_for_inversion on 1 glaciers\n",
      "2020-06-04 13:46:25: oggm.workflow: Execute entity task mass_conservation_inversion on 1 glaciers\n",
      "2020-06-04 13:46:25: oggm.workflow: Execute entity task filter_inversion_output on 1 glaciers\n"
     ]
    }
   ],
   "source": [
    "gdirs = workflow.init_glacier_directories(rgidf) #,from_prepro_level=1, prepro_border=10) \n",
    "workflow.execute_entity_task(tasks.define_glacier_region, gdirs)\n",
    "workflow.gis_prepro_tasks(gdirs)\n",
    "workflow.climate_tasks(gdirs)\n",
    "workflow.inversion_tasks(gdirs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['RGI60-13.52594'], dtype=object)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rgidf.RGIId.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from rgitools.funcs import compute_intersects\n",
    "# new_intersects = compute_intersects(gl)\n",
    "# new_intersects\n",
    "divides= utils.get_rgi_intersects_entities([gl])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPkAAAD4CAYAAAAuE0A1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3xT1f/H8dcnabpbSmkLLbS0gGwBoWxURhkiQ/GrIoo4EVAEFcGNoqIgKG5FUMGFMgTEwZ4OVoGyl2wos3SPNDm/Pxr81VogbZOmTc/z8eiDpPeecz8B3rm5N/eeI0opNE1zXwZXF6BpmnPpkGuam9Mh1zQ3p0OuaW5Oh1zT3JyHqwsoipCQEBUdHe3qMjStzNm8efM5pVRoYcvKVcijo6PZtGmTq8vQtDJHRI5cbpn+uK5pbk6HXNPcnA65prk5HXJNc3M65Jrm5nTINc3N6ZBrmpvTIdfy5OTABx/A8uWurkRzMB3yiuzIERg2DCIjUV5eMHw4xMWRW7UqTJ/u6uo0B9Ehr4i++AKio1HR0fDxx1iPH2e/CIt9K7HdYMJw5gw89BDm4GBYtcrV1WolVK4ua9UcYMAA1HffAbBTDEyvXo/vWvbDM6Y5BpMXAFUunGDC3FfpcuE4qlMnsp96Cu9Jk1xZtVYCUp6Gf4qNjVX62vUS6NoVli3jrBjoHNuH5PYDMHj5Xnb1xkcS+P6HF/G1Wki9914CZ8woxWK1ohCRzUqp2MKW6Y/rFYHVCjfcAMuWcdhgoFnPJ0nt/NAVAw6wo2YTrh82gxSjiYCZM0latqyUCtYcSYfc3eXkQIsWsHYtBwwG2vUZjalxR7ubX/AL4qE7xgGQdvfdTipScyYdcne2aRNUrQpbt7LVYOT6fi/iWa9DkbvZGHUtiX6ViThzhszUVCcUqjmTDrk7slrhySdRrVqhLl7kKy8fbrrrTbxqtyx2l3vCamEENvzyi+Pq1EqFPrvuTnJz4bnnUB98gGRmkinC/TUb8+ctL+Dl7V+iritlp6GA7IAAx9SqlRodcneQmwtDh6JmzkRycshF+N7Hj2c7DcbQuDMGkRJvokpGMhagWo0aJa9XK1U65OXd8uXQrx+kpJAtwky/SrzS8WE8Gl6P0WAscfdKWclJPEC1lHOkAJGRkSWvWStVOuTl2WuvoV58EYCvvH14uucovOu0wuSAPbdSVjK2LyMnfj6BPiaUNZe1QN/KlUvct1a6dMjLq/vuQ82YQbYIferEsqfXGHw8vR3Stfn8MTKWfUCNIG8+mfst7dq1Y83YsZwweTqkf6106ZCXR88+i5oxg/Mi3NCuP+ntBzjkuPuSzB1L6dG2CV9/9RVGY95H/hvHjeNGh21BK036stbyJiMDFRBAplI0vX4g5rZ3OHwTWcd3cn7OWPz9A4ipfQ01oyJp3bwZQ1q1JKhzZ3DgG4rmGPqyVnfy8suI1cqHlas6JeAA3jUaUX3kHPxuf5OTtXqxNjOCrdN+ICgujpPDhjllm5rz6JCXJ1Yr6sMPyQWmtL3H6ZvzqFQV75pN8G/chVNdhgCQ/OOPTt+u5lg65OXJ/v2orCxme/ni0fD6Ut30keDqpJq8qXP6NNYLF0p121rJ6JCXJ/XqMeyaejzZ+WHEAd+BF9X77fpjApIaNiz1bWvFd9WQi4i3iGwQkW0islNEXimwfJSIKBEJKUpbEXlLRPaISIKI/CgiQY55Se7r2LFjzDx5Cs/GnV2y/alt/semsBiqnD7NyYcfdkkNWtHZsyfPBjorpZoCzYAeItIGQEQiga7A0aK2BZYCjZVSTYB9wLPFfxkVw9KlS8lMvYg1O8NlNQwYOJlcEQI//zzvRhitzLtqyFWeNNtTk+3n0vdu7wCj8z23u61SaolSKte27C9AXxR9FX379qX7TT3JiF/oshpyPDz57Zq2+FutXGhZ/LvatNJj1zG5iBhFZCtwBliqlFovIn2AE0qpbUVtW8hqDwC/Xqb9YBHZJCKbzp49a0+5bqtKlSrMbNSQG+N/IjfNdSe/Hu87miO+QQTHx5MTFgZr1risFu3q7Aq5UsqilGpG3t62lYg0AZ4HXipG28b5l4vI80Au8M1l2k9VSsUqpWJDQwudY71CCZkzh+8z0vDcttRlNVgNHnR/ZBqrQ2tiOnsWdeON5NaqBS++CPrMe5lTpLPrSqmLwCqgLxADbBORw+QFOF5EqtnRtsel34nIIKAXcLcqT5feuZAhOhoTEOXi70WyPL0Z9MCH3Nx/PAe9/TEeOpR3w0yVKpg7d4YM15030P7NnrProZfOfIuIDxAHbFFKhSmlopVS0cBxoLlSKtGOtntsz3sAY4A+Sin9P8Ie77+PWrWKDBH2RF7r6moA2FWzCXEjZlF/6Je8Wq8DiR6emFauJLdyZdSsWa4uT8O+PXk4sFJEEoCN5B1XL7rcyiISISK/2NH2AyAAWCoiW0Xkk2K/Ckc7eRKefho6dIDoaAgPh7p1oW1beOmlvOWlKS0N2rZFPf44uUCfem0wnz1M5t+bUapsnOHODgxh+i3P0PapeXzYsCOSkwN33UV6v36uLq3C0zeoQF6IZs+GhQth7VrU+fMIeV8D5Np+TIAR/vm9NSAA4yOPwCuvgO+VhzYukTlzUAMHIllZ7BMDt7a9g9SWt3Ds3f7UbXgtieeT8e06HK/qDZxXQzHUOHeU+TOeICQ3m4zWrfFdsgQCA11dltu60g0qFTPk33wDixfDli2ow4chLY1L91WZgQSDkR/Ca/Ntk5uwVKuDMTAUZc7G72Iid8Yv4uZjCTRPT877PlAEa4MGGB99FAYPBg8H3b178CDcfTdq/XoUMCmgMlP+Nw7PsBiUUpz56G4O7d/Ln3/+ycD7HySg93NlLugBWWn88unDRGalYhVBXX89xgULIEhf9+RoOuT5DRqUNxYaeXvkFOCgwcifQVWZE92c3Q1uwCui3lUvGzXmZDLs5ync+fdGqufmIIDVaIRbb8Xw6acQHFy8+lasgFdfRdnmIDssBvo3605il8GI0fTPaumzn2H+zE9o37498+bN475hI6l8z7uIR9kb2GHwko95eNtiQq25WIxGDAMGIF9+CQZ9VbWjXCnkFW/QiKAgBHjZvwqfdhuGMfwaPPz/P5D2jq1i8fTh/Vuf5X2gyulDPLrsY24/sYeAOXNQc+eiWrXC8PTTcOuthf9nvngRduyAXbvg999hwwbU4cNIVhYAZ0QYE9WIpd0fx1Q5goJ3cIvBiNlsJiUlhTW//4E5LRmrOQtjGQz51G5DmdptKI8t/pBhW3/D96uvyFi2DN/SPrdRQVW8PXlUFOrYMW5ocwfHbrzXMYVdYrUyYMVnPLH1V0Isuf98WkAkL+hK/fNTMLQKuACs9/LhvTptSWh9G56hNS+7qcxF45n6+mg2bo7n/Vm/4td+IJ5Vazn29TiBslhYMH0IzZJOcfi994gePtzVJbkF/XH9kieegClTWG704L4hX2D0d96ghNF/b2bY6i+JTjlDYK4Zk9VCrghmEXLFQJKHB4lefhz2DWJ9cCR/1m2PZ/V6GH0r2dV/6rqv6Brlwab4LaS1uBfvyMZXb1RGxJw/zoppQzhSOZjoC+ddXY5b0B/XL/n8c3KA+7o96tSAAxyu1YLRtVrYvb5PEfv3va4Xy1d8DL4R+EbUK2Jr1zpUpQYXvfyomnQBq9WKQR+bO1XFCfny5aiUFDYaPTBcG+fqakrM6FcZv97PubqMYrvgG0RMdjoHDhygbt26ri7HrVWMt9D161E9eqCAide0RvRAhC7nZ87EDPg68xoDDXD3kH/3HcTEoNq0gdxcHguNYkvPJ11dlQZYEYyAj09RD1S0onK/kFut8PHHeWfRBwxAHT7MHjFwf3QTFt09CYPJy9UVVmjhKWf5dN5rVE27gBE48tNPri7J7bnXMfnUqahHHvnnq6u/DEYea9Gbs+0HYPDydcN3tPLDLyuDT+aPp8ORrQiQBsyqEsL1bdu6ujS3514h378fAY4AN7frT0rLWzB4++twu1i3fX/w/sKJeFlySRTh2xaxdPjwAx5s2VKfHykF7hXyV17BOnkyfmIg7fp7dLhLiVJWso4kkH1sJwZLJqaaLTBVyZv9tNfOlXyyZgZW4L2wanT49WeebN68+BtLS0PVqAGpqeDpCf7+SFAQVK0KXbvCww9DRIRjXpibcK+Q+/pC3boE792LJT0Jo5+egdPZrDlZpC+ZQmDOOfrf3IO01FTmzX+PHAARrJ5eLPP0Iv722xn5+XQ8PUt42a2/P3/cdBNhs2ZRKTsbv6wsvM6dw3jgAPL776iXX0b5+iL16yO1auXdItykSd5twlFRjnjJ5Y57hRwwXHMN7N3LdXvWkdCit6vLcXuZa6bTJiaYH+cu/yfA06f9dz1HXpkwu2pVZrS9nUo3DEJZLVgzUvC/cJzbt/zCTce20zg9Gf/4eIiP/1c7ZTRCaChSrx60aJG3x8/KyhuyascOOHAg7/Hff4MbTdHsfpe1/v47qkMHjnr6cOMTs0unsArs7CeDSNi8ntq1a5fK9l59fTxvvfcx/r2fwzOk8D2zUgpr0glqHUmg8cm91L1wgsapZ6mXkUyYJRcT/Ofegfx21m9Ao927nFK/s1Ssy1rbt8fcujU1169nzvRH+d/97+tbGp3IOyiUAwcOlErIv/rqaya+8z6V7nwTj4D/zOXxDxHBGFyDI8E1OHJdT37Ot0zlmvE9uYfr964j6GIiWQYPUkzerI+oh09ACPHzx5N03r2up3e/kAOe69aRGhND7PEjTP1mNIMHTnJ1SW7LWrMVv/y2mO7du3Ps2DGeePoZKgUGMG7si1SvXt1h21mxYgVDh48g8LZXrxjwqxEPE5lR17Ik6r9j5Pmn5oVbytGnW3u45y7Ow4OAI0fICg6m68k93Lhzpasrclue4XX55ttZ3DlgIA0aN2XlSWHWis1MeffdK7a7cOECu3btIikpiT/++IM1a9Zw8jL3l2dmZnL/Q4/g23noFW+/LSmLbaAQHfLywmDAe+VKFDD157cJS67YEzM4i090M4ydHmXZ+UAq9Z9IQLu78A0M5po6dQpdPykpiZq16hBRI4pWHToSXj2S3gMe5LYHH6dO/YZ0jOvOtm3/P1/HtM+/oFJQZdIrxeBzTWunvpbgjOS8B24Wcrf8uP6PJk2wTJiA55gxrJz6EBM7P8iMFn1cXZXb8YluBtHNAFAWMxkHN9Kly0cApKSksHTpUry8vMjJyeHVNyaSFlCTqre/g8i/9zE+5iy2b1/GDZ3juK1fP9q2imX69On41WmFf/cRTn8dD2/Mm3t9W736tHP61kqPe4ccMI0eTer+/fhOm8Yry6byzMrPORoUwfhOD7C6dqEnI7USyEk8QET1Gv+ciJs7dy7DR7+Ab2gkiAFrZAv8m3T/T8ABDCZvApr3Iqd6A+b/vZ1F8XNRUhnPJqUzF3vTxP0oIHL0qFLZXmlx+5ADBHz2GTmvvsqRbt2pvHMHdc8f5bO5r9Dy0a9I9tMjhzqSwSeQ06dOkpiYSLVq1bBarfhGXYtvV/uHefKsWhvPqqXzldy/2D6lN4iMLP1tO5H7HpMX4FmtGrUStlHZYiF17Fg8lWLmjJGuLsvtmIKr49WsF81atOTo0aNERUVBWvk4H7Kg4Y0IkPvEE64uxaEqTMjzC3z5ZbLq1KFJ6jnuXz3D1eW4Hb/Wt5MT1YpJb08hPDwcswtnYC2Kz1reQprRxDVr15K5bp2ry3GYChlyAO/Vq7EYjbz012y67NZT7zqad4NOzPphNjExMWQmncWane7qkq4q18OTR257EQFO3X67q8txmAobciIiMNq+Ypvw8zuursbtmEKjyfUOYu7cuTRtHkvW8fJxmejvMc3J8PAkLPF03gAkbqDihhyQ66+Htm2pYjFT+8B6V5fjVkQErw7389iIJ/A0Qu7FxKs3KiPWxjTHH8Xe8eNdXYpDVOiQAxiGD0eAe3/X0+w6mldEPbza3cu6NavBYnZ1OXYb2WsUv4XVps+331OebuC6nAofcm6/HavBwL2J+/np88eoXo72OOWB37VxVL17Ir4NO7q6FLsl7V7JnWSRfDFJh9wteHhg/vVXzgQE0PjsYVZNfZhe+lp3h/Ku0fBf882VdZZ9axj33NOcOHbELSZ+KP+vwAG8unWjakoK58ePx6AU7y+azM/Thuq9egVlDamF2WzGaLzyzLblhQ55PiHPPot1507OhIbS8Pwx1n36ED9PH0anAxtcXZpWiuTCEaKjo11dhsPokBdgatiQqmfOcOaDDzjt50fDc0f5Yu449r7Vl6+/e46mJ/e4ukTNibJP7sUj/Rx9+/Z1dSkO437DPznY6aVLOTtyJFF79hBotaKAZE9f1tZqzt7QaM76VmZx3TZctHM2Uq3syjq6nfSl7zFlwms8+MADri6nSEo0dbGIeANrAC/ybmiZo5Qam2/5KOAtIFQpdc7etiISDHwPRAOHgTuUUklXqsUVIc/v0BdfkPvCi9Q8eYL8Y44q4LxPAAerRHE8MIw0L1+qZFzEoGCPbZCD7dXqsLJOK5fUrV1eZFIiJyqFknFiN6mLJjBj+lT+97//ubqsIitpyAXwU0qliYgJWAeMUEr9JSKRwDSgPtCikJBfqe1E4IJS6k0ReQaorJQac6VaXB3ySyxmMzs++ojT69bheeoUkdt3EJmSfMUBAhWwNTSGJ/uO4VCVGqVYrVaQwZrLo3/MZlD8IoIzU7i1fhuWHt/BnO9n0bNnT1eXVywlGshR5b0LpNmemmw/l94Z3gFGAwuK0bYv0NH2eAawCrhiyMsKo8lE0xEjYMT/D2RgNpvZunEjKevWkXrkCLnVq1PFYiFt+3awWrluxQquO3uIFdOGcNo3iDlNuvJp69tI9fZ34SupeB7c8CNjVs/A05qLBVjmYaJJ99Z8OOoHatRwzzdfu47JRcQIbAbqAB8qpcaISB+gi1JqhIgcBmIL7skv19b2+4tKqaB86yUppf4z2LWIDAYGA0RFRbU4cuRIMV5m2XDs/fcxvzKOmufPYSTv3e6MbxCrascy6fqBnA2o4uoS3VZI2gW+/v5F6p87Qg7wQ3Awhjfe4LZBg/DyKv+TYJbo43qBjoKAH4ERwGdAN6VU8pVCXkjb4UqpHfaGPL+y8nG9pDLPn2f/0KEE/vobNdJS8QCswObwujx/0+PsC412cYXu5YGN83lu5RcYlYWtBiMbXn2Vh58Z4xYXulxypZAX6VUqpS6S97G6LxADbLMFvAYQLyLV7Gjbw/ar0yISbiswHDhTlFrKM58qVWjyww9Ep6Zgychgx+OPc8Hbm9hT+1j8+WMsnTqYXrv07a8lZrUyffbLvLhiGlZlYUJ0LQL37eWR5551q4BfzVVfqYiE2vbCiIgPeTPebFFKhSmlopVS0cBxoLlSKtGOtpe+aF4IDLI9HsRljuvdnZePD43ffZeQzExOTJ/OqUqVqJN0kg9+msiOyf0Yvu5bt7nlsTSJsrLqs8F0+XsTJ8TAxy++yOiD+0ttppeyxJ63s3BgpYgkABuBpUqpRZdbWUQiROQXO9q+CXQVkf1AV9vzCq3GAw8QcfEi5zdsYE+jRnjl5vDU79+SMOV24vb96eryyhelOGbyZoHJxO5ff2HEuHEVau+dn74YpgyzJCdztHt3aq5fjwDLo5ry0J2v6mmf7JB9ci+pP41n2a+LaNOmjavLcTqHHZNrpctYqRIxf/1Fdnw8yf7+xB3dxvJPHsAjN8fVpZVZSinSd60m7afX+XbG5xUi4FdTIYZkLu98rrsOn+RkLjRpQu2dO1n8yYPc+sCHpPgGurq0MsGSlUb20e1YT+zAcjyBSj4mlq5cTrNmzVxdWpmgP66XM6ktWhAQH48FWFq7NcP6PYvVUHHeqy3pSaTvWI7HxWOQlYw5+SyZSWdp3qoNvbrH0blzJ1q1alXhjr8d9j25q+mQ5znx9NMETJlCYG4ua8Jqce/977m6pFKhlJXkb56kZ8e29OgWR7Vq1QgPD6dBgwZ4eFScN7rC6JC7I6uVjOrV8UlMZHdwDUb1HsWuaoVPMuguMo9sI2DLN+zfvYO82yK0S/SJN3dkMOC7fTvJVarQ8MJxvp8xkuQ/ZmHNznB1ZU6Te2QLA++6Uwe8iHTIy7OQEILOnePMkiU8eW0TmnufJemr4WQedNNPO9ZcfHy8XV1FuaND7gbCunbls4RtLF/yGz/O+hrDX5+TtngKucmnXV2aQxnCGzB7/k+uLqPc0SF3M3FxcRzYs4s72jfk1GePYMlIdnVJDuNbuxU7d+4kISHB1aWUKzrkbsjPz4+PPniXhk2vw3zuqKvLcRjxMOF3w/307H0LSUlXHERIy0eH3I01a3ItOWcPu7oMh/Jr1Jms8Gbcent/rPrGHbvokLux9q1bYkwqv4NsXI7fDfex7eAJ3plSMa4PKCkdcjcWEhKCITvt6iuWM2Iw4tttJGNfGcf27dtdXU6Zp0Puxho0aEDW6UNuMZ9XQabK4fjccD9db7qZgwcPurqcMk2H3I01aNCA4Er+ZOz7w9WlOIVfo85YGvehdbsO7NpVPuY/dwUdcjdmMBh4d/JE0lZ8isotP1MHF4Vfs5sgtj+dunZHX/JcOB1yN9e7d29u7NCetBWfoJR7no32vzYOa4v+dIzrzgcffkhWVparSypTdMjdnIgwe9Y3RBkvkrbsI6xm9wyAX8OOBN42jpc/mElYeHUGDx1GfHy8q8sqE3TIK4CAgABWL19Ch5r+XJz5OBkHN7q6JKfwDI3Gv+9LBNw5kbk7U7ihS3cWLlzo6rJcTt9qWsEsXbqU+x4cTFblGHxveBCjX97Q98qSS+ahzeSc/htDbhbKOxC/BtfjERjm4oqLL3XLL3QLTeO7r2e4uhSn07eaav/o2rUr+3bv4K6OTbkw4zFS/5xF6toZnJ/2EBFHljCkQxSj+8Zyax0PLn7zJOk7lru65GLzqtGI3377lYwM97391h56T16B7du3jxdffpWIiGoMfvABGjRo8K/le/fupWXb9vj3exXPkCgXVVkyab9Opm+ra5j6yUcYjUZXl+M0emQYrdg6dbuJHX7X4VevvatLKRZLVhrpP0+kZpCJObO+cdvJFfTHda1YTp8+zaYN6/Eqx8NKGb39Cej3MscqNaZJ81hGPvlUhbuDTYdcu6wnRo3Bq2FnPCpVdXUpJSJiwD/2VirfPYWvVm6n/Y2dsVgsri6r1OiQa4XKzc3lx3lz8Ynt5+pSHMYjMAT/bo9zOtvEnQPuISenYkxSoUOuFSohIQGfymEYfSu5uhSHEhH8bh7Diu1H6Nn7Fsxm97zcNz8dcq1QycnJePgFXX3Fcsjg6U3AzWOIP3qBW2+/0+336DrkWqFyc3NB3Pe/hxg9CLj5Gf48eI76jZu69SWw7vuvqJWIwWAAN72h5RLxMOHfczSp9XtxfccudOgUx1uTJnH8+HFXl+ZQOuRaofz8/LCas11dhtOJCH4NOxJ8/0fsq9yaCbOWU79REzZudJ/r+3XItUJFRUWRce642961VpDRJxDfuu3wj3sMz9jbmPL+R64uyWF0yLVCRURE0KZNW9J3rXF1KaXOK7Ixa9audXUZDqNDrl3WC888jXnz3ArxsT0/z6q1OH/hgtscm+uQa5fVuXNnOnVoS9qaL8g6vpuUTQtIWTODlN+/Jev4LpTVPa8aEzHgV6cVr70xwdWlOMRVQy4i3iKyQUS2ichOEXmlwPJRIqJEJKSQtpEislJEdtvajsi3rJmI/CUiW0Vkk4i0csxL0hxp2qcfEWE+SUD8TPrVMvDUzU25r3UEvhu/4OKMYaRtXuiWw0p5dxjEzK++cospma56F5rkzRPrp5RKExETsA4YoZT6S0QigWlAfaCFUupcgbbhQLhSKl5EAoDNwC1KqV0isgR4Ryn1q4j0BEYrpTpeqRZ9F1rZsn79eoYNH8n+I8fxaNAFvybd3OoKufTdq7H8MYN1q1fSsGFDV5dzRSW6C03luTRCv8n2c+md4R1gdL7nBdueUkrF2x6nAruB6pcWA4G2x5WAk1d/KVpZ0rp1a9b/sZZFs7/hpig4N/0R0he+Rsa+P11dmkP4NbgRY9uBdLixE3v27HF1OcVm1/3kImIkby9cB/hQKTVGRPoAXZRSI0TkMBBbcE9eoI9oYA3QWCmVIiINgMWAkPdm004p9Z85fURkMDAYICoqqsWRI+437Y+7SEpKYvHixTz59DOYr7sdv4YdXV2SQ6RtX4op4UfiN/xFeHi4q8spVInvJ1dKWZRSzYAaQCsRaQI8D7xkZwH+wFxgpFIqxfbrocATSqlI4Alg+mW2PVUpFauUig0NDbVnc5qLmEwmVq1ZR0Z6GmRedHU5DuN/bVdyat1AXI+b8y73LWeKdHZdKXURWAX0BWKAbba9eA0gXkSqFWxjO46fC3yjlJqXb9Eg4NLz2YA+8VaOmc1mevXtxw/rduF7y1h8m/d1dUkO5d/mTk5lCi+Pe9XVpRSZPWfXQ0UkyPbYB4gDtiilwpRS0UqpaOA40FwplVigrZC3h96tlHq7QNcngRttjzsD+0v0SjSXUUrx0CNDSTiRTECPkXiG1SLvn959iAh+3Z7gnfc+YNGiRa4up0js2ZOHAytFJAHYCCxVSl32VYpIhIj8YnvaHhgIdLZ9VbbVdiYd4GFgsohsA8ZjO+7Wyp9Jk99mwZLV+N80CjG472CJRv/KBPZ+lrvuuZe9e/e6uhy76YEctRL56aefuGvQgwTdOQGPSuV3jPaiSE9YgufOhfy5djVRUWVjFFs9kKPmFH/99RcDBt5HQK9nKkzAAfyadCO7Tmdi27Rj//6yf5SpQ64Vy65du+hxc298uz6OV0Q9V5dT6vxb9kM1u43rYlvz8ivj2Ldvn6tLuiwdcq3IDh8+zA2d4zC1G4RP7UI/IVYIfk26UemON/hg0Qaat2rH7XfdQ1k8/NUh1+yWnJzM5MmTad6yNdKkL36NOrm6JJczVamBf5ehBD84lV8XLymTe3QPVxeglQ/Z2dl0jOvO32fT8On+NN41Gly9UQViMHnhE1yNM2fOUK9e2Tp80SHX7PLw0Ec5muFB0J1vIG48wGNJGHwqcUn5KakAABR9SURBVObMGVeX8R/6X0uzy4L5C/C58WEd8CvI9Q3m77//dnUZ/6H/xTS75JpzEJOXq8so0zyirmPal19htZat++t1yDW7iEiZHaJZWXJJ2/QjGXOfI+2n8WSfcs3JL++Y5pxJMzNr1iyXbP9ydMg1u4SFVyc35ayry/gPpRRpv06mbs5BZr4/gbFD7yL7lzdJWzm11MemE4MRrw73MWz4CI4ePVqq274SHXLNLrVq1SI36ZSry/iPjN1rqJRzlpVLf6NHjx4Mf+wxDu7bQ7vqJlLmvIAlM+XqnTiQd2RjDHVvZNxr40t1u1eiQ67ZpW3L5ljOla2TSum715C5ZhrffzMTL6//P19QuXJlFsydzX233UTKnBexZCSXal1+rW7jux9ml5nRZHTINbu0ad0Kj/OHXF0GALmp50ldNAGf7fNYuXQxLVu2/M86IsI7kyfx0IDbSJ03tlT36EafQHyu682QRx8vE/Og65BrdmnUqBGZZ117nGnNziB1/RySvh7Jg706sHdnQqEBv0REmDTxTQbd0YfUeS9jyUwttVr9Ym9h24FjzJ49u9S2eTk65JpdoqKiyE5Ldtm0SVnHd3Fh5nCur5LBxj/XMfHN8fj4+Fy1nYgwZfIk7r71JtJ/eavUri0XownPVnfy3EuvuPwrNR1yzS4Gg4HQahFYXHCGPTf5NKmL3uS7Lz9j/pzvizw8sojw7tuTCPU0k7G79KZ98q4Vy8UcYc6cOaW2zcLokGt2CwwKwpqVXqrbVFYL6b+9zYvPjaF3797F7sfDw4OP3nuHnA2zSu2rNRHB1PJOxrww1qV7cx1yzW7RUVHknCvdIbHTN82nYVQYo0eNKnFfXbp0oVvHDmSsm+GAyuzjXasFyTnCjz/+WGrbLEiHXLPbiEeHYN3xW6kd11oyU8jcNI+vv5yOwVD8/6pKKdauXcsLL77Elq1byU0svdFcRARjk5uZ/N6HpbbNgvRdaJrdunXrRmRoJU5uXoh/rPOHXLakniOsajgxMTFFbms2m0lISGDt2rW8//FUzqdmYohpjce1dxBYo5ETqr0837rt2PbZdP7++29q1apVqtsGHXKtCAwGA/PnfE/T5rGYajQi99AmMhN+wyuyEX5xwzF4ejt0ex6BYRw+uJeeffoRUiWYBnVrExERgZeXF2fOnOHI0WMYjAY8TSZycy0cPXGSYydOkngqkWOH/8a3SjU8qtXD0GwAlaKbuWyYaPHwxKfBjXz62TQmvFH6V8Lp0Vq1Ips58ysGDbqXsPDqLFv8K7ffdQ9JDW7Fp/blv7MurpxzRzGfOYQ1Ox2VnIhHTgpiMWPxDsTiXRnIOzknYsDgF4TRrzJG/2BMlcMxePk5vJ7iyjl3lKwFL5N44hienp4O7/9Ko7XqPblWZAMH3oPZnENcXBw1a9akXr16HDx/gk7ZGWyrdg1HgyMcti3PkCg8Q5w77LH53DEyd60goP3diNE5kfAMiSInOJLvv/+egQMHOmUbl6P35FrxffIJ1lGjsKan/7O3mBlWi5F1YkEpDL5B+LXoU6ZnU7FkJJM8+zkyz50g8LqeBHR23hwfGQc2ELxnPru3b3X434ked11zjqwsDOnpHAfW+AewsG5dLt4Sx6ibmzK6T3NUwgKyDm91dZWXZclIJmXuiwwZNIDTp0/jc2YHaTtXOm17PrVjOXUuidLeUemP61rxjRzJ2bvvpnpQENEm038WR0ZGMuzpF/CqMRlDGRtV5p+A39ufN8e/hojw84J5tLuhE96R1+IRGOLwbYoYMNVpx2fTv7jiNfeOpvfkWomEhoZiKiTgAAMGDKBBrSgy9/9VylVdmSUrjdQfxzL4njv+CThA06ZNGfn4Y2St+9Jp2/ZuchPfzprFxYulN7WzDrnmNCLC48Mewbx1gctubCnImp1B2vxXGNivFxPfHP+fY+Nnx4zGmrib7MQDTtm+R2AI3jEt+HTqVKf0Xxgdcs2p+vfvT7f2saT/9rbLZxexZqeTuuBVbo3rwHtT3i705Je/vz/jxr5I1prpKCeNaefRoAsffTqt1P4+dMg1pzIYDHwz80siPLNJ/+Nbl9VhyUwlZd5Y/tetPdOmfnLFs9tDhwyhdlgA6RvnOaUWr8jGJKVmsHHjRqf0X5AOueZ0JpOJ5Yt/wefYetI2/1Tq27ekJZEy5wUG3daTqR9/dNXr4I1GI/N++I7s+AXknHX8DTkigunaHoyfMMnhfRdGh1wrFWFhYaxdtRzT7l9IX1t6d4GZzx8n+YdnGDn4XqZMnmT399NRUVGMe3ksWas/c8rHar8m3Vi6fDm7du1yeN8F6ZBrpSY6OprdOxLwOrGJ1MVTyE1x7JRCBmsu/bYvp+3hrVjN2aRuWkDyD88y6fWxvDL2pSJfgPL48McI9jCT+bfjv9c2ePnh3aQHb709xeF9F6S/J9dKVVBQENs2b2T8G2/w2Yzn8WzVH9/6HYr9Pbpnbg6DN8zj9oSlRCafwYBiSWAoC5SZtu3a8t4fa2nUqHh3nRmNRsa99ByPvfAGqlYLh08R5XVNexYvecuhfRZGX9aqucySJUt47c232Lx5M361mpMbUgev6g3wDItBDMbLN7RauXXnSh7ZMI+6545iQGEFDosQH1aN8/fdS4+hQ6lZs2aJa7RYLFzbPJaztW7Cr177EveXn7LkcvbTQRzYs4uIiJJd71+iG1RExBtYA3jZ1p+jlBqbb/ko4C0gVCl1rkDbSGAmUA2wAlOVUu/mWz4ceAzIBX5WSo0u4mvTyrFu3brRrVs3Dh06xMqVK1mxZh3rfv+YxBPH8b+mFaYmN+NVvf4/67c7vJXnVn5O/bNH8FAWFHBMhNURkchzz9Dnnnv4X2CgQ2s0Go2MfmIET0/6DBwccjF64F/zWtasWUP//v0d2ve/tnO1PbnkHcj4KaXSRMQErANGKKX+soV4GlAfaFFIyMOBcKVUvIgEAJuBW5RSu0SkE/A8cLNSKltEwpRSVzxI03vyiuH8+fN8OWMm4157HY8GXQhtFMeab8cQlp6EAs4g/BUUxKkhQ+j3xEjCwsKcWs+FCxeoWesaAm57Fc/Qkn86yC9l43z6RMOX0z8rUT8lukFF5UmzPTXZfi69M7wDjM73vGDbU0qpeNvjVGA3UN22eCjwplIq27a87E3srLlElSpVeOrJJ9ixbQu96gWQ+c2TXMxOY6PByORB92E5cZy+SRcY8sZ4pwccIDg4mIlvvk7GkndQuTkO7duzam3W/fGnQ/ssyK4zCSJiFJGtwBlgqVJqvYj0AU4opbbZ2Uc0cB2w3varusD1IrJeRFaLSOldsa+VC5GRkcz4fBrxRw9zcukSmmSkM+rLL0p8/FocQx55hHbNGpH++zcO7derRkMOH9iH2Wx2aL/52RVypZRFKdUMqAG0EpEm5H3Ufsme9iLiD8wFRiqlLs1X4wFUBtoATwM/SCHfcYjIYBHZJCKbzp4te7Nqas5XpUoVOnbs+K/5zkqbiDDzi2lk71mJ+fxxx/VrMGLy8iY52XnztRXpOwGl1EVgFdAXiAG2ichh8sIfLyLVCraxHcfPBb5RSuW/TvA4MM92OLCBvBNz/7m/Tyk1VSkVq5SKDQ0NLUq5muZQISEhPDVyBNkbHTv1kW+VCA4dct48c1cNuYiEikiQ7bEPEAdsUUqFKaWilVLR5AW2uVIqsUBbAaYDu5VSbxfoej7Q2bZeXcATOIemlWFPPjES89GtmM8dc1ifElqLVatXO6y/guzZk4cDK0UkAdhI3jH5osutLCIRIvKL7Wl7YCDQWUS22n562pZ9DtQSkR3ALGCQKk9f2msVUmBgIHcPGEDWvnUO69MQ1ZxZc+Y7rL+C9MUwmlZE27Zto32nrgTf/4lDRryxmrM5/dFAzp89jb+/f7H60GO8aZoDNW3alBs6tCdji2PuqDOYvAiIqM2ffzrnqzQdck0rhklvvk5m/EKHjXhjrVqflaucc1yuQ65pxdCwYUPatm1D+s5VDunPM6YFc+YvdEhfBemQa1oxPTn8Udi/yiF9eUXU4/ixoxw+fNgh/eWnQ65pxdS5c2csF0875Os0MZrwadiJqZ9Nc0Bl/6ZDrmnF5O3tzdAhg8nessAh/XnW78hnn3+BxWJxSH+X6JBrWgmMHvUUliObHTKEs1e1Olg8A1m+fLkDKvt/OuSaVgKVK1fmpReex+yoS13rdODzGV85pi8bHXJNK6EhjwzGcnqfQ47N/Rp2ZOHChQ6dYUWHXNNKyMfHh8EPP0R2ws8l7svoF4Rfjfr89ttvDqgsjw65pjnAyMeHk7F7DZastKuvfBXWmq34YZ7jrmXXIdc0BwgPD6db9+5kbv3l6itfhXfUtaxZvcZh473rkGuag0x4/VUyt/yEJT2pRP14VI7ALCY2bNjgkLp0yDXNQerVq0fv3r3I2L6sRP2ICB4NOvHeh584pC4dck1zoJdfeI6sLQuxZqeXqB+fhp1YMP9HcnNzS1yTDrmmOVD9+vWJi+tCesKSEvXjERCCV0gNh1wYo0OuaQ720nPPkLN1EcpasstTVWRz5i+87CBMdtMh1zQHa9GiBdXDq5J5sGSjGHlFX8dvS0t2fA865JrmFBNeH0fu5pJd6uoZFsPJY0dJSUm5+spXoEOuaU5w8803k3X+BJb04l+eKkYTgVH1WbasZHtzHXJNcwIPDw+6de9B5valJerHUqMFCxaV7AIbHXJNc5JJb44nI35+iS519YpsxNrfSzbAow65pjlJnTp16NatOxlbfy12Hx4BVTiTeLJEl7jqkGuaE01683Us238m8/DWYrU3+AaRlZFOTk7xZ1PVIdc0J6pTpw5zvv+OjGUfYM3OKHJ7EcEnsDKJiYlXX/kydMg1zcni4uK4pddNZPxRvBFffEIi2b17d7G3r0OuaaXgnUkTyd77O+YLJ4rcVjxMJbqGXYdc00pBSEgIwx8bRvaWok+goJQVo9FY7G3rkGtaKXls2FAy9qzDkpFctIYKDIbiR1WHXNNKSUREBHcPuIvMDUW73FWZM4s92ynokGtaqXpt3Mtk7l5FbvJpu9tYszPx8/Mr9jZ1yDWtFFWtWpVhQx4ha7N9AzUqi5m0cyeJiYkp9jZ1yDWtlD05cgTpu1bZNXpMzum/qVEzhkqVKhV7ezrkmlbKwsPD6dHjJtLjrz4gRPbxnXS6oUOJtqdDrmkuMGnCeLK2/HTV780Nx7fSs0e3Em1Lh1zTXKB27dq8NWE8aYvewGrOKnQd88VEshIP0rNnzxJt66ohFxFvEdkgIttEZKeIvFJg+SgRUSISUkjbSBFZKSK7bW1HFLLOZdtrmjsbNnQo3W9sR/qKT/9zl5lSiqw103h61FN4e3uXaDv27Mmzgc5KqaZAM6CHiLSBvBADXYGjl2mbCzyllGoAtAEeFZGGlxba0V7T3Nq0Tz+mmvUcaUvexXzhBJb0JHJOHyRt8RQivMw8M2Z0ibdx1ZCrPJfuejfZfi697bwDjM73vGDbU0qpeNvjVGA3UD3fKldsr2nuLjAwkD/XruLO9g3I+WkcyV+PxGPNB9zRoSF/rVuNl5dXibfhYc9KImIENgN1gA+VUutFpA9wQim1TUTs6SMauA5Yb3tuV3sRGQwMBoiKirKnXE0rVwIDA/nog/f46IP3nNK/XSFXSlmAZiISBPwoIk2A5wG7TvuJiD8wFxiplEoREV972yulpgJTAWJjY/UeX9OKqEhn15VSF4FVQF8gBtgmIoeBGkC8iFQr2EZETOQF/Bul1Dzbr2vb217TtJK56p5cREIBs1Lqooj4AHHABKVUWL51DgOxSqlzBdoKMB3YrZR6+9LvlVLbgau21zSt5OzZk4cDK0UkAdgILFVKXfZSHRGJEJFLY8i2BwYCnUVkq+2nZF/6aZpWJFfdkyulEsg7YXaldaLzPT4J9LQ9Xgdc9axc/vaapjmWvuJN09ycDrmmuTkdck1zc1KSmRlKm4icBY64YNMhQHk686/rda6yWG9NpVRoYQvKVchdRUQ2KaViXV2HvXS9zlXe6tUf1zXNzemQa5qb0yG3z1RXF1BEul7nKlf16mNyTXNzek+uaW5Oh1zT3FyFC7mIPGEbb26HiHxnG8MuWESWish+25+VL9O2h4jsFZEDIvJMvt/b1b4M1fuWiOwRkQQR+dE2TkCZrTffcoePB+isekVkuG3ZThGZ6Kh6i0UpVWF+yBt66hDgY3v+A3AfMBF4xva7Z8i7lbZgWyNwEKgFeALbgIa2ZVdtX8bq7QZ42B5PKOv12pZHAovJuxgqpCzXC3QClgFetudhrvx/X+H25OTdeecjIh6AL3CSvEEwZtiWzwBuKaRdK+CAUupvpVQOMMvWDjvbl5l6lVJLlFKXJrz+i7xBO8psvTbOGg/QGfUOBd5USmUDKKXOOLjmIqlQIVdKnQAmkTc67CkgWSm1BKiqlDplW+cU+Qa0yKc6cCzf8+P8/6CU9rQvS/Xm9wDwa1muN/94gI6o09n1AnWB60VkvYisFpGWjqy7qCpUyG3HVpeGrooA/ETkHnubF/I7p37/6Ox6ReR58obN/qYkdebrz+H15hsP8CVH1PivDTrv79cDqEzeMORPAz+IPaOdOkmFCjl5Q1cdUkqdVUqZgXlAO+C0iIQD2P4s7OPVcfKOCy+pQd5HO+xsX5bqRUQGAb2Au5XtwLGM1uvM8QCd9fd7HJin8mwArOTd1OISFS3kR4E2IuJre2ftQt5Y8AuBQbZ1BgELCmm7EbhGRGJExBPob2uHne3LTL0i0gMYA/RRSmU4qFan1KuU2q6UClNKRau8EYSOA82VUollsV7bsvlAZwARqUveiTnX3bXmyrN+rvgBXgH2ADuArwAvoAqwHNhv+zPYtm4E8Eu+tj2BfeSdVX0+3+8LbV+G6z1A3vHkVtvPJ2W53gL9H8ZBZ9ed+PfrCXxt6zOevBmIXPZ/Xl/WqmlurqJ9XNe0CkeHXNPcnA65prk5HXJNc3M65Jrm5nTINc3N6ZBrmpv7P/oYFV47H8eOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f, ax = plt.subplots()\n",
    "rgidf.plot(ax=ax, edgecolor='k');  #for_intersects\n",
    "divides.plot(ax=ax, edgecolor='r');  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdir=gdirs[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_new = xr.open_dataset(gdir.get_filepath('gridded_data')).load()\n",
    "ds_new.close() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### calculate the length in meters of the divides plotted in red above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "shape=gdir.read_shapefile('intersects')\n",
    "#     len(shape.geometry)\n",
    "ds_new['divides_len']=shape.geometry.length"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "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>RGIId_1</th>\n",
       "      <th>RGIId_2</th>\n",
       "      <th>min_x</th>\n",
       "      <th>max_x</th>\n",
       "      <th>min_y</th>\n",
       "      <th>max_y</th>\n",
       "      <th>geometry</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>RGI60-13.52589</td>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>-4178.261806</td>\n",
       "      <td>-1582.388811</td>\n",
       "      <td>3.794838e+06</td>\n",
       "      <td>3.796791e+06</td>\n",
       "      <td>LINESTRING (-4178.262 3796790.794, -4169.071 3...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>RGI60-13.52593</td>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>-347.775523</td>\n",
       "      <td>-162.668086</td>\n",
       "      <td>3.791703e+06</td>\n",
       "      <td>3.791739e+06</td>\n",
       "      <td>LINESTRING (-347.776 3791702.737, -345.728 379...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>RGI60-13.52595</td>\n",
       "      <td>29.685256</td>\n",
       "      <td>297.439355</td>\n",
       "      <td>3.796691e+06</td>\n",
       "      <td>3.796871e+06</td>\n",
       "      <td>LINESTRING (294.956 3796871.224, 297.439 37968...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>RGI60-13.52589</td>\n",
       "      <td>-4169.071378</td>\n",
       "      <td>-1585.868390</td>\n",
       "      <td>3.794848e+06</td>\n",
       "      <td>3.796792e+06</td>\n",
       "      <td>LINESTRING (-1591.287 3794864.811, -1585.868 3...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>RGI60-13.52609</td>\n",
       "      <td>1303.832237</td>\n",
       "      <td>1697.370107</td>\n",
       "      <td>3.795211e+06</td>\n",
       "      <td>3.795509e+06</td>\n",
       "      <td>LINESTRING (1697.370 3795508.861, 1688.872 379...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>RGI60-13.52609</td>\n",
       "      <td>1827.243298</td>\n",
       "      <td>1943.640691</td>\n",
       "      <td>3.795926e+06</td>\n",
       "      <td>3.796584e+06</td>\n",
       "      <td>LINESTRING (1830.580 3796582.760, 1840.798 379...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>RGI60-13.52593</td>\n",
       "      <td>-339.519930</td>\n",
       "      <td>-161.587891</td>\n",
       "      <td>3.791708e+06</td>\n",
       "      <td>3.791741e+06</td>\n",
       "      <td>LINESTRING (-161.588 3791741.128, -178.504 379...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>RGI60-13.52600</td>\n",
       "      <td>753.712337</td>\n",
       "      <td>1215.382870</td>\n",
       "      <td>3.791054e+06</td>\n",
       "      <td>3.791248e+06</td>\n",
       "      <td>LINESTRING (1215.383 3791247.776, 1206.703 379...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>RGI60-13.53522</td>\n",
       "      <td>-855.525793</td>\n",
       "      <td>-538.542660</td>\n",
       "      <td>3.798163e+06</td>\n",
       "      <td>3.798727e+06</td>\n",
       "      <td>LINESTRING (-855.526 3798721.669, -846.239 379...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>RGI60-13.53522</td>\n",
       "      <td>-540.017538</td>\n",
       "      <td>297.439355</td>\n",
       "      <td>3.796856e+06</td>\n",
       "      <td>3.797979e+06</td>\n",
       "      <td>LINESTRING (-540.018 3797979.386, -536.675 379...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>RGI60-13.53522</td>\n",
       "      <td>890.382786</td>\n",
       "      <td>1845.376371</td>\n",
       "      <td>3.796563e+06</td>\n",
       "      <td>3.796700e+06</td>\n",
       "      <td>LINESTRING (890.383 3796688.811, 891.370 37966...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>RGI60-13.53507</td>\n",
       "      <td>-3786.346141</td>\n",
       "      <td>-1784.542609</td>\n",
       "      <td>3.797148e+06</td>\n",
       "      <td>3.797588e+06</td>\n",
       "      <td>LINESTRING (-3786.247 3797252.297, -3786.346 3...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>RGI60-13.53507</td>\n",
       "      <td>-1829.435605</td>\n",
       "      <td>-837.294206</td>\n",
       "      <td>3.797851e+06</td>\n",
       "      <td>3.798727e+06</td>\n",
       "      <td>LINESTRING (-1829.436 3797850.827, -1824.524 3...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>RGI60-13.53498</td>\n",
       "      <td>-4169.071378</td>\n",
       "      <td>-3777.171849</td>\n",
       "      <td>3.796780e+06</td>\n",
       "      <td>3.797266e+06</td>\n",
       "      <td>LINESTRING (-4165.020 3796780.298, -4169.071 3...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>RGI60-13.52595</td>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>38.266351</td>\n",
       "      <td>300.002738</td>\n",
       "      <td>3.796685e+06</td>\n",
       "      <td>3.796861e+06</td>\n",
       "      <td>LINESTRING (47.607 3796684.678, 47.794 3796686...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>RGI60-13.52600</td>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>762.481142</td>\n",
       "      <td>1209.513032</td>\n",
       "      <td>3.791050e+06</td>\n",
       "      <td>3.791244e+06</td>\n",
       "      <td>LINESTRING (773.434 3791050.256, 762.481 37910...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>RGI60-13.52609</td>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>1310.299730</td>\n",
       "      <td>1693.075798</td>\n",
       "      <td>3.795209e+06</td>\n",
       "      <td>3.795505e+06</td>\n",
       "      <td>LINESTRING (1314.390 3795209.240, 1310.300 379...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>RGI60-13.52609</td>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>1827.243298</td>\n",
       "      <td>1949.182307</td>\n",
       "      <td>3.795937e+06</td>\n",
       "      <td>3.796589e+06</td>\n",
       "      <td>LINESTRING (1949.182 3795936.598, 1937.332 379...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>RGI60-13.53498</td>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>-4178.261806</td>\n",
       "      <td>-3785.604773</td>\n",
       "      <td>3.796791e+06</td>\n",
       "      <td>3.797275e+06</td>\n",
       "      <td>LINESTRING (-3786.483 3797274.669, -3786.346 3...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>RGI60-13.53507</td>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>-3786.483275</td>\n",
       "      <td>-1784.542609</td>\n",
       "      <td>3.797148e+06</td>\n",
       "      <td>3.797587e+06</td>\n",
       "      <td>LINESTRING (-1816.619 3797587.395, -1809.775 3...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>RGI60-13.53507</td>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>-1835.183548</td>\n",
       "      <td>-846.238682</td>\n",
       "      <td>3.797862e+06</td>\n",
       "      <td>3.798737e+06</td>\n",
       "      <td>LINESTRING (-852.041 3798736.865, -851.810 379...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>RGI60-13.53522</td>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>-852.041416</td>\n",
       "      <td>-538.542660</td>\n",
       "      <td>3.798157e+06</td>\n",
       "      <td>3.798737e+06</td>\n",
       "      <td>LINESTRING (-553.796 3798157.346, -561.682 379...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>RGI60-13.53522</td>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>-533.565104</td>\n",
       "      <td>300.002738</td>\n",
       "      <td>3.796850e+06</td>\n",
       "      <td>3.797971e+06</td>\n",
       "      <td>LINESTRING (300.003 3796849.809, 297.439 37968...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>RGI60-13.53522</td>\n",
       "      <td>RGI60-13.52594</td>\n",
       "      <td>882.198518</td>\n",
       "      <td>1848.732855</td>\n",
       "      <td>3.796563e+06</td>\n",
       "      <td>3.796701e+06</td>\n",
       "      <td>LINESTRING (1848.733 3796589.270, 1841.827 379...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           RGIId_1         RGIId_2        min_x        max_x         min_y  \\\n",
       "0   RGI60-13.52589  RGI60-13.52594 -4178.261806 -1582.388811  3.794838e+06   \n",
       "1   RGI60-13.52593  RGI60-13.52594  -347.775523  -162.668086  3.791703e+06   \n",
       "2   RGI60-13.52594  RGI60-13.52595    29.685256   297.439355  3.796691e+06   \n",
       "3   RGI60-13.52594  RGI60-13.52589 -4169.071378 -1585.868390  3.794848e+06   \n",
       "4   RGI60-13.52594  RGI60-13.52609  1303.832237  1697.370107  3.795211e+06   \n",
       "5   RGI60-13.52594  RGI60-13.52609  1827.243298  1943.640691  3.795926e+06   \n",
       "6   RGI60-13.52594  RGI60-13.52593  -339.519930  -161.587891  3.791708e+06   \n",
       "7   RGI60-13.52594  RGI60-13.52600   753.712337  1215.382870  3.791054e+06   \n",
       "8   RGI60-13.52594  RGI60-13.53522  -855.525793  -538.542660  3.798163e+06   \n",
       "9   RGI60-13.52594  RGI60-13.53522  -540.017538   297.439355  3.796856e+06   \n",
       "10  RGI60-13.52594  RGI60-13.53522   890.382786  1845.376371  3.796563e+06   \n",
       "11  RGI60-13.52594  RGI60-13.53507 -3786.346141 -1784.542609  3.797148e+06   \n",
       "12  RGI60-13.52594  RGI60-13.53507 -1829.435605  -837.294206  3.797851e+06   \n",
       "13  RGI60-13.52594  RGI60-13.53498 -4169.071378 -3777.171849  3.796780e+06   \n",
       "14  RGI60-13.52595  RGI60-13.52594    38.266351   300.002738  3.796685e+06   \n",
       "15  RGI60-13.52600  RGI60-13.52594   762.481142  1209.513032  3.791050e+06   \n",
       "16  RGI60-13.52609  RGI60-13.52594  1310.299730  1693.075798  3.795209e+06   \n",
       "17  RGI60-13.52609  RGI60-13.52594  1827.243298  1949.182307  3.795937e+06   \n",
       "18  RGI60-13.53498  RGI60-13.52594 -4178.261806 -3785.604773  3.796791e+06   \n",
       "19  RGI60-13.53507  RGI60-13.52594 -3786.483275 -1784.542609  3.797148e+06   \n",
       "20  RGI60-13.53507  RGI60-13.52594 -1835.183548  -846.238682  3.797862e+06   \n",
       "21  RGI60-13.53522  RGI60-13.52594  -852.041416  -538.542660  3.798157e+06   \n",
       "22  RGI60-13.53522  RGI60-13.52594  -533.565104   300.002738  3.796850e+06   \n",
       "23  RGI60-13.53522  RGI60-13.52594   882.198518  1848.732855  3.796563e+06   \n",
       "\n",
       "           max_y                                           geometry  \n",
       "0   3.796791e+06  LINESTRING (-4178.262 3796790.794, -4169.071 3...  \n",
       "1   3.791739e+06  LINESTRING (-347.776 3791702.737, -345.728 379...  \n",
       "2   3.796871e+06  LINESTRING (294.956 3796871.224, 297.439 37968...  \n",
       "3   3.796792e+06  LINESTRING (-1591.287 3794864.811, -1585.868 3...  \n",
       "4   3.795509e+06  LINESTRING (1697.370 3795508.861, 1688.872 379...  \n",
       "5   3.796584e+06  LINESTRING (1830.580 3796582.760, 1840.798 379...  \n",
       "6   3.791741e+06  LINESTRING (-161.588 3791741.128, -178.504 379...  \n",
       "7   3.791248e+06  LINESTRING (1215.383 3791247.776, 1206.703 379...  \n",
       "8   3.798727e+06  LINESTRING (-855.526 3798721.669, -846.239 379...  \n",
       "9   3.797979e+06  LINESTRING (-540.018 3797979.386, -536.675 379...  \n",
       "10  3.796700e+06  LINESTRING (890.383 3796688.811, 891.370 37966...  \n",
       "11  3.797588e+06  LINESTRING (-3786.247 3797252.297, -3786.346 3...  \n",
       "12  3.798727e+06  LINESTRING (-1829.436 3797850.827, -1824.524 3...  \n",
       "13  3.797266e+06  LINESTRING (-4165.020 3796780.298, -4169.071 3...  \n",
       "14  3.796861e+06  LINESTRING (47.607 3796684.678, 47.794 3796686...  \n",
       "15  3.791244e+06  LINESTRING (773.434 3791050.256, 762.481 37910...  \n",
       "16  3.795505e+06  LINESTRING (1314.390 3795209.240, 1310.300 379...  \n",
       "17  3.796589e+06  LINESTRING (1949.182 3795936.598, 1937.332 379...  \n",
       "18  3.797275e+06  LINESTRING (-3786.483 3797274.669, -3786.346 3...  \n",
       "19  3.797587e+06  LINESTRING (-1816.619 3797587.395, -1809.775 3...  \n",
       "20  3.798737e+06  LINESTRING (-852.041 3798736.865, -851.810 379...  \n",
       "21  3.798737e+06  LINESTRING (-553.796 3798157.346, -561.682 379...  \n",
       "22  3.797971e+06  LINESTRING (300.003 3796849.809, 297.439 37968...  \n",
       "23  3.796701e+06  LINESTRING (1848.733 3796589.270, 1841.827 379...  "
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### calculate the length in meters of the outlined plotted in black above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "shape=gdir.read_shapefile('outlines')\n",
    "#     len(shape.geometry)\n",
    "ds_new['perimeter']=shape.geometry.length[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<xarray.DataArray 'perimeter' ()>\n",
       " array(38652.19334194),\n",
       " <xarray.DataArray 'divides_len' ()>\n",
       " array(26657.7669086))"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds_new.perimeter, ds_new.divides_len.sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}