From 35f4de708f3420dda9d419454fc0db9179da3c8b Mon Sep 17 00:00:00 2001 From: Brendon Gory <79238147+BGory@users.noreply.github.com> Date: Wed, 28 Feb 2024 11:07:47 -0500 Subject: [PATCH] Add files via upload --- icedyno/preprocess/preprocess_netcdf.ipynb | 289 +++++++++++++++++++++ 1 file changed, 289 insertions(+) create mode 100644 icedyno/preprocess/preprocess_netcdf.ipynb diff --git a/icedyno/preprocess/preprocess_netcdf.ipynb b/icedyno/preprocess/preprocess_netcdf.ipynb new file mode 100644 index 0000000..6a9295d --- /dev/null +++ b/icedyno/preprocess/preprocess_netcdf.ipynb @@ -0,0 +1,289 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "c8a671ab", + "metadata": {}, + "source": [ + "## Import libraries" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "1c3ccd02", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.colors import ListedColormap\n", + "import numpy as np\n", + "\n", + "import requests\n", + "from bs4 import BeautifulSoup\n", + "\n", + "import h5py\n", + "\n", + "years = range(2014, 2025)" + ] + }, + { + "cell_type": "markdown", + "id": "b9cce684", + "metadata": {}, + "source": [ + "## Download netcdf files" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "03faece8", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Downloading files for 2014\n", + "Total downloaded files for 2014: [30]\n", + "\n", + "Downloading files for 2015\n", + "Downloading files for 2015: [50]\n", + "Downloading files for 2015: [100]\n", + "Downloading files for 2015: [150]\n", + "Downloading files for 2015: [200]\n", + "Downloading files for 2015: [250]\n", + "Downloading files for 2015: [300]\n", + "Downloading files for 2015: [350]\n", + "Total downloaded files for 2015: [364]\n", + "\n", + "Downloading files for 2016\n", + "Downloading files for 2016: [50]\n", + "Downloading files for 2016: [100]\n", + "Downloading files for 2016: [150]\n", + "Downloading files for 2016: [200]\n", + "Downloading files for 2016: [250]\n", + "Downloading files for 2016: [300]\n", + "Downloading files for 2016: [350]\n", + "Total downloaded files for 2016: [366]\n", + "\n", + "Downloading files for 2017\n", + "Downloading files for 2017: [50]\n", + "Downloading files for 2017: [100]\n", + "Downloading files for 2017: [150]\n", + "Downloading files for 2017: [200]\n", + "Downloading files for 2017: [250]\n", + "Downloading files for 2017: [300]\n", + "Downloading files for 2017: [350]\n", + "Total downloaded files for 2017: [363]\n", + "\n", + "Downloading files for 2018\n", + "Downloading files for 2018: [50]\n", + "Downloading files for 2018: [100]\n", + "Downloading files for 2018: [150]\n", + "Downloading files for 2018: [200]\n", + "Downloading files for 2018: [250]\n", + "Downloading files for 2018: [300]\n", + "Downloading files for 2018: [350]\n", + "Total downloaded files for 2018: [365]\n", + "\n", + "Downloading files for 2019\n", + "Downloading files for 2019: [50]\n", + "Downloading files for 2019: [100]\n", + "Downloading files for 2019: [150]\n", + "Downloading files for 2019: [200]\n", + "Downloading files for 2019: [250]\n", + "Downloading files for 2019: [300]\n", + "Downloading files for 2019: [350]\n", + "Total downloaded files for 2019: [365]\n", + "\n", + "Downloading files for 2020\n", + "Downloading files for 2020: [50]\n", + "Downloading files for 2020: [100]\n", + "Downloading files for 2020: [150]\n", + "Downloading files for 2020: [200]\n", + "Downloading files for 2020: [250]\n", + "Downloading files for 2020: [300]\n", + "Downloading files for 2020: [350]\n", + "Total downloaded files for 2020: [366]\n", + "\n", + "Downloading files for 2021\n", + "Downloading files for 2021: [50]\n", + "Downloading files for 2021: [100]\n", + "Downloading files for 2021: [150]\n", + "Downloading files for 2021: [200]\n", + "Downloading files for 2021: [250]\n", + "Downloading files for 2021: [300]\n", + "Downloading files for 2021: [350]\n", + "Total downloaded files for 2021: [363]\n", + "\n", + "Downloading files for 2022\n", + "Downloading files for 2022: [50]\n", + "Downloading files for 2022: [100]\n", + "Downloading files for 2022: [150]\n", + "Downloading files for 2022: [200]\n", + "Downloading files for 2022: [250]\n", + "Downloading files for 2022: [300]\n", + "Downloading files for 2022: [350]\n", + "Total downloaded files for 2022: [365]\n", + "\n", + "Downloading files for 2023\n", + "Downloading files for 2023: [50]\n", + "Downloading files for 2023: [100]\n", + "Downloading files for 2023: [150]\n", + "Downloading files for 2023: [200]\n", + "Downloading files for 2023: [250]\n", + "Downloading files for 2023: [300]\n", + "Downloading files for 2023: [350]\n", + "Total downloaded files for 2023: [364]\n", + "\n", + "Downloading files for 2024\n", + "Downloading files for 2024: [50]\n", + "Total downloaded files for 2024: [56]\n", + "\n" + ] + } + ], + "source": [ + "root_dir = \"D:/IceDyno/netcdf\"\n", + "for yr in years:\n", + " print(f\"Downloading files for {yr}\")\n", + " save_dir = f\"{root_dir}/{yr}/\"\n", + " os.makedirs(save_dir, exist_ok=True)\n", + "\n", + " html_page = f\"https://noaadata.apps.nsidc.org/NOAA/G02186/netcdf/1km/{yr}/\"\n", + " response = requests.get(html_page)\n", + "\n", + " try:\n", + " if response.status_code == 200:\n", + " soup = BeautifulSoup(response.text, \"html.parser\")\n", + " links = soup.find_all(\"a\")\n", + " file_cnt = 0\n", + " for link in links:\n", + " href = link.get(\"href\")\n", + " if \".nc\" in href:\n", + " download_url = f\"{html_page}/{link['href']}\"\n", + " response = requests.get(download_url, stream=True)\n", + " response.raise_for_status()\n", + "\n", + " file_name = f\"{save_dir}{os.path.basename(download_url)}\"\n", + "\n", + " with open(file_name, \"wb\") as file:\n", + " for chunk in response.iter_content(chunk_size=512):\n", + " file.write(chunk)\n", + "\n", + " file_cnt += 1\n", + " if file_cnt % 50 == 0:\n", + " print(f\"Downloading files for {yr}: [{file_cnt}]\")\n", + "\n", + " except requests.exceptions.RequestException as e:\n", + " print(f\"Error: {e}\")\n", + "\n", + " print(f\"Total downloaded files for {yr}: [{file_cnt}]\\n\")" + ] + }, + { + "cell_type": "markdown", + "id": "04b105e1", + "metadata": {}, + "source": [ + "## Show example" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "23d3cd7c", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Unique values in HDF file [0 1 2 3 4 5]\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAp0AAAGeCAYAAADWqboNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAABrr0lEQVR4nO3dd3gUdf4H8PcESC90SIMgIRCkRlpEktAuAirtsIA/iCLkRKRIUe5AilhAmp14IHCKynmUQ0RQOBJAsBAIoHSkhQRRkASQFnZ+f8RZdjdbZndndmd236/nyQPZnf3ubMnOez/fMoIoiiKIiIiIiFQU4O0dICIiIiLfx9BJRERERKpj6CQiIiIi1TF0EhEREZHqGDqJiIiISHUMnURERESkOoZOIiIiIlIdQycRERERqY6hk4iIiIhUx9BJqkpISMCPP/7o7d2Q7cqVK8jMzETNmjVRs2ZN2bcrKipCZmYmGjdujBYtWuDhhx/GxYsXjdcfPXoU9957L5KSktCuXTscOHDAeN0rr7yCxo0bIyAgAOvWrTNrNysrC3FxcWjVqhVatWqFCRMmuP8giYiIvIChk3xeWVmZ7G2rVKmCiRMnYtOmTU7dR6VKlTBlyhQcPnwY+/btQ/369fHCCy8Yr8/Ozsbw4cNx5MgRTJw4EUOHDjVe17VrV6xfvx5paWlW237hhRdQUFCAgoICvP76607tFxERkVYwdJJHzJs3D23btkXr1q3Rrl07fPfdd8brBEHArFmz0L59ezRo0ABLliwxXmdZKW3Tpg1yc3NltTl37lxkZGRg0qRJaN68OXbu3Gm8PicnB4888kiF/QwKCkLXrl1RtWpVpx5fnTp1cN999xl/b9++PX7++WcAwPnz57F79248/vjjAID+/fvjxIkTOHnypHHbhg0bOnV/REREesPQSR7xf//3f/jhhx+wZ88evPnmm2aVPgAIDg7Gd999h/Xr12PUqFGyqpOO2rxx4wZyc3Px+uuvY9SoUXjnnXeM173zzjsYOXKkMg/Owu3bt/HOO+/gwQcfBACcOXMGMTExqFy5MoDyQFyvXj2cPn1aVnvz5s1DixYt8MADD6CgoECVfSYiIlJbZW/vAPmHPXv24OWXX8aFCxdQuXJlHDhwADdv3kRgYCAAYNCgQQCA5ORkVK5cGefOnUNcXJxbbT755JPGbR9//HFMnToV58+fx8GDByEIAjp16qT44xRFESNGjEDVqlXx7LPPGi8XBKHCdnK8/PLLiI6ORkBAAFavXo0ePXrg6NGjCA8PV3S/iYiI1MbQSaozGAzo378/cnNzcc8996C0tBRRUVFmATE4ONi4faVKlYyVzsqVK+P27dvG665fvw4AuHnzpsM2TYNZSEgIhgwZgkWLFmHPnj2qVTlHjRqFM2fOYM2aNQgIKO9IiI+PR2FhIcrKylC5cmWIoogzZ86gXr16DtuLjY01/r9v37544YUXcPjwYdxzzz2q7D8REZFa2L1OHnHr1i3Ex8cDAN566y3Zt2vYsKFxrOb333+Pw4cPAygPn862+cwzz+C9995DXl6esbLqrK5du+L777+3et2oUaNw7NgxrF692hh8AaB27dpo3bo1PvroIwDAypUrkZCQgISEBIf3V1hYaPz/t99+iwsXLiAxMdGlfSciIvImVjpJVWVlZQgNDcWMGTPQrl071KtXDw899JDs27/88ssYMmQIFi9ejJSUFNx9990AgMjISKfblJYeSkpKQmhoqM3tUlJSUFxcjN9//x1xcXHo3LkzPvzwQ9y+fRt79+612u3/zTff4K233kKTJk3Qvn17AECDBg2wevVqAOUTl7KysvDKK68gMjISy5YtM9721VdfxTvvvINff/0VWVlZCA4Oxp49e1CrVi1kZWXhl19+QaVKlRASEoLPPvsMUVFRsp8/IiIirRBEuYPLiJxUXFyMJk2a4Ny5cwgJCfH27uDKlSto0qQJtm3bhgYNGjh9+927d+Pdd9/FokWLVNg7IiIi38budVLFvHnzkJGRgTlz5mgicC5cuBBNmjTBiBEjXAqcQHkFlIGTiIjINax0EhEREZHqWOkkIiIilwmCgDVr1nh7NyrIyspCnz59vL0bZIKhk4iIdO38+fPIzs5GvXr1EBQUhLp16yIzM9PsLGRasWfPHgwYMAB16tRBcHAwkpKSMGzYMBw5csTbu6Yahj+SMHQSEZGu9e/fH3v37sWyZctw5MgRrF27FhkZGbh48aK3d83MunXr0KFDB9y4cQPLly/HwYMH8eGHHyIqKgpTpkzx6r7dunXLq/dP/oGhk4iIdOvSpUvYvn07Zs2ahc6dO6N+/fpo164dJk2ahF69ehm3EwQBixYtQt++fREaGopGjRph7dq1Zm3l5eWhXbt2CAoKQnR0NF544QXjiSo+//xzVK1aFQaDAQBQUFAAQRAwYcIE4+2zs7Px2GOPWd3PP/74A0888QR69uyJtWvXolu3bmjQoAHat2+POXPmICcnR9Z+5OTkIDY21rgfkoceeghDhgwx/v7555/jnnvuQXBwMO666y5Mnz7d7PTCgiBg4cKF6N27N8LCwjBz5kxZtzt69CjS0tIQHByMpk2b4uuvv5bxKtk3b948NG/eHGFhYYiPj8eIESNw5coV4/VLly5F1apVsXHjRiQnJyM8PBz3338/iouLjdvcvn0bzz33HKpWrYoaNWpg4sSJss/8Rh4kEhERKeDatWtiSUmJIj+XLl2qcNn169cr3OetW7fE8PBwccyYMVavlwAQ4+LixI8//lg8evSoOGrUKDE8PFy8cOGCKIqiWFhYKIaGhoojRowQDx48KK5evVqsWbOmOHXqVFEURfHSpUtiQECAuGvXLlEURXHBggVizZo1xbZt2xrvIykpSXzvvfes3v+qVatEAOKOHTvsPoeO9uPChQtiYGCguGnTJuNtLl68KAYGBoobN24URVEUN2zYIEZGRopLly4Vjx8/Ln711VdiQkKCOG3aNLPno3bt2uLixYvF48ePiydPnnR4u9u3b4vNmjUTMzIyxD179oh5eXli69atRQDi6tWrbT6mIUOGiL1797Z5/fz588X//e9/4s8//yxu3rxZbNy4sfj0008br1+yZIlYpUoVsVu3buIPP/wg5ufni8nJyeLAgQON28yaNUuMiooS//Of/4gHDhwQhw4dKkZERNi9X/I82bPXi4QiNbMvERGpLEaMUa3t69evo35ICM4r1F54eLhZtQsApk6dimnTplXYduXKlRg2bBiuXbuGlJQUpKen49FHH0WLFi2M2wiCgMmTJ+Oll14CAFy9ehURERFYv3497r//fvzjH//AypUrcfDgQQiCAAB499138fzzz6OkpAQBAQG45557MHDgQIwbNw59+/ZF27ZtMX36dPz222+4evUqoqOjcfDgQTRp0qTCPs6ePRvPP/88Ll68iGrVqtl83HL2o3fv3qhZsyYWL14MAHj//fcxdepUFBYWolKlSkhLS0OPHj0wadIkY7sfffQRJk6ciKKiIuPzMWbMGMyfP9+4jaPbffXVV+jZsydOnjxpPEnGhg0b0KNHD6xevdrmuM2srCxcunRJ9mSjzz77DE8//TR+++03AOWVzieeeALHjh1Dw4YNjc/JjBkzcO7cOQBATEwMRo8ejeeffx5A+YlJGjRogHvuuUeTk5z8Fc9IREREbrt58ybOA/gBQISbbV0G0PbKFZw5cwaRkZHGy4OCgqxu379/f/Tq1Qvbtm3Dzp07sWHDBsyePRuLFi1CVlaWcTvTEBoWFoaIiAicP18ekw8ePIjU1FRj0AOAjh074sqVKygsLES9evWQkZGB3NxcPPfcc9i2bRtmzpyJlStXYvv27bh06RLq1KljNXACkN3VK2c/Bg0ahOHDh+Pdd99FUFAQli9fjkcffRSVKlUCAOTn5+OHH37Ayy+/bGzj9u3buH79Ov744w/jGdnatGljdt+Obnfw4EHUq1fP7Kxsqampsh6XPVu2bMErr7yCAwcOoLS0FGVlZbh+/TquXr2KsLAwAEBoaKgxcAJAdHS08bUrKSlBcXGx2b5UrlwZbdq0YRe7xjB0EhGRYiLgfuiUREZGmoVOe4KDg9G9e3d0794dL774Ip566ilMnTrVLHRWqVLF7DaCIBjHRoqiaBb0pMuk7QAgIyMDixcvxt69exEQEICmTZsiPT0deXl5+P3335Genm5z/5KSkgAAhw4dshvU5OzHgw8+CIPBgC+++AJt27bFtm3bMG/ePOP2BoMB06dPR79+/aw+TxIp0Mm9nbUAZ7mvzjp16hR69uyJv/3tb3jppZdQvXp1bN++HUOHDjWb3GTttWOg1B9OJCIiIp/TtGlTXL161antd+zYYRZkduzYgYiICMTGxgIo736+fPkyFixYgPT0dAiCgPT0dOTm5iI3N9du6PzLX/6CmjVrYvbs2Vavv3Tpkuz9CAkJQb9+/bB8+XJ88sknSEpKwj333GPcPiUlBYcPH0ZiYmKFn4AA24d9R7dr2rQpTp8+beyiB+D2slS7du1CWVkZ5s6diw4dOiApKcmsfTmioqIQHR2Nb7/91nhZWVkZ8vPz3do3Uh4rnUREpFsXLlzAgAED8OSTT6JFixaIiIjArl27MHv2bPTu3Vt2OyNGjMCCBQvw7LPPYuTIkTh8+DCmTp2K5557zhjUoqKi0KpVK3z00Ud44403AJQH0QEDBuDWrVvIyMiw2X5YWBgWLVqEAQMG4KGHHsKoUaOQmJiI3377Df/+979x+vRpfPrpp7L2AwAGDRqEBx98ED/99BMef/xxs/t68cUX8cADDyA+Ph4DBgxAQEAA9u3bh/379xtnqVvj6HbdunVD48aNMXjwYMydOxelpaX4xz/+Iev5LSkpQUFBgdll1atXR8OGDVFWVoa33noLDz74IL755hssXLhQVpumRo8ejddeew2NGjVCcnIy5s2bZwzypB0MnUR/ih1ue5LF2fc5kY5Ii8LDw9G+fXvMnz8fx48fx61btxAfH49hw4bh73//u+x2YmNjsX79ekyYMAEtW7ZE9erVMXToUEyePNlsu86dO2P37t3GgFmtWjU0bdoURUVFSE5OtnsfvXv3xo4dO/Dqq69i4MCBKC0tRXx8PLp06WIMg3L3o0uXLqhevToOHz6MgQMHml2XmZmJdevWYcaMGZg9ezaqVKmCJk2a4KmnnrK7f45uFxAQgNWrV2Po0KFo164dEhIS8Oabb+L+++93+Pzm5uaidevWZpcNGTIES5cuxbx58zBr1ixMmjQJaWlpePXVVzF48GCHbZoaN24ciouLkZWVhYCAADz55JPo27cvSkpKnGqH1MXZ6+Qz7IVGOYQc98Ymidnmf0oMqqQ1as5eLy0tRVRUFA5BmYlETVBeHZM7ppOItI+VTtIFdwOlJ1iG1ricituYBlOGUiIi8icMnaQp3gyXYrbodrXTEdP2TUOpFEYZRImIyFcxdJLXaLF6KYU/tcOnJen+GESJiMhXMXSSR2gxYGodgygREfkShk5Shd5Dpie62l1hGUQZQomISC8YOkkReg+ZemUthDKAEhGRFjF0kst8PWhqtdppj5AjMIASEZEmMXSSU3w9aFry1sQiJVgGUIDd8ERE5D0MneSQvwVNX2TaDc8ASkRE3sDQSVYxaJrTY1e7LZYBlOGTiIg8gaGTjBg0/Q/HgBIRkacEeHsHyPtih8cwcDrgK1VOe4QcAXFiLN8LRESkClY6/RSDhXN8qXvdEVY/iYhIDax0+hlWNckZrH4SEZFSGDr9BMOme/ylymkLwycREbmL3es+jAFBOf7UvW6Padd7oXDWuztDRES6wtDpgxg2lcfAWVGcGAuA4ZOIiORh6PQhDJvkDVL45KQjIiKyh2M6fQDHa5IWcNwnERHZw0qnjvHg7hnsWneONO6T3e5ERGSKlU6dYuBUhpAjGH9sXU+uYdWTiIhMMXTqDLvSlWMaKMVs0Yt74rvY5U5ERBJ2r+sED9rKcbV6yWWT3BMnxnKyERGRH2OlU+NY2XSfFBQddaNbXsdKqLKk55DvZyIi/8RKp4bx4KwcRwFSyBFQKABxYsVQarq9J6udhq1AQJpH7kp1Qo5g9jzGDo9hxZOIyM8wdGoQw6bnidmiWeB0pbKpdEiU25Zew6n0Pmf4JCLyD4IoirKOrkUCDwxqY9hUlmVFUukucrkVT7Wro+UVWtWaV4Sc557hU30xonqfMaWlpYiKisIhABFutnUZQBMAJSUliIyMdH/niEgTWOnUCAZOZWhtoo/SQdfTVU13Aq2zj51d7kREvo2h08sYNj3DVyYCWQucWqxyuvp8s8udiMh3MXR6EQOnf9Ba9VViKxia7q8rgVaJgM+qJxGR72Ho9AKGTXUV/pmZpMCkVpXTtF13z2jkbre5M7d39Hy4MwZVyeeaVU8iIt/C0OlhDJzqkgJnTJ5n71cKW6ZhzbAVEJZXvNzabQPcrIbKCZzOBEJrj8fdNl3FqicRkW/g7HUPYuBUl7TWJlCxW9iTYzrtrQlqLcTZ2kZOJRWoWOX0lfGrlhg83cfZ60TkTax0egDDpvqEHAGGrQDyvL9mpWmV0HJRdNPrHd3e8jJbwdPbj9dT2N1ORKRvDJ0qY+BUnxQ4tRa+lF6f01qXt69WNe1hdzsRkT6xe10lDJue4cwC7VpkrRLqajuAdh+nWhg+ncPudSLypgBv74AvYuD0HL2HLCX3X+/PhSv4t0ZEpB/sXlcYD4Lqc6XLWqmKolb58mNzhN3tRET6wEqnghg41SVNzCGyxL89IiLtY6VTATzgqYtBk+Tg7HYiIm1jpdNNDJzqYWWTXMG/SSIibWLodAMPbupReqkhfx7z6I/4t0lEpD0MnS7iQU09jgKnYav1y20tqk7+iX+jRETawjGdLuDBzLsC0myHSYZMMsWZ7URE2sFKp5MYONUlt1udYz3L8XlwjH+zRETawNDpBB681CUFKFvd57a292es7MoTOzyGf79ERF7G0CkDD1jqMw2QWjuHOvkO/h0TEXkPx3Q6wIOUsiyrk2K2qNmKpbV9Jf3jOE8iIu8QRFGUdSQtEvzvQ5qBUzlqBktvhkHpcZnug7XLSHv8MXjGiOp9ppWWliIqKgqHAES42dZlAE0AlJSUIDIy0v2dIyJNYPe6DQycylG7kunNSqnlGqCm+6LW4vZarQzrDf/GiYg8i6HTCh6MlONPAcnWY5XCp7UQ6srzwwqqcvi3TkTkORzTaYEHIeU4E6isBSkhRzBe7ivh1dbjtHUdqY9jPImIPINjOk0wcCrHmZBo2AoIyx2/DR21qVZos3W/tu6PE5D0yR+CJ8d0EpE3sdL5JwZOZThbkRSzRci9hadnutu7L3tBkiFTn1jxJCJSF8d0goHTVygZSK21JU0aYqj0XfwsICJSj9+HTh5klONKldNZagc+e2GT/AM/E4iI1MHudVKEJ7u9XQmAppOSrF1n63eGTf/ErnYiIuX5dehkRUMZpudM1/IpLC1niduqahIREZHy/DZ0MnAqQy/nTDcNmgybJAernUREyvLL0MnAqQx3utS9FfIYLskZDJ5ERMrxu4lEDJxE5Ax+ZhARKcOvQicPHsqxrHIatnppR1TkK2dBIvfxs4OIyH1+Ezp50FCOPwROgF3xZI6fIURE7vGL0MmDhboC0rQ9icgRVjRJLn6WEBG5zi9CJ2mPVoKetB9a2R8iIiJf5fOhk5UJ5flKt7OcoMkwSpb4mUJE5BqfDp08OKhHieCplUBnb0yqrwRsUhY/W4iInOez63TyoED2mAbeonQAzJbkJK7haV30eiAyzL02Sq8C6KnI7hCRhvhkpZOBk5wRJ7KiSa7hZw0RkXw+GTpJP2ydltJTfHW5JyIiIq3xudDJygM5Q89LPZE28DOHiEgenwqd/PD3HK1MAlICu9bJXfzsISJyzGdCJz/0SS5fCsykHfwMIiKyz2dCJ+mL5VhKNYOgadsMnERERN7hE6GTFQb1qBXSrI2ldOa+pAlI4iB5E5FsTVhi1zopiZ9FRES26X6dTn7Ie56nqpJAxVAoDhIQkFZeKbV2znfp9mK2yKomeQXX7yQisk73oZPUo4XQZrkPUmneVti09bu99lntJCIiUp+uu9dZ5VQfAxmR8/jZRERUkW5DJz/UPcPVCqKe+OJjIu/jZxQRkTndhk5Sl78FMW+fGYmIiMjX6TJ0soLgPb7a3c7TYZIa+FlFnjJt2jS0atXK27tBZJfuQic/xD3DVrjUWzVQbpiUJibp7fGR9vEzS33nz59HdnY26tWrh6CgINStWxeZmZnYuXOnt3dNtoyMDIwZM8bbu2EmISEBCxYs8PZukA/h7HXyaZaz3KWlluwxXXaJiLSvf//+uHXrFpYtW4a77roLv/zyCzZv3oyLFy96e9cquHXrFqpUqeLt3SDyCl1VOlkx0Aa9dEVb209HgdMUx3mSUvjZpZ5Lly5h+/btmDVrFjp37oz69eujXbt2mDRpEnr16mXcrqSkBMOHD0ft2rURGRmJLl26YO/evcbrjx8/jt69e6NOnToIDw9H27ZtsWnTJof3/95776Fhw4YIDAxE48aN8eGHH5pdLwgCFi5ciN69eyMsLAwzZ8506XE+//zzSEpKQmhoKO666y5MmTIFt27dsrn9iRMnkJiYiKeffhoGgwE3b97ExIkTERsbi7CwMLRv3x65ublO7cOlS5cwfPhw1KlTB8HBwWjWrBnWrVtnvH7Hjh1IS0tDSEgI4uPjMWrUKFy9etWlx0u+STehkx/a2uFMcPMmvewn+Qd+hjmvtLTU7OfGjRsVtgkPD0d4eDjWrFlj9XoAEEURvXr1wrlz57B+/Xrk5+cjJSUFXbt2NVZDr1y5gp49e2LTpk3Ys2cPMjMz8eCDD+L06dM292/16tUYPXo0xo0bhx9//BHZ2dl44oknsGXLFrPtpk6dit69e2P//v148sknXXouIiIisHTpUhw4cABvvPEG/vnPf2L+/PlWt/3xxx/RsWNHDBgwAO+99x4CAgLwxBNP4JtvvsGnn36Kffv2YcCAAbj//vtx9OhRWfdvMBjQo0cP7NixAx999BEOHDiA1157DZUqVQIA7N+/H5mZmejXrx/27duHFStWYPv27Rg5cqRLj5d8kyCKoqw+xCLBu2fY4Ae2Z5lW+KRuZn+u+rGrnZTg7TMVxYjqfY6WlpYiKioKJeuByDA327oKRPWsePnUqVMxbdq0CpevXLkSw4YNw7Vr15CSkoL09HQ8+uijaNGiBQDgf//7H/r27Yvz588jKCjIeLvExERMnDgRw4cPt7ofd999N55++mmbwaljx464++678f777xsve/jhh3H16lV88cUXAMornWPGjLEZECUZGRlo1aqV7DGUr7/+OlasWIFdu3YBKJ9ItGbNGrz33nt44IEHMGnSJIwfPx5AeRW3UaNGKCwsREzMnfdAt27d0K5dO7zyyitW7yMhIQFjxozBmDFj8NVXX6FHjx44ePAgkpKSKmw7ePBghISEICcnx3jZ9u3bkZ6ejqtXryI4OFjW4yLfposxnQycROQLeIpM55w5cwaRkZHG300Do6n+/fujV69e2LZtG3bu3IkNGzZg9uzZWLRoEbKyspCfn48rV66gRo0aZre7du0ajh8/DgC4evUqpk+fjnXr1qGoqAhlZWW4du2a3UrnwYMHKwTWjh074o033jC7rE2bNk49bmv+85//YMGCBTh27BiuXLmCsrIys+cGAE6fPo1u3bph5syZGDt2rPHy3bt3QxTFCmHxxo0bFZ4TWwoKChAXF2c1cAJAfn4+jh07huXLlxsvE0URBoMBJ06cQHJystyHSj5MF6GTPMufK5rWsMpJ5B2RkZEVgpUtwcHB6N69O7p3744XX3wRTz31FKZOnYqsrCwYDAZER0dbHcNYtWpVAMCECROwceNGzJkzB4mJiQgJCcFf//pX3Lx50+79CoL556UoihUuCwtzr/T77bff4tFHH8X06dORmZmJqKgofPrpp5g7d67ZdrVq1UJMTAw+/fRTDB061PjcGQwGVKpUCfn5+cbucEl4eLisfQgJCbF7vcFgQHZ2NkaNGlXhunr16sm6D/J9mg+drHJ6ni90pxcKQByzImkQq52e0bRpU6xZswYAkJKSgnPnzqFy5cpISEiwuv22bduQlZWFvn37Aigf43ny5Em795GcnIzt27dj8ODBxst27NiheFXvm2++Qf369fGPf/zDeNmpU6cqbBcSEoJ169ahZ8+eyMzMxFdffYWIiAi0bt0at2/fxvnz59GpUyeX9qFFixYoLCzEkSNHrFY7U1JS8NNPPyExMdGl9sk/6GYiEZEzlAyceg7fpE38Mq2cCxcuoEuXLvjoo4+wb98+nDhxAp999hlmz56N3r17Aygfu5iamoo+ffpg48aNOHnyJHbs2IHJkycbx0QmJiZi1apVKCgowN69ezFw4EAYDAa79z1hwgQsXboUCxcuxNGjRzFv3jysWrXKOJbSWb/++isKCgrMfs6dO4fExEScPn0an376KY4fP44333wTq1evttpGWFgYvvjiC1SuXBk9evTAlStXkJSUhEGDBmHw4MFYtWoVTpw4gR9++AGzZs3C+vXrZe1beno60tLS0L9/f3z99dc4ceIEvvzyS2zYsAFA+ez6nTt34plnnkFBQQGOHj2KtWvX4tlnn3XpuSDfpOnQyQ9m8jZ2rRNpW3h4ONq3b4/58+cjLS0NzZo1w5QpUzBs2DC8/fbbAMq7wNevX4+0tDQ8+eSTSEpKwqOPPoqTJ0+iTp06AID58+ejWrVquPfee/Hggw8iMzMTKSkpdu+7T58+eOONN/D666/j7rvvRk5ODpYsWYKMjAyXHsvHH3+M1q1bm/1Iyy2NHTsWI0eORKtWrbBjxw5MmTLF7nPy5ZdfQhRF9OzZE1evXsWSJUswePBgjBs3Do0bN8ZDDz2E7777DvHx8bL3b+XKlWjbti0ee+wxNG3aFBMnTsTt27cBlFdC8/LycPToUXTq1AmtW7fGlClTEB0d7dJzQb5Js7PXGTi9T8gRzEKXP1b8xGyxwvNApBRPd7PrbfZ6SUmJ7DGdRKR9mq50kvdIAdPa0kn+xNrzQERERM7TZOhkldP7rAVMfw9e/v74SXn8rCMif6LJ0Enaw8BFRERE7tBc6OQ3f+1h4LyDzwUpjZ95ROQvNBc6iTzBsLV8Lc9Cofz/zmDwJCIicp6mQie/8WuLL0wckoKlPUXpzodPIUdg+CTF8LOPiPyBpkInaZej4KY3AWlATJ75IvJF6ebB09kKKBEREdmmmdDJb/rapvdTSloLzQFpd06XKT0+KXj6Wsgm7eNnIBH5Os2ETvINhq3arBBahmZrwdI0eJr+a3obIiIico0mQie/4euTtRAWkFb+o5WAJgVJy/2xFSjjxDvd7jF55tsEpFm/DZFS+FlIRL5ME6GT9MkyhMm9zhssQ6Zp5bNQML9e2ndHjyEgjTPZiYiI5PJ66OQ3e+3ytUBlOkNdyTGbljPZhRwB4iDfeu7Ic/iZSES+yuuhk7RJS4HT3W5syxnqpoFT6kKPE92/H+k5E7NFzVV6iYiIvI2hkzTPlQBnGSClcZrSj/R7QJp7M/Mt78ey4qml8E76wWonEfkir4ZOfrCSI0pM1jFsvVPhlLrYrQVZVyYKWbuNZdBk8CQiImKlk1SkRGCUU+U0vR8pWJqGTMtJRPbalNqyFUyt3aettk3P6MTgSURE/s5roZNVTt9nGfZc5Si8mt6PtaWOTNm7DnBu5rpcUvhk8CRn8DOSiHwNK51klRLnXXcU8ORyFPCsradpK3zKrV46+39TYrZofP5Mn0dfOJc9ERGRq7wSOvkN3j94ewa36SQh01Nd2tve8v+W3eyW/5cbJBk4yRX8rCQiX8JKJ+meFASVXHtTqmJaO7uS6ak+hRyBgZKIiEgGhk6ySS9hSjqPuukpL92dxGSr0ildZxpGpeBp2q1uC8d1apflcldERKQsj4dOdheRM1xZtqgo3fwyV8aWmlY6LffD8jpr7IUWaf1OruOpDdLrIC13pbXqNT8zichXVPb2DhDZI3dcqBQsDVuBAJOKp1SVdMQySDqqdFpeZ3o2ItN/Hd2PtduSumwFfb4mRETq8mjo5Dd2UkqhUB40i9LvnFkI+HP5pD9zgtzAalnFtDZxSG5bloHGMrTYaofhRptMX09vvkaxw2Nw9v0ir90/EZESOKaT7JJ7oHVmDKUSi8YD5QHOsuvc1gx1ufdprTvdHXK6zxk4Pced4Qz2hkRwqAQRkWPsXidFOLM8kpJLKUmVTdMKpWW10tHZhey1a0lqy5k2LccIMmR6h5Kh0FZbWhsPSkSkJR6rdLJrXb+8cRB1VGm0tiC8La6EXGsVU1tBU4nZ8qQuX6hC8jOUiPSO3eski6eDp9zzozu6rZJh0FbXu5zJSux+9R5PP+98rYmIrGPoJNk8GTztnW6yKP3Oj+VtLAOps1VOuUs0WQZQVjq1RRzk/SWpvH3/RERa45HQyW4h3+Fq8HQ2lNk69SRwp2vdtIvd2XGbUkC1tl9yLwPuBN+idOtnLiJ1WYY6Keh5+xSsppQMnvwsJSI9Y6WTnOZK8FQiBFhbSzNONK82OrofueMvTdu0d3pN03O721rP07RNTjJRjhTm9LDQvtb3j4jIEzh7nVwiZoseP4iahjjTiT7WgqatCT/O3ldAGqzOYjdtU874U1KWXgMcF50nIn+meqWT3UG+y96B05kuantsje10tn1379sWe7PZTW8vLBchLGfQUIJeA6cpdx4DP1OJSK/YvU5usRU85VYfHbE2WcgWpZczktNVL+cyhk2yxhfCMxGRMxg6yW2OFj53NvSZbi91o1ueR9109nihoPxambbaMr3c2phNSwycymJQIyLSL4ZOUoS9rnZnK5xyqqSWwbQoXV4IlHv/jsZpOmpbztqdRK6GaHaxE5EeqTqRiB+M/kfpSlShUL40krUAVygAASazx00pVfWUOyPeEieKKI9VTiIifWOlkzTN8nSXkqL0iut0agW71JWnduDU0vuHiMhXMXSSplk7+w9wp/pp68xD7nZtW2tXzsQhBk73WAuX9gKnUmGRQyGIiNSnWuhk17p/cjUE2DudpLQ4u+VZiky7vpWuVMkZH8rqmPIshyU4qnAqGRb19nryM5aI9IaVTlKUqyHAcla6Kcvxmtbuy7Tq6e6+WJK7XiernJ6hVjj0RrWTY3+JyJ8wdJLmSUsiWTsdpbsVSEeLu8sNnAwPypAzdtOXusI5OYqI/AlDJylGjeBl2HpnwpDlpCKpe91a97vp9Y7al0MaW2qtPQZOZQg5gtlz6YmxnJ5g6wsTEZG/USV0cqyR/1I6gJmGPVuThWLyzKuS0lmMLLvfrbVt6zpr20n34UuVNi2RGzgB/bwGjk7d6i5+1hKRmrZt24bHH38cqampOHv2LADgww8/xPbt211qj5VOUo3aB1xTAWnlYdO0Mmp6nfSvo1nu9hZ1t3Udu0iV5+nqsZrjRAPSyscl6yUoExEBwMqVK5GZmYmQkBDs2bMHN27cAABcvnwZr7zyikttMnSSahxVGl1hGg5szWo33c5amLBcBslaODatrlq24eh3e6SAyqCqLQyERETmZs6ciYULF+Kf//wnqlSpYrz83nvvxe7du11qU9UzEhG5w1FglLrVpctszXJ3pn3Ty6XxogGi49s4YhoyxUECsBXAIIEz3m3wp1BuOZaViEgLDh8+jLS0it/IIyMjcenSJZfaVLzSyTFG/kuJoGAa6orSzX+AOzPZJbbGajoarym3shUnVqyaOlMVE3KECs8Lz8tOSuJnLhGpITo6GseOHatw+fbt23HXXXe51Ca710lxYrbocuXGtHs8Ju/O+ExpTJxpNVMaw+ksa13o9raTAq/p/snhaPY1q5wk8afKLhHpQ3Z2NkaPHo3vvvsOgiCgqKgIy5cvx/jx4zFixAiX2mT3OilGzBYV6yo0DXeFFtcVpQOGPAB59kOgrZnmlmczsnW9RAq6rE56hqsBjCsLaENUR0CIdK8NsVSZfSEi102cOBElJSXo3Lkzrl+/jrS0NAQFBWH8+PEYOXKkS20KoijKSghFQpGsBtnVQ4CylRtrYcJWV3ehUHFsp5zbW3bZ2wsw1qqrUtVSzuNmldM2f6/4ufqF7ez78j6fY0T1Pp9LS0sRFRUFlCgUOqOAkpISREa62RgRueWPP/7AgQMHYDAY0LRpU4SHh7vcFiudpBqp8ukOe9VKa+IsJv0UpQPIcy5E8jzrnqXVoMnKKREREBoaijp16kAQBLcCJ6DwmE5frHJq9YCodUrNxrV3mkpbTK+PsdIFL/fUls5sD8h/rzDI3KHlvy9Pv07u/M344mcvEXlXWVkZpkyZgqioKCQkJKB+/fqIiorC5MmTcevWLZfa5EQiO7imonvcHd9pa3F5aSF4R7exdp2jAClnAXlrhBxBdjjl8jjl+HdFRKRdI0eOxPvvv4/Zs2djz5492LNnD2bPno3Fixfj2WefdalNdq/bwXDgPneChbXxlhLLbnR7tzPtJrW23JKt+3HmvOxytzVsBRi1tE/JrnV20xORHn3yySf49NNP0aNHD+NlLVq0QL169fDoo49i4cKFTrfJSiepRonQXijcWaez0CKtyalcOmJrzU9n2hWzRdkTgxg+9IGvExH5u+DgYCQkJFS4PCEhAYGBgS61ydBJHqPGZBzL86rLnSxkbfF4e6fOtEeq5soJKoat/tOtbO9x+stzANyphFt+abL8nYhIS5555hm89NJLxnOuA8CNGzfw8ssvu7xkErvXncDT1bnHlepRnFh+cLY2IcjW0kaWl0nbqTkD3ZlJRL7axe5PQdJZAWnla8uavQdNVlXg5woRac2ePXuwefNmxMXFoWXLlgCAvXv34ubNm+jatSv69etn3HbVqlWy2lQsdPrD7EkeGJzn7rJJhq3m5z6vcF1axf9bYznuUtre2r8u7aOTgdqXvsA4+/r6Yzi1PL2rdKYtJb+ExA6Pkb1eJxGRI1WrVkX//v3NLouPj3erTVY6neBLQUHvrAVIe9tamzRk2Z3ubOB0Z9yfYSuAQYIuq1zOhEbLvxm9BE6lJ/+YtSVWvE5I42cLEWnLkiVLFG+TYzqdwIOC59mbwS5dbqvb3db2XOjdda6ERiFHMP7ohZ7W6CQiUsO0adNw6tQpRdtk6CRdsFahlHsb6XbWrnc0AUlO2+7Qa9hwNriL2aLxh8zxSxARadHnn3+Ohg0bomvXrvj4449x/fp1t9tk6CSfYq3b3BrLyUXOHPiVrIKJg/RT/QPuhGRHFWi98+Tj4vJMRKRF+fn52L17N1q0aIGxY8ciOjoaTz/9NH744QeX22ToJFWp2aUqVSjlVj/lLqfkShvOsNWGaTe0XrqifTUw+erjIiJyRosWLTB//nycPXsWH3zwAc6ePYuOHTuiefPmeOONN1BSUuJUewydpHmOqpH2usetdbGbjvm0171uerkrXfBy2GtTjeDprUBrep/sYici0heDwYCbN2/ixo0bEEUR1atXx3vvvYf4+HisWLFCdjsMnaRp9gKnvdNWOprNbo219T4dteUu0wCmVhhTOmgyNGqXPyxdR0Sek5+fj5EjRyI6Ohpjx45F69atcfDgQeTl5eHQoUOYOnUqRo0aJbs9RUInP+hIDXK6vW1NEJLbhq2uebUqm5b34QlSSNRLl72/8dVxsUSkT5UqVcL58+fRokULdOjQASdOnMDixYtx5swZvPbaa0hMTDRuO3jwYPz666+y22alk1ThyYBjGRztjfG07DL3JGv7aUmNKqKnw6ajx6X0/ngjtCm59FZRujLtEBEpQRTLP68HDBiAkydP4osvvkCfPn1QqVKlCtvWqlULBoNBdtsMnaQ4KVR4q7LmKFha64a311Wv9H55q7Jla9kiy9dJzutmLxxbPudqd8fL/fKgxPMuhU0lK+FxHK1ARBo0ZcoUxMbGKtomz0hEbpNCiuUpL909Baar7FUTHQUPpUOhvdNrCsutpw2lQ5q918He5Ursh5bGfzobEk1fO9NTVxIR+bqNGzciKirK7jYPPfSQ0+0ydJJbLMOJ6fhBPYwhtDWeU+ngadm+mmHM2vPuSve26ZcJR+3bui89M6uWu/CQnJmIxlPsEpGWDBkyxO71giDg9u3bTrfL0EluMa2iKVnZ9EQXtCdmp1ujt3DhzGvq6LFp/YuI5WlV3T0HuzO3ZfAkIq04d+4cateurXi7HNNJbjENEWff13agsMWTYyxtdak71Yad4GYruMhddF7roVBt1qrcrk4aciWs6u3kAETkewRBxZO6qNYy+TRrB0Y9j3mz1gXubhhVK8zaq4ZZDnVwtnLGSpt5WLR38gC1MXgSkTdIs9fVwO51cpq1sX5CjqD7UwfaWoDelQko3mZr0Xk1K52W98MA63lCjgC8f9bbu0FEOjZkyBCEhISo0jYrneQUW5NLfJncReoddcMq0bXuLsvqp1rrZ+r1/aGFLwzu0OvzTkTasWTJEkRERKjSNiud5DS1D2yFAgCNdtW7O7FEKyyDp6uB0xdCjumSSEXpQOGfl3P9TCIiZbHSSbLZm6SiJC0f7C0Xk5c7yUQLVU4t0GJIDUi7Mx45Jq/8/afl9yARkV4xdJIs9sbomXbZ2hoX6QwtdnE6On2lPVoPnJywYj5ZSIvvPyIiX8DQSYqxd7B2J7RphStL5/h64HT29loPuL4wdIKISEnHjh3Dxo0bce3aNQDuzW5n6CRZ5HaLFqXLC5jOXq5FetpXS0qtBanF7nJ3aTV4yn29YofHqLwnROQPLly4gG7duiEpKQk9e/ZEcXExAOCpp57CuHHjXGqToZPscjaYSGPj7FUFrXVjuroAtxYFpJVXONWucpouJG75o1XeCqmuvre08p7kElTalJWVhT59+ijWniAIWLNmjUfvk8iWsWPHonLlyjh9+jRCQ0ONlz/yyCPYsGGDS226HTr5rdq3KXmgM518Y3mZr1D7vOpaYutxajn0OsvbVU/pS4S/vKdc5StBrLi4GD169AAAnDx5EoIgoKCgwCv7Ivc5PX/+PLKzs1GvXj0EBQWhbt26yMzMxM6dO9XfSVLVV199hVmzZiEuLs7s8kaNGuHUqVMutcklk0gxUmVPHGQ7dCgx0UjLPBUOxEGC7G+M7qytKneJKMtgZBo81VoLlMqZPvfG55cLxOtS3bp1vb0LTuvfvz9u3bqFZcuW4a677sIvv/yCzZs34+LFi97eNXLT1atXzSqckt9++w1BQUEutcnudVKcsFw0djGT8hyd/clWyHOl+11upc9aoJV7Gk5PfenQ85cby9fM8rXU+rAKb5s3bx6aN2+OsLAwxMfHY8SIEbhy5Yrx+qVLl6Jq1arYuHEjkpOTER4ejvvvv984hg0Abt++jeeeew5Vq1ZFjRo1MHHiRLsTKkRRRK1atbBy5UrjZa1atULt2rWNv+/cuRNVqlQx7otp93qDBg0AAK1bt4YgCMjIyDBrf86cOYiOjkaNGjXwzDPP4NatW8brfv/9dwwePBjVqlVDaGgoevTogaNHjxqvnzZtGlq1amXW3oIFC5CQkGC8ftmyZfjvf/8LQRAgCAJyc3MrPMZLly5h+/btmDVrFjp37oz69eujXbt2mDRpEnr16mXcThAELFq0CH379kVoaCgaNWqEtWvXmrWVl5eHdu3aISgoCNHR0XjhhRdQVlYGAPj8889RtWpVGAwGAEBBQQEEQcCECROMt8/OzsZjjz1WYR/JdWlpafjXv/5l/F0QBBgMBrz++uvo3LmzS20ydJIq2B2oDnvBwlqokl4HKQC6+7qY3odp4LG8T3vLa1lSYkhCoWA/VHrr/OlK0sOYXaWVlpaa/dy4ccOldgICAvDmm2/ixx9/xLJly/C///0PEydONNvmjz/+wJw5c/Dhhx9i69atOH36NMaPH2+8fu7cufjggw+wePFibN++HRcvXsTq1att3qcgCEhLSzOGtd9//x0HDhzArVu3cODAAQBAbm4u7rnnHoSHh1e4/ffffw8A2LRpE4qLi7Fq1SrjdVu2bMHx48exZcsWLFu2DEuXLsXSpUuN12dlZWHXrl1Yu3Ytdu7cCVEU0bNnT7Ngas/48ePx8MMPG4N3cXEx7r333grbhYeHIzw8HGvWrHH42kyfPh0PP/ww9u3bh549e2LQoEHGaujZs2fRs2dPtG3bFnv37sV7772HxYsXY+bMmQDKw8/ly5exZ88eAOUBtWbNmsjLu3MWkdzcXKSnp8t6fCTP66+/jpycHPTo0QM3b97ExIkT0axZM2zduhWzZs1yqU2GTlKVZcXTFw7+tgSkVawsKslR4AxIK//XWghUal8sXztX2jUNpZb/d1Wc6LvvK38WHx+PqKgo48+rr77qUjtjxoxB586d0aBBA3Tp0gUvvfQS/v3vf5ttc+vWLSxcuBBt2rRBSkoKRo4cic2bNxuvX7BgASZNmoT+/fsjOTkZCxcuRFRUlN37zcjIMIbOrVu3omXLlujSpYvxstzc3AoVTEmtWrUAADVq1EDdunVRvXp143XVqlXD22+/jSZNmuCBBx5Ar169jPt69OhRrF27FosWLUKnTp3QsmVLLF++HGfPnnU4SUkSHh6OkJAQ4xjNunXrIjAwsMJ2lStXxtKlS7Fs2TJUrVoVHTt2xN///nfs27evwrZZWVl47LHHkJiYiFdeeQVXr141But3330X8fHxxsfUp08fTJ8+HXPnzoXBYEBUVBRatWpl9ryNHTsWe/fuxeXLl3Hu3DkcOXLE5nNJrmnatCn27duHdu3aoXv37rh69Sr69euHPXv2oGHDhi61ydBJqpKChC+HTYlUaZOqUYat9se3yiEOEiAOclzdMn1uCwVU6HpVm6vhk8ieM2fOoKSkxPgzadIkl9rZsmULunfvjtjYWERERGDw4MG4cOECrl69atwmNDTU7EAaHR2N8+fPAwBKSkpQXFyM1NRU4/WVK1dGmzZt7N5vRkYGfvrpJ/z222/Iy8tDRkYGMjIykJeXh7KyMuzYscOl6tzdd9+NSpUqWd3XgwcPonLlymjfvr3x+ho1aqBx48Y4ePCg0/flSP/+/VFUVIS1a9ciMzMTubm5SElJMau8AkCLFi2M/w8LC0NERITZPqempkIQ7nyOdOzYEVeuXEFhYfmJaaUAL4oitm3bht69e6NZs2bYvn07tmzZgjp16qBJkyaKPz5/V7duXUyfPh3r1q3D+vXrMXPmTERHR7vcntuh8+z7Re42QX7AMmD4Wgg1XQbKdPknZ8e1Wo7Rk/M8ma4KIJ1D3DQAK0mpbnoiOSIjI81+XJm8cOrUKfTs2RPNmjXDypUrkZ+fj3feeQcAzLqbq1SpYnY7QRDcWgQbAJo1a4YaNWogLy/PGDrT09ORl5eHH374AdeuXcN9993ndLvW9lUa72hrn0VRNIa6gICACtvJ7Xq3Jjg4GN27d8eLL76IHTt2ICsrC1OnTnVqn00Dp+njkC7PyMjAtm3bsHfvXgQEBKBp06bG55Jd6+pYsmQJPvvsswqXf/bZZ1i2bJlLbbLSSRWoVRmzNubOWqDSYxi1NZ7QmefSsjLp7MSXgLQ75wy3fA5NQ6JSgdGyHXuPVe7zwDBLStu1axfKysowd+5cdOjQAUlJSSgqcq5YEhUVhejoaHz77bfGy8rKypCfn2/3dtK4zv/+97/48ccf0alTJzRv3tzYlZ+SkoKIiAirt5W6s2/fvu3UvjZt2hRlZWX47rvvjJdduHABR44cQXJyMoDyrvtz586ZBU/LpZkCAwOdvm/TfTCtIsvZfseOHWb7s2PHDkRERCA2NhbAnXGdCxYsQHp6OgRBQHp6OnJzcxk6VfLaa6+hZs2aFS6vXbs2XnnlFZfaZOgkj7MMU3oMmXLJDY6WoUzuUkXAnbGc0n3ZC5zW7ssZ1saLunpbV7ch+/w1uJeUlKCgoMDs5/Tp02jYsCHKysrw1ltv4eeff8aHH36IhQsXOt3+6NGj8dprr2H16tU4dOgQRowYgUuXLjm8XUZGBj7++GO0aNECkZGRxiC6fPlyu2MQa9eujZCQEGzYsAG//PILSkpKZO1no0aN0Lt3bwwbNgzbt2/H3r178fjjjyM2Nha9e/c27tOvv/6K2bNn4/jx43jnnXfw5ZdfmrWTkJCAffv24fDhw/jtt9+sVkIvXLiALl264KOPPsK+fftw4sQJfPbZZ5g9e7bxvuQYMWIEzpw5g2effRaHDh3Cf//7X0ydOhXPPfccAgLKY4o0rvOjjz4yPm9paWnYvXs3x3Oq5NSpU8ZVFEzVr18fp0+fdqlNhk6qQI2DljS20dFZiuxtoydSt7iry0Y5WhLJ8sd0e8vrTCkR6uwFT3fb99fARO7Lzc1F69atzX5efPFFtGrVCvPmzcOsWbPQrFkzLF++3KUJSePGjcPgwYORlZWF1NRUREREoG/fvg5v17lzZ9y+fdssFKWnp+P27dt2q3OVK1fGm2++iZycHMTExDgV4pYsWYJ77rkHDzzwAFJTUyGKItavX2/s4k5OTsa7776Ld955By1btsT3339vNlMfAIYNG4bGjRujTZs2qFWrFr755psK9xMeHo727dtj/vz5SEtLQ7NmzTBlyhQMGzYMb7/9tuz9jY2Nxfr16/H999+jZcuW+Nvf/oahQ4di8uTJZttZPpfVqlVD06ZNUatWLWMVl5RTu3Ztq5PC9u7dixo1arjUpiDKHLRSJNjujuBZicge08k0cscn+gK5gdOZoOZuKFOykmgv0DpznSU9VDsLhTtDGbTWXqFge2H4GFG9z+rS0tLy2dwlgBDpXltiKYCo8uplZKSbjRGRSyZOnIh///vfWLJkCdLSyg/eeXl5ePLJJ/HXv/4Vc+bMcbpNnpGIVCOFB6mc7kuBUimOApZhq3KL7KuxhJPlWFFb9+GLFcxCoXzSlhLDQ2LyHG9DRORJM2fOxKlTp9C1a1dUrlweFw0GAwYPHswxnaQ91mas+wulH6vlrHZXFgp3Nfi5+mXB1ZCrp9nxRenKfJnyp78NItKHwMBArFixAocOHcLy5cuxatUqHD9+HB988IHVdVvlYOgkTfHXg2/hn/nM2iQraXiC5ax2R8siWTt1oivsvSbW2tRLYHRHTN6d6mQRJ80SkQ9LSkrCgAED8MADD6B+/fputcXudVKNOEhw6luNP3a/S13Slksdmc5Et5xkJc1Wt2Ta3W1rYXhnZsW7Qo2lmbQoIO3OFwV2jRORryosLMTatWtx+vRp3Lx50+y6efPmOd0eQydphmXg0islziVuLRhKi78b8u4sAg9YH0tpa3ylGoFTHCQoNu5UT+LE8uDpr9V5IvJtmzdvxkMPPYQGDRrg8OHDaNasGU6ePAlRFJGSkuJSm+xeJ1VIZ9Mh+UzDdqEgbz1Ta7PEnR3r6S61Xmc9VEpZ5SQiXzVp0iSMGzcOP/74I4KDg7Fy5UqcOXMG6enpGDBggEttMnSSKtypVnozrLpy30pU+Syfrzix4r5I20hjCLUU6r0RdrVAS68BEZGSDh48iCFDhgAoXzf22rVrCA8Px4wZMzBr1iyX2mToJFX448FYWC5CWO7azGvTsZumZxeSSL9LYz9NK2xCjmD1Np4mjXH0hfDp7eeSiMjbwsLCcOPGDQBATEwMjh8/brzut99+c6lNjukkTfPkOE9Xg7I4qHwogZAmKNIlbK3CaXlZUToA0WQtVA2EfMvFzS3X8XSFvbU/1SSNndXC80pE5A0dOnTAN998g6ZNm6JXr14YN24c9u/fj1WrVqFDhw4utcnQSapwNyy4crB3J6C6elvL/ZQeszthS8wWjcskWXsepAksSp8VRw3OBk/L589bFVOlFn0nItKrefPm4cqVKwCAadOm4cqVK1ixYgUSExMxf/58l9pk9zrpgpwA4G411JXbS93aprc1bJXXxWwvoJo+Xml5HtP7iMkzD5yFGunNtvZY5Ha3m27j7S56pRZ9JyLSkzfffBPXr18HUD6Os3nz5gCA0NBQvPvuu9i3bx9WrVrl8nqdioTOs+/bPi87+S9b4xJduS1wZ9yjZQC1XMfSFUoEDGvjKu0FJ8NW8/PS27tNnGgehCyfAy1VPPUw69warQR3IiJvee6551BaWgoAaNCgAX799VdF22f3OqnOVlByl7RIuruB0dZi6+4wXcTdNMuYdjdL14uDBMDKuM0Ki8L/2a0ODWY6y6ApDa+QG0C1MPFImpwlPdd6VCic9fYuEJGOxcTEYOXKlejZsydEUURhYaGx8mmpXr16TrfP0EmqcSdk2rqtFOak610NjNYCnlrBE8vNL7c8naWc50naJk5U/6xC/sr0PUVE5I8mT56MZ599FiNHjoQgCGjbtm2FbURRhCAIuH37ttPtM3SSKiwrV1Koc/eAbmvtSiUo0U1vi1T1M30OXH0upJnr7vBEcHWmm91bs9SJiOiO4cOH47HHHsOpU6fQokULbNq0CTVq1FCsfYZOUpS1c31b/qtk2LFVoXS3cqnkUk1S9zkGKXOWJmn8pjvPpRqBU053uhLLKBERkXoiIiKQnJyMDz74AMnJyYiOjlasbc5eJ8V4q1Jl6zzl9oKV5WxzW5NIHIUzaxObTK9TYwa05ZmJPEHMFiv8mF4nvfaW/xIRkf5UqlQJf/vb32yO53QVK52kCG+HDNOAZ7nckDXSBCQpNBZaud5aO7ZCpL2AqcYEKqDijHVXKp+21vqUU4201SXu7fcCERG5r3nz5vj555/RoEEDxdpkpZNUYznpB1B3HKG9qqMjpqeVlNqy9n/T+7EMmtZu48lJKd6YAGMvnFoLn/a2Z7c7EZF2vPzyyxg/fjzWrVuH4uJilJaWmv24QrFK59n3ixA7PEap5khH7FW2jEHCynqU3mRveSJbTJdBsryNv84oV7KqyQqpfNbeb+V/a1wzmYiUcf/99wMAHnroIQjCnc9nzl4nr7IVuKwtF+QLLCuatrri9cDVrnWAIdGbrL3HeJIOIlLSli1bFG+ToZPcIuQINsdoWIYXPYUxa+Qs++TJx6hWdVWJ88eTevy1qk5EnpWervxsVYZOUkWFM/EsFyH+eYmeD5hqrhPqLDWeR8uZ6LaCp60qp+XQAy6RpDw9//0QkX5s3Wr/AJeW5vyHEUMnqUJYXjFomM4A95UDpy89JikcOjqFpa3AWSgAyDN/Hhg4iYj0KSMjo8JlpmM7OaaTPMZW8JA1O1ljk4rskbPep6OZ7Hpg+rq5EjiBiuNDGTiV5wtfbohIH37//Xez32/duoU9e/ZgypQpePnll11qU9HQyRns/suZgKHUKTGt8fRB2dGkIj2wHArhbFjU42O2R8uPx9Z+FQpnPbsjROTzoqKiKlzWvXt3BAUFYezYscjPz3e6TVY6yWWuVrLEbBEBOYIqB3dvt6fVsGKL9BoIaSbLYTj5uirxmLUyE17LgZOISAtq1aqFw4cPu3Rbhk5yibtdp2K2CAxSJ3j6g0KhfEF7d587LTz3WgmcgDaeDyIiLdi3b5/Z76Ioori4GK+99hpatmzpUpsMneQ0pWYkS93RegientxH00XobZ5X/s+n39ZpLJ25H7mcCYfS+8PR8ktaCpxERHRHq1atIAgCRNH887tDhw744IMPXGqToZOcptQEEW8uN+QsT48TtRcITS+3PH2nKcs2pOe7KL08qNpq390vFZa3dTS5zNvBUw9feoiIPO3EiRNmvwcEBKBWrVoIDg52uU3FQycnExG5T24IciYsGbdVoEoqlx5msDvzHGoxoGpuEtHoEoiBke61cbMUQMVJDETkOfXr11e8TVsnkyEinfNGOHI2ZHq7yuksrQVOIiKlfffdd/jyyy/NLvvXv/6FBg0aoHbt2hg+fDhu3LjhUtsMnUR+ypkqp9xwKOQILgdJPVRFiYh83bRp08wmEe3fvx9Dhw5Ft27d8MILL+Dzzz/Hq6++6lLbHNNJHsNzeuuTKyFSzskDLLeRxpLqrfpJRORLCgoK8NJLLxl///TTT9G+fXv885//BADEx8dj6tSpmDZtmtNtqxI6Oa6TTAl/rskJaeKQjs5I5K/UCn6O2mXwdE55iC/y9m4QkQ/5/fffUadOHePveXl5uP/++42/t23bFmfOnHGpbXavk2qEHAEi1+L0Cm+sDCD3Ph1tp9SSXP7g7PsMnESkrDp16hhnrt+8eRO7d+9Gamqq8frLly+jSpUqLrXN0EmqkMImwMDpDbaec2ldVE/epzPbidmiS6fhJCIiZdx///144YUXsG3bNkyaNAmhoaHo1KmT8fp9+/ahYcOGLrXNMZ2kONGk+1wKGJYBlEHBswr/fEmkulgh1F8uyV32utkdLaBPRESumTlzJvr164f09HSEh4dj2bJlCAwMNF7/wQcf4C9/+YtLbQui5VLzNhQJznfjcFynb5HT7WkrcNoLBgwO6il0MDTSmeDpiXU9gYoTzUyDJ98rd7iyPmeMqN5ncmlpKaKiooCsEkCJdTqXRqGkpASRkW62RUQuKSkpQXh4OCpVqmR2+cWLFxEeHm4WROVi9zrJJjdwBqRZr2jaCkAMEeqQU0229ZpY3taw1f7Zj9Rk+r7je4WIyDOioqIqBE4AqF69ukuBE2DoJIUIOQIC0gBh+Z2AYBlctN6d62uK0l2/rWW4M/0ioTZrXeqcWEREpH+qhk7OrPQPxiWRbGB1yjtMK5P2Ar9ptVMrY225ZBIRke9hpZPcJgUV08AiVaYYOJXjbCAMSCsPm1LgNP2/vdtY42hsqCe4MrPdV2nufOtERDIwdJIipO5XrrGoHMuQ6c0AH5NnP3iqUSGVc1YjIiLSD4ZOcoq1IBCQdicIcNKHOXfWxfTk8+eokilVTT2N3exERL5D9dDJcZ2+z1rgpDv0Er4LhfIfV0Iy11/1HHatE5FesdJJTtNb9cmb4xE9FTiVDHtF6Y6fM1v3586MeVs4m52IyDcwdJJLhBzB+OMPPFHBM2x1PiBL+yUFRVfDorNshWlPdsEzeBIR6YtHQie72H2Hvx7oPVGxdHbcpOXZeWLyHC/gbrm93vnb+5Fd60SkZ6x0klv0ctDXwpI/SrMMkEXp8sNxnKifsaaA/SEdXEqJiEgfGDrJ59mqHvrSpJeANNuVS3tnhtJTtdPRUA4GTyIibfNY6GQXu+/Qa2XJstqpxqQXb7JVubT3OPVU7QT8O3iya52I9I6VTnKZng7wUnXPNHj6y7ngHT1Oa9dr+bnx1OQ1X6qEExFpAUMn+Q0pSDFMVCQ9NzF52g6cEk8ET29UgaUVDCyr8nr6gkdEZItHQye72PXN3WWS1JrM40yIjBPLu5stb8MgWh449dbd7ksKhfL3prXgz89OIvIFrHSSLEpUlrRSQbM2c9vZ8Z32QqpeA6yvBE69VgXjRP2tKkBE5AyPh05+Y9c/rR3UlThIy5nFbVqptXefzuyP1J2q16BqSSuTzLSwD0rxpcdCRP6NlU5yijMHQD0dLOUERUeVWtPgKGcogbS9tUlOeuTp11vpmexaDf78ok5EvoKhkxxyp2tdK5Uvd8kJJHLO9mPYeucnIO3ObfTepWr5GnvqNVcyeFob60tERMrxSujkN3f9UGqWsN6Dp7Oh0HJ7qYopBU1r7elpoXZf5M2JVNZmrJdfzrU5ich3sNJJXqW3ypKtLnB7j6NQkDeJSs/VTk+tnanmfWvt+df7FzUiIkteC52sdmqfr66F6A5b4VEPj0NvAd+fxIkVK938jCQiX8NKJznNm1UtPbJX5fRkEFQ7GGv9faHV0C3tl+nrwyonEfkir4ZOfpPXNmsHPiFH4AFRQXqokLrKk+8TOYFXq8+1VveLiEhprHQSeZGt0x7qkbWQqYUvKKb74Mlqpzvrr/ILORH5IoZOsksLocEfaOVsTXom5xStnqwqxuSVL8PkS4v/ExG5w+uhk9/otU+pShEDbEUBab4bOKUQqLVF410lra8qd1vptZXCp1xcJomIfJXXQyfpg2lw0PqEEdIG6T3jjfeLp+7TVgg1rahK4ZPVTiLyd5oInax26od0MHXmoC5VvPQaVr0dFrx9/+7w5hmpTN9vlvvgynNqbVF/Z7rr5VQ7WeUkIl+midBJ+mEaIuQESVvXG7beacs0EGjxtJmOgoXaoVAPs5u19ppJbL3/vPGc2hpGoecvFUREztBM6GS1U/uc7S61t52w3Dvn6laDHkKhmhy9dt6ucHv7/h2R3j+schKRr6vs7R0gfRKzRfuh0sp1zgRLR+2Td8l9LbX0GnpzjCkREWmo0gmw2qkXUiBU+iCu52onmdPSGF49vK9Y5SQif8BKJ7lFbkXSlQO/O6FWWrKGlFEomI9JVONLhxpsLViv5X0mIvJVmqp0Aqx26oWjwOFMyFRjdjsDp7Ji8mxfp8VKohYnpNnCKicR+QvNhU7SD9OQ6M4YTmtrgOp9mSVfYxnitRzotLxvRET+TJOhk9VOfbFXVZJbceLi86QEvQVOVjmJyJ9oMnSSvjAk6pOr60PKDXaeXH/Sme50vl+JiLxDs6GT1U7tkw7ySlWXtFSl4oLd1jnzGjlzvnF3OBM2tRQ4WeUkIn+j2dAJMHjqgSeCojcCoD9MRHL2Mdp6rb0R5KydzcoeLYVNIiJ/xSWTSPN8OQDqZWknVwKnvRnvSu+H3rDKSUT+SNOVToDVTn/h7a7PQsF+RVWtaqueA6ej65R+bK4ETm+/r6zxleCsNWfOnMHQoUMRExODwMBA1K9fH6NHj8aFCxe8vWu4cuUKqlSpghUrVphd/sgjj0AQBBw/ftzs8oYNG+Lvf/+7rLYTEhKwYMECpXaVSFWaD53k+7QQCmLy7IckRwHKV8eAeiMgWXabu7rmphbeV9bwi7Tyfv75Z7Rp0wZHjhzBJ598gmPHjmHhwoXYvHkzUlNTcfHiRa/uX3h4ONq0aYMtW7aYXZ6Xl4f4+HizywsLC/Hzzz+jc+fOHt3HmzdvevT+yD/pInTyQ9p3WQaDQi/lBGercr4aMiVyg56tYOfq66hE2ATu7JfWXid2q6vjmWeeQWBgIL766iukp6ejXr166NGjBzZt2oSzZ8/iH//4h3HbhIQEvPTSSxg4cCDCw8MRExODt956y6y9kpISDB8+HLVr10ZkZCS6dOmCvXv3Gq+fNm0aWrVqhQ8//BAJCQmIiorCo48+isuXL9vcx86dOyM3N9f4+8GDB3Ht2jWMGDHC7PItW7agSpUq6NixI44fP47evXujTp06CA8PR9u2bbFp0ybjthkZGTh16hTGjh0LQRAgCHf+8Hbs2IG0tDSEhIQgPj4eo0aNwtWrV82eh5kzZyIrKwtRUVEYNmyYU885kSt0ETrJ99jq9lRjHKCz5AQVPXSLu8qZoGcaDE1vFzvcubCo5BmETN9XnppBT+ooLS01+7lx40aFbS5evIiNGzdixIgRCAkJMbuubt26GDRoEFasWAFRvPP+ev3119GiRQvs3r0bkyZNwtixY/H1118DAERRRK9evXDu3DmsX78e+fn5SElJQdeuXc0qpsePH8eaNWuwbt06rFu3Dnl5eXjttddsPpbOnTvj8OHDKC4uBlAeLjt16oQuXbpUCJ3t27dHaGgorly5gp49e2LTpk3Ys2cPMjMz8eCDD+L06dMAgFWrViEuLg4zZsxAcXGxse39+/cjMzMT/fr1w759+7BixQps374dI0eONNun119/Hc2aNUN+fj6mTJki5yUhcotuQierndqiVtelFsKctA/OVMm0sN/ucif4mQZPOUtpyd3OGda+yGjhS4yEVU7nxcfHIyoqyvjz6quvVtjm6NGjEEURycnJVttITk7G77//jl9//dV4WceOHfHCCy8gKSkJzz77LP76179i/vz5AMpD3/79+/HZZ5+hTZs2aNSoEebMmYOqVaviP//5j7ENg8GApUuXolmzZujUqRP+7//+D5s3b7b5WDp27IgqVaoYA2Zubi7S09ORkpKCkpISHD161Hi51LXesmVLZGdno3nz5mjUqBFmzpyJu+66C2vXrgUAVK9eHZUqVUJERATq1q2LunXrAigPkwMHDsSYMWPQqFEj3HvvvXjzzTfxr3/9C9evXzfuU5cuXTB+/HgkJiYiMTHR4etB5C5dzV4/+34RYofHeHs3CMoEBcPW8kqUFAy0Fty0tj9qUXPcpqO21ahumtLKa1j+OPnF2VlnzpxBZGSk8fegoCCn25AqnKZdz6mpqWbbpKamGifj5Ofn48qVK6hRo4bZNteuXTOb8JOQkICIiAjj79HR0Th//rzN/QgNDUW7du2Qm5uLxx57DHl5eZgwYQIqV66Mjh07Ijc3F0FBQThx4gS6dOkCALh69SqmT5+OdevWoaioCGVlZbh27Zqx0mlLfn4+jh07huXLl5s9DwaDASdOnDAG9DZt2thth0hpugqd5FsC0gDDnxN4pKqiXpYQ8iVCjqDbGdVanSxkiT01romMjDQLndYkJiZCEAQcOHAAffr0qXD9oUOHUK1aNdSsWdNuO1IoNRgMiI6ONuvyllStWtX4/ypVqlS4vcFgsHsfnTt3xooVK/DTTz/h2rVrSElJAQCkp6djy5YtCAwMRHBwMDp06AAAmDBhAjZu3Ig5c+YgMTERISEh+Otf/+pw0o/BYEB2djZGjRpV4bp69eoZ/x8WFma3HSKl6S50stqpb1JIELNFCDmCMWBa/uuPvBW49Ro4gTv7ruXwyW51ddWoUQPdu3fHu+++i7Fjx5qN6zx37hyWL1+OwYMHm1U6v/32W7M2vv32WzRp0gQAkJKSgnPnzqFy5cpISEhQdF87d+6MmTNn4uOPP8Z9992HSpUqASgPnW+99RaCgoKQmpqK4OBgAMC2bduQlZWFvn37AihfeunkyZNmbQYGBuL27dtml6WkpOCnn35ilzlpjm7GdJpi1UD/tBwSlGBaubV3vSlvBW5feC20FpydPWMSueftt9/GjRs3kJmZia1bt+LMmTPYsGEDunfvjtjYWLz88stm23/zzTeYPXs2jhw5gnfeeQefffYZRo8eDQDo1q0bUlNT0adPH2zcuBEnT57Ejh07MHnyZOzatcut/bz33nsRFBSEt956C+npd2a5tW3bFiUlJVi5cqXZUkmJiYlYtWoVCgoKsHfvXgwcOLBCNTUhIQFbt27F2bNn8dtvvwEAnn/+eezcuRPPPPMMCgoKcPToUaxduxbPPvusW/tP5C5dhk7ybZaBTGvL3jhDGjogPQZpKSF/ruiqxZshzzRkmu4DvyB7RqNGjbBr1y40bNgQjzzyCBo2bIjhw4ejc+fO2LlzJ6pXr262/bhx45Cfn4/WrVvjpZdewty5c5GZmQmgvJt8/fr1SEtLw5NPPomkpCQ8+uijOHnyJOrUqePWfkpd55cvX0ZGRobx8ipVqiA1NRWXL182C53z589HtWrVcO+99+LBBx9EZmamsUteMmPGDJw8eRINGzZErVq1AAAtWrRAXl4ejh49ik6dOqF169aYMmUKoqOj3dp/IncJouk6EnYUCdr78GQ3u/7IraqZdjXrYZynFCq1vp/W6KEa5+q4U7WruPb2SYuBM0ZU7zOztLQUUVFRQFYJEGh/HKZDN0uBpVEoKSlxOKbTWQkJCRgzZgzGjBmjaLtE5JjuxnSa4vhO/ZHGcjqit/AWkPZnFVP7+c2MHgIn4Pp+unLOeHfalWgxcBIReZuuQyf5B70E0Dh95DeC/DCql1BORKQHug+drHb6Dz10s8uhlcfBQFWREs8Jq5zaZjn7m4g8xycmEvFD3n/oeVKRRCuB09apSMl1/CwiIrLNJ0In6Yur1aSANG0ENi1yNowzbBIRkaf5TOhkhUE/GHiU52oYZxe7cvgZRERkn8+EToAf+lrELlzyB/zsISJyzKdCJ8APfz1gdc07+Lyrg585RETy6H72OhHJx+BJRETe4nOVToCVBy3i+ofki/hZQ0Qkn0+GToAHA61hyJRPrWWhOLZWWfyMISJyjs+GToAHBa1jECW94mcLEZHzBFEUZR35iwT9fsjyjEXeI+QIDsOlFipwapwlSCtnHrLEsO8ePQfOGFG9z8LS0lJERUXhEA4hAhFutXUZl9EETVBSUoLIyEiF9pCIvM2nK53kfXoJOFoMh0RERL7EL0KnnisT/kDMFj0aTj11Kk2tBlktVJb1ip8lRESu84vQCfBgoQeeCp9aDYOkbfwMISJyj9+EToAHDb3wdOVTazxViSX5+NlBROQ+vwqdAA8eeuKvwZOVWG3hZwYRkTL8LnQCPIjoid6qnnrZV47rlIefFUREyvHb02Cefb+ISynpiJgt6iYoyQ2eenk8/oqBk4hIWX5Z6ZTwoKIcIUcw/qjFNMxJ/9dLZZGIiMjf+XXoJHV4I3jqMXzqcZ/9Bb+QEhEpz2+71yXsZnePp7uI7d2ft7vgXQmRlrdhl7v3MXASEamDlU7wIKM0tSp4UiCz177eq4d633+942cBEZF6GDr/xIONfqg9dtRZSnfv63W4gN7xM4CISF0MnSZ40HGeJ8OR3K5oXwlsDJ+ew799IiL1MXRa4MHHPd4eV+mLGD7Vxb95IiLPYOi0ggch1wk5gqoBydrsdW/ydKXX9MfT9++L+LdOROQ5DJ028GBUkbW1OC2rmp4IQaz8lXM3ePp7RZp/40REnsXQaQcPSndYCyhCjgDDVi/sjB1am2RE2sS/bSIiz/P7dTod8cd1PK0tTWStmiZmiwiQsYyRr/LHx+wLGDiJiLyDlU4Zzr5f5DcHKleqhFoJX1paqJ60yV/+jomItIih0wn+csAyHTNpL1hpJXSxS12ZsyH5Mn/64khEpFUMnU7y9QOXtOSRnLP/yLlebZ4Im5Yzxj19/+QeX/+bJSLSC4ZOF/jiQUwKT9KSR6YBS6uLsHsj8Fk+Zm8/B2SfL/6tEhHpFUOni3ztYOZo+R2tdWF7al9sTaDS4rJNWtsfb/O1v1EiIr1j6HSDLxzUrIVJy3U4HYUZTuDRDgbPcr7wt0lE5GsYOt2k54Ob3HGbjrZh0NEWf3899Pw3SUTkyxg6FaDHmbGunFVIC2GGVU7l+Npzqce/QyIif8LQqSA9H/BsnXHI1dsSeZKe//aIiPwFQ6fC9HDwsxUSrVU8tVDddIa3T8uplQCu1clOatDD3xwRETF0qkIPB0ElZ6lLa3tqQUBaxcsKHeyarXDGBde1Tw9/a0REVI7nXleJdDDU4nnb5QQjObPWnW3TGwoFICbP+dtp9fHIZfr62VpzVc+PkWGTiEh/WOlUmZ4OjnIXhbdFK9VOU3Gi9eqnmt3P3n4eTBf6t+QL3e56+psiIqI7WOn0AC1WPe0FIznnXdcrPYctf8ewSUSkb6x0epBWDprWus6VCGOeGNupZmhUqm0GW+Vp5W+HiIhcx9DpYVo4eJqO9VO6q9UTgcvV+/CHMGjv7FJ6pYW/GSIich+7173Ak93tzk4I0gtnhwD44nPg6xg2iYh8CyudXuSJg6o7YcvV23qyuibn9JyuzML3hQqhnjFwEhH5HoZOL/PFU/d5uqqo1BqbvkpPAdoX/x6IiKgcQ6dGuHOg1VKo8Na+aPl88SQPwyYRkW9j6NQQVnncw4Cpz6ov3/dERP6BE4k0SO5EI2mSkLRUkTfDhWmF05v74s6yTVqqGLtD6yFTwqBJRORfWOnUMLkHZW+FJek87VoLa0qFLq09Ll/CwElE5H9Y6dQ4R1VPbwZOX+Jrj0erGDaJiPwXQ6dO2Auf3uhOdbS0kLe7+93libMr+ROGTSIiYve6zpgevL09htNRKNNTaNPjBBw94CQhIiKSMHTqkOmBXE4wUiP8yblfPYU26TlS+rSg/oxhk4iITLF7Xcc8eTpNaxx1QZsGOS3SUyVWTxg2iYjIGlY6fYCcLkxvBiwthjtr+6TVcKwX7EonIiJ7GDp9iK2DvpphSm7bWgmetsaiMnC6jmGTiIjkYPe6D7LV7a5Gd7dWwiR5HoMmERE5g6HTh5mGgljcGX9pGhTdCaDOBk5vLqOkdDhW6jnUGwZNIiJyFUOnnzj7fhHw/lmz6qc7Xc16qnAqsa+22vCXwMmwSURE7mLo9DOmXe+mgclTIVKLi8Y72h9/DpwMm0REpBSGTj9l1vVuEkCl6qe1QKVUMNVi8LRFTxVdpTBoEhGRGhg6yWr10zIYqjUmUsvh095j1vJ+u4phk4iI1MTQSUaWE49MqRVAvR0+bd2vv3SpM2gSEZGnMHSSVZbd777IXwMngyYREXkDQyc5ZG38py+OddTTElDOYtAkIiJvY+gkpxjHf0L54Gka4pQKdHKrlpbbGbYCAWm2t7d1mZYwaBIRkZYwdJJLpHU/gfLqp7sB1NrC9e7y9W5yaxg0iYhIqxg6yW2WARRQbkF2JQOiM5OhpCqn1gMqQyYREekFQycpyhiCVAihzlDiNJVaDJwMmUREpFcMnaQqyxAKeC6IOhq36QxvTRpiyCQiIl/B0EkeJzeIuhPy3F1X1NOnCGW4JCIiX8fQSZpgLYgC7q8VqkRgVLLCyXBJRET+iqGTdEFOWLMMpp4eR8pASUREZBtDJ/mMCqHPSve9/HYYIImIiJTE0El+gVVIIiIi7wrw9g4QERHp2bRp09CqVSu722RlZaFPnz4e2R+lZWRkYMyYMd7eDc2Q83o7Q8/vDWcxdBIRke6dOXMGQ4cORUxMDAIDA1G/fn2MHj0aFy5ccKqdkydPQhAEFBQUyL7N+PHjsXnzZif32HkZGRkQBAGffvqp2eULFixAQkKC2+3n5uZCEARcunTJ7bbUZCsEL126FFWrVlX9/j31evsihk4iItK1n3/+GW3atMGRI0fwySef4NixY1i4cCE2b96M1NRUXLx4UdX7Dw8PR40aNVS9D0lwcDAmT56MW7duKdqu0u1p7f6U5MnX21W3b9+GwWDw9m5UwNBJRES69swzzyAwMBBfffUV0tPTUa9ePfTo0QObNm3C2bNn8Y9//MO4rSAIWLNmjdntq1atiqVLlwIAGjRoAABo3bo1BEFARkYGgPIqYLt27RAWFoaqVauiY8eOOHXqFICK3a23b9/Gc889h6pVq6JGjRqYOHEiRNF86TVRFDF79mzcddddCAkJQcuWLfGf//zH4WN97LHHUFJSgn/+8592t3vvvffQsGFDBAYGonHjxvjwww/NrhcEAQsXLkTv3r0RFhaGp556Cp07dwYAVKtWDYIgICsry7i9wWDAxIkTUb16ddStWxfTpk0za6+kpATDhw9H7dq1ERkZiS5dumDv3r3G66Xn6IMPPsBdd92FoKAgiKIIQRCwaNEi9O3bF6GhoWjUqBHWrl3r8HmQ6/PPP8c999yD4OBg3HXXXZg+fTrKysrMnoecnBw88MADCA0NRXJyMnbu3Iljx44hIyMDYWFhSE1NxfHjxys8Fomc90ZOTg7i4+MRGhqKAQMGWK0mz5kzB9HR0ahRowaeeeYZs2B+8+ZNTJw4EbGxsQgLC0P79u2Rm5trvF6q8q5btw5NmzZFUFAQTp065fB2nsbQSUREunXx4kVs3LgRI0aMQEhIiNl1devWxaBBg7BixYoKoc+W77//HgCwadMmFBcXY9WqVSgrK0OfPn2Qnp6Offv2YefOnRg+fDgEwfqybHPnzsUHH3yAxYsXY/v27bh48SJWr15tts3kyZOxZMkSvPfee/jpp58wduxYPP7448jLy7O7f5GRkfj73/+OGTNm4OrVq1a3Wb16NUaPHo1x48bhxx9/RHZ2Np544gls2bLFbLupU6eid+/e2L9/P2bMmIGVK1cCAA4fPozi4mK88cYbxm2XLVuGsLAwfPfdd5g9ezZmzJiBr7/+GkB5gO7VqxfOnTuH9evXIz8/HykpKejatatZlfnYsWP497//jZUrV5oNX5g+fToefvhh7Nu3Dz179sSgQYMUqU5v3LgRjz/+OEaNGoUDBw4gJycHS5cuxcsvv2y23UsvvYTBgwejoKAATZo0wcCBA5GdnY1JkyZh165dAICRI0davQ857w3pcX/++efYsGEDCgoK8Mwzz5i1s2XLFhw/fhxbtmzBsmXLsHTpUuMXIQB44okn8M033+DTTz/Fvn37MGDAANx///04evSocZs//vgDr776KhYtWoSffvoJtWvXlnU7jxKJiIjcVFJSIgIQf8AP4iEccuvnB/wgAhDPnDkjlpSUGH+uX79e4X6//fZbEYC4evVqq/s1b948EYD4yy+/iKIoWt02KipKXLJkiSiKonjixAkRgLhnzx7j9RcuXBABiLm5uVbvY+rUqWLLli2Nv0dHR4uvvfaa8fdbt26JcXFxYu/evUVRFMUrV66IwcHB4o4dO8zaGTp0qPjYY49ZvQ9RFMX09HRx9OjR4vXr18X69euLM2bMEEVRFOfPny/Wr1/fuN29994rDhs2zOy2AwYMEHv27Gn8HYA4ZswYs222bNkiAhB///33Cvd73333mV3Wtm1b8fnnnxdFURQ3b94sRkZGVnh9GjZsKObk5IiiWP4cValSRTx//rzZNgDEyZMnG3+/cuWKKAiC+OWXX9p9HqpUqSKGhYWZ/QQFBYlRUVHG7Tp16iS+8sorZrf98MMPxejoaJv3v3PnThGAuHjxYuNln3zyiRgcHGz83fT1lvPeqFSpknjmzBnjZV9++aUYEBAgFhcXi6IoikOGDBHr168vlpWVGbcZMGCA+Mgjj4iiKIrHjh0TBUEQz549a9Z2165dxUmTJomiKIpLliwRAYgFBQXG6+XcztO4ZBIREbktMDAQdevWRdtzbRVpLzw8HPHx8WaXTZ06tUK3riPinxVOW1VJOapXr46srCxkZmaie/fu6NatGx5++GFER0dX2LakpATFxcVITU01Xla5cmW0adPGuC8HDhzA9evX0b17d7Pb3rx5E61bt3a4P0FBQZgxYwZGjhyJp59+usL1Bw8exPDhw80u69ixo1nlEgDatGnj8L4kLVq0MPs9Ojoa58+fBwDk5+fjypUrFcY5Xrt2zaxbun79+qhVq5bdtsPCwhAREWFs25ZBgwaZDZsAgFWrVuGVV14x/p6fn48ffvjBrLJ5+/ZtXL9+HX/88QdCQ0Mr3H+dOnUAAM2bNze77Pr16ygtLUVkZKTZfcp5b9SrVw9xcXHG31NTU2EwGHD48GHUrVsXAHD33XejUqVKxm2io6Oxf/9+AMDu3bshiiKSkpLM7vvGjRtmz3lgYKDZY5F7O09i6CQiIrcFBwfjxIkTuHnzpiLtiX+O9zMVFBRUYbvExEQIgoADBw5YXXbm0KFDqFatGmrWrAmgPHyKFl3tcia1LFmyBKNGjcKGDRuwYsUKTJ48GV9//TU6dOjgxKMqJ03w+OKLLxAbG2t2nbXHaM3jjz+OOXPmYObMmVZnrls+d9aez7CwMNn7XKVKlQrtS4/DYDAgOjra6lhB09nktu7PXtu2REVFITEx0eyy2rVrm/1uMBgwffp09OvXr8Ltg4ODrd6/9BxZu8zWPjn73pDaM309HD2/lSpVQn5+vlkwBcq/nElCQkLM2pR7O09i6CQiIkUEBwebHcw9oUaNGujevTveffddjB071mxc57lz57B8+XIMHjzYeDCuVasWiouLjdscPXoUf/zxh/H3wMBAAOUVMUutW7dG69atMWnSJKSmpuLjjz+uECyioqIQHR2Nb7/9FmlpaQDKx/1J4xwBGCd6nD59Gunp6S497oCAALz66qvo169fhWpncnIytm/fjsGDBxsv27FjB5KTk+22ae+x25OSkoJz586hcuXKiizdpJSUlBQcPny4QjhVg733xunTp1FUVISYmPIz4+3cuRMBAQEVKpD22r59+zbOnz+PTp06ObVPrtxOTQydRESka2+//TbuvfdeZGZmYubMmWjQoAF++uknTJgwAbGxsWbdq126dMHbb7+NDh06wGAw4PnnnzerMtWuXRshISHYsGED4uLiEBwcjIsXL+L999/HQw89hJiYGBw+fBhHjhwxC3WmRo8ejddeew2NGjVCcnIy5s2bZzZbOSIiAuPHj8fYsWNhMBhw3333obS0FDt27EB4eDiGDBki63H36tUL7du3R05OjrFbGAAmTJiAhx9+2DiZ5/PPP8eqVauwadMmu+3Vr18fgiBg3bp16NmzJ0JCQmRVxLp164bU1FT06dMHs2bNQuPGjVFUVIT169ejT58+TnXjK+nFF1/EAw88gPj4eAwYMAABAQHYt28f9u/fj5kzZypyHydOnHD43ggODsaQIUMwZ84clJaWYtSoUXj44YeNXeuOJCUlYdCgQRg8eDDmzp2L1q1b47fffsP//vc/NG/eHD179lT0dmri7HUiItK1Ro0aYdeuXWjYsCEeeeQRNGzYEMOHD0fnzp2xc+dOVK9e3bjt3LlzER8fj7S0NAwcOBDjx483ju0Dysdfvvnmm8jJyUFMTAx69+6N0NBQHDp0CP3790dSUhKGDx+OkSNHIjs72+r+jBs3DoMHD0ZWVhZSU1MRERGBvn37mm3z0ksv4cUXX8Srr76K5ORkZGZm4vPPPzcu2STXrFmzcP36dbPL+vTpgzfeeAOvv/467r77buTk5GDJkiXG5Z9siY2NxfTp0/HCCy+gTp06NmdsWxIEAevXr0daWhqefPJJJCUl4dFHH8XJkyfNwrCnZWZmYt26dfj666/Rtm1bdOjQAfPmzUP9+vUVuw85743ExET069cPPXv2xF/+8hc0a9YM7777rlP3s2TJEgwePBjjxo1D48aN8dBDD+G7776rMO5ZqdupRRAtB7cQERERkdumTZuGNWvWOHWGK1/GSicRERERqY6hk4iIiIhUx+51IiIiIlIdK51EREREpDqGTiIiIiJSHUMnEREREamOoZOIiIiIVMfQSURERESqY+gkIiIiItUxdBIRERGR6hg6iYiIiEh1DJ1EREREpLr/B13/QNz18jzWAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Total HDF keys: ['_nc4_non_coord_y', 'polar_stereographic', 'sea_ice_extent', 'time', 'x', 'y']\n" + ] + } + ], + "source": [ + "# Define a custom colormap with five colors\n", + "# (Colors used from IMS product map)\n", + "\"\"\"\n", + "0 (outside Northern Hemisphere).\n", + "1 (open water)\n", + "2 (land without snow)\n", + "3 (sea or lake ice)\n", + "4 (snow covered land)\n", + "\"\"\"\n", + "colors = [\"#E500E5\", \"#0066FF\", \"#01FF00\", \"#FFC100\", \"#E50000\"]\n", + "cmap = ListedColormap(colors, name=\"custom_colormap\", N=len(colors))\n", + "\n", + "with h5py.File(f\"{root_dir}/2015/masie_all_r00_v01_2015001_1km.nc\", 'r') as file:\n", + " sie = file['sea_ice_extent']\n", + " print(f\"Unique values in HDF file {np.unique(sie)}\")\n", + "\n", + "\n", + " # Plot\n", + " plt.imshow(sie[0], cmap=cmap, vmin=0, vmax=len(colors) - 1)\n", + " plt.title('January 1, 2015', fontsize=8)\n", + " plt.axis(\"off\")\n", + "\n", + " cbar = plt.colorbar(ticks=np.arange(len(colors)))\n", + " cbar.ax.set_yticklabels(\n", + " [\n", + " \"Outside Northern Hemisphere\",\n", + " \"Open Water\",\n", + " \"Land without Snow\",\n", + " \"Sea or Lake Ice\",\n", + " \"Snow Covered Land\",\n", + " ]\n", + " )\n", + " cbar.set_label(\"Surface Type\")\n", + "\n", + " plt.show()\n", + "\n", + " keys = list(file.keys())\n", + " print(f\"Total HDF keys: {keys}\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}