{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Compute smoother centerlines" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook describes how to compute centerlines with OGGM and write them to disk. It is meant for users who are mostly interested in the centerlines, not so much the rest of the model.\n", "\n", "We use an example of a user-provided glacier inventory and DEM (thanks to [Liss Andreassen](https://www.nve.no/hydrology/our-researchers/liss-marie-andreassen/) for providing the data)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import rioxarray as rioxr\n", "import geopandas as gpd\n", "import matplotlib.pyplot as plt\n", "\n", "from oggm import cfg, utils, workflow, tasks" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data preparation" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "fpath_inventory = 'zip://BreFlate20182019.zip/BreFlate20182019.shp'\n", "fpath_dem = 'norgedtm10nov21.tif'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Read the data" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "inventory = gpd.read_file(fpath_inventory)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this inventory, a few geometries have a topological error (the figure of eight where the outlines touch):" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
objektTypebreIDbrenavnbreTypeareal_km2forstedatadatakildemetodesceneIDhovedbreID...vassdragsOFTEMAGlobalIDcreated_uscreated_dalast_editelast_edi_1Shape_STArShape_STLegeometry
8None3161None00.022019-08-27NoneNone470...None0NoneFS_BRE2021-11-12FS_BRE2021-11-122.248363e+041021.901894POLYGON ((100216.542 6912509.442, 100236.493 6...
92None1956None00.182019-08-27NoneNone470...None0NoneFS_BRE2021-11-12FS_BRE2021-11-121.826024e+054047.790603POLYGON ((92498.356 6900920.848, 92500.210 690...
159None2456None00.192019-08-27NoneNone470...None0NoneFS_BRE2021-11-12FS_BRE2021-11-121.932463e+054669.023657POLYGON ((91962.885 6883230.562, 91965.661 688...
387None2605Liabreen00.622019-08-27NoneNone470...None0NoneFS_BRE2021-11-12FS_BRE2021-11-126.171481e+056595.974453POLYGON ((120594.079 6853553.546, 120595.925 6...
523None2656Fannaråkbreen02.752019-08-27NoneNone470...None0NoneFS_BRE2021-11-12FS_BRE2021-11-122.747803e+0615241.345864POLYGON ((124543.417 6842414.252, 124543.816 6...
..................................................................
6643Breflate2290None011.412019-08-27NoneNone457...None0NoneFS_BRE2021-11-12FS_BRE2021-11-121.141095e+0724954.300295POLYGON ((63155.410 6869314.958, 63181.449 686...
6649Breflate2306Flatebreen00.402019-08-27NoneNone452059...None0NoneFS_BRE2021-11-12FS_BRE2021-11-124.046787e+055277.768078POLYGON ((63964.380 6864753.248, 63965.304 686...
6719Breflate2266Lodalsbreen08.762019-08-27NoneNone451...None0NoneFS_BRE2021-11-12FS_BRE2021-11-128.758818e+0642808.850408POLYGON ((91018.242 6874020.521, 90986.386 687...
6721Breflate2237None00.592019-08-27NoneNone450...None0NoneFS_BRE2021-11-12FS_BRE2021-11-125.943676e+056680.173218POLYGON ((88863.476 6878889.301, 88873.452 687...
6725Breflate6762Storefonnbreen00.932019-08-27NoneNone451...None0NoneFS_BRE2021-11-12FS_BRE2021-11-129.274839e+059173.311529POLYGON ((89582.043 6864493.669, 89582.967 686...
\n", "

67 rows × 22 columns

\n", "
" ], "text/plain": [ " objektType breID brenavn breType areal_km2 forstedata \\\n", "8 None 3161 None 0 0.02 2019-08-27 \n", "92 None 1956 None 0 0.18 2019-08-27 \n", "159 None 2456 None 0 0.19 2019-08-27 \n", "387 None 2605 Liabreen 0 0.62 2019-08-27 \n", "523 None 2656 Fannaråkbreen 0 2.75 2019-08-27 \n", "... ... ... ... ... ... ... \n", "6643 Breflate 2290 None 0 11.41 2019-08-27 \n", "6649 Breflate 2306 Flatebreen 0 0.40 2019-08-27 \n", "6719 Breflate 2266 Lodalsbreen 0 8.76 2019-08-27 \n", "6721 Breflate 2237 None 0 0.59 2019-08-27 \n", "6725 Breflate 6762 Storefonnbreen 0 0.93 2019-08-27 \n", "\n", " datakilde metode sceneID hovedbreID ... vassdragsO FTEMA GlobalID \\\n", "8 None None 47 0 ... None 0 None \n", "92 None None 47 0 ... None 0 None \n", "159 None None 47 0 ... None 0 None \n", "387 None None 47 0 ... None 0 None \n", "523 None None 47 0 ... None 0 None \n", "... ... ... ... ... ... ... ... ... \n", "6643 None None 45 7 ... None 0 None \n", "6649 None None 45 2059 ... None 0 None \n", "6719 None None 45 1 ... None 0 None \n", "6721 None None 45 0 ... None 0 None \n", "6725 None None 45 1 ... None 0 None \n", "\n", " created_us created_da last_edite last_edi_1 Shape_STAr \\\n", "8 FS_BRE 2021-11-12 FS_BRE 2021-11-12 2.248363e+04 \n", "92 FS_BRE 2021-11-12 FS_BRE 2021-11-12 1.826024e+05 \n", "159 FS_BRE 2021-11-12 FS_BRE 2021-11-12 1.932463e+05 \n", "387 FS_BRE 2021-11-12 FS_BRE 2021-11-12 6.171481e+05 \n", "523 FS_BRE 2021-11-12 FS_BRE 2021-11-12 2.747803e+06 \n", "... ... ... ... ... ... \n", "6643 FS_BRE 2021-11-12 FS_BRE 2021-11-12 1.141095e+07 \n", "6649 FS_BRE 2021-11-12 FS_BRE 2021-11-12 4.046787e+05 \n", "6719 FS_BRE 2021-11-12 FS_BRE 2021-11-12 8.758818e+06 \n", "6721 FS_BRE 2021-11-12 FS_BRE 2021-11-12 5.943676e+05 \n", "6725 FS_BRE 2021-11-12 FS_BRE 2021-11-12 9.274839e+05 \n", "\n", " Shape_STLe geometry \n", "8 1021.901894 POLYGON ((100216.542 6912509.442, 100236.493 6... \n", "92 4047.790603 POLYGON ((92498.356 6900920.848, 92500.210 690... \n", "159 4669.023657 POLYGON ((91962.885 6883230.562, 91965.661 688... \n", "387 6595.974453 POLYGON ((120594.079 6853553.546, 120595.925 6... \n", "523 15241.345864 POLYGON ((124543.417 6842414.252, 124543.816 6... \n", "... ... ... \n", "6643 24954.300295 POLYGON ((63155.410 6869314.958, 63181.449 686... \n", "6649 5277.768078 POLYGON ((63964.380 6864753.248, 63965.304 686... \n", "6719 42808.850408 POLYGON ((91018.242 6874020.521, 90986.386 687... \n", "6721 6680.173218 POLYGON ((88863.476 6878889.301, 88873.452 687... \n", "6725 9173.311529 POLYGON ((89582.043 6864493.669, 89582.967 686... \n", "\n", "[67 rows x 22 columns]" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "inventory.loc[~inventory.is_valid]" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAN4AAAEDCAYAAABasLryAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAcJUlEQVR4nO2de5ScZX3HP7+9zGz2kmQ3u0mXXAjQRElUaFxAtCISNWqPpqC08VQLisVDUXs9B2Jty7FNi7TipT1F4wU5rQEjSMWqYKCIViUQIAJJCAkLJEtCdnYTkr1kZ3Z2fv1j3iUzu7MzszPz3mZ+n3PmvO887/O872+fne/7PM/vuYmqYhiGt9T5bYBh1CImPMPwAROeYfiACc8wfMCEZxg+YMIzDB8IjfBE5FMisldEdonITTPE+TMRedqJ8+cZ4Zc7YSkR6ckIf6eIPCYiTznHS4qwQ0Rkk4g8KyJ7ROTTFfkDjZqiwW8DpiIiFwNXquqVGWFvB9YDb1DVuIgszJHudcCfAOcDCeBeEfmRqu4DngYuA742JdkA8D5VPeSkvw9YXMDEK4GlwGtVNZXLFsMoRFhKvGuAG1U1DqCq/TninA08rKqjqpoEHgIudeLvUdW9UxOo6hOqesj5ugtoEpEogIi8S0R+LSKPi8j3RKQ1w5bPqWoqjy2GkZewCG8l8FYR2S4iD4nIeTniPA1cJCILRKQZeC/pkqlYPgA84ZSoncBngXeo6hpgB/CXTryzgD8UkR0i8hMRWVHyX2XULIGpaorIdiAKtAIdIrLTuXQdaTvbgTcB5wFbReRMzRjvpqp7ROTzwDZgGPgNkCzy2auBzwPvcoLeBKwCfikiABHg1861KDCmqj0ichnwLeCtpfzNRu0SGOGp6gUwYxvvL4DvO0J7RERSQCcQm3KPbwLfdNL8E9BX6LkisgS4G/hjVX1uMhjYpqofypGkD7jLOb8buLW4v9AwThGWquZ/A5cAiMhK0iXQwNRIk44OEVlG2plye76bish84EfARlX9Zcalh4G3iMhvO/Ganedm2QK8DXi2lD/IqHFUNVAf4GLg21PCIsB/kW7HPQ5c4oSfBvw4I94vgN2kq5lrM8IvJV1SxYEjwH1O+GeBEWBnxmehc+0S4FHgSefzfid8PmmxPkW6+nmO33lmn/B9RNWmBRmG14SlqmkYVUUgnCudnZ26fPlyv80wjII89thjA6raVe59AiG85cuXs2PHDr/NMIyCiMiLlbiPVTUNwwdMeIbhAyY8w/ABE55h+IAJzzB8wIRnGD5gwjMMHzDhGYYPBKID3ZjOYy8e5asP9dLVFqWrNUqnc+xqi7KwLcqiuU1EGuy9GVZMeAHluf4Rtu0+MuP1LR+/gDf/dqeHFhmVxF6ZASU2HM97vast6pElhhuY8AJKbMiEV82Y8AJKPuE11gvz5jR6aI1RaUx4ASWf8LpaoziLMBkhxYQXUAbytPE6rZoZesyr6RNf+OleHno29moXwasf53t/gRLPCDcmPJ949sgQT/YdLyntkaExbn/kAJ0Zou1sjRBtqK+wlYZbmPB8YmA4UXLap186wcbvPzUtfG5TA11tUf7zqgs4bf6ccswzXMaE5xOFugtK4cRYkhNjSVqb7N8adMy54gOq6orwAKINdbRFTXhBx4TnAyOJCU6OT7hy707raggFJjwfcKu0AxvREhZMeD6Qr4+uXEx44cCE5wNW4hkmPB9wVXjWuR4KTHg+4KbwbDhZODDh+YCVeIYJzwfMuWIUFJ6IvEZEdmZ8TojIn4tIh4hsE5F9zrE9I81GEdkvIntFZJ27f0L4KDS7vBwWmvBCQUHhqepeVT1XVc8F3giMkt77+3rgAVVdATzgfEdEVgEbgNXAu4H/EBEbvZuBq208q2qGgtlWNdcCz6nqi8B64DYn/Dbg953z9cAdqhpX1eeB/cD5FbC1Kkil1LWqZmu0gTkRe8eFgdkKbwNwu3O+SFUPAzjHhU74YuBgRpo+JywLEblaRHaIyI5YLDZLM8LL8ZPjjE+4s/21te/CQ9HCE5EI8H7ge4Wi5gib9ktT1c2q2qOqPV1dZW+wGRpcdaxYNTM0zKbEew/wuKpOLvZ4RES6AZxjvxPeByzNSLcEOFSuodWCjVoxYHbC+xCnqpkA9wBXOOdXAD/ICN8gIlEROQNYATxSrqHVgpseTRNeeChq4paINAPvBD6REXwjsFVErgIOAJcDqOouEdkK7AaSwLWq6s4cmBDirkcz4tq9jcpSlPBUdRRYMCVskLSXM1f8TcCmsq2rQqyqaYCNXPEcq2oaYMLzHHfHaTa5dm+jspjwPMZN4bW32LLuYcFWxfEYN/vx3nrTg3Q0R7IWx53cV+/cZfM5b3mHa882ZocJz0OSEykGR0pfT7MQqjA4kmBwJMEzLw9lXbvyzctNeAHCqpoecnQ0gbozWqwg5ngJFiY8D3GzfVcIE16wMOF5iAnPmMSE5yG+Cs8GUAcKE56HuNl5XgibmR4sTHgeMjDknkczHyLQ0WLjOIOECc9D/CrxOpojNNTbvzpI2H/DQ2JDY7481xwrwcOE5yF+OVdMeMHDhOchvgnPPJqBw4TnEWPjE5wYS/rybCvxgocJzyPcHKNZCFtrM3iY8DzCRq0YmZjwPMKEZ2RiwvMIE56RiQnPI2ycppGJCc8j3Jx5no+GOmHeHFsSImiY8CrESDzJibFxdIaZrn6VeJ2tUerqcq2qb/iJLf1QIb7/xEv87X8/TbShjq62KJ2t0ay1T/a8fMIXu6x9F0xMeBViskSLJ1P0HTtJ37GTPluUxoQXTKyqWSH8asMVwhwrwcSEVyH89Frmw0q8YGLCqxBBFZ5tZBJMTHgVIqjC62qzZd2DSFHCE5H5InKniDwjIntE5EIRuUFEXhKRnc7nvRnxN4rIfhHZKyLr3DM/GKiqr+up5MOqmsGkWK/ml4F7VfWDzpbMzcA64Iuq+q+ZEUVkFem90lcDpwH3i8jKat4j78RYkkQy5bcZOTHhBZOCJZ6IzAUuAr4JoKoJVX0lT5L1wB2qGlfV54H9wPkVsDWwBNWjCSa8oFJMVfNMIAbcKiJPiMg3RKTFufZJEXlSRL4lIu1O2GLgYEb6PicsCxG5WkR2iMiOWCxWzt/gO0Ft381prKclUu+3GUYOihFeA7AGuEVVfwcYAa4HbgHOAs4FDgNfcOLnGp80bRyVqm5W1R5V7enq6irB9OAQVOEtMI9mYCmmjdcH9Knqduf7ncD1qnpkMoKIfB34n4z4SzPSLwEOVcDWwBJU4fUdO8nKz/7k1PC1zGFsbVHetrKL0xe0FL6RUXEKCk9VXxaRgyLyGlXdS3rf890i0q2qh51olwJPO+f3AFtE5GbSzpUVwCMu2B4YgurRBBifUA4fH+Pw8elLC97yR2tMeD5RrFfzU8B3HI9mL/BR4Csici7pauQLwCcAVHWXiGwFdgNJ4Npq9mgCDAS0xCuEOV78oyjhqepOoGdK8EfyxN8EbCrdrHAR5BIvHyY8/7CRKxUgqG28QtjqY/5hwqsAYRReS6SelqjNCvMLE16ZTKTU1zUzS6XTqpm+YsIrk2OjCSZSPm1sXgY2T89fTHhlEuThYvkwx4q/mPDKJIztOzDh+Y0Jr0xCKzyravqKCa9Mwio8c674iwmvTMIqPCvx/MWEVyY2asUoBRNemZhX0ygFE16ZhLWqOZpIMprwZ4daw1aSLpuwCu8dN/8cSA8dm5yfNzlvb+3Zi3jbynBPTg46JrwySCRTHBsd99uMshhJTDAyOMoLg6Ovhi2a22TCcxmrapbB4Eg4S7tCmMfTfazEK4P25gh3XXMhsaH4qc9wnNhQgthwnN7+YYbi4WtHmePFfUx4ZdDUWM8bT++Y8fqnbn+CH/4mfMvNmPDcx6qaLmJLQhgzYcJzkbB2rne02LKAbmPCc5EwdjV0tERorLefhdtYDrtEPDnB8ZPh62owj6Y3mPBcYmA4fMtBgLXvvMKE5xJhrGaCbWTpFSY8lzCPppEPE55LhNWjacLzBhOeS4S1qmnC8wYTnkuEVXi2urQ32JAxlwir8J7rH2benEa62qIsaIkSabB3sxuY8FwirDPTb/jh7qzv7c2Nr87Te+uKLq65+CyfLKsuinqdich8EblTRJ4RkT0icqGIdIjINhHZ5xzbM+JvFJH9IrJXRNa5Z35wCatzZSrHRsfZ1z/Mr54bZF//kN/mVA3FlnhfBu5V1Q86e+Q1A58BHlDVG0XketLbM18nIquADcBq0htT3i8iK6ttj7wt2w+w5/CJrJnbp2ZyR0Jb1cyHOV4qR0Hhichc4CLgSgBVTQAJEVkPXOxEuw34GXAdsB64Q1XjwPMish84H/h1hW33lQf39rNt95HCEasIG05WOYqpap4JxIBbReQJEfmGiLQAiya3YnaOC534i4GDGen7nLAsRORqEdkhIjtisVhZf4QfVGOJVggr8SpHMcJrANYAt6jq7wAjpKuVMyE5wqZtp6Oqm1W1R1V7urrCt75HTQrPSryKUYzw+oA+Vd3ufL+TtBCPiEg3gHPsz4i/NCP9EiB807DzoKqh9VqWg5V4laOg8FT1ZeCgiLzGCVoL7AbuAa5wwq4AfuCc3wNsEJGoiJwBrAAeqajVPjMUTxJPpvw2w3PabYJsxSjWq/kp4DuOR7MX+Chp0W4VkauAA8DlAKq6S0S2khZnEri22jyatVjNBDhv0/0saImc8uJmeHPXnN7OmmXthW9iAEUKT1V3Aj05Lq2dIf4mYFPpZvnPx2/bQZ2kd9XJ/IF1tkbpOzZa+AZViGp6nuHAcIJnXs7u07vm4rNMeLPARq7kIJVSHnq2n/GJ8G2x7BfmeJkdNhAvB8dPjpvoZok5XmaHCS8HteixLBcT3uww4eWgVp0n5WDTiWaHCS8H1TLA2UusxJsdJrwcWIk3OyINdcxtMj/dbLDcyoEJb3akUspHv/1ozv69yS6YuU0NiOQaTVibmPByYFXN2ZFMKT/bm3+ge6ShLqcgu9qifOi8pTTU2OrVJrwcWIlXeRLJFC+9cpKXXjmZFR6pr+PDFyzzySr/qK3XTJGY8Lyjqy1ak1VQE14OrB/POzpr1BtqwptCciLF4Eg49z0II7U61MyEN4WjownURot5RldbbU41MuFNwdp33mIlngGY8LymVke8WHfCFEx43nLrr17gl/sHs/r3ulqjdLZFWdgW5bT5c/w20RVMeFOwznNv6Y2N0BsbyXltSfsc/u+6Szy2yBusqjmFgSHzaAaFaq6GmvCmYCVecKhmx4sJbwqxoTG/TTAcrMSrIcy5EhyqWXg151z5xi96ebh3MOf0la7WJvpNeIHBhFdFPHHgFe7f0184ouE71dzGqznhWVUyPPzs2RivjI5nzd1b0BqhsQrm7tWe8MxrGRq2bD/Alu0HpoV3tER428ouvviH53pvVIUI/6tjlliJF36OjiQ4mQj3rgA1JbzRRJLheNJvM4wKEHbHS00Jz0alVA9hF17VtfE+eMuviCdTObsLbIJr9VATwhORF4AhYAJIqmqPiNwA/AnpbZoBPqOqP3bibwSucuJ/WlXvq7DdOVFVdh06wcnxcNf/jcKEvathNiXe21V1YErYF1X1XzMDRGQVsAFYDZwG3C8iK73YI28kMWGiqxHCXuK50cZbD9yhqnFVfR7YD5zvwnOmYR7L2uHYaILYUJyJVDjX6Si2xFPgpyKiwNdUdbMT/kkR+WNgB/BXqnoMWAw8nJG2zwnLQkSuBq4GWLasMusqmvBqhytvfRSAOoGOlqmTaCNZ7fuFTgf8vDmNgVlKsFjhvUVVD4nIQmCbiDwD3AL8A2lR/gPwBeBjQK6/bNpryRHvZoCenp6KvLZMeLVHStPLMQ4Mx9lzOH/cSH0dna0Rvv2x81m5qM0bA2egqKqmqh5yjv3A3cD5qnpEVSdUNQV8nVPVyT5gaUbyJcChShibSimpPFULm9Jj5CMxkeLQ8TFaov478wtaICItQJ2qDjnn7wI+JyLdqjr5jrkUeNo5vwfYIiI3k3aurAAeqYSxvQPDrPvSL1jQEpm2PkdXa5SHewcr8RijylnQ4v+SgsVIfxFwt1M3bgC2qOq9IvKfInIu6WrkC8AnAFR1l4hsBXYDSeDaSnk0Y0MJJlJK/1Dcpu8YJTG3qYGmxnq/zSgsPFXtBc7JEf6RPGk2AZvKM206NsDZKJegdEOEasiYOU+McjHhlYAJzyiXoOzVbsIzaoqUKr2xYYbGxlEfN8nw3686C2z7LKNcfvzUy/z4qZcBaGqsOzW73elw3/jes2n1oLshVMKzEs+oJGPjKQ4ePcnBo6d2qf3796325NnhqmpaiWe4SHtzI5EGbyQRGuFNpJRBE57hIl46XkIjvKMjCUI6EN0ICV52NYRGeOZYMdzGhJcDc6wYbuPlrPbQeDVNeEa5LGmfwzUXn8XAUILY8BixoXj6Mxyn/0ScTg9LvPAIz6qaRpmc0dnCH11wes5rquqpDyE0whuwEs8ok3xtOBGh3sPJ6eFp41mJZ5RJkFYmC4/wrMQzyiQoMxPAhGfUECa8ErCqplEuQapqhsK5kkimeGV03G8zKs6n167g8jcuod9xaw8Mn3JvZ30fihNPpvw2N/BctmYxv/f67qxugszzhXOb/DbxVUIhvMGR6iztlrTPYWlHM0s7mvPGU1WG4kkGMn5EDz4T467H+zyyNBy87rR5rD17kd9mFEUohFet7bti2xwiwtymRuY2NXJmVysA/Sfi3PW4m9aFjyC14QoRijZe1QqvjDaHtXmnY8KrMFUrvDJ+KNWaJ+Vgwqsw1TgzQSS9l3epVGOelIsJr8JU49u9ozlCY33p2V+NeVIO0YY62gKwNHuxhMLSamzP1NUJP3nq8KnFdtqis1rTvxaF99UPr2FJe/P07oLhOCiB2QmoGMIhvCr8kcWG4lzznWy3ZHOkPmsL6c7W7O2nJs87WiI1ua30a39rLss7W/w2oyKY8ALEaGKCFwdHeXFwNG+8BS0RHvvsO3K++WNDcXYefIXe2IhHVntHmNpwhQiF8D63/nUcPn4y64eWnswY5+XjYzW3/XJna5T5zRHmN0dYkWOft3/+yR6+9lCvD5a5R3OkPhDba1WKUPwlF63smvHalu0H+MzdT3lojf8UevNXYw2hmko7KNKrKSIviMhTIrJTRHY4YR0isk1E9jnH9oz4G0Vkv4jsFZF1bhkP1fkjK0RNCi9AA5wrwWxKvLer6kDG9+uBB1T1RhG53vl+nYisAjYAq0lvTHm/iKys1B55U4kN194usGEV3lW/ewZvWDJv+gDmoTgDwwkGR+LMtJ1BtZV45VQ11wMXO+e3AT8DrnPC71DVOPC8iOwnvU3zr8t41owE9UfmJp2t+Tveg9q5fuGZC3jHqpkHMScnUhwdTUwTZGwozmu7/d2zvNIUKzwFfioiCnxNVTcDiya3YlbVwyKy0Im7GHg4I22fE5aFiFwNXA2wbNmyEs2HgeHac6vne/tPpJSjAe1qKFRqNdTXsbCtiYVtwZm+4xbFCu8tqnrIEdc2EXkmT9xcvZjTKhCOeDcD9PT0lLy+Uy2WeF2tM/8wB0figV1xu9qqi+VQlPBU9ZBz7BeRu0lXHY+ISLdT2nUD/U70PmBpRvIlwKEK2pxpV20KL88PuDnSwM1/cM60Pr7JSbXHfJxQvKBAFbmWKCg8EWkB6lR1yDl/F/A54B7gCuBG5/gDJ8k9wBYRuZm0c2UF8IgLtjORUi5bs3jazO2x8eqerZ1PeK3RBi5bs2TG64lkisGR6TPcp/aPxobiDMeTFbN5fnMj0Yb6it0v7BRT4i0C7nbGwTUAW1T1XhF5FNgqIlcBB4DLAVR1l4hsBXYDSeBatzyaDfV1bLr09VlhqspwPPlqozz9GTv1PaMEOHw8fB7R+jph/pzGktNHGuronjeH7nlzCsYdTSSzhJhrSYrJ8ESBpSmqrTugXAoKT1V7gXNyhA8Ca2dIswnYVLZ1JSAitDU10tbUyBl5xvWpKqv+7r7QjXrpbI1QV+fNYODmSAPLFjSwbEFxS1NkeyOzhXna/MJCryVCMXLFDUYSE6ETHQTTQZG5NMVZztIURn5CMR/PDcK6JLxV2aqDmi3xgjrHr7Fe+MRFZ00b3TEwHCeZUk93LTXco3aFF9ASb2FbE3+97jXTwlMp5fjJcSZmGlNlhAoTXsCYaY+2ujqhvYw1WoxgUbttvIBWNbusk7kmqNkSb1X3XD6wZonTaZxuTw0O+z/cKoheS6Py1Kzw3vP6bt7z+u6ssMkBxlP7oAamODoOvXKSkYQ7XRHmtawNalZ4uaivk1cXFDq7e+Z4N/90L1/53/2u2GAlXm1gwiuBcroizlk6nw+uWeyUponsUnUobt0FNYIJrwTK8Yi+7rS5fOTC5TmvqarvbUzDG0x4JRArY/JtvqqkiFAfnjVZjTKo2e6EcihnuJm14Qww4c2aciffmtfSABPerDlxMkliovSJtlbiGWBtvFmjKFe+efmpSaGON3KoyNnaJjwDTHizZn5zhBvev3pa+Nj4xAzrRZ7qiO8/Yd0FRhoTXoVoaqxnaUczSzvyz9Y2DLA2nmH4ggnPMHzAhGcYPmDCMwwfMOEZhg+Y8AzDB0x4huEDJjzD8AETnmH4gGgA1mkUkRjwYpHRO4GBgrG8x+yaHWG163RV7Sr3IYEQ3mwQkR2q2uO3HVMxu2ZHrdtlVU3D8AETnmH4QBiFt9lvA2bA7JodNW1X6Np4hlENhLHEM4zQY8IzDD9QVU8+wLeAfuDpjLAOYBuwzzm2Z1zbCOwH9gLrMsI3AQeB4Sn3vxKIATudz8czrl3hPGMfcEWl7QKagR8BzwC7gBsz4keB7zpptgPLA2KXb/nlhN8L/Max66tAvd/5VcCukvJrRj14KLyLgDVTMuYm4Hrn/Hrg8875KuePjwJnAM9lZMCbgG5yC+/fczy3A+h1ju3OeXsl7SL9A3+7EycC/AJ4j/P9T4GvOucbgO8GxC7f8su5Ntc5CnAXsMHv/CpgV0n5NdPHs6qmqv4cODoleD1wm3N+G/D7GeF3qGpcVZ8n/WY637nPw6p6eBaPXgdsU9WjqnqM9Jvv3ZW0S1VHVfVB534J4HFgSY573QmsFREJgF2+5ZdznxNOnAbSLwXNcS9P86uAXTOR166Z8LuNt2hSRM5xoRO+mHR1cpI+J6wQHxCRJ0XkThFZWsa9SrZLROYD7wMemJpGVZPAcWBBAOwCn/NLRO4jXT0cIi2yrDR+5dcMdkHl8st34c1Erh0ECr15fki6PfAG4H5OvelKuVdJdolIA3A78BVV7S2Qxm+7fM8vVV1HutkQBS4pkMZvuyqaX34L74iIdAM4x34nvA9YmhFvCXAo341UdVBVJ9dW/zrwxlLvVYZdm4F9qvqljLBX0zgCmEe6SuSrXQHJL1R1DLiHdPUvK42P+TXNrgrnl3fOFachupzsxu+/kN34vck5X01247cXp/GbkXaqc6U74/xS4OGMxu/zpBu+7c55R6XtAv6RdGO8bsq9ryXbWbA1IHb5ll9A6+TzSbelvgt80u/8KmBXyfmVUwseiu524DAwTvotcRXpuvsDpN2wD2QaDPwNaW/TXhxPnBN+k5M+5RxvcML/mbQL+DfAg8BrM9J8jHQDej/w0UrbRfotp8AepribgSbge86zHwHODIhdfubXIuBR4EnHhn8DGgKQX/nsKim/ZvrYkDHD8AG/23iGUZOY8AzDB0x4huEDJjzD8AETnmH4gAnPMHzAhGcYPvD/fB2UCIJFwfMAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "inventory.loc[~inventory.is_valid].iloc[[0]].plot();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's correct them:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "inventory.loc[~inventory.is_valid, 'geometry'] = inventory.loc[~inventory.is_valid].buffer(0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A final preparation step is to convert the format of the inventory to a file which resembles the RGI (see [use_your_own_inventory.ipynb](use_your_own_inventory.ipynb)):" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# We keep the original ID for later reference\n", "gdf = utils.cook_rgidf(inventory, o1_region='08', assign_column_values={'breID':'breID'})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute the centerlines" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We use the standard OGGM procedure for this:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2022-01-05 21:19:34: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.\n", "2022-01-05 21:19:34: oggm.cfg: Multiprocessing switched OFF according to the parameter file.\n", "2022-01-05 21:19:34: oggm.cfg: Multiprocessing: using all available processors (N=32)\n", "100% of 22.3 MiB |######################| Elapsed Time: 0:00:02 Time: 0:00:02\n", "2022-01-05 21:19:38: oggm.cfg: PARAMS['use_rgi_area'] changed from `True` to `False`.\n", "2022-01-05 21:19:38: oggm.cfg: PARAMS['use_intersects'] changed from `True` to `False`.\n", "2022-01-05 21:19:38: oggm.cfg: PARAMS['border'] changed from `40` to `10`.\n", "2022-01-05 21:19:38: oggm.cfg: PARAMS['d1'] changed from `14.0` to `7`.\n", "2022-01-05 21:19:38: oggm.cfg: PARAMS['d2'] changed from `10.0` to `5`.\n", "2022-01-05 21:19:38: oggm.cfg: PARAMS['dmax'] changed from `200.0` to `100`.\n" ] } ], "source": [ "cfg.initialize(logging_level='WARNING')\n", "\n", "# Parameters\n", "cfg.PARAMS['use_multiprocessing'] = False # True is often a good idea but not here \n", "cfg.PARAMS['continue_on_error'] = True # Yeah, shit happens\n", "cfg.PARAMS['use_rgi_area'] = False # this is required for user-defined inventories\n", "cfg.PARAMS['use_intersects'] = False # we don't care about intersects for centerlines\n", "cfg.PARAMS['border'] = 10 # no need to make a large map\n", "\n", "# Optional: change the grid resolution\n", "# E.g. fixed grid spacing\n", "# cfg.PARAMS['grid_dx_method'] = 'fixed'\n", "# cfg.PARAMS['fixed_dx'] = 10\n", "# Or variable but twice higher than default \n", "cfg.PARAMS['grid_dx_method'] = 'square'\n", "cfg.PARAMS['d1'] = 7 # (default is 14)\n", "cfg.PARAMS['d2'] = 5 # (default is 10)\n", "cfg.PARAMS['dmax'] = 100 # (default is 200)\n", "\n", "# Tell OGGM to use our user DEM (important!)\n", "cfg.PATHS['dem_file'] = fpath_dem\n", "\n", "# Where to work\n", "cfg.PATHS['working_dir'] = 'wd_2xdefres'" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now the workflow: " ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2022-01-05 22:10:42: oggm.workflow: Execute entity tasks [glacier_masks] on 6744 glaciers\n", "2022-01-05 22:10:42: oggm.workflow: WARNING: you are trying to run an entity task on 6744 glaciers with multiprocessing turned off. OGGM will run faster with multiprocessing turned on.\n", "2022-01-05 22:15:32: oggm.core.gis: InvalidDEMError occurred during task glacier_masks on RGI60-08.01696: (RGI60-08.01696) min equal max in the masked DEM.\n", "2022-01-05 22:32:52: oggm.workflow: Execute entity tasks [compute_centerlines] on 6744 glaciers\n", "2022-01-05 22:32:52: oggm.workflow: WARNING: you are trying to run an entity task on 6744 glaciers with multiprocessing turned off. OGGM will run faster with multiprocessing turned on.\n" ] } ], "source": [ "# gdirs = workflow.init_glacier_directories(gdf)\n", "\n", "# workflow.execute_entity_task(tasks.define_glacier_region, gdirs, source='USER'); # Use the user DEM\n", "\n", "workflow.execute_entity_task(tasks.glacier_masks, gdirs);\n", "workflow.execute_entity_task(tasks.compute_centerlines, gdirs);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Note: the default in OGGM is to use a grid size of varying resolution for each glacier. I think it makes sense in many cases, but you may prefer to use the native resolution of your DEM. You can do so by commenting / un-commenting the options above.**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Write the data to a shapefile: optional smoothing" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The relevant task is \"write_centerlines_to_shape\", which writes everything to a shapefile:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2022-01-05 22:46:43: oggm.utils: Applying global task write_centerlines_to_shape on 6744 glaciers\n", "2022-01-05 22:46:43: oggm.utils: write_centerlines_to_shape on outputs/Norway_Centerlines_default.shp ...\n", "2022-01-05 22:46:43: oggm.workflow: Execute entity tasks [get_centerline_lonlat] on 6744 glaciers\n", "2022-01-05 22:46:43: oggm.workflow: WARNING: you are trying to run an entity task on 6744 glaciers with multiprocessing turned off. OGGM will run faster with multiprocessing turned on.\n" ] } ], "source": [ "from oggm.utils import write_centerlines_to_shape\n", "\n", "write_centerlines_to_shape(gdirs, # The glaciers to process\n", " path='outputs/Norway_Centerlines_default.shp', # The output file\n", " to_tar=True, # set to True to put everything into one single tar file\n", " to_crs=inventory.crs, # Write into the projection of the original inventory\n", " keep_main_only=True, # Write only the main flowline and discard the tributaries\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's have a look at the output:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
RGIIDSEGMENT_IDLE_SEGMENTMAINgeometrybreID
0RGI60-08.000010322.01LINESTRING (118270.857 6914371.200, 118278.618...1855
1RGI60-08.000020358.01LINESTRING (118573.658 6913783.719, 118579.630...3156
2RGI60-08.000030314.01LINESTRING (117411.558 6913657.999, 117412.239...3157
3RGI60-08.000040332.01LINESTRING (117804.659 6913611.777, 117809.950...1859
4RGI60-08.000050215.01LINESTRING (138401.165 6909220.080, 138407.784...3158
\n", "
" ], "text/plain": [ " RGIID SEGMENT_ID LE_SEGMENT MAIN \\\n", "0 RGI60-08.00001 0 322.0 1 \n", "1 RGI60-08.00002 0 358.0 1 \n", "2 RGI60-08.00003 0 314.0 1 \n", "3 RGI60-08.00004 0 332.0 1 \n", "4 RGI60-08.00005 0 215.0 1 \n", "\n", " geometry breID \n", "0 LINESTRING (118270.857 6914371.200, 118278.618... 1855 \n", "1 LINESTRING (118573.658 6913783.719, 118579.630... 3156 \n", "2 LINESTRING (117411.558 6913657.999, 117412.239... 3157 \n", "3 LINESTRING (117804.659 6913611.777, 117809.950... 1859 \n", "4 LINESTRING (138401.165 6909220.080, 138407.784... 3158 " ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cls_default = gpd.read_file('tar://outputs/Norway_Centerlines_default.tar.gz/Norway_Centerlines_default.shp')\n", "cls_default['breID'] = gdf['breID'] # This only works this way because we have one centerline per glacier!\n", "\n", "cls_default.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`LE_SEGMENT` is the length of the centerline in meters. The RGI \"IDs\" are fake (OGGM needs them) but the breID are real. Lets use them as index for the file:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "cls_default = cls_default.set_index('breID')\n", "orig_inventory = inventory.set_index('breID')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can plot an example:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS4AAAEDCAYAAACLcumrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAvK0lEQVR4nO3de3RV9Z338ff33HK/J4iExHAJKNaKErzXC4gCnWpttdV2HK2tiDiu4uN0Ou0z07X6zKxnxj6OD856ahmq1dra0mLVtlZuKniroAFR7oa73K8hgZDbOd/nj70Dh3BCcuAk5+yT72utLHZ+e/9Ofj+SfPLbv73P/omqYowxXuJLdgOMMSZeFlzGGM+x4DLGeI4FlzHGcyy4jDGeY8FljPGctAkuEfmFiOwVkVU9PP5rIrJGRFaLyG96u33GmMSRdLmPS0SuBY4Az6vq57o5thr4PTBOVQ+JyABV3dsX7TTGnL20GXGp6tvAwegyERkmIvNEZJmIvCMi57u77gd+qqqH3LoWWsZ4SNoEVxdmAQ+r6hjgH4Cn3PIRwAgReU9ElojIxKS10BgTt0CyG9BbRCQXuAqYIyIdxRnuvwGgGrgeGAy8IyKfU9X6Pm6mMeYMpG1w4Ywm61V1dIx924ElqtoGbBaR9ThB9mEfts8Yc4bS9lRRVRtwQukOAHFc7O5+BbjBLS/FOXXclIx2GmPilzbBJSK/Bd4HRorIdhH5NvBN4Nsi8jGwGrjVPXw+cEBE1gCLgO+p6oFktNsYE7+0uR3CGNN/pM2IyxjTf6TF5HxpaalWVVUluxnGmARbtmzZflUt61zebXCJyEjgd1FFQ4EfqeqMGMeOBZYAX1fVF0WkAngeGAhEgFmq+qR77MXATCAX2AJ8U1UbRGQC8B9ACGjFmX9683RtrKqqora2truuGGM8RkS2xirv9lRRVder6mj3toIxQBPwcowv4Acew5n47tAOPKqqFwBXAA+JyCh339PAP6nqRe7rfc8t3w98yS2/B/hV990zxvQn8c5xjQc2qmqsFHwY+ANw/O0zqrpLVZe7243AWqDc3T0SeNvdXgh81T3uI1Xd6ZavBjJFpOPGUWOMiTu47gR+27lQRMqB23BO/WISkSrgEmCpW7QKuMXdvgOoiFHtq8BHqtoSZzuNMWmsx8ElIiGcoJkTY/cM4PuqGu6ibi7OaGy6e2MowH04p47LgDyc+azoOhfinHo+0MVrThGRWhGp3bdvX0+7YYxJA/FcVZwELFfVPTH21QCz3fcElgKTRaRdVV8RkSBOaL2gqi91VFDVdcBNACIyAvhixz4RGYwz7/V3qroxVmNUdRbOm6ipqamxm9GM6UfiCa67iHGaCKCqQzq2ReQ54FU3tAR4Blirqk9E1+l4BpaI+IB/xj3NFJFC4C/AD1T1vTjaZ4zpJ3p0qigi2cAE4KWosqkiMrWbqlcDdwPjRGSF+zHZ3XeXiHwKrAN2As+65X8PDAf+JarOgJ53yRiT7tLiLT81NTVq93EZk35EZJmq1nQuT4s7542p3VXP0baY14Z6bERxDufmZiaoRaY3WXCZtHCouY3G1vazeo2WcFaCWmN6m73J2hjjOTbiMp6yen8j2w43nVLe3B4569dOg+nefsOCy3hKWzjCsQSEVCzhiCWXV1hwmaSIpODwpl17JxBN4llwmaR4beMeWsOpFV424vIOm5w3xtWegqNAE5uNuEyfaWoLE3bDIRUzwkZc3mHBZfrM0p2HONTcluxmdMmCyzvsVNEY4zk24jK9qr65jZawc7WuLZziV+2cxzIZD7DgMr1q/cEj7GhsTnYzesRiyzvsVNEY4zkWXMa4bMTlHRZcxnSw5PIMCy5jXJZb3mGT8ybhdh9pZl+Ts2jT4ZbUvW+rM7Ho8gwLLpNw+4+1UnfoaLKbET/LLc+wU0WPUlXCESUd1gxIFZZb3tHtiEtERgK/iyoaCvxIVWfEOHYssAT4uqq+KCIVwPPAQCACzFLVJ91jL8ZZkiwX2AJ8U1UbRKQEeBEYCzynqn9/5t1Lfe2RCHuOttAWUdojSnsk4v578uex9iswIDtE0Of8/ak5txC/Lzm/fnuOtrCl3nnA3+FW75weGm/qNrhUdT0wGkBE/MAOnMVaT+LuewyYH1XcDjyqqstFJA9YJiILVXUN8DTwD6r6lojcB3wP+Beg2f33c+5HWtt1pIUPd9Wfcf29TScWAD9lKZQ+1NQWZscRb9xo2hUbcXlHvKeK44GNqro1xr6HcVas3ttRoKq7VHW5u90IrAXK3d0jgbfd7YXAV93jjqrquzgBlpaa28PsPtrMugONrDvQmOzmGJfYW348I97J+TuJsZq1iJQDtwHjcE7xTiEiVcAlwFK3aBVwC/BH4A6gIp6GiMgUYApAZWVlPFV7zab6o8dPl7rS0ouPHo7H5vomNtcnbgK9JdXfh2jSSo+DS0RCOEHzgxi7ZwDfV9VwrL9aIpKLMxqbrqoNbvF9wH+JyI+APwGtp1Q8DVWdBcwCZ0HYeOr2lvaIUt9ydktk9ZWWcNgzbe0rNt7yjnhGXJOA5aq6J8a+GmC2G1qlwGQRaVfVV0QkiBNaL6jqSx0VVHUdcBOAiIwAvniGfUi4NfsbaW4PE1bn6p1PhIBP8ItQnBWkPC/2+nsdk+TGmN4VT3DdRYzTRABVHdKxLSLPAa+6oSXAM8BaVX0iuo6IDFDVvSLiA/4Z5wpj0vx1+0EOHHMGfW2neaBcZTiry+BK9unSaxtj/U2JLWy3UZzCpri8o0dDBBHJBiYAL0WVTRWRqd1UvRq4GxgnIivcj8nuvrtE5FNgHbATeDbqtbcATwD3ish2ERnV0w6dqY5bDk4XWuBc9m9uP3Wp99ZwhLqDR3qreT3SFtWH7j7sYZ/Gy3o04lLVJqCkU1nMEZKq3hu1/S5dTB2493M92cW+qp6062y1R/T4MlnOXVHdawlHWLztAMMLcwj4BVU40trO3qaWbkPPpLau3vKjqnbFMcX067f8fLz3MFsPH4u7XlNbmE/2NXR/oPGUWNm060gzBRkBsoP9+lcl5dhssjFdONYeZtnuw8luhonBgssY187G5hNTB6os21VPq92flpIsuIxxNbS2U3fwKK3hCJ8ePHrS26lMarETd2OirN7fyOr99jasVGcjLmOM51hwGWM8x4LLGOM5NsdlTAzleZlkB/wABOw9qCmn3wfXjgNNbNpzBJ8IIs4zmTpuRPSJcy91dLkAmSE/I8vz8dnd1GlraGE2ZdkZyW6G6UK/D64Nu4/wxw+2x13v8hEl3HFVpYWXMUnQr4OrIi+LB64Zwj2XV7KxvomjrWEUJeLec6goqqBK1I2JULvxAG98sgefCLdfWWHvYzOmj/Xr4BqQk8GAHOd04Eg4wv5jPbvhcPKlg1CFN1fuwSfwlSssvIzpS/06uM6UiPDFMYNQVRat2otPhC9fPtjCy5g+YsF1hkSEv6kpJ6Lw1uq9+HzCLWPLLbw8bFhRNue4E/L5GcEkt8acjgXXWRBxwiocUd5avRcR+FKNhZdXFWQEGZibmexmmB6w4DpLIsJtlw9GVVnsnjZ+ccwgCy9jepEFVwKICLddUUE44kzY1x9tpTg3dDy8Ou4D8xF1P5hAZtBPzfASQgG7wdGYeFhwJYhPhNuvqsDvE5bW7ScScW+l6Kbeyq313Dd+GEELL2N6zIIrgXwifPXKCr565clr26qeCLHo7Y82HeR3723j2Tc3cd/4oQT8Fl7G9ES3vykiMjJqhZ4VItIgItO7OHasiIRF5Hb38woRWSQia0VktYh8N+rYi0XkfRFZKSJ/FpH8qH0/EJENIrJeRG5OQD+TSkTw+QS/Twj4fQQDPkIBH5ePKOWOqypZt6OB597cRLs9bbPPjSjO4dbqgdxaPZDz8mMvO2dST7fBparrVXW0qo4GxgBNwMudjxMRP/AYMD+quB14VFUvAK4AHopaauxp4J9U9SL39b7nvs4o4E7gQmAi8JT72mnpypGl3H5lBWu2N/D84s0WXn1McP6g+H1iF1Q8JN5zk/HARlXdGmPfwzgrVu/tKFDVXaq63N1uBNYC5e7ukcDb7vZC4Kvu9q3AbFVtUdXNwAbgsjjb6SlXnV/GV66oYNW2w/xq8WbCtsyZMacVb3DdSYzVrEWkHLiN06xGLSJVwCXAUrdoFXCLu30H0DExVA58FlV1OyfCLvr1pohIrYjU7tu3L75epKBrLijjy5cPZuW2w/z6LQsvY06nx8ElIiGcoJkTY/cM4PuqeuoSz07dXJzR2HRV7ViQ8D6cU8dlQB7Q8UbBWOP1U36LVXWWqtaoak1ZWVlPu5HSrh01gFsvK+fjLfX85u0tFl7GdCGeq4qTgOWquifGvhpgtjtHUApMFpF2VX1FRII4ofWCqr7UUUFV1wE3AYjICOCL7q7tnBh9AQwGdsbRTk+77sJzCEfg1dodiMA3vlCFz2dzL4kU8vvICznTptnBtJ0+TWvxBNddxDhNBFDVIR3bIvIc8KobWgI8A6xV1Sei64jIAFXdKyI+4J85cZr5J+A3IvIEMAioBj6Io52eN+6ic4io8tqynfhEuPOa8yy8Euic7BBjBxUluxnmLPQouEQkG5gAPBBVNhVAVbuc1wKuBu4GVorICrfsh6r6GnCXiDzklr0EPOu+3moR+T2wBueq5ENdnYKmsxs/P5BIRJn30S58PuFrV9tDC43pIKren0epqanR2trauOsdONbKkdZ2AD49eITG1tTLx3kf7WTBit32xNUEqsjL9PyI6+CxVhrdn92uVORnef7nRUSWqWpN5/J+fef8lsNNbD18LNnNOK2bR59LJKK8/ske/O6d+Xa/kdnWcIxN9U2nPWZwXlbsS11poF8HlxeICJMuHUTEfeKqiPCVK+yhhaZ/s+DygI4nrkYiyuLVe/H54MuXWXjFIzPgo6ogG4CCDO/82De2tLO98dSzgkPNbd3WXXfgCJ2v6WQF/FQVZieqeUnjne9gPycifGms88TVt9c4z/2yJ672XHbAz6jSvGQ3I26Nre2sPXDkjOquP3hqveLMoAWX6Vsiwq2XlRNR54mrPhH+psYeWmj6n34XXDuPNLOzsRlwrip6TfQTVxet2oPfB5MutfDqa83tYVbta+z1r9PUltgr3UfawtTuqj+l/NzcDMrzvPN0jH4XXA0tbWxrSO0rid2JfuLq65/sOT6Bb/pOe0Q9+XPUGo7EbHdW0E+5h86k+11wpYuOJ66qKgs/3o3PJ9w8+txkN8uYPmHB5WE+Ee64upKIKvM/2oVPYMLFFl69RVV5+7ODAITT4MbtaNsON7G/yZk6GVGcw7kpvtqRBZfH+UT4+tXnEVGYu3wXPhHGf35gspuVtrw4L9oTx9ojHGt3+tYSTv25LguuNODzCXddcx4RVf6ybCd+n3D9585JdrOM6TUWXGnC5xO+8YUqNKL86cMd+ES49sIByW5WUoX8PsadVwpwyo2Ypmsr9zawdr9zxfSyQUWUZIWS3KJTWXClEb9P+OZ1Q4joZl75YDsi8IVR/Te8BHve1ploiyht7kMsIyk6l2frYaUZv0+4+/ohfK6ygJeXbue9dd5/rLVJHlUnvCKq/PaDbWzY2/v3rvWEjbjSkN8n/N31Q/jlos384f3P8Ilw5cjSZDfLeNC7252rqG+v3ssrH2znrssq+PevfD7JrbIRV9oK+H3cc8MQLhicz5y/bmPpp/uT3STjUYtX7eGVD7YzYdQAfnzL55LdHMCCK60F/D7uvWEoI8vz+f172/iw7kCym9QncoJ+coN+cmx+66y9uXI3f/pwBxdXFfL41y4mFEiNyLBTxTQXDPj41rih/OKNjcx+dysiUDO8JNnN6lU3VpXht8uIZ+31j3fz2vKdXDKkiG9cW0XQnxqhBTbi6hdCAR/3jR/G8HPz+O27W1m+6WCym2RS3MIVu3ht+U4uHeqEVqr9Ieg2uERkpIisiPpoEJHpXRw7VkTCInK7+3mFiCwSkbUislpEvht17GgRWeK+Zq2IXOaWh0TkWRFZKSIfi8j1CelpP+eE11CGnpPLC29vYcXmQ8lukmfsa2ph99Fm9hxtSXZTep2qMu+jncz9aBc1w4r5xhdSL7SgB6eKqroeGA0gIn5gB/By5+PcfY8B86OK24FHVXW5iOQBy0RkoaquAX4C/FhV54rIZPfz64H73a97kYgMAOaKyFhVjZx5Nw1ARtDPd24cxs8XbuTXb21GBC6u8vaiEX2hdtdhjrWn3kIqiXK4qZWdB4/hE2H9jgYWr97LZdUlfO2qypRdFi/eOa7xwEZV3Rpj38M4C7+O7ShQ1V3ALne7UUTWAuU4S48pkO8eWsCJRV9HAW+4dfaKSD3OgrP9am3F3pIR9POdCcOYtWADv1q8GblB+Px5hclulkmSz/YfZeb8DRyLWuHKCytKxRtcdxJjUVgRKQduA8YRFVydjqkCLgGWukXTgfki8jjOKetVbvnHwK0iMhtnResx7r8WXAmSGfQzZcJwZi6o4/lFm7h33FA+V1mY7GadlfK8TEI+Z+YjhX/fUsrWfUf57wUbyA75ueeGIQT9PgJ+H4NLslL+wZQ9Di4RCQG3AD+IsXsG8H1VDcfqsIjk4ozGpqtqg1v8IPCIqv5BRL6Gs+L1jcAvgAuAWmAr8FecU87OrzkFmAJQWVnZ024YV2bIzwM3VTNzfh2/XLSZb40byqiKgmQ364xdUJJLfkYw2c1ISR9tOshf1zv38flEjgf71n1Hyc0M8ODEaopzM5LYwvjFc1VxErBcVffE2FcDzBaRLcDtwFMi8mUAEQnihNYLqvpSVJ17cFawBpgDXAagqu2q+oiqjlbVW4FCoK7zF1TVWapao6o1ZWVlcXTDdMgK+XngpuGcW5TJs29uYt2Ohu4rGU9Z8ul+fv3WFhqa2lCF9nCE1vYIrW0Rhg3MZdrEEZ4LLYjvVPEuYpwmAqjqkI5tEXkOeFVVXxFn+PUMsFZVn+hUbSdwHbAY5xSzzq2fjbPC9lERmQC0u5P5phdkZwR44OZqZs6r4xdvbOTb44cxsjy/+4opoDI/i/yQ8yOc4U/szaab6o9y1J33aYt487rQX9ft48X3P+P88nzuHTc0ZW4eTYQeBZcbJhOAB6LKpgKo6szTVL0auBtYKSIr3LIfquprOFcPnxSRANCMe9oHDMCZ+4rgXMG8u8e9MWckxw2vn82r45k3NnL/jcOpHpT6DyAvz8vstSd1bm9oZr+HHxr43tp9/GHJZ4wanM89NwwlmEahBT0MLlVtAko6lcUMLFW9N2r7XbpYBNzdNyZG+RZgZE/aZRLHmesYzlNz63j69Q3cf9Nwhg9M/fAyp3pnzV5eXrqdCysKuOeGIQRS6I73REm/HnVjUG4mlw8q5PJBhZRlp94D0pIpNzPoTNTmZfD0wo1s2n1mC5Ga5HlrtRNaF1Wmb2hBPwyu/Iwg5XlZlOdl2UPmYsjLCvLgzdUU5gT5+cINbN5j4eUVi1bu4Y8fbOfiqkL+7oahaRta0A+Dy3QvP9sZeeVnB5m1cANb9h5NdpNMN974ZDd/rt3B6CFF/O11Q1LybTqJJJqij2aNR01NjdbW1sZd72hbO61hp/8f7T5MfUtbopvmafVHW/np3DqONrcx9eZqKstykt2kk1xZXpTQyfm1+xvZdcRZ5byxNeyZJcgWfryLuct3cenQIu7qxfcW5ocC7rP7hXFVffNgShFZpqo1ncv79YgrJxigKDNIUWaQQJr/hToThTkhpk2sJjsjwH8v2MBn+5uS3aRe1dQWpr6lnfqWdk+E1vE3RC/fxZg+eEN0Q2u7+/+T/D/w/Tq4TPeKckNMm1RNZsjPzPl1bD+Q3uHlFU5o7WLBit2MHV7MXdecl7JviO4NFlymW8W5GUybWE1m0AmvHQdTI7w+2FnPn+t28+e63TSkwCigr6gqry3fycKPd3N5dQlf72ehBRZcpodK8jJ4cGI1Ib+PmfM2sOvQsWQ3ibDqSUtp9Qeqyqu1O3jjkz1cMaKUO65O7ac49BYLLtNjpfkZPDipGr9f+Nm8OnbXJz+8OrRFlNZwhNZwhHguOIWj6qXqGoIdVJ3Ffhet2svV55dy+1UVSQutjv+ztnBy3g7Vr68qRnt72wFPv8WjL+093MxP534KwLSJIzinsHfednOmbq0e2ONJ6nUHGlmzP/XvVVNVXlm6nXfW7uOaC8q47fLBKfHomZygn5uH9t6iw3ZV0STMgIJMpk0cAcCzb25Mcmv6h5fd0LruwgEpE1rJZKv8uHJDftoizn9HQ0s73h+H9q5zCjO5dGgx769PvfUa61va8PfwF/tYe+o/+eFAYwvvrt3HVeeXcsvY8n4fWmDBddylAwuPb8/buJemNH7GeCKl4u/QW9vSa/3Ijnmk4QPzLLRcdqpozlh+VpCWtghvrd6b7KaktTSYhk44G3HFcE5OBi3hCKDsPJL+S1KdqWsvHMC2/Uf54wfb8Ql8YVTvTdL2Zx0X0PrZrVqnZcEVwyUDnWevqyovf7o7ya1JXX6f8LfXDSES2cTLS7fjE+HqC+wx2ol2fMRlp4nH2amiOSt+n3D39UO4sKKAPyz5LCUn673ueG4ltRWpxYLLnLWA38c9Nwxh1OB85vx1G0s+tfBKpI4Rl50qnmDBZRLCCa+hnF+ez5z3tvFhXXpd2Uum4zeJ26nicRZcJmGCAR/fGjeUEYPymP3uVmo3WHglgp0qnqrb4BKRkSKyIuqjQUSmd3HsWBEJi8jt7ucVIrJIRNaKyGoR+W7UsaNFZIn7mrUicplbHhSRX4rISrderAVoTYoKBnx8a/wwhp+bx2/f3cqyjQeT3STP6xhxpeKAqyUc4aM9h/loz2E2HOq7J+V2G1yqut5dnHU0zqo8TcDLnY8TET/wGDA/qrgdeFRVLwCuAB4SkVHuvp8AP3Zf90fu5wB3ABmqepH79R4Qkar4u2aSJRTw8e0bhzH0nFx+884WPtpk4XU2th9w3swe8KXeCVJ7RNlc38Tm+iZ2u0+P7Qvx/k+MBzaq6tYY+x7GWbH6+N2IqrpLVZe7243AWqC8YzfQsfJoAc4CsR3lOe56i1lAK2BLLHtMKODjOzcOY8iAXF54ewsfbzmU7CZ50ocbDvDyks8YNjCXIeek1qOzkyne4LqTGKtZi0g5cBvQ5eKw7qjpEmCpWzQd+D8i8hnwONBxSvgicBTYBWwDHlfVU/5ki8gU9xSzdt++fXF2o+euKi/iqvIiLj2noNe+RrrKCPr5zoRhVJbl8KvFm1m5tT7ZTfKUZRsPMvudrVQPyuP+CcPTetWeePX4f0JEQsAtwJwYu2cA31fVmG/wE5FcnNHYdFXtGD09CDyiqhXAI8AzbvllQBgYBAwBHhWRoZ1fU1VnqWqNqtaUlfXOTY8iwsDcTAbmZlJqazCekcygnykThlNRmsPzizezelt9spvkGW+u3M2g4izuGz+MUJqtRH224vnfmAQsV9U9MfbVALNFZAtwO/CUiHwZnMl2nNB6QVVfiqpzD9Dx+RycwAL4BjBPVdtUdS/wnvv6xqMyQ36m3DScQcVZPLdoM2s+O5zsJnlCOKKU5WdYaMUQz//IXcQ4TQRQ1SGqWqWqVTinetNU9RVx3sr+DLBWVZ/oVG0ncJ27PQ6oc7e3AePEkYMzqb8ujnb2iuygnwlVZUyoKuPSgXbaGK+skJ+pNw3n3KIsnn1zE2u3W3h1R5XUvJTYhQPHWlmweS8LNu9lzf7GXv1aPQouEckGJnBihISITBWRqd1UvRq4GyeIOm6nmOzuux/4TxH5GPjfwBS3/KdALrAK+BB4VlU/6WmHeotPhLyMAHkZAbIDtgL2mcjKCPDAzcMZWJjJs29uYv0Ou+ZyOqreuls+rHCkNcyR1jDNvfycsx69yVpVm4CSTmUxJ+JV9d6o7Xfp4r45d9+YGOVHcG6JMGkoJyPA1JureWpeHc+8sZHv3DiMEYPyu6/YDylqz9/qgp08mz6XkxngwYnVlOZl8MzrG9mwq3dPK7xK1e6W74oFl0mKXDe8ivMyePr1jWzcbeHVmaqnprj6lAWXSZq8rCAP3lxNYU6Qny/cyOY9qb/aTl9S9e6pYjgS4VhbmGNtYdoiiZ/vsuAySZWfHWTaxBEUZAeZtWADW/ZaeHXw8hObP2tsZu6mvczdtJdNhxK/8rkFl0m6/OwgD06sJi/LCa+t+/ruzbqpzGtXFfuSBZdJCYU5IaZNqiYnM8B/L9jAZ/v7d3ht399EU2s7QXubT0z2v2JSRmFOiGkTR5Ad8jNz/ga270/8KYYXbNt3lJ/NryMvK8h1F9oCJLFYcJmUUpQbYtrEajKDfmYuqGPHwf4VXlv3HWXmgg1khfw8NMm56mpOZcF1BjICPirzs6jMz6IoM5js5qSd4rwMpk2qJuT3MXNeHTsPHkt2k/rE5j1HmDm/jpyMAA9NGkFxbnqEVn1LG1sPN7H1cBNH29oT8poWXGegICNIzbmF1JxbSFVBdrKbk5ZK8jJ4cFI1fr+PmfPr2H0ovcNr054jzFqwgfysIA9NqqYoN32eRrKjsZlluw+zbPdhDh1rS8hrWnCZlFWWn8m0idX4RHhqXh176tMzvDbsbmTWgg0U5ASZNqmawpz0Ca3eYsF1lgozA5xfksv5JbmUZNlpY6INKMjkwYnViMBT8+rYe7jvHg/cF+p2NvLzBRsoci9MFKThc9/KskPHf0fyMhKzBrUF11kqygwxqjSPUaV5lGSl3w9dKjin0AkvVSe89jWkR3h9urOBn7++gRJ3Ti8/Oz3/8A3Izjj+O1KQkZg+WnAZTxhYmMWDE6sJh5Wfza1jf0NLspt0VtbtaODp1zdSlu+EVp6N1uNiwWU849yiLKZOHE5rOMJT8z7lQKM3w2vNZ4f5xRsb3dPgEeSm4ZXp0qwQYwYWMGZgAQN74eqoBZfxlPLibKbeXE1LW4Sn5tVx8Ii3wmv1tnqefXMTA93T39zMxMz5pJrcUIDzCrI5ryA7YaeH0Sy4jOcMLnHC61hrmKfm1nHoSGuym9Qjq7bW89yizQwqzmLqzdXkJGiiuj+y4EqgIQXZXFtRzLUVxZTaRH2vqijN5oGbhtPU0s5T8+qoP5ra4fXJlkM8t2gT5SVZTL1pONlpGFqFGcHjP//Vxb27BqQFVwLlhAKUZmdQmp1Bhq3M0uvOK8thyk3VHGlu42fz6mhoSszNjYm2YvMhnl+8mcrSHKbeVE1WGoYWQMgvx3/+80K920f77TKeVjUghykThnO4qY2n5n1K/dFW2sMRwhEloopqcp9qtXzTQX791mbOG5DDlJuHkxmyhVYSQbr7xorISOB3UUVDgR+p6owYx44FlgBfV9UXRaQCeB4YCESAWar6pHvsaJyVrzOBdpwlzT4QkW8C34t62c8Dl6rqiq7aWFNTo7W1tafvaR9rcX95AN7fcZDDLYl5j5aJbdPuI8xauIHWGKvLCM4jkEXE+ZeobZET+6PL6dgffYzgkxOPUz6pbtQxx+shbNt/lCEDcrl/wjAygukXWtkBP9dWOuvo+AUyErwClogsU9VT1lXtdjynquuB0e6L+IEdwMsxvoAfeAyYH1XcDjyqqstFJA9YJiILVXUN8BPgx6o6112y7CfA9ar6AvCC+5oXAX88XWilqgy/D9zvoc+jj9/1kqEDc3l48gjWbm9AUVSdJ4iqutsKEfePtKri/k0hEtHjTxpVVbeOe+xJdZxtBTSqTsStoETXOfG1L68u4dbLBqdlaIET0NlJ6Fu8J6LjgY2qujXGvodxVqwe21GgqruAXe52o4isBcqBNTjf2451qQpwFojtrMtFaL1E6H61Fi8/pjdVlJdkU15ib3rvD+INrjuJESQiUg7chrMi9djO+91jqoBLgKVu0XRgvog8jjPXdlWMal8Hbu3i9abgLiJbWVkZRxf63vXnlZ52/+HmNt7Yur+PWmOM9/V4cl5EQsAtwJwYu2cA31fVcBd1c3FGY9NVtWP54geBR1S1AngEeKZTncuBJlVdFes1VXWWqtaoak1ZWVlPu2GMSQPxjLgmActVdU+MfTXAbHcppVJgsoi0q+orIhLECa0XVPWlqDr3AN91t+cAT3d6zZiju3TkEyHXnSdojyjN4d5dvtyYs9Xx85qVpLm7eIKry/kmVR3SsS0izwGvuqElOCOptar6RKdqO4HrgMU4p5h1Ua/hA+4Aro2jfZ6VlxHgpqHOs8V3NjazZOehJLfImNPr+HlNlh4Fl4hkAxOAB6LKpgKo6szTVL0auBtYKSIr3LIfquprwP3AkyISAJpx56tc1wLbVXVTD/thjOlHehRcqtoElHQqixlYqnpv1Pa7dHFBzd03pot9i4EretK2dJMR8DEg23k3fVN7O0daY04bGtOvped7DzysJCvENRXFANQdPMLKfY1JbpExqcfe8mOM8RwbcaWwgowgQ9xVhA4ca6Wh1d42ZJKnIi+TgM/X/d3UfcCCK4UNyMlgQI4z3/XxnsMWXCapLizLIzuYGpFhp4rGGM+x4DLGeI4FlzHGcyy4jDGekxozbcaYlHRBSS757qOmQ/7UGedYcBljulSaHaIsO/HrIp6t1IlQY4zpIQsuY4znWHAZYzzHgssY4zk2OW+MOcnFA/KpyM8CIOBLgTcmxmDBZYw5id8nKXXrQyyp3TpjjInBgssY4zkWXMYYz7E5LmMMGX4fmQFnHBPypf54ptsWishIEVkR9dEgItO7OHasiIRF5Hb38woRWSQia0VktYh8N+rY0SKyxH3NWhG5LGrf50XkfbfOShHJTEBfjTFdqCrMZnxVGeOryhiUl/q/bt2OuFR1PTAaQET8wA7g5c7HufseA+ZHFbcDj6rqchHJA5aJyEJVXQP8BPixqs4Vkcnu59e7y5X9GrhbVT8WkRKg7Ww6aYxJL/GeKo4HNqrq1hj7HsZZsXpsR4Gq7gJ2uduNIrIWKAfWAArku4cW4CwQC3AT8ImqfuzWOxBnG40xPZAV8FGUGQIgP+StWaN4W3snMVazFpFy4DacFanHdt7vHlMFXAIsdYumA/NF5HGcU9ar3PIRgIrIfKAMmK2qP4nxelNwF5GtrKyMsxvGmNKsEGMHFSW7GWekx7NwIhICbgHmxNg9A/i+qsZcvVREcnFGY9NVtcEtfhB4RFUrgEeAZ9zyAHAN8E3339tEZHzn11TVWapao6o1ZWVlPe2GMSYNxDPimgQsV9U9MfbVALNFBKAUmCwi7ar6iogEcULrBVV9KarOPUDHZP0c4Gl3ezvwlqruBxCR14BLgTfiaKsxxuUXYWhh9inlBZnBJLQmMeIJrruIcZoIoKpDOrZF5DngVTe0BGcktVZVn+hUbSdwHbAY5xSzzi2fD/yjiGQDre4x/zeOdhpjogR8wkUD8rs/0EN6FFxuiEwAHogqmwqgqjNPU/Vq4G5gpYiscMt+qKqvAfcDT7pXEZtx56tU9ZCIPAF8iDOB/5qq/iWeThlj0luPgktVm4CSTmUxA0tV743afpcu1r11943pYt+vcW6JMMacoc+X5SOSuk94OBveugZqjOmxIYXZ+NMwtMDeq2iM8SAbcRmTRkafk092wA9Amg62AAsuYzzNJ1CYEaQoM0hRVohBuZlpOafVmQWXMR6SG/RTnBWiKDNIcVaIgowAPkn/oOrMgsuYFBUQoSjLCahiN6gyUvyRyn3FgsuYFJEX8lOcGXKCKitIfiiA9MPRVE9YcBmTBEGfHD/dK84KUpwZSvkFKlKJBZcxfSA/FKAoK0hJphNUeTaaOisWXMYkWNAnJ81LFWcGCdpoKqEsuIw5S/mhAMVZIUrcU77ckN9GU73MgsuYOIT84k6gOyFVlBUk6IHFJdKNBZcxp1GQETh+uleSFSInaKOpVGDBZYwrw+87PpIqznLuRg/YaColWXCZfkmAgoygE1TuiMpGU95hwWX6hQy/7/jkeXFWiMLMYL94T1+6suAyaUeAwszgidsRsoJkB2w0lU4suIznZQZ8x+elSjKd0VS6PkDPOCy4jKd0PMYl+gbP7KA/2c0yfazb4BKRkcDvooqGAj9S1Rkxjh0LLAG+rqovikgF8DwwEIgAs1T1SffY0cBMIBNoB6ap6gfuwrFrgfXuyy5R1aln1DvjeVkB30m3IxRk2GjK9CC4VHU9MBpARPzADuDlzse5+x7DWV6sQzvwqKouF5E8YJmILFTVNcBPgB+r6lwRmex+fr1bb6Oqjj7TThlv8gnOG4+jbvDMstGUiSHeU8XxOKGyNca+h3EWfh3bUaCqu4Bd7najiKwFyoE1OEuPdSz2VoCzzqLpR7KD/pPez1eYGeyXD8Uz8Ys3uO4kxqKwIlIO3IazsOvYzvvdY6qAS4ClbtF0YL6IPI6zaMdVUYcPEZGPgAbgn1X1nRivNwV3LcbKyso4u+E9FflZFLorD284dJTDLe1JblF8/AJFUSOp4qwgmQEbTZkz0+PgEpEQcAvwgxi7ZwDfV9VwrEvOIpKLMxqbrqoNbvGDwCOq+gcR+RrOitc34ozQKlX1gIiMAV4RkQuj6gGgqrOAWQA1NTXa0354lXNZPwTAjsbmXgsuwVmHz+8T/CLOtjifn7Qt3R9zfJ9PyAr4bTRlEiaeEdckYLmq7omxrwaY7YZWKTBZRNpV9RURCeKE1guq+lJUnXuA77rbc4CnAVS1BWhxt5eJyEZgBFAbR1vTWtDvIzvoPyk8ugqW48ecJmyi91m4GC+IJ7juIsZpIoCqDunYFpHngFfd0BKckdRaVX2iU7WdwHXAYpxTzDq3fhlw0B29DQWqgU1xtDPtjT23MNlNMCapehRcIpINTAAeiCqbCqCqM09T9WrgbmCliKxwy36oqq8B9wNPikgAaMadrwKuBf6XiLQDYWCqqh7scY+MMWlPVL0/PVRTU6O1tXYmaUy6EZFlqlrTudye2WGM8RwLLmOM51hwGWM8x4LLGOM5FlzGGM+x4DLGeI4FlzHGc9LiPi4R2QfEemJFqisF9ie7Eb2sP/QR+kc/k9HH81S1rHNhWgSXV4lIbayb69JJf+gj9I9+plIf7VTRGOM5FlzGGM+x4EquWcluQB/oD32E/tHPlOmjzXEZYzzHRlzGGM+x4DLGeI4F1xkSEb+IfCQir7qfjxaRJSKyQkRqReSyqGN/ICIbRGS9iNwcVT5GRFa6+/7LfWIsIpIhIr9zy5e6C4101LlHROrcj3tSvZ8iki0ifxGRdSKyWkT+I+r4lOhnIr6XUfv/JCKr0rGPIhISkVki8qn7/fxq0vqoqvZxBh/A/wB+g/OYaoAFwCR3ezKw2N0eBXwMZABDgI2A3933AXAlzhoVc6PqTwNmutt3Ar9zt4txHmNdDBS520Wp3E8gG7jBPSYEvJNq/UzE99Ld/xX3dVZFlaVNH4EfA//mbvuA0mT10UZcZ0BEBgNfxF3gw9XVOpG3ArNVtUVVNwMbgMtE5FwgX1XfV+e7/Dzw5ag6v3S3XwTGu6Oxm4GFqnpQVQ8BC4GJvdFHSEw/VbVJVRcBqGorsBwYHFUnqf1MRB/d18nFCYd/6/Ql0qaPwH3AvwOoakRV90fV6dM+xruuonHMAP4RyIsqm07sdSLLgSVRx213y9rc7c7lHXU+A1DVdhE5DJREl8eo0xtmcPb9PE5ECoEvAU9G1Ul2P2eQmD7+K/CfQFOn10+LPrrfO4B/FZHrcUZif6/Oql993kcbccVJRP4G2Kuqyzrt6lgnsgJ4BGd1I3BOAzvT05SfaZ2ESmA/O14vgLNK1H+p6qZu6vRJPxPVRxEZDQxX1ZdjfZlYdU5TnlAJ/D4GcEbK76nqpcD7wOPd1Om1Plpwxe9q4BYR2QLMBsaJyK9x1onsWDdyDieG19uBiqj6g3GG5ds5ccoUXX5SHfcXvgA4eJrX6g2J6meHWUCdqs6IKkt2PxPVxyuBMe7rvAuMEJHFnet4vI8HcEaTL0fVubRznT7rY29NBvaHD+B6Tkx2rgWud7fHA8vc7Qs5ebJzEycmOz8EruDE5Pxkt/whTp7s/L2emOzcjDPRWeRuF3ugn/+Gsyiwr9Prpkw/z7aPUa9TxcmT82nTR9zgc7fvBeYkq49J/+X38kenH4RrgGXuN30pMCbquP+JMyewHvdKjlteA6xy9/0/TryTIRPnL9oGnCuPQ6Pq3OeWbwC+ler9xPkrq+4vygr34zup1s+z/V5G7a/i5OBKmz4C5wFvA58AbwCVyeqjveXHGOM5NsdljPEcCy5jjOdYcBljPMeCyxjjORZcxhjPseAyxniOBZcxxnP+PyUis4Pj4NRkAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "sel_breID = 1189 # 5570\n", "\n", "f, ax = plt.subplots(figsize=(9, 4))\n", "orig_inventory.loc[[sel_breID]].plot(ax=ax, facecolor='lightblue');\n", "cls_default.loc[[sel_breID]].plot(ax=ax);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What can we see?\n", "\n", "- the centerline does not end *exactly* at the glacier outline\n", "- the line seems \"crooked\", it has sudden turns\n", "\n", "Both effects are due to the algorithm we use to compute the centerlines ([Kienholz et al., (2014)](https://tc.copernicus.org/articles/8/503/2014/)),\n", "which works on the underlying glacier grid. Each vertice (point) in the line corresponds to the center of the grid point.\n", "\n", "**We have implemented new options in OGGM v1.6, which allow to circumvent these limitations**. We illustrate them here:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2022-01-05 22:55:51: oggm.utils: Applying global task write_centerlines_to_shape on 6744 glaciers\n", "2022-01-05 22:55:51: oggm.utils: write_centerlines_to_shape on outputs/Norway_Centerlines_smooth.shp ...\n", "2022-01-05 22:55:51: oggm.workflow: Execute entity tasks [get_centerline_lonlat] on 6744 glaciers\n", "2022-01-05 22:55:51: oggm.workflow: WARNING: you are trying to run an entity task on 6744 glaciers with multiprocessing turned off. OGGM will run faster with multiprocessing turned on.\n" ] } ], "source": [ "write_centerlines_to_shape(gdirs, # The glaciers to process\n", " path='outputs/Norway_Centerlines_smooth.shp', # The output file\n", " to_tar=True, # set to True to put everything into one single tar file\n", " to_crs=inventory.crs, # Write into the projection of the original inventory\n", " keep_main_only=True, # Write only the main flowline and discard the tributaries\n", " ensure_exterior_match=True, # NEW! Ensure that the lines are touching the outlines\n", " simplify_line=0.5, # NEW! this option reduces the number of vertices along the line\n", " corner_cutting=5, # NEW! this then augments the number of vertices again\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `simplify_line` and `corner_cutting` options are cosmetic and subjective. The former will simplify the line, by making it look less edgy but also less precise, while the latter then \"smoothes\" it. Users may try different combinations to see their effect (see the [documentation](https://docs.oggm.org/en/latest/generated/oggm.global_tasks.write_centerlines_to_shape.html))." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "cls_smooth = gpd.read_file('tar://outputs/Norway_Centerlines_smooth.tar.gz/Norway_Centerlines_smooth.shp')\n", "cls_smooth['breID'] = gdf['breID']\n", "cls_smooth = cls_smooth.set_index('breID')" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS4AAAEDCAYAAACLcumrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA1WElEQVR4nO3deXhV1dn38e99hpzMCQlBMCQEFFBxAAzWaq0DooBWa9VWa33so1VRa4WqdbZV0apVi76tInWqQ0Vp0VaLIo9DHaECDsgkM4QZAiFzcs653z/2Dh7DAXIg4Qy5P9eVi521h6xFkl/WXntYoqoYY0wy8cS7AsYYEysLLmNM0rHgMsYkHQsuY0zSseAyxiQdCy5jTNJJmeASkadEZIOIfNXG7X8sIvNEZK6I/K2j62eMaT+SKvdxicj3gRrgWVU9dDfb9gVeBk5S1S0i0k1VN+yLehpj9l7K9LhU9X2gMrJMRA4QkTdFZJaIfCAiB7mrLgX+rKpb3H0ttIxJIikTXDsxAbhaVY8ErgMedcv7Af1E5CMRmS4iw+NWQ2NMzHzxrkBHEZFs4Bhgkoi0FAfcf31AX+AEoCfwgYgcqqpb93E1jTF7IGWDC6c3uVVVB0ZZVwFMV9VmYJmILMQJsk/3Yf2MMXsoZU8VVXUbTiidCyCOI9zVrwInuuVdcU4dl8ajnsaY2KVMcInIi8AnQH8RqRCRS4ALgEtE5AtgLnCmu/lUYLOIzAPeBa5X1c3xqLcxJnYpczuEMabzSJkelzGm80iJwfmuXbtqWVlZvKthjGlns2bN2qSqRa3LdxtcItIfeCmiqA9wu6qOi7LtEGA68BNV/buIlADPAt2BMDBBVR92tz0CGA9kA8uBC1R1m4gMA+4F0oAmnPGnd3ZVx7KyMmbOnLm7phhjkoyIrIhWvttTRVVdqKoD3dsKjgTqgFeifAEvcB/OwHeLIHCtqh4MHA1cJSKHuOueAG5U1cPc413vlm8CfuCWXwQ8t/vmGWM6k1jHuIYCS1Q1WgpeDfwD2P74jKquVdXZ7nI1MB8odlf3B953l6cBZ7vbfaaqa9zyuUC6iLTcOGqMMTEH13nAi60LRaQYOAvn1C8qESkDBgEz3KKvgDPc5XOBkii7nQ18pqqNMdbTGJPC2hxcIpKGEzSToqweB9ygqqGd7JuN0xsb7d4YCnAxzqnjLCAHZzwrcp8BOKeel+/kmJeJyEwRmblx48a2NsMYkwJiuao4ApitquujrCsHJrrPBHYFRopIUFVfFRE/Tmi9oKqTW3ZQ1QXAKQAi0g84rWWdiPTEGff6H1VdEq0yqjoB5yFqysvL7WY0YzqRWILrfKKcJgKoau+WZRF5BnjdDS0BngTmq+pDkfu0vANLRDzArbinmSKSD/wbuElVP4qhfsaYTqJNp4oikgkMAyZHlI0SkVG72fVY4ELgJBH53P0Y6a47X0S+BhYAa4Cn3fJfAgcCt0Xs063tTTLGpLqUeOSnvLxc7T4uY1KPiMxS1fLW5Slx57wxM9dupbY56rWhNutXkEWP7PR2qpHpSBZcJiVsaWimuim4V8doDGW0U21MR7OHrI0xScd6XCapzN1Uzcqquh3KG4LhvT52Cgz3dhoWXCapNIfC1LdDSEUTCltyJQsLLhMX4QTs3gS1YwLRtD8LLhMXU5aspymUWOFlPa7kYYPzxriCCdgLNNFZj8vsM3XNIUJuOCRiRliPK3lYcJl9ZsaaLWxpaI53NXbKgit52KmiMSbpWI/LdKitDc00hpyrdc2hBL9q57yWySQBCy7ToRZW1rC6uiHe1WgTi63kYaeKxpikY8FljMt6XMnDgsuYFpZcScOCyxiX5VbysMF50+7W1TSwsc6ZtKmqMXHv22pNLLqShgWXaXeb6ptYtKU23tWIneVW0rBTxSSlqoTCSirMGZAoLLeSx257XCLSH3gpoqgPcLuqjouy7RBgOvATVf27iJQAzwLdgTAwQVUfdrc9AmdKsmxgOXCBqm4TkULg78AQ4BlV/eWeNy/xBcNh1tc20hxWgmElGA67/37782jrFeiWmYbf4/z9Ke+Rj9cTn1+/9bWNLN/qvOCvqil5Tg9NctptcKnqQmAggIh4gdU4k7V+i7vuPmBqRHEQuFZVZ4tIDjBLRKap6jzgCeA6Vf2PiFwMXA/cBjS4/x7qfqS0tTWNfLp26x7vv6HumwnAd5gKZR+qaw6xuiY5bjTdGetxJY9YTxWHAktUdUWUdVfjzFi9oaVAVdeq6mx3uRqYDxS7q/sD77vL04Cz3e1qVfVDnABLSQ3BEOtqG1iwuZoFm6vjXR3jEnvkJ2nEOjh/HlFmsxaRYuAs4CScU7wdiEgZMAiY4RZ9BZwB/BM4FyiJpSIichlwGUBpaWksu3aYpVtrt58u7UxjB756OBbLttaxbGv7DaA3JvpziCaltDm4RCQNJ2huirJ6HHCDqoai/dUSkWyc3thoVd3mFl8MPCIitwP/App22HEXVHUCMAGcCWFj2bejBMPK1sa9myJrX2kMhZKmrvuK9beSRyw9rhHAbFVdH2VdOTDRDa2uwEgRCarqqyLixwmtF1R1cssOqroAOAVARPoBp+1hG9rdvE3VNARDhNS5eucRwecRvCIUZPgpzok+/17LILkxpmPFElznE+U0EUBVe7csi8gzwOtuaAnwJDBfVR+K3EdEuqnqBhHxALfiXGGMm48rKtlc73T6mnfxQrnSUMZOgyvep0tTlkT7mxJdyG6j2IENcSWPNnURRCQTGAZMjigbJSKjdrPrscCFwEki8rn7MdJdd76IfA0sANYAT0cceznwEPBzEakQkUPa2qA91XLLwa5CC5zL/g3BHad6bwqFWVRZ01HVa5PmiDbs7sNe9mmSWZt6XKpaBxS2KovaQ1LVn0csf8hOhg7c+7ke3sm6srbUa28Fw7p9miznrqjdawyFeW/lZg7Mz8LnFVShpinIhrrG3YaeSWw7e+RHVe2KY4Lp1I/8fLGhihVV9THvV9cc4suN23a/oUkq0bJpbU0DeQEfmf5O/auScGw02ZidqA+GmLWuKt7VMFFYcBnjWlPd8M3QgSqz1m6lye5PS0gWXMa4tjUFWVRZS1MozNeVtd96nMokFjtxNybC3E3VzN1kj2ElOutxGWOSjgWXMSbpWHAZY5KOjXEZE0VxTjqZPi8APnsGNeF06uDyzPwvnllz2VTTBB6P8yEewj4vwbR0mrJzqc/tQm1hN9SfhojzGEB6mpf+xbl47G7qlNUnP5OizEC8q2F2olMHl+/tt+j/xr/pv5vtQuJhZc5+zCsoY+Z+/ZnVrT+DD+7OuceUWngZEweSCpMtlJeX68yZM2Peb/2mrWzeUk1DUzPLt9ZR39iMahiag3jr6/FVV5G2eSPpayvIXr6Y3K/n4WuopyE9k3+VfIcVw3/EaUMPsefYUtBxJQXW40oAIjJLVXd4K3mn7nHt1zWf/brmA1CzcjOb6ne84TAE1Lof64NBsuZ8Rt5bUzjnk/eoHT+dd5b9jL6/OM/Cy5h9qFMHV8x8PmoHDaF20BD8K5cRuP8+znz9L0zfWEHujdciNohrzD5hv2l7qLm0N9Xj/sSsI0/m6Blv0Hjn3WizTcuVzA7okskxxV04prgLuQF/vKtjdsGCay+Iz0fGrTfy0XdPZ9Bn75Iz5mp8ayriXS2zh/ICfrpnp9M9O52A1341Epl9d/aSeDx0uWEMr/zgcgrWrqT3Ly9mv8cfIW1VtBncjDHtwca42oGI0O+SHzO+50EMmPIiw958jYIpr7KlYD/WFR9AdV4hzYEMVvc+hPWl/Zz7wQTS/V7KDywkzWd/P4yJhQVXO/GIMPzUw3klr4AX5izl2FVfcMSGRfRZPJ8DG7bhV+c99Z8V9eWxw3/Iqpz9AJizYisXDz0Av4WXMW3Wqe/jivT+Tm6HaA+qitTUkP/OVIpefg5PYwNrL7iEaQefwEufVHBQcS4XD+2Dz8ZV4mpw9zzK8jLjXQ0TYWf3ce32N0VE+kfM0PO5iGwTkdE72XaIiIRE5Bz38xIReVdE5ovIXBG5JmLbI0TkExGZIyKviUhuxLqbRGSxiCwUkVP3qMUJREQgJ4etZ57D0j8/Te2gIRQ/M56zJ43jgoEFLFi9jWfeWUrQ3ra5z/UryOLMvt05s293euVGn3bOJJ7dBpeqLlTVgao6EDgSqANeab2diHiB+4CpEcVB4FpVPRg4GrgqYqqxJ4AbVfUw93jXu8c5BDgPGAAMBx51j50SQvkFVNx8F+suu5rsz2fy4/G3cGlvZV7FNp59b5mF1z4mCF6P82E3ESePWM9NhgJLVDXaJbOrcWas3tBSoKprVXW2u1wNzAeK3dX9gffd5WnA2e7ymcBEVW1U1WXAYuCoGOuZ2ETYctpZrLj7j0hjE2c+djuj8zbw1coqnntvGSGb5syYXYo1uM4jymzWIlIMnMUuZqMWkTJgEDDDLfoKOMNdPhcocZeLgVURu1bwTdhFHu8yEZkpIjM3btwYWysSRP1BA1j+4KM0FZdwyrN/4Nbgl8xZWcXz/7HwMmZX2hxcIpKGEzSToqweB9ygqjtO8ezsm43TGxutqi0TEl6Mc+o4C8gBWkbGo/XXd/gtVtUJqlququVFRUVtbUbCCRYWseKeP1L9nWM59vVnuW/j23y5rJK/vb/cwsuYnYjldogRwGxVXR9lXTkw0R0j6AqMFJGgqr4qIn6c0HpBVSe37KCqC4BTAESkH3Cau6qCb3pfAD2BNTHUM+loegarf/Nbmp8ez+Gv/YNHDq3k1+Ef8qLAT48rw+OxsZf2lOb1kJPmDJtm+lNm+LRTiSW4zifKaSKAqvZuWRaRZ4DX3dAS4Elgvqo+FLmPiHRT1Q0i4gFu5ZvTzH8BfxORh4D9gb7Af2OoZ3Lyetnwi6sIFnXjwKce49Haaq4Jn89EEc77Xi8Lr3a0X2YaQ/bvEu9qmL3QplNFEckEhgGTI8pGicio3ex6LHAhcFLE7RQj3XXni8jXwAKcHtXTAKo6F3gZmAe8CVy1s1PQVFR55rms/vUtFK/6msdmPcHiuSt4+eOV2ycqNcZ08htQN9c3UdMUBODryhqqmxInH7Nmf0rPe39LVVYeo4+8hN6D+tobV9tJSU560ve4KuubqHZ/dnemJDcj6X9e7EWCUSyvqmNFVX28qxFV7eAhrLjrAUrvuJE/TX+Ma8O/YLIIZ3+3xO43MqzcVs/SrXW73KZnTkb0S10pwJ4xSWAN/Q9hxT3jyPDCuE/Gs37G50yeXkEq9JKN2RuduseVDBrL+rD89w9Tevt1PPjJ49wUuoRXPfDDo3pazysG6T7P9ucQ8wLJ82Nf3RikonrHs4ItDbt/aeWCzTW0vqaT4fNSlp/8z2Mmz3ewE2vevycrfv8Ipbdfy73T/8Lt4Wb+JcdyxpBiC682yvR5OaRrTryrEbPqpiDzN9fs0b4LK3fcryDdnxLBZaeKSSJY1I0V9zxMuEcP7prxFDXvvM/rM9fYaaPplDpdj2tNTQNrqhsA56piMgl1KWDl3X+k5PbruWPG04wNB3nDM5QRg/e3ntc+1hAM8dXG6g7/OnXN7Xulu6Y5xMy1W3co75EdoDgned6O0emCa1tjMyu3JeaVxLYI5eaxcuyDlPzuBm799DnuDYd4U05lxOD94121TiUY1qT8OWoKhaPWO8PvpTiJzqTtVDEJhbNzWHXnH2jofzA3zXyBxn+/ydTP18a7WsbsM52ux5UqwplZrPrtfZTcdRM3zHqBBzTMNDmNYUf0iHfVUpaq8v6qSgBCKTa2uLKqjk11ztBJv4IsemSnx7lGu2Y9riSmGRmsuv331B06kOtmv0jjq6/z9pfr4l2tlLa5vonN9U1sbcPtCMmkPhje3rbGJHiZpQVXktP0DCpuu5u6wwZx7eyJNE7+J+99Fe0FHsakDjtVTAEaSKfi1rvpefetjPnsZf4IvC9n8f0B3eJdtbhK83o4qVdXgB1uxDQ7N2fDNuZvcq6YHrV/Fwoz0uJcox1ZjytFaCBAxS1jqTtiMGM+e5naSZP5YN6G3e+YwgTnfVuZfi/pPnvvVls1h5X6YJj6YDhh30piwZVCWsKrduCRjP5sEtWTXuGjBcn5WmuTGFQhrEpYlRf/u5LFGzr+3rW2sOBKMRoIsPrmu6g9YjBjZr/M1pde4ZOFm+JdLZOkPqyo5NWv1/HgH17g5n98wZMfLot3lQALrpSkgQCrbxlL7eGD+PXsl9j00ivM+NrCy+yZysef5vSn7mZM03zuOOPQeFcHsOBKWRoIsPrWsdQeNpDrZk9kw8RX+XTR5nhXa5/I8nvJ9nvJsvfJ77Utj/6FY6c8x5d9j+R/7ryKNF9iREZi1MJ0CA2ks/rWsdQdcjjXzX6RNRP/yczFqR9eJ5cVcUqfbpzgXlE0e0CVqocf5ZipL/JlvyF4fv970tID8a7VdhZcKc65z+se6g46lOtn/Y2VL73G7KWV8a6WSWDS1ET4jjs5+p2/M/uQY/Dcczdef2LdObXb4BKR/hETXXwuIttEZPROth0iIiEROcf9vERE3hWR+SIyV0Suidh2oIhMd485U0SOcsvTRORpEZkjIl+IyAnt0tJOTDMyWH3776nvdxA3znyeJRNf5/NlW+JdraSxsa6RdbUNrK9tjHdVOlxg0QK6XvULBnz2H949+gcE7roj4UIL2nADqqouBAYCiIgXWA280no7d919wNSI4iBwrarOFpEcYJaITFPVecD9wB2q+oY788/9wAnApe7XPUxEugFviMgQVU385xASWDgzk9W/u4+et1/PzZ8+x1iPBznvNI4oS+5JI/aFmWurqA8mzkQq7a2qromqhUs5aOo/KP3kHSrTc3ju7DEM/tnpCTstXqxROhRYoqoroqy7Gmfi1yEtBaq6FljrLleLyHygGGfqMQVy3U3z+GbS10OAt919NojIVpwJZ1N/bsUOFs7MouJ391Ny27XcMuNZ7hAvct4IDu+VH++qmTjwbqui4f0P8b8+hWFrFxD0eJl8wHF8ffpPOOPEgxJ6hqBYg+s8okwKKyLFwFnASUQEV6ttyoBBwAy3aDQwVUQewDllPcYt/wI4U0Qm4sxofaT7rwVXOwhnZ7Pqzj9QcsuvuX3GM/zW48Xzk1M4tDQ/3lXbK8U56aR5nJGPBP59iztPbQ15b79J3n/eJn3J14gqlZl5LBp+DhUnnUZOYVfOLMxI+BdTtjm4RCQNOAO4KcrqccANqhqK1mARycbpjY1W1W1u8RXAGFX9h4j8GGfG65OBp4CDgZnACuBjnFPO1se8DLgMoLS0tK3NMEA4J5eKsQ9ScvMYfjf9KW7zePH8+GQOKcmLd9X22MGF2eQG/PGuRkL6bGkln361hqFfTmPoF9PIaGpgZVEvPjzkVBaWHMIpFwylIDeDZHoVZSw9rhHAbFWN9uqBcmCiG1pdgZEiElTVV0XEjxNaL6jq5Ih9LgJaBusnAU8AqGoQGNOykYh8DCxq/QVVdQIwAZwJYWNoh8F5k+qqsQ9QetNo7vjkSW7xePH8+CQOKs7d/c4maUz/ehOLX5nGzV/+g6KaSj4vO4Ipg0eyqrCUzHQvZx9dSpfsxHuIendiCa7ziXKaCKCqvVuWReQZ4HU3tASnJzVfVR9qtdsa4HjgPZxTzEXu/pk4M2zXisgwIOgO5pt2FsovYNXYBym5aTRjP/4LN3u8yDnH0z9Jwqs0N4PcNOdHOOBt35tNl26tpdad2bw5nJzXhWbMWU3hk+P532Uf0VBcyoqbbiNw6BGcFe+KtYM2BZcbJsOAyyPKRgGo6vhd7HoscCEwR0Q+d8tuVtUpOFcPHxYRH9CAe9oHdMMZ+wrjXMG8sM2tMTELFhaxauyDlN40mrEfTeBGjxfP2cfRd//EfwF5cU56h72ps2JbA5uSbDKVSF98PI9jx/+evlWr2Xj6j9h80WVoWvL1rHamTcGlqnVAYauyqIGlqj+PWP6QnUwC7q47Mkr5cqB/W+pl2kewW/ft4XX3h49zo8eLnH0MB3ZP/PAyO/p6ygcMf/p+AoRZfuOd1H/3e/GuUrvrdHfO75+dznf2z+c7++dTlJk6f4H2VnOPYlaNfZDMgJe7P3yc1175L0vX7dlEpCZ+Vk38FyP/cifNmVmsfPDRlAwt6ITBlRvwU5yTQXFOBpn2EO63NPUsZdVdD5LjCXP3B4/xyj8/Zdl6C69kseGJ5zn5xYdZ3b03mx4eT7i0V7yr1GE6XXCZXWvs1ZuKO++nizYy9oPxTPrXbJZvqI13tcxuVP75Lxz/2lMs7H04tQ89Avn58a5Sh+rUwXVQYTYn9urKib26km/3AG3XcEA/Kn57L0VN1dz1wXhefO0zVm5M/fCav6mad5Zv5J3lG9mSRLP4bPl/4zn2rRf5ov9RhO6/H09Gx1yw+Lhii/v/E/93u3Xq4Mry++iS7qdLuh9fgj6TFS/1Bw2g4ta76VFXyR0fTOC51+ewalNdvKvVoeqaQ2xtDLK1MZgU8yaqKmvHP8Mx//cyXxx0NL67x+JN67g/wNuagu7/T/xDvVMHl9m1usMHsfrmO+m1bR23fziBZ17/iorNqR1eyUJVWfjSFE544znmHzAQ39g78STgWxw6igWX2aXawUex5vrb6Fu5kls+foKn/j2X1ZWJEV7/XbOV1xat47VF69iWAL2AfUVV+XjabE6d9CfW7dcLvefuThVaYMFl2qD6u8exZvRNHLJhCTd+9DRP/ns+a7fUx7tahFRpDjsfnYWqMmXGCo5/fhz4fNSMvQdJz4h3tfY5Cy7TJtuOH8raX17HEesWcN3Hz/CXf89n3db4h1eL5rDSFArTFAqjMYxPhSL2S9Q5BFuoKv/6dDXdJ79I36rVbBr9G0LduselLi3/Z82h+DwO1bn6l2avVJ08Amlu4sjxDzP6k2cZ57mIy0cezH75HXMVKxb/WfnNu/TP7NsdbxuvtSzaUsO8TYl/r5qq8uqMCpZ+Oo/rF71D1XEnUfvd4+JWn9cXO+9ayPJ7ObXPvp8x3XpcJiZbR5zJul9cxXcqvuRXHz/HX//v63hXqVN4ZUYFH8zfyA0r3oJAgPWXXBnvKsWV9bhc2WlemsPOf8e2xiCJfdIQX1t+cDYSDvPdpx6j6b2/wln3Qzu/nWFvbG1sxtvGF+HVBxP/zQ+bqxv5cP5Gzsmtou+iz9jws0sIdSmId7XiyoLLNbh7/vblN5dsoC6F3zHeHirPPJd5yys58Z2XqBp3L2tG35gw4RV52pgKWsaRTp39JsGcXCpP/1GcaxR/Flxmj80/6YcsW1fNxe9PAUio8EolqtC9djPFc2ey6dyfoRmd7ypiaxZcUeyXFaAxFAaUNTWpPyXVnvr+gG48f8qPeArc8FLWjL7JwqudqSqnrpiBImwdfnq8q5MQLLiiGNTdefe6qvLK1+viXJvE5fUIPzu+N8+Gz9oeXhIKs/rXN4PPfrTai4aVE1d9xvqDjiBYWBTv6iQEu6po9orXI1x4Qm/mDT2LJwacTu5H71H8wFho7jx3sne07BVL2K9+C6sGHxvvqiQMCy6z13xeDxed2Jv5J57B44eeQe4n79Pzvt8hzcn76uNE0vXLTwkjrBuwwwuDOy0LLtMunPDqw4LjT+fPh59Fzqef0PPu25DGhnhXLel1WTCHpXn705SbH++qJAwLLtNu/D4P/3tSHxZ9bzh/HPRjsj6fScmdN+GpS4yHspNSKETe0oXMLewdffKGTmq3wSUi/UXk84iPbSIyeifbDhGRkIic435eIiLvish8EZkrItdEbDtQRKa7x5wpIke55X4R+auIzHH3izYBrUlQfp+H/x16AMuOHsr9R/6UjHlzKP3t9XhqquNdtaSUtnolvqZGvs4vScgZuhtDYT5bX8Vn66tYvGXfvWxyt8GlqgtVdaCqDsSZlacOeKX1diLiBe4DpkYUB4FrVfVg4GjgKhE5xF13P3CHe9zb3c8BzgUCqnqY+/UuF5Gy2Jtm4iXN5+GSkw9gVflxjB3yP6QtWUSvW8bg3VoZ76olncCqFQAsz+2Oz5N4J0jBsLJsax3LttaxrmbfDQvE+j8xFFiiqiuirLsaZ8bqDS0FqrpWVWe7y9XAfKC4ZTXQMvNoHs4EsS3lWe58ixlAE7AtxnqaOEvzefjFyQewfuDR3H70xfjWrKbsxmvwbbDbS2KxYY4ziXvggDJ675cV59okjliD6zyizGYtIsXAWcBOJ4d1e02DgBlu0WjgDyKyCngAaDkl/DtQC6wFVgIPqOoOf6pF5DL3FHPmxo0bY2xG2x1T3IVjirsweL+8DvsaqSrg9/KLYQdQOWAQNx59KWzdStmNvyJt5fJ4Vy0pzFpSyfp5S6nOyOGikQPweROvxxUvbf6fEJE04AxgUpTV44AbVDXqA34iko3TGxutqi29pyuAMapaAowBnnTLjwJCwP5Ab+BaEenT+piqOkFVy1W1vKioY27KExG6Z6fTPTudrjYH4x5J93u5bNiB1PQfwK+PGUWwOUivm0eTvnBevKuW8N6Zs44+jZV4SktI81loRYrlf2MEMFtV10dZVw5MFJHlwDnAoyLyQ3AG23FC6wVVnRyxz0VAy+eTcAIL4KfAm6rarKobgI/c45sklZ7m5bJTDqSp9wH88rtX0hDIpNdt15E1a8bud+7EQqEwPbeupam0LN5VSTixBNf5RDlNBFDV3qpapqplOKd6V6rqqyIiOD2p+ar6UKvd1gDHu8snAYvc5ZXASeLIwhnUXxBDPTtEpt/LsLIihpUVMbi7nTbGKiPNy6hTDsTTsyejjrqcbd32p2TsLeS9/Wa8q5awulZvJquhhoYD+sW7Km2yub6Jt5Zt4K1lG5i3qWOvIrcpuEQkExjGNz0kRGSUiIzaza7HAhfiBFHL7RQj3XWXAg+KyBfAPcBlbvmfgWzgK+BT4GlV/bKtDeooHhFyAj5yAj4yffYQ8Z7ICPi4/NQDSe/ejcsGXcLGfoex/yP3U/jy884rEMy39F2/BHCmiksGIYWaphA1TSEaOvg9Z216ElZV64DCVmVRB+JV9ecRyx9C9Pvm3HU7PMOgqjU4t0SYFJQV8DHq1L48+uYiLj/4pzySm0/JC0/h37SBdZdfY2+WiHDo2gXUZmTT2Kt3vKuScGzEz+xzWek+rhjel/y8LK7q+QMWDT+HLlNfp+TuW+wu+xahEEesnseissMgAe/fijf7HzFxke2GV0FuOtdlH8OXF1xJ1mcz6XXTr/Bt3LD7A6S4rC9mkdNUx7x+9mB1NBZcJm5yMvxccWpf8rP83NZ4IJ9e/Vv8G9bT+/orSV84P97Vi6u8t6dSk5bJ4j6Hx7sqeyQUDlPfHKK+OURzuP3Huyy4TFzlZvq5cng/8jL93LM6h49/cx/htAC9bhlN7n/ejnf14sK3eRO5n7zP+72HEPL5412dPbKquoE3lm7gjaUbWLql/U//LbhM3OVm+rlieF9yMvw8NLeJD278A/X9Dqb4obsp+usECHWuiUsKXvsHqPJm3+/jScAHqxOBBZdJCPlZaVw5oi9Z6T7+9MlGPv7V79gy/Ad0nTyRkrG34KnuHI+rerdU0mXKP1k35DiWBfLx22M+Udn/ikkY+VlpXDm8H5lpXh57ewUzz72ctVeMIevL2fS+9goCSxfHu4odruiFpyDYzNiCY8nJ8HP8gH0/S3QysOAyCaVLdhpXDu9Lut/L+LcWMfeooSy/ZxwSbKbshl+S939vxLuKHSZjwVy6TJvCawccx9bC7lw1oi8FOYF4VyshWXDtgYDPQ2luBqW5GXRJT87B00RWkBPgyhF9SfN6GP/mIpYW9WbZQ49Tf9Ch7P///kCPh+9D6uvjXc12JQ31dH3oXjZkduG1wadx1Yh+FGSnRmhtbWxmRVUdK6rqqG0OtssxLbj2QF7AT3mPfMp75FOWlxnv6qSkwpwAV4zoi9frYfzURazWdFb+7j42/uR/yHv3LXpfO4rAsiXxrma7yXrkj2SuX8OTx17AJWccRpfs1HkbyerqBmatq2LWuiq21LfP7E8WXCZhFeWmc+XwvnhEePTNRayvbmLTT3/OyjsfwFNXR9l1V1Lwz0nQAfcJ7UvNL71MyUf/x2uHn8KJF59OflbqhFZHseDaS/npPg4qzOagwmwKM+y0sb11y0vniuF9EYFH31zEhqoG6g4fxLJHnqD2yKPY76nHKL39Ovzrk/PNqtumvsOAFyfwec8BdL/hGvJS8L1vRZlp239HcgLtM1GwBdde6pKexiFdczikaw6FGan3Q5cI9st3wkvVCa+N2xoI5eZRcdOdrLn6etIXf03vX11C/hv/TKre15a3P2Dg+HtZWdCT4F1jyc1Oj3eVOkS3zMD235G8QPv8cbfgMkmhe34GVwzvSyikPPbGIjZtawQRqk4ewdJHnqC+/yH0GP8wvW4ZkxSvht72+puU/+kuNuR1Y+u9D5BVYO94i4UFl0kaPbpkMGr4gTSFwjz65tdsrm4EINitO6vuuJ81V19PYOVy+oy+lKK/TkjMN000NyN/+hPf+cv9rOxawsY/PEJ6t67xrlW765qRxpHd8ziyex7dO+DqqAWXSSrFBZmMOrUvjc1hHn1zEZU1Tni19L6WPPpXqk4YRtfJEzngyv8h/61/J8wjQxkL5tLjmss5aNpk3u93DFUP/T8yigriXa0OkZ3mo1deJr3yMtvt9DCSBZdJOj0LnfCqbwrx6BuL2FLTtH1dKC+ftb/6Dcv+8GeaunWnx58fpM+vfkHOB+/GbfwrsHQRxff+jrIbria0uZI/D7ucrLG/IzPHbqXZU6Ip8Mrc8vJynTlzZryrQW1TkPqg89d93qYaNtU37WYPszdWbKzl8amLyEr3c9WIvjveRqBKzvQPKHr+aQIVK2jsWcrmH/6Ybd8figY69uZOT10d2TM+JH/aFLLmfklzIINJvb/H9KNH8PORh5LRTlfXEkl+wM/h3XIACPi85KTtfRtFZJaq7jBZjgVXB5mxZgurq/fdzL6d1fINtTz+1iJyM/xcNaIfuZlRTktCIXI/+g+Fr0wkfeligjm5bDvhZKpOPIWGPn1pl7ntQyECq5aTOecLsj/7lMwvZ+Npbqapew/mHnUyvw/1o3D/rlx2yoGkp6Xm66m7ZabxvZLC3W8YAwuufcyCa99Ztr6Gx99aTH6Wn1Gn9iU73YeIIOJMeCAtwaRK5ldf0OWNf5I942M8wWaa9utB7eAh1A44goYD+tK8X4/o771XRRoa8FZX4du6Bf/mjfjXryNtTQWBFctIX74ET4Pz/W7qUUxN+XfYdszxfJTWg799uJJe3bK4dNiBpPtTM7QgwYJLRPoDL0UU9QFuV9VxUbYdAkwHfqKqfxeREuBZoDsQBiao6sPutgNxZr5OB4I4U5r9V0QuAK6POOzhwGBV/XxndUzE4GoMhQmFnf/bT1ZXUtXYPs9omeiWrqthwrTFNEWZXcYJL74VZrlNdRy9dg5Hr/6KARsXkxF0TuvDItT4Mwl6vYiCR8P4w82kB5vwRPldqQ5ksTa/OysLe7K8axmLuh/IltxC9+sIKzfV0rtbNpcOO4BACoZWps/L90udsPKKc4rYntqlxyUiXmA18B1VXRFl3TSgAXjKDa4eQA9VnS0iOcAs4IeqOk9E3gL+qKpvuFOW/UZVT2h1zMOAf6rqDjNZR0rE4Ir07opNbGlon2e0zM6t3lzH/IptKIoqKKDqLiuE3Z91VcX9m0I4rEgoSLeNFXRbv5L8rRtIr6/FEwoCQkg8BH0+mnwBGtLSqU3PpiYjh6rsLmzJLaA2kA2qKO7xFcIRX7sgO40zj+qZkqEFkOX3cmqfjnv1zs6CK9bRs6HAktah5boaZ8bqIS0FqroWWOsuV4vIfKAYmIfzvc11N83DmSC2tZ1OQptMhJ3M0RYh+U/Y46+4MJPiwj29UmdTgCWTWIPrPKIEiYgUA2fhzEg9pPV6d5syYBDQMu/6aGCqiDyAc1vGMVF2+wlw5k6OdxnuJLKlpaUxNGHfO6HXrm8wrGpo5u0Vm/ZRbYxJfm2+j0tE0oAzgElRVo8DblDVqHf6iUg2Tm9stKq2vIP3CmCMqpYAY4AnW+3zHaBOVb+KdkxVnaCq5apaXlRU1NZmGGNSQCw9rhHAbFVdH2VdOTDRvXrTFRgpIkFVfVVE/Dih9YKqTo7Y5yLgGnd5EvBEq2NG7d2lIo8I2e4YSDCsNISS50Fh0zm1/LxmxGnsLpbg2ul4k6puHyAQkWeA193QEpye1HxVfajVbmuA44H3cE4xF0UcwwOcC3w/hvolrZyAj1PcAc411Q1MX7MlzjUyZtdO6cAB+bZoU3CJSCYwDLg8omwUgKqO38WuxwIXAnNE5HO37GZVnQJcCjwsIj6cK5GXRez3faBCVZe2sR3GmE6kTcGlqnVAYauyqIGlqj+PWP6QnVxQc9dFnV9cVd8Djm5L3VJNwOehW6bzOEpdMEhNU2I8IGxMIkm9B6aSXGFGGt8rcd4YsKiyhjkbq+NcI2MSj70dwhiTdKzHlcDyAn56u7MIba5vYluTPTZk4qckJx2fx7P7u6n3AQuuBNYtK0C3LGe864v1VRZcJq4GFOWQ6U+MyLBTRWNM0rHgMsYkHQsuY0zSseAyxiSdxBhpM8YkpIMLs8l134+f5k2cfo4FlzFmp7pmplGU2bETi+yJxIlQY4xpIwsuY0zSseAyxiQdCy5jTNKxwXljzLcc0S2XktwMAHyeBHgwMQoLLmPMt3g9klC3PkST2LUzxpgoLLiMMUnHgssYk3RsjMsYQ8DrId3n9GPSPInfn9ltDUWkv4h8HvGxTURG72TbISISEpFz3M9LRORdEZkvInNF5JqIbQeKyHT3mDNF5KiIdYeLyCfuPnNEJL0d2mqM2Ymy/EyGlhUxtKyI/XMS/9dttz0uVV0IDAQQES+wGnil9XbuuvuAqRHFQeBaVZ0tIjnALBGZpqrzgPuBO1T1DREZ6X5+gjtd2fPAhar6hYgUAs1700hjTGqJ9VRxKLBEVVdEWXc1zozVQ1oKVHUtsNZdrhaR+UAxMA9QINfdNA9ngliAU4AvVfULd7/NMdbRGNMGGT4PXdLTAMhNS65Ro1hrex5RZrMWkWLgLJwZqYe0Xu9uUwYMAma4RaOBqSLyAM4p6zFueT9ARWQqUARMVNX7oxzvMtxJZEtLS2NshjGma0YaQ/bvEu9q7JE2j8KJSBpwBjApyupxwA2qGnX2UhHJxumNjVbVbW7xFcAYVS0BxgBPuuU+4HvABe6/Z4nI0NbHVNUJqlququVFRUVtbYYxJgXE0uMaAcxW1fVR1pUDE0UEoCswUkSCqvqqiPhxQusFVZ0csc9FQMtg/STgCXe5AviPqm4CEJEpwGDg7RjqaoxxeUXok5+5Q3leuj8OtWkfsQTX+UQ5TQRQ1d4tyyLyDPC6G1qC05Oar6oPtdptDXA88B7OKeYit3wq8BsRyQSa3G3+GEM9jTERfB7hsG65u98wibQpuNwQGQZcHlE2CkBVx+9i12OBC4E5IvK5W3azqk4BLgUedq8iNuCOV6nqFhF5CPgUZwB/iqr+O5ZGGWNSW5uCS1XrgMJWZVEDS1V/HrH8ITuZ99Zdd+RO1j2Pc0uEMWYPHV6Ui0jivuFhbyTXNVBjTJv1zs/Em4KhBfasojEmCVmPy5gUMnC/XDJ9XgBStLMFWHAZk9Q8AvkBP13S/XTJSGP/7PSUHNNqzYLLmCSS7fdSkJFGl3Q/BRlp5AV8eCT1g6o1Cy5jEpRPhC4ZTkAVuEEVSPBXKu8rFlzGJIicNC8F6WlOUGX4yU3zIZ2wN9UWFlzGxIHfI9tP9woy/BSkpyX8BBWJxILLmH0gN81Hlww/helOUOVYb2qvWHAZ0878HvnWuFRBuh+/9abalQWXMXspN81HQUYahe4pX3aa13pTHcyCy5gYpHnFHUB3QqpLhh9/EkwukWosuIzZhbyAb/vpXmFGGll+600lAgsuY1wBr2d7T6ogw7kb3We9qYRkwWU6JQHyAn4nqNwelfWmkocFl+kUAl7P9sHzgow08tP9neKZvlRlwWVSjgD56f5vbkfI8JPps95UKrHgMkkv3efZPi5VmO70plL1BXrGYcFlkkrLa1wib/DM9HvjXS2zj+02uESkP/BSRFEf4HZVHRdl2yHAdOAnqvp3ESkBngW6A2Fggqo+7G47EBgPpANB4EpV/a87cex8YKF72OmqOmqPWmeSXobP863bEfIC1psybQguVV0IDAQQES+wGnil9XbuuvtwphdrEQSuVdXZIpIDzBKRaao6D7gfuENV3xCRke7nJ7j7LVHVgXvaKJOcPILz4HHEDZ4Z1psyUcR6qjgUJ1RWRFl3Nc7Er0NaClR1LbDWXa4WkflAMTAPZ+qxlsne8nDmWTSdSKbf+63n+fLT/Z3ypXgmdrEG13lEmRRWRIqBs3Amdh3Ser27TRkwCJjhFo0GporIAziTdhwTsXlvEfkM2AbcqqofRDneZbhzMZaWlsbYjORTkptBvjvz8OIttVQ1BuNco9h4BbpE9KQKMvyk+6w3ZfZMm4NLRNKAM4CboqweB9ygqqFol5xFJBunNzZaVbe5xVcAY1T1HyLyY5wZr0/G6aGVqupmETkSeFVEBkTsB4CqTgAmAJSXl2tb25GsnMv6aQCsrm7osOASnHn4vB7BK+Isi/P5t5Zl99tsX+cRMnxe602ZdhNLj2sEMFtV10dZVw5MdEOrKzBSRIKq+qqI+HFC6wVVnRyxz0XANe7yJOAJAFVtBBrd5VkisgToB8yMoa4pze/1kOn3fis8dhYs27fZRdhErrNwMckgluA6nyiniQCq2rtlWUSeAV53Q0twelLzVfWhVrutAY4H3sM5xVzk7l8EVLq9tz5AX2BpDPVMeUN65Me7CsbEVZuCS0QygWHA5RFlowBUdfwudj0WuBCYIyKfu2U3q+oU4FLgYRHxAQ2441XA94E7RSQIhIBRqlrZ5hYZY1KeqCb/8FB5ebnOnGlnksakGhGZparlrcvtnR3GmKRjwWWMSToWXMaYpGPBZYxJOhZcxpikY8FljEk6FlzGmKSTEvdxichGINobKxJdV2BTvCvRwTpDG6FztDMebeylqkWtC1MiuJKViMyMdnNdKukMbYTO0c5EaqOdKhpjko4FlzEm6VhwxdeEeFdgH+gMbYTO0c6EaaONcRljko71uIwxSceCyxiTdCy49pCIeEXkMxF53f18oIhMF5HPRWSmiBwVse1NIrJYRBaKyKkR5UeKyBx33SPuG2MRkYCIvOSWz3AnGmnZ5yIRWeR+XJTo7RSRTBH5t4gsEJG5InJvxPYJ0c72+F5GrP+XiHyVim0UkTQRmSAiX7vfz7Pj1kZVtY89+AB+DfwN5zXVAG8BI9zlkcB77vIhwBdAAOgNLAG87rr/At/FmaPijYj9rwTGu8vnAS+5ywU4r7EuALq4y10SuZ1AJnCiu00a8EGitbM9vpfu+h+5x/kqoixl2gjcAYx1lz1A13i10Xpce0BEegKn4U7w4drZPJFnAhNVtVFVlwGLgaNEpAeQq6qfqPNdfhb4YcQ+f3WX/w4MdXtjpwLTVLVSVbcA04DhHdFGaJ92qmqdqr4LoKpNwGygZ8Q+cW1ne7TRPU42TjiMbfUlUqaNwMXA7wFUNayqmyL22adtjHVeReMYB/wGyIkoG030eSKLgekR21W4Zc3ucuvyln1WAahqUESqgMLI8ij7dIRx7H07txORfOAHwMMR+8S7neNonzbeBTwI1LU6fkq00f3eAdwlIifg9MR+qc6sX/u8jdbjipGInA5sUNVZrVa1zBNZAozBmd0InNPA1nQX5Xu6T7tqx3a2HM+HM0vUI6q6dDf77JN2tlcbRWQgcKCqvhLty0TbZxfl7aodv48+nJ7yR6o6GPgEeGA3+3RYGy24YncscIaILAcmAieJyPM480S2zBs5iW+61xVAScT+PXG65RV8c8oUWf6tfdxf+DygchfH6gjt1c4WE4BFqjouoize7WyvNn4XONI9zodAPxF5r/U+Sd7GzTi9yVci9hncep991saOGgzsDB/ACXwz2DkfOMFdHgrMcpcH8O3BzqV8M9j5KXA03wzOj3TLr+Lbg50v6zeDnctwBjq7uMsFSdDOsTiTAntaHTdh2rm3bYw4ThnfHpxPmTbiBp+7/HNgUrzaGPdf/mT+aPWD8D1glvtNnwEcGbHdLThjAgtxr+S45eXAV+66P/HNkwzpOH/RFuNceewTsc/Fbvli4H8TvZ04f2XV/UX53P34RaK1c2+/lxHry/h2cKVMG4FewPvAl8DbQGm82miP/Bhjko6NcRljko4FlzEm6VhwGWOSjgWXMSbpWHAZY5KOBZcxJulYcBljks7/BwWxWD7QJuZ0AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "sel_breID = 1189\n", "\n", "f, ax = plt.subplots(figsize=(9, 4))\n", "orig_inventory.loc[[sel_breID]].plot(ax=ax, facecolor='lightblue');\n", "cls_default.loc[[sel_breID]].plot(ax=ax);\n", "cls_smooth.loc[[sel_breID]].plot(ax=ax, color='C3');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Final remarks" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "While the centerline algorithm is quite robust, the results will vary as a function of the resolution of the underlying grid, and the smoothing options. After trying a little, it seems difficult to find a setting which works \"best\" in all circumstances, and we encourage users to try several options and see what they prefer." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What's next?\n", "\n", "- return to the [OGGM documentation](https://docs.oggm.org)\n", "- back to the [table of contents](welcome.ipynb)" ] } ], "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.8.8" }, "metadata": { "interpreter": { "hash": "705f036afebab14ba3958dfbf5720c1e1e37a03d5afe33574ff09620abf8737d" } } }, "nbformat": 4, "nbformat_minor": 4 }