diff --git a/quickstarts/stacsearch_visualizations.ipynb b/quickstarts/stacsearch_visualizations.ipynb new file mode 100644 index 0000000..70f7780 --- /dev/null +++ b/quickstarts/stacsearch_visualizations.ipynb @@ -0,0 +1,907 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "311ba37a-5e26-45fb-819c-6fe354aeb583", + "metadata": {}, + "source": [ + "# Common Interactions with STAC Search Results \n", + "Here are some common functionalities that allow to get information directly from a STAC search (list of pystac.item).\n", + "- plot a time series for available images\n", + "- plot the bboxes\n", + "- plot a heatmap of bboxes (for many bboxes that overlap)" + ] + }, + { + "cell_type": "markdown", + "id": "918c1980-5ad7-4dc9-b417-97b638e8bec1", + "metadata": {}, + "source": [ + "## STAC Search" + ] + }, + { + "cell_type": "code", + "execution_count": 37, + "id": "4740e183-7654-4004-bb17-b5cb200270d5", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found 2483 Scenes\n" + ] + } + ], + "source": [ + "from pystac_client import Client as pystacclient\n", + "\n", + "catalog_url = 'https://stac.terrabyte.lrz.de/public/api'\n", + "catalog = pystacclient.open(catalog_url)\n", + "\n", + "collection = 'sentinel-1-slc'\n", + "start = '2018-01-01T00:00:00Z'\n", + "end = '2025-12-31T23:59:59Z'\n", + "\n", + "bbox_larsen = [-66.16915256, -69.69490353, -59.94131864, -65.3833819]\n", + "\n", + "res = catalog.search(collections=[collection],\n", + " datetime=[start, end],\n", + " bbox = bbox_larsen,\n", + " limit=500)\n", + "\n", + "items = list(res.items())\n", + "print(f'Found {len(items)} Scenes') " + ] + }, + { + "cell_type": "markdown", + "id": "503eda24-1577-4419-bc77-b11b8e98df8f", + "metadata": {}, + "source": [ + "## Timeseries\n", + "After a STAC search it is important to get an overview of the dates that are within the date limits set in the serach. \n", + "A histogram helps to understand if the temporal distribution of the data is as expected. " + ] + }, + { + "cell_type": "markdown", + "id": "9fc11ab8-9b6a-428d-8f16-85d5d56a45c3", + "metadata": {}, + "source": [ + "This function allows to create a histogram that is grouped by any combination of day, week, month, year." + ] + }, + { + "cell_type": "code", + "execution_count": 42, + "id": "65bba08a-5d6d-4b90-b950-8e86369fdc3f", + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "\n", + "def plot_date_histogram(items, group_by=(\"year\", \"month\"), figsize=(12, 5)):\n", + " \"\"\"\n", + " Plot a histogram of datetime entries aggregated by year/month/week/day.\n", + "\n", + " Parameters:\n", + " dates: list of datetime\n", + " group_by: tuple or list\n", + " Any combination of: \"year\", \"month\", \"week\", \"day\"\n", + " Example: (\"year\",) or (\"year\", \"month\") or (\"year\", \"week\")\n", + " figsize: tuple\n", + " Output figure size\n", + "\n", + " Returns:\n", + " (fig, ax)\n", + " \"\"\"\n", + "\n", + " dates = []\n", + " for item in items:\n", + " dates.append(item.datetime)\n", + "\n", + " print(min(dates))\n", + " print(max(dates))\n", + "\n", + " \n", + " valid_parts = {\"year\", \"month\", \"week\", \"day\"}\n", + " group_by = tuple(group_by)\n", + " if not all(p in valid_parts for p in group_by):\n", + " raise ValueError(f\"group_by must use only {valid_parts}\")\n", + "\n", + " df = pd.DataFrame({\"date\": dates})\n", + "\n", + " # Extract numeric fields for correct sorting\n", + " df[\"year\"] = df[\"date\"].dt.year\n", + " df[\"month\"] = df[\"date\"].dt.month\n", + " df[\"week\"] = df[\"date\"].dt.isocalendar().week.astype(int)\n", + " df[\"day\"] = df[\"date\"].dt.day\n", + "\n", + " # Build a numeric tuple for sorting\n", + " df[\"sort_key\"] = df.apply(lambda r: tuple(r[p] for p in group_by), axis=1)\n", + "\n", + " # Build the string label for display\n", + " df[\"label\"] = df.apply(\n", + " lambda r: \"-\".join(f\"{r[p]:02d}\" if p != \"year\" else str(r[p]) for p in group_by),\n", + " axis=1\n", + " )\n", + "\n", + " # Count by sort_key\n", + " grouped = df.groupby([\"sort_key\", \"label\"]).size().reset_index(name=\"count\")\n", + "\n", + " # Sort by the numeric tuple\n", + " grouped = grouped.sort_values(\"sort_key\")\n", + "\n", + " # Extract sorted labels and values\n", + " labels = grouped[\"label\"].tolist()\n", + " values = grouped[\"count\"].tolist()\n", + "\n", + " # Plot\n", + " fig, ax = plt.subplots(figsize=figsize)\n", + " ax.bar(labels, values)\n", + " plt.xticks(rotation=90)\n", + " ax.set_xlabel(\" / \".join(group_by).title())\n", + " ax.set_ylabel(\"Count\")\n", + " ax.set_title(f\"Histogram Grouped by: {', '.join(group_by).title()}\")\n", + " plt.tight_layout()\n", + "\n", + " return fig, ax\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "076e475e-5e48-4a61-bd8f-11e1c3499a86", + "metadata": {}, + "source": [ + "The histogram shows the distribution of the available time steps in your search." + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "id": "1612a06a-60eb-42c2-a790-6464fd6b2a2c", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2019-01-01 08:18:00.864852+00:00\n", + "2025-12-09 08:03:36.304436+00:00\n" + ] + }, + { + "data": { + "text/plain": [ + "(
,\n", + " )" + ] + }, + "execution_count": 43, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAABKUAAAHqCAYAAADVi/1VAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWv1JREFUeJzt3XlcFvX6//H3IAiCguYCkgpquG+VZmIJuWZmHa1Trmm222bWMc1TormknkxPLm1qtpidU+qxzTRzydRSyxY129xS0DJDMyOFz++Pvtw/7wAFhJkP3K/n43E/Hs3MNTPXNfeHe26vZuZ2jDFGAAAAAAAAgIuCvE4AAAAAAAAAgYemFAAAAAAAAFxHUwoAAAAAAACuoykFAAAAAAAA19GUAgAAAAAAgOtoSgEAAAAAAMB1NKUAAAAAAADgOppSAAAAAAAAcB1NKQAAAAAAALiOphQAAEXs+eefl+M42rRpU67Lr7zySsXHx/vNi4+P18CBAwu0n3Xr1iklJUW//PJL4RItZdauXavevXurVq1aCg0NVUREhBo3bqz7779fX331ldfpuSI5OVnJyclnjHMcR3fddVfxJ/R/fvnlF9WoUUOtW7dWZmZmjuVr165VmTJlNGLECNdyKojk5GQ5jqM6derIGJNj+Zo1a+Q4jhzH0fPPP1+sucycOTPXfaxatUqO4+i1114r1v0DAFCUaEoBAGCBRYsW6eGHHy7QOuvWrdPo0aNpSkn65z//qUsvvVS7d+/WP//5Ty1dulSLFy/WoEGDtHz5cjVs2DDXZgjcUbFiRc2ZM0cff/yxJk6c6Lfst99+04033qjGjRtr9OjRHmV4ZhUqVNDOnTv1/vvv51g2Z84cRUZGupJHXk0pAABKomCvEwAAANL555/vdQoFduLECTmOo+Bgb79OvPLKKxo3bpxuv/12zZw5U47j+JZ16tRJQ4cO1cyZM8+4nd9++03h4eHFmWpA69y5s+644w6NHj1a3bt3V9OmTSVJw4cP1+7du7Vx40aVLVu22PM4fvy4wsLC/MZJftSqVUsVKlTQnDlz1KFDB9/8o0eP6r///a/69u2rZ599tqjTBQCgVONKKQAALPDX2/eysrI0duxY1a9fX+XKlVPFihXVrFkzTZs2TZKUkpKif/zjH5Kk2rVr+24dWrVqlW/9SZMmqUGDBgoNDVW1atV0ww036IcffvDbrzFG48ePV1xcnMLCwtSyZUstX748x21g2bcGvfjii7r//vt17rnnKjQ0VN9++61+/PFHDR48WI0aNVL58uVVrVo1tW/fXh988IHfvnbt2iXHcTR58mRNnDhR8fHxKleunJKTk/X111/rxIkTGj58uGJjYxUVFaUePXro4MGDZzx2Y8eOVZUqVfTEE0/k2mhwHEd33nmnypQp45uXnJysJk2aaM2aNUpMTFR4eLgGDRokSdqzZ4/69eunatWqKTQ0VA0bNtTjjz+urKysHMcj+3j/tcZTr2QZOHCgypcvr61bt6pDhw6KiIhQ1apVddddd+m3337L8X7MnDlTLVq0ULly5VSpUiVde+21+v7773PETZo0yfe+XXDBBXrnnXfOeKz+6umnn1a9evUUGhqqRo0aacGCBX61BAcHa8KECTnWy75d7b///W+B9jd58mTVrFlTAwYM0IkTJ7RmzRpNnz5dKSkpat68uSTp1VdfVZs2bRQREaHy5curS5cu+vTTT/22s2nTJvXq1cs3huLj49W7d2/t3r3bLy77Vtply5Zp0KBBqlq1qsLDw5WRkVGgvLMNGjRICxcu9Ls6MfuY9erVK9d11q5dqw4dOqhChQoKDw9XYmKi3nrrrVzzXLlype644w5VqVJFlStXVs+ePbV//35fXHx8vLZu3arVq1f7/ub/eivwiRMnNHLkSMXGxioyMlIdO3bUjh07ClUvAADFjaYUAADFJDMzUydPnszxyu2ZNH81adIkpaSkqHfv3nrrrbf06quv6qabbvL9Y/jmm2/W3XffLUlauHCh1q9fr/Xr1+uCCy6QJN1xxx168MEH1alTJy1ZskSPPvqoli5dqsTERP3000++/YwcOVIjR47U5Zdfrv/973+6/fbbdfPNN+vrr7/ONa8RI0Zoz549euqpp/TGG2+oWrVq+vnnnyVJo0aN0ltvvaW5c+eqTp06Sk5OztG0kaQZM2boww8/1IwZM/Tcc8/pq6++Uvfu3XXTTTfpxx9/1Jw5czRp0iS99957uvnmm097nPbv369t27apU6dOCgsLO+NxPVVqaqr69eunPn366O2339bgwYP1448/KjExUcuWLdOjjz6qJUuWqGPHjnrggQfO6hlMJ06c0BVXXKEOHTpo8eLFuuuuu/T000/r+uuv94u77bbbNGTIEHXs2FGLFy/WzJkztXXrViUmJurAgQO+uNGjR/ve38WLF+uOO+7QLbfcUqDmw5IlS/Tvf/9bY8aM0Wuvvaa4uDj17t3b90yi+Ph4XXXVVXrqqady3Po4ffp0xcbGqkePHr4GXUpKyhn3GRERoXnz5umzzz7TQw89pBtvvFEXXXSRHnzwQUnS+PHj1bt3bzVq1Ej/+c9/9OKLL+ro0aO69NJLtW3bNt92du3apfr162vq1Kl69913NXHiRKWmpqpVq1Z+4zvboEGDFBISohdffFGvvfaaQkJC8n2cTtWrVy+VKVNGr7zyim/e7Nmzde211+Z6+97q1avVvn17paena/bs2XrllVdUoUIFde/eXa+++mqO+JtvvlkhISGaP3++Jk2apFWrVqlfv36+5YsWLVKdOnV0/vnn+/7mFy1a5LeNhx56SLt379Zzzz2nZ555Rt988426d+/O7asAADsZAABQpObOnWsknfYVFxfnt05cXJwZMGCAb/rKK680LVq0OO1+Jk+ebCSZnTt3+s3fvn27kWQGDx7sN/+jjz4yksxDDz1kjDHm559/NqGhoeb666/3i1u/fr2RZJKSknzzVq5caSSZdu3anbH+kydPmhMnTpgOHTqYHj16+Obv3LnTSDLNmzc3mZmZvvlTp041ksxVV13lt50hQ4YYSSY9PT3PfW3YsMFIMsOHD88zj+xXVlaWb1lSUpKRZFasWOG3zvDhw40k89FHH/nNv+OOO4zjOGbHjh1+x2PlypV+cdk1zp071zdvwIABRpKZNm2aX+y4ceOMJLN27VpjzP8/7o8//rhf3N69e025cuXMsGHDjDHGHD582ISFhfkdW2OM+fDDD3O8b3mRZMqVK2fS0tJ8806ePGkaNGhgzjvvPN+87DoXLVrkm7dv3z4THBxsRo8ebYwxZtWqVaZMmTK+6fwYNmyYL4fsY7pnzx4THBxs7r77br/Yo0ePmpiYGHPdddflub2TJ0+aX3/91URERPgd5+y/xRtuuCHfueUmKSnJNG7c2Bjz5/vZsmVLY4wxW7duNZLMqlWrzMaNG3O89xdffLGpVq2aOXr0qF+uTZo0MTVq1PCNyew8//o3O2nSJCPJpKam+uY1btw41/c4+7264oor/Ob/5z//MZLM+vXrz+oYAABQHLhSCgCAYvLCCy9o48aNOV6XXHLJGde96KKL9Nlnn2nw4MF69913deTIkXzvd+XKlZKU49f8LrroIjVs2FArVqyQJG3YsEEZGRm67rrr/OIuvvjiHLcEZbvmmmtynf/UU0/pggsuUFhYmIKDgxUSEqIVK1Zo+/btOWKvuOIKBQX9/68gDRs2lCR169bNLy57/p49e/Ko9PQqV66skJAQ3+v111/3W16pUiW1b9/eb97777+vRo0a6aKLLvKbP3DgQBljcn3IdX717dvXb7pPnz6S/v/79eabb8pxHPXr18/vyrqYmBg1b97cd9XZ+vXr9fvvv+fYXmJiouLi4vKdT4cOHRQdHe2bLlOmjK6//np9++23vts8k5OT1bx5c82YMcMX99RTT8lxHN16662SpKSkJJ08eVKPPPJIvvc9ZswYSVK/fv1Ur149SdK7776rkydP6oYbbvCrPywsTElJSX5X3f3666968MEHdd555yk4OFjBwcEqX768jh07luuYy2vcFsagQYO0adMmffHFF5o9e7bq1q2rdu3a5Yg7duyYPvroI1177bUqX768b36ZMmXUv39//fDDDzmubLvqqqv8pps1ayZJOW5LPJ2i2AYAAG7hQecAABSThg0bqmXLljnmR0VFae/evaddd8SIEYqIiNBLL72kp556SmXKlFG7du00ceLEXLd5qkOHDkmSqlevnmNZbGys7x+n2XGnNiay5TYvr21OmTJF999/v26//XY9+uijqlKlisqUKaOHH3441wbBOeec4zed/XDrvOb//vvvueYiSTVr1pSU+z+4V61apZMnT2rz5s26/fbb81XLoUOHcm3IxcbG+pYXRnBwsCpXruw3LyYmxm+bBw4ckDEmz2Nfp04dv/js9XPbZn6cbv1Dhw6pRo0akqR77rlHN998s3bs2KE6dero2Wef1bXXXlugff1VaGioJPk92Dz79sRWrVrlus6pjcw+ffpoxYoVevjhh9WqVStFRkbKcRxdccUVOn78eI51c3uvC6tdu3ZKSEjQ008/rf/85z8aMmRIrs8yO3z4sIwxef4dSjnH01/HSPZxyq2mvBTFNgAAcAtNKQAALBQcHKyhQ4dq6NCh+uWXX/Tee+/poYceUpcuXbR3797T/kpc9j9KU1NTfY2FbPv371eVKlX84k59VlG2tLS0XJszuf3j+6WXXlJycrJmzZrlN//o0aOnL7IIxMbGqnHjxlq+fLl+//13v+dKtWjRQtKfV9XkJrdaKleurNTU1Bzzsx82nX3ssvfz1wdm5/Y8I0k6efKkDh065NcwSEtL8+0ze9uO4+iDDz7wNRJOlT0vOz57/VPl9b7lJq/1T92H9GcD6MEHH9SMGTN08cUXKy0tTXfeeWe+9lEQ2cc2+/lWeUlPT9ebb76pUaNGafjw4b75GRkZvueb/VVBf2nvTG688Ub985//lOM4GjBgQK4xlSpVUlBQUL7GEwAAgYrb9wAAsFzFihV17bXX6s4779TPP/+sXbt2Scr7CojsW9Jeeuklv/kbN27U9u3bfT9n37p1a4WGhuZ44PKGDRsKdKuP4zg5miiff/651q9fn+9tnI2RI0fqp59+0tChQ/P1EPnT6dChg7Zt26ZPPvnEb/4LL7wgx3F02WWXSZKv8fP555/7xS1ZsiTPbb/88st+0/Pnz5ck368cXnnllTLGaN++fWrZsmWOV9OmTSX9eXtlWFhYju2tW7euQO/bihUr/BqSmZmZevXVV1W3bl2/ZmZYWJhuvfVWzZs3T1OmTFGLFi3Utm3bfO8nv7p06aLg4GB99913udaffYWg4zgyxuQYc88995xrD/MeMGCAunfvrn/84x8699xzc42JiIhQ69attXDhQr+/0aysLL300kuqUaOG79bFgggNDeWqJwBAqcGVUgAAWKh79+5q0qSJWrZsqapVq2r37t2aOnWq4uLilJCQIEm+JsW0adM0YMAAhYSEqH79+qpfv75uvfVWPfnkkwoKClLXrl21a9cuPfzww6pZs6buu+8+SX/eLjd06FBNmDBBlSpVUo8ePfTDDz9o9OjRql69ut/tUqdz5ZVX6tFHH9WoUaOUlJSkHTt2aMyYMapdu7ZOnjxZPAfoFL1799bWrVs1btw4ffbZZxo4cKASEhKUlZWlvXv36sUXX5QkVahQ4Yzbuu+++/TCCy+oW7duGjNmjOLi4vTWW29p5syZuuOOO3xNhJiYGHXs2NF37OLi4rRixQotXLgw1+2WLVtWjz/+uH799Ve1atVK69at09ixY9W1a1ffM8batm2rW2+9VTfeeKM2bdqkdu3aKSIiQqmpqVq7dq2aNm2qO+64Q5UqVdIDDzygsWPH6uabb9bf//537d27VykpKQW6pa5KlSpq3769Hn74YUVERGjmzJn66quvtGDBghyxgwcP1qRJk7R582Y999xzfstWr16tDh066JFHHinQc6X+Kj4+XmPGjNHIkSP1/fff6/LLL1elSpV04MABffzxx4qIiNDo0aMVGRmpdu3aafLkyapSpYri4+O1evVqzZ49WxUrVizQPh3HyfG8qvyIjY3V4sWLzxg3YcIEderUSZdddpkeeOABlS1bVjNnztSXX36pV155pVBXcDVt2lQLFizQq6++qjp16igsLMz3WQAAQElDUwoAAAtddtllev311/Xcc8/pyJEjiomJUadOnfTwww/7fs4+OTlZI0aM0Lx58/Tss88qKytLK1eu9N1KV7duXc2ePVszZsxQVFSULr/8ck2YMMHv1qxx48YpIiJCTz31lObOnasGDRpo1qxZGjlyZL7/gT9y5Ej99ttvmj17tiZNmqRGjRrpqaee0qJFiwr8j/3CGjt2rLp06aIZM2ZozJgxOnDggEJCQhQfH6+kpCRNnDhRF1544Rm3U7VqVa1bt04jRozQiBEjdOTIEdWpU0eTJk3S0KFD/WJffPFF3X333XrwwQeVmZmp7t2765VXXsn1mV8hISF68803dc8992js2LEqV66cbrnlFk2ePNkv7umnn9bFF1+sp59+WjNnzlRWVpZiY2PVtm1bv4evjxkzxtdIevHFF9WgQQM99dRT+te//pXvY3bVVVepcePG+uc//6k9e/aobt26evnll3X99dfniD333HN1ySWX6PPPP/c9oD2bMUaZmZnKysrK977zMmLECDVq1EjTpk3TK6+8ooyMDMXExKhVq1Z+zwWbP3++7r33Xg0bNkwnT55U27ZttXz58hwPyz+d7Ns6i/J5U3+VlJSk999/X6NGjdLAgQOVlZWl5s2ba8mSJbryyisLtc3Ro0crNTVVt9xyi44ePaq4uDjf1ZMAAJQ0jjnb69wBAECpsnPnTjVo0ECjRo3SQw895HU6Jd7AgQP12muv5flsq5Lg4MGDiouL0913361JkyZ5nU6RePvtt3XllVfqs88+40ojAAA8wpVSAAAEsM8++0yvvPKKEhMTFRkZqR07dmjSpEmKjIzUTTfd5HV68NgPP/yg77//XpMnT1ZQUJDuvfder1MqMitXrlSvXr1oSAEA4CGaUgAABLCIiAht2rRJs2fP1i+//KKoqCglJydr3Lhxio6O9jo9eOy5557TmDFjFB8fr5dffjnPh3qXRH+9dRIAALiP2/cAAAAAAADguvz9rA4AAAAAAABQhGhKAQAAAAAAwHU0pQAAAAAAAOC6Uv+g86ysLO3fv18VKlSQ4zhepwMAAAAAAFCqGWN09OhRxcbGKigo7+uhSn1Tav/+/apZs6bXaQAAAAAAAASUvXv3qkaNGnkuL/VNqQoVKkj680BERkZ6nA0AAAAAAEDpduTIEdWsWdPXk8lLqW9KZd+yFxkZSVMKAAAAAADAJWd6jBIPOgcAAAAAAIDraEoBAAAAAADAdTSlAAAAAAAA4DqaUgAAAAAAAHAdTSkAAAAAAAC4jqYUAAAAAAAAXEdTCgAAAAAAAK6jKQUAAAAAAADX0ZQCAAAAAACA62hKAQAAAAAAwHU0pQAAAAAAAOA6mlIAAAAAAABwHU0pAAAAAAAAuI6mFAAAAAAAAFxHUwoAAAAAAACuC/Y6AQAAAAAAAC/FD38rX3G7HutWzJkEFs+vlNq3b5/69eunypUrKzw8XC1atNDmzZt9y40xSklJUWxsrMqVK6fk5GRt3brVw4wBAAAAAABwtjxtSh0+fFht27ZVSEiI3nnnHW3btk2PP/64Klas6IuZNGmSpkyZounTp2vjxo2KiYlRp06ddPToUe8SBwAAAAAAwFnx9Pa9iRMnqmbNmpo7d65vXnx8vO+/jTGaOnWqRo4cqZ49e0qS5s2bp+joaM2fP1+33Xab2ykDAAAAAACgCHh6pdSSJUvUsmVL/f3vf1e1atV0/vnn69lnn/Ut37lzp9LS0tS5c2ffvNDQUCUlJWndunW5bjMjI0NHjhzxewEAAAAAAMAunjalvv/+e82aNUsJCQl69913dfvtt+uee+7RCy+8IElKS0uTJEVHR/utFx0d7Vv2VxMmTFBUVJTvVbNmzeItAgAAAAAAAAXmaVMqKytLF1xwgcaPH6/zzz9ft912m2655RbNmjXLL85xHL9pY0yOedlGjBih9PR032vv3r3Flj8AAAAAAAAKx9OmVPXq1dWoUSO/eQ0bNtSePXskSTExMZKU46qogwcP5rh6KltoaKgiIyP9XgAAAAAAALCLp02ptm3baseOHX7zvv76a8XFxUmSateurZiYGC1fvty3/I8//tDq1auVmJjoaq4AAAAAAAAoOp7++t59992nxMREjR8/Xtddd50+/vhjPfPMM3rmmWck/Xnb3pAhQzR+/HglJCQoISFB48ePV3h4uPr06eNl6gAAAAAAADgLnjalWrVqpUWLFmnEiBEaM2aMateuralTp6pv376+mGHDhun48eMaPHiwDh8+rNatW2vZsmWqUKGCh5kDAAAAAADgbDjGGON1EsXpyJEjioqKUnp6Os+XAgAAAAAAOcQPfytfcbse61bMmZQO+e3FePpMKQAAAAAAAAQmmlIAAAAAAABwHU0pAAAAAAAAuI6mFAAAAAAAAFxHUwoAAAAAAACuC/Y6AQAAAAAAgJKEX+srGlwpBQAAAAAAANfRlAIAAAAAAIDraEoBAAAAAADAdTSlAAAAAAAA4DoedA4AAAAAAEqV/DyInIeQe48rpQAAAAAAAOA6mlIAAAAAAABwHU0pAAAAAAAAuI6mFAAAAAAAAFxHUwoAAAAAAACuoykFAAAAAAAA19GUAgAAAAAAgOtoSgEAAAAAAMB1NKUAAAAAAADgOppSAAAAAAAAcB1NKQAAAAAAALiOphQAAAAAAABcR1MKAAAAAAAArqMpBQAAAAAAANfRlAIAAAAAAIDraEoBAAAAAADAdTSlAAAAAAAA4DqaUgAAAAAAAHAdTSkAAAAAAAC4jqYUAAAAAAAAXEdTCgAAAAAAAK6jKQUAAAAAAADX0ZQCAAAAAACA62hKAQAAAAAAwHU0pQAAAAAAAOA6mlIAAAAAAABwHU0pAAAAAAAAuI6mFAAAAAAAAFxHUwoAAAAAAACuoykFAAAAAAAA19GUAgAAAAAAgOtoSgEAAAAAAMB1NKUAAAAAAADgOppSAAAAAAAAcJ2nTamUlBQ5juP3iomJ8S03xiglJUWxsbEqV66ckpOTtXXrVg8zBgAAAAAAQFHw/Eqpxo0bKzU11ff64osvfMsmTZqkKVOmaPr06dq4caNiYmLUqVMnHT161MOMAQAAAAAAcLY8b0oFBwcrJibG96pataqkP6+Smjp1qkaOHKmePXuqSZMmmjdvnn777TfNnz/f46wBAAAAAABwNjxvSn3zzTeKjY1V7dq11atXL33//feSpJ07dyotLU2dO3f2xYaGhiopKUnr1q3zKl0AAAAAAAAUgWAvd966dWu98MILqlevng4cOKCxY8cqMTFRW7duVVpamiQpOjrab53o6Gjt3r07z21mZGQoIyPDN33kyJHiSR4AAAAAAACF5mlTqmvXrr7/btq0qdq0aaO6detq3rx5uvjiiyVJjuP4rWOMyTHvVBMmTNDo0aOLJ2EAAAAAAAAUCc9v3ztVRESEmjZtqm+++cb3K3zZV0xlO3jwYI6rp041YsQIpaen+1579+4t1pwBAAAAAABQcFY1pTIyMrR9+3ZVr15dtWvXVkxMjJYvX+5b/scff2j16tVKTEzMcxuhoaGKjIz0ewEAAAAAAMAunt6+98ADD6h79+6qVauWDh48qLFjx+rIkSMaMGCAHMfRkCFDNH78eCUkJCghIUHjx49XeHi4+vTp42XaAAAAAAAAOEueNqV++OEH9e7dWz/99JOqVq2qiy++WBs2bFBcXJwkadiwYTp+/LgGDx6sw4cPq3Xr1lq2bJkqVKjgZdoAAAAAAAA4S542pRYsWHDa5Y7jKCUlRSkpKe4kBAAAAAAAAFdY9UwpAAAAAAAABAaaUgAAAAAAAHAdTSkAAAAAAAC4jqYUAAAAAAAAXEdTCgAAAAAAAK6jKQUAAAAAAADX0ZQCAAAAAACA62hKAQAAAAAAwHU0pQAAAAAAAOA6mlIAAAAAAABwHU0pAAAAAAAAuC7Y6wQAAAAAAADyEj/8rXzF7XqsWzFngqLGlVIAAAAAAABwHU0pAAAAAAAAuI6mFAAAAAAAAFxHUwoAAAAAAACuoykFAAAAAAAA19GUAgAAAAAAgOtoSgEAAAAAAMB1NKUAAAAAAADgOppSAAAAAAAAcB1NKQAAAAAAALiOphQAAAAAAABcR1MKAAAAAAAArqMpBQAAAAAAANfRlAIAAAAAAIDraEoBAAAAAADAdTSlAAAAAAAA4DqaUgAAAAAAAHAdTSkAAAAAAAC4jqYUAAAAAAAAXEdTCgAAAAAAAK6jKQUAAAAAAADXBXudAAAAAAAACBzxw9/KV9yux7oVcybwGldKAQAAAAAAwHU0pQAAAAAAAOA6mlIAAAAAAABwHU0pAAAAAAAAuI6mFAAAAAAAAFxHUwoAAAAAAACuoykFAAAAAAAA19GUAgAAAAAAgOtoSgEAAAAAAMB1NKUAAAAAAADgOppSAAAAAAAAcF2w1wnALvHD38pX3K7HuhVzJijp8jOWGEcAAAAAELisuVJqwoQJchxHQ4YM8c0zxiglJUWxsbEqV66ckpOTtXXrVu+SBAAAAAAAQJGwoim1ceNGPfPMM2rWrJnf/EmTJmnKlCmaPn26Nm7cqJiYGHXq1ElHjx71KFMAAAAAAAAUBc+bUr/++qv69u2rZ599VpUqVfLNN8Zo6tSpGjlypHr27KkmTZpo3rx5+u233zR//nwPMwYAAAAAAMDZ8rwpdeedd6pbt27q2LGj3/ydO3cqLS1NnTt39s0LDQ1VUlKS1q1bl+f2MjIydOTIEb8XAAAAAAAA7OLpg84XLFigTz75RBs3bsyxLC0tTZIUHR3tNz86Olq7d+/Oc5sTJkzQ6NGjizZRlGqB+HD3QKzZDTzcveQJxL+FQKwZAAAAdvLsSqm9e/fq3nvv1UsvvaSwsLA84xzH8Zs2xuSYd6oRI0YoPT3d99q7d2+R5QwAAAAAAICi4dmVUps3b9bBgwd14YUX+uZlZmZqzZo1mj59unbs2CHpzyumqlev7os5ePBgjqunThUaGqrQ0NDiSxwAAAAAAABnzbMrpTp06KAvvvhCW7Zs8b1atmypvn37asuWLapTp45iYmK0fPly3zp//PGHVq9ercTERK/SBgAAAAAAQBHw7EqpChUqqEmTJn7zIiIiVLlyZd/8IUOGaPz48UpISFBCQoLGjx+v8PBw9enTx4uUAQAAAAAAUEQ8fdD5mQwbNkzHjx/X4MGDdfjwYbVu3VrLli1ThQoVvE4NAAAAAAAAZ8GqptSqVav8ph3HUUpKilJSUjzJBwAAAAAAAMXDs2dKAQAAAAAAIHDRlAIAAAAAAIDraEoBAAAAAADAdTSlAAAAAAAA4DqaUgAAAAAAAHCdVb++B5RW8cPfOmPMrse6Fev2z3YftikNNRf3uChupeE9cENJf58BAACA4sKVUgAAAAAAAHAdTSkAAAAAAAC4jqYUAAAAAAAAXEdTCgAAAAAAAK7jQeclTGl4sHBpqAFFr7jHRUG3b+M4LQ01FFQg1lxQtv3tAAAAAPnFlVIAAAAAAABwHU0pAAAAAAAAuI6mFAAAAAAAAFxHUwoAAAAAAACuoykFAAAAAAAA1/Hreyh13PilqPzso7i3f7b7AAAAAADAS1wpBQAAAAAAANfRlAIAAAAAAIDraEoBAAAAAADAdYVqStWpU0eHDh3KMf+XX35RnTp1zjopAAAAAAAAlG6Fakrt2rVLmZmZOeZnZGRo3759Z50UAAAAAAAASrcC/frekiVLfP/97rvvKioqyjedmZmpFStWKD4+vsiSAyR+iQ4AAAAAgNKoQE2pv/3tb5Ikx3E0YMAAv2UhISGKj4/X448/XmTJAQAAAAAAoHQqUFMqKytLklS7dm1t3LhRVapUKZakAAAAAAAAULoVqCmVbefOnUWdBwAAAAAAAAJIoZpSkrRixQqtWLFCBw8e9F1BlW3OnDlnnRgAAAAAAABKr0I1pUaPHq0xY8aoZcuWql69uhzHKeq8AAAISPn5cQd+2AEAAAClQaGaUk899ZSef/559e/fv6jzAQAAAAAAQAAIKsxKf/zxhxITE4s6FwAAAAAAAASIQjWlbr75Zs2fP7+ocwEAAAAAAECAKNTte7///rueeeYZvffee2rWrJlCQkL8lk+ZMqVIkgMAAAAAAEDpVKim1Oeff64WLVpIkr788ku/ZTz0HAAAAAAAAGdSqKbUypUrizoPFJP8/IqTdHa/5MQvRZV+bowjeI/PCwAAAABuKtQzpQAAAAAAAICzUagrpS677LLT3qb3/vvvFzohAAAAAAAAlH6FakplP08q24kTJ7RlyxZ9+eWXGjBgQFHkBQAAAAAAgFKsUE2pJ554Itf5KSkp+vXXX88qIQAAAAAAAJR+RfpMqX79+mnOnDlFuUkAAAAAAACUQoW6Uiov69evV1hYWFFuEgAAlCD8WicAoLTh3AYUn0I1pXr27Ok3bYxRamqqNm3apIcffrhIEgMAAAAAAEDpVaimVFRUlN90UFCQ6tevrzFjxqhz585FkhgAAAAAAABKr0I1pebOnVvUeQAAAAAAACCAnNUzpTZv3qzt27fLcRw1atRI559/flHlBQAAAAAAgFKsUE2pgwcPqlevXlq1apUqVqwoY4zS09N12WWXacGCBapatWpR5wkAgOd40CkAAEBOfEdCYQUVZqW7775bR44c0datW/Xzzz/r8OHD+vLLL3XkyBHdc889+d7OrFmz1KxZM0VGRioyMlJt2rTRO++841tujFFKSopiY2NVrlw5JScna+vWrYVJGQAAAAAAABYpVFNq6dKlmjVrlho2bOib16hRI82YMcOvqXQmNWrU0GOPPaZNmzZp06ZNat++va6++mpf42nSpEmaMmWKpk+fro0bNyomJkadOnXS0aNHC5M2AAAAAAAALFGoplRWVpZCQkJyzA8JCVFWVla+t9O9e3ddccUVqlevnurVq6dx48apfPny2rBhg4wxmjp1qkaOHKmePXuqSZMmmjdvnn777TfNnz+/MGkDAAAAAADAEoVqSrVv31733nuv9u/f75u3b98+3XffferQoUOhEsnMzNSCBQt07NgxtWnTRjt37lRaWpo6d+7siwkNDVVSUpLWrVtXqH0AAAAAAADADoV60Pn06dN19dVXKz4+XjVr1pTjONqzZ4+aNm2ql156qUDb+uKLL9SmTRv9/vvvKl++vBYtWqRGjRr5Gk/R0dF+8dHR0dq9e3ee28vIyFBGRoZv+siRIwXKBwAAAAAAAMWvUE2pmjVr6pNPPtHy5cv11VdfyRijRo0aqWPHjgXeVv369bVlyxb98ssvev311zVgwACtXr3at9xxHL94Y0yOeaeaMGGCRo8eXeA8AAAAAAAAv6YH9xTo9r33339fjRo18l191KlTJ919992655571KpVKzVu3FgffPBBgRIoW7aszjvvPLVs2VITJkxQ8+bNNW3aNMXExEiS0tLS/OIPHjyY4+qpU40YMULp6em+1969ewuUDwAAAAAAAIpfgZpSU6dO1S233KLIyMgcy6KionTbbbdpypQpZ5WQMUYZGRmqXbu2YmJitHz5ct+yP/74Q6tXr1ZiYmKe64eGhioyMtLvBQAAAAAAALsUqCn12Wef6fLLL89zeefOnbV58+Z8b++hhx7SBx98oF27dumLL77QyJEjtWrVKvXt21eO42jIkCEaP368Fi1apC+//FIDBw5UeHi4+vTpU5C0AQAAAAAAYJkCPVPqwIEDCgkJyXtjwcH68ccfC7S9/v37KzU1VVFRUWrWrJmWLl2qTp06SZKGDRum48ePa/DgwTp8+LBat26tZcuWqUKFCgVJGwAAAAAAAJYpUFPq3HPP1RdffKHzzjsv1+Wff/65qlevnu/tzZ49+7TLHcdRSkqKUlJSCpImAAAAAAAALFegptQVV1yhRx55RF27dlVYWJjfsuPHj2vUqFG68sorizRBAAAAAABQePn5NT1+SQ9eKFBT6p///KcWLlyoevXq6a677lL9+vXlOI62b9+uGTNmKDMzUyNHjiyuXAEAAAAAAFBKFKgpFR0drXXr1umOO+7QiBEjZIyR9Odtdl26dNHMmTMVHR1dLIkCAAAAAACg9ChQU0qS4uLi9Pbbb+vw4cP69ttvZYxRQkKCKlWqVBz5AQAAAAAAoBQqcFMqW6VKldSqVauizAUAAAAAAAABIsjrBAAAAAAAABB4aEoBAAAAAADAdTSlAAAAAAAA4DqaUgAAAAAAAHAdTSkAAAAAAAC4jqYUAAAAAAAAXBfsdQIAAAAAAJQW8cPfylfcrse6FXMmgP24UgoAAAAAAACuoykFAAAAAAAA19GUAgAAAAAAgOtoSgEAAAAAAMB1NKUAAAAAAADgOppSAAAAAAAAcB1NKQAAAAAAALiOphQAAAAAAABcR1MKAAAAAAAArqMpBQAAAAAAANfRlAIAAAAAAIDraEoBAAAAAADAdTSlAAAAAAAA4DqaUgAAAAAAAHAdTSkAAAAAAAC4jqYUAAAAAAAAXBfsdQIAAAAAgNIpfvhb+Yrb9Vi3Ys4EgI24UgoAAAAAAACuoykFAAAAAAAA19GUAgAAAAAAgOtoSgEAAAAAAMB1NKUAAAAAAADgOn59DwAAAAACFL+OB8BLXCkFAAAAAAAA19GUAgAAAAAAgOtoSgEAAAAAAMB1NKUAAAAAAADgOh50DgAAUMoU94OL3Xgwcn72wYOX7VIaxkVp+NsJRHxeACUXV0oBAAAAAADAdTSlAAAAAAAA4DqaUgAAAAAAAHAdTSkAAAAAAAC4jqYUAAAAAAAAXMev7wEAgBKlNPx6Fb8UhdwwLpAbxgWA0szTK6UmTJigVq1aqUKFCqpWrZr+9re/aceOHX4xxhilpKQoNjZW5cqVU3JysrZu3epRxgAAAAAAACgKnjalVq9erTvvvFMbNmzQ8uXLdfLkSXXu3FnHjh3zxUyaNElTpkzR9OnTtXHjRsXExKhTp046evSoh5kDAAAAAADgbHh6+97SpUv9pufOnatq1app8+bNateunYwxmjp1qkaOHKmePXtKkubNm6fo6GjNnz9ft912mxdpAwAAAAAA4CxZ9aDz9PR0SdI555wjSdq5c6fS0tLUuXNnX0xoaKiSkpK0bt26XLeRkZGhI0eO+L0AAAAAAABgF2sedG6M0dChQ3XJJZeoSZMmkqS0tDRJUnR0tF9sdHS0du/enet2JkyYoNGjRxdvsgAAoNQqDQ9SD0TF/b6VhnFR0BpKQ81AbhjbgD2suVLqrrvu0ueff65XXnklxzLHcfymjTE55mUbMWKE0tPTfa+9e/cWS74AAAAAAAAoPCuulLr77ru1ZMkSrVmzRjVq1PDNj4mJkfTnFVPVq1f3zT948GCOq6eyhYaGKjQ0tHgTBgAAAAAAwFnx9EopY4zuuusuLVy4UO+//75q167tt7x27dqKiYnR8uXLffP++OMPrV69WomJiW6nCwAAAAAAgCLi6ZVSd955p+bPn6///e9/qlChgu8ZUlFRUSpXrpwcx9GQIUM0fvx4JSQkKCEhQePHj1d4eLj69OnjZeoAAAAAAAA4C542pWbNmiVJSk5O9ps/d+5cDRw4UJI0bNgwHT9+XIMHD9bhw4fVunVrLVu2TBUqVHA5WwAAAAAAABQVT5tSxpgzxjiOo5SUFKWkpBR/QgAAoNThV5bgBcYdgPzi8wKBzJpf3wMAAAAAAEDgoCkFAAAAAAAA19GUAgAAAAAAgOtoSgEAAAAAAMB1NKUAAAAAAADgOk9/fQ8AACA/vzpk8y8OlYZfTSoNNQQi/na834eNfzslfVwgf3ifUVpwpRQAAAAAAABcR1MKAAAAAAAArqMpBQAAAAAAANfRlAIAAAAAAIDreNA5AACA5Ur6A21tfBg0ADvxeYHSirGdO66UAgAAAAAAgOtoSgEAAAAAAMB1NKUAAAAAAADgOppSAAAAAAAAcB1NKQAAAAAAALiOX98DAABAwOFXkIDC4W8HQFHiSikAAAAAAAC4jqYUAAAAAAAAXEdTCgAAAAAAAK6jKQUAAAAAAADX0ZQCAAAAAACA62hKAQAAAAAAwHU0pQAAAAAAAOA6mlIAAAAAAABwHU0pAAAAAAAAuI6mFAAAAAAAAFxHUwoAAAAAAACuoykFAAAAAAAA19GUAgAAAAAAgOtoSgEAAAAAAMB1NKUAAAAAAADgOppSAAAAAAAAcB1NKQAAAAAAALiOphQAAAAAAABcR1MKAAAAAAAArqMpBQAAAAAAANfRlAIAAAAAAIDraEoBAAAAAADAdTSlAAAAAAAA4DqaUgAAAAAAAHAdTSkAAAAAAAC4jqYUAAAAAAAAXEdTCgAAAAAAAK6jKQUAAAAAAADX0ZQCAAAAAACA6zxtSq1Zs0bdu3dXbGysHMfR4sWL/ZYbY5SSkqLY2FiVK1dOycnJ2rp1qzfJAgAAAAAAoMh42pQ6duyYmjdvrunTp+e6fNKkSZoyZYqmT5+ujRs3KiYmRp06ddLRo0ddzhQAAAAAAABFKdjLnXft2lVdu3bNdZkxRlOnTtXIkSPVs2dPSdK8efMUHR2t+fPn67bbbnMzVQAAAAAAABQha58ptXPnTqWlpalz586+eaGhoUpKStK6devyXC8jI0NHjhzxewEAAAAAAMAu1jal0tLSJEnR0dF+86Ojo33LcjNhwgRFRUX5XjVr1izWPAEAAAAAAFBw1jalsjmO4zdtjMkx71QjRoxQenq677V3797iThEAAAAAAAAF5OkzpU4nJiZG0p9XTFWvXt03/+DBgzmunjpVaGioQkNDiz0/AAAAAAAAFJ61V0rVrl1bMTExWr58uW/eH3/8odWrVysxMdHDzAAAAAAAAHC2PL1S6tdff9W3337rm965c6e2bNmic845R7Vq1dKQIUM0fvx4JSQkKCEhQePHj1d4eLj69OnjYdYAAAAAAAA4W542pTZt2qTLLrvMNz106FBJ0oABA/T8889r2LBhOn78uAYPHqzDhw+rdevWWrZsmSpUqOBVygAAAAAAACgCnjalkpOTZYzJc7njOEpJSVFKSop7SQEAAAAAAKDYWftMKQAAAAAAAJReNKUAAAAAAADgOppSAAAAAAAAcB1NKQAAAAAAALiOphQAAAAAAABcR1MKAAAAAAAArqMpBQAAAAAAANfRlAIAAAAAAIDraEoBAAAAAADAdTSlAAAAAAAA4DqaUgAAAAAAAHAdTSkAAAAAAAC4jqYUAAAAAAAAXEdTCgAAAAAAAK6jKQUAAAAAAADX0ZQCAAAAAACA62hKAQAAAAAAwHU0pQAAAAAAAOA6mlIAAAAAAABwHU0pAAAAAAAAuI6mFAAAAAAAAFxHUwoAAAAAAACuoykFAAAAAAAA19GUAgAAAAAAgOtoSgEAAAAAAMB1NKUAAAAAAADgOppSAAAAAAAAcB1NKQAAAAAAALiOphQAAAAAAABcR1MKAAAAAAAArqMpBQAAAAAAANfRlAIAAAAAAIDraEoBAAAAAADAdTSlAAAAAAAA4DqaUgAAAAAAAHAdTSkAAAAAAAC4jqYUAAAAAAAAXEdTCgAAAAAAAK6jKQUAAAAAAADX0ZQCAAAAAACA62hKAQAAAAAAwHU0pQAAAAAAAOA6mlIAAAAAAABwHU0pAAAAAAAAuI6mFAAAAAAAAFxXIppSM2fOVO3atRUWFqYLL7xQH3zwgdcpAQAAAAAA4CxY35R69dVXNWTIEI0cOVKffvqpLr30UnXt2lV79uzxOjUAAAAAAAAUkvVNqSlTpuimm27SzTffrIYNG2rq1KmqWbOmZs2a5XVqAAAAAAAAKCSrm1J//PGHNm/erM6dO/vN79y5s9atW+dRVgAAAAAAADhbwV4ncDo//fSTMjMzFR0d7Tc/OjpaaWlpua6TkZGhjIwM33R6erok6ciRI8WXqIuyMn7LV1x2vcUdn991Chp/NjmV1Jo5RnbVYEvNHCO7arClZo6RXTXYUjPHyK4abKmZY2RXDbbUzDGyqwZbauYY2VdDSZZdhzHmtHGOOVOEh/bv369zzz1X69atU5s2bXzzx40bpxdffFFfffVVjnVSUlI0evRoN9MEAAAAAADAX+zdu1c1atTIc7nVV0pVqVJFZcqUyXFV1MGDB3NcPZVtxIgRGjp0qG86KytLP//8sypXrizHcYo1Xy8cOXJENWvW1N69exUZGVni4m3MiZq9j7cxJ2qmZi/ibcyJmr2PtzEnavY+3sacqNn7eBtzombv423MiZrzV0NJYozR0aNHFRsbe9o4q5tSZcuW1YUXXqjly5erR48evvnLly/X1Vdfnes6oaGhCg0N9ZtXsWLF4kzTCpGRkQUayLbF25gTNXsfb2NO1Fz08TbmZFu8jTlRs/fxNuZEzd7H25gTNXsfb2NO1Ox9vI05UXPpEhUVdcYYq5tSkjR06FD1799fLVu2VJs2bfTMM89oz549uv32271ODQAAAAAAAIVkfVPq+uuv16FDhzRmzBilpqaqSZMmevvttxUXF+d1agAAAAAAACgk65tSkjR48GANHjzY6zSsFBoaqlGjRuW4ZbGkxNuYEzV7H29jTtRc9PE25mRbvI05UbP38TbmRM3ex9uYEzV7H29jTtTsfbyNOVFz4LL61/cAAAAAAABQOgV5nQAAAAAAAAACD00pAAAAAAAAuI6mFAAAAAAAAFxHUwoAAAAAAACuoymFIrdq1SodP3682LafkZGh7777ThkZGcW2jwMHDigtLe20MZmZmTpw4IB++umnfG83e52DBw8qMzPzbNMs9YpzLNkyjqSCjyXGUcGU9HEkFc9nEuOoYBhHp49nHOUP35FOH884yh/G0ZnXYSzlT0k/t/Fdu5QwKDW2bdtmateu7Tdvy5Yt5tFHHzUzZswwP/74o9+y9PR0c+ONN/rNe/bZZ80NN9xg5syZY4wxZsGCBaZBgwamdu3a5pFHHslXHiEhIWbbtm25LtuxY4fJysryTX/wwQfm6quvNo0aNTIdOnQwixcv9oufO3euWb9+vTHGmOPHj5ubbrrJlClTxgQFBZng4GBz2223md9//90X36RJEzNmzBizZ8+efOV66NAh07NnT1OrVi0zePBgc/LkSXPTTTcZx3FMUFCQadOmjdm/f7/fOm+++aa59NJLTWhoqAkKCjJBQUEmKirK9OvXz+zevTvX/SxcuNAkJiaasmXL+tYpW7asSUxMNIsWLcpXrtmK4n0u7u0XxTgyJu+xVBrGkTEFH0slfRwVZh/F+Zlk2zgyxp3PpJI0jopqHyV5HBlj57mtuMeRMaXz3FZU48gYO89tXo4jYwLnOxLftXMqSec2W8aRMSXnOxLftUs3mlKlyJYtW0xQUJBv+t133zVly5Y1jRs3NrVq1TJVqlQx77//vm95WlqaX/wTTzxhIiIiTM+ePU316tXN2LFjTeXKlc3YsWPNmDFjTFRUlHn66ad98eeff36uL8dxTMOGDX3TpwoKCjIHDhwwxhizcuVKExQUZLp3727GjRtnrrnmGhMUFGSWLl3qiz/vvPPMxo0bjTHGPPDAAyY+Pt4sXLjQbN++3SxevNjUq1fP/OMf//DFO45jKleubMqUKWO6dOliXnvtNXPixIk8j9mNN95omjRpYp588kmTlJRk/va3v5lmzZqZtWvXmnXr1plWrVqZG264wRf/wgsvmAoVKpghQ4aY4cOHm+joaDN8+HAza9Ysk5SUZKpUqWK+/vprv3089dRTpmzZsub22283ixYtMuvWrTMffvihWbRokbn99ttNaGioeeaZZ07/5p7ibN/n4t5+QceRMQUfSyV9HBlT8LFU0sdRYfZR3J9Jto0jY4r/M6mkjaOi2EdJH0fG2HduK+5xZEzJP7cV9zgyxr5zm9fjyJjS9x2J79p8185t+3zX5rt2aUNTqgS57777Tvvq16+f30Bu06aNeeihh4wxxmRlZZlJkyaZ8uXLm3feeccYk3PgN2jQwLz88svGGGM++eQTExwcbJ577jnf8jlz5pgLL7zQNx0cHGwuv/xyk5KS4nuNGjXKBAUFmcGDB/vmncpxHN8HXIcOHczgwYP9lg8fPty0a9fONx0aGurrYterV8+Xe7bVq1ebWrVq+W1/3759ZtGiRaZ79+4mODjYVK1a1dx///25/l+A6tWrmw8//NB3PBzHMcuWLfMtX7t2rTn33HP9jtGCBQt80xs3bjQ1atTw/Z+E66+/3vTo0cNvH3Xr1vU7jn81e/ZsU6dOHd90cb/Pto0jYwo+lkr6OMo+TgUZSyV9HBVmH8X9mWTbODKm+D+TbBtHbuyjpI+j7H3YdG4r7nFkTMk/txX3ODLGvnNbcY8jY+w7t9k2joyx79zGd23vx5Ex9p3b+K6dv+9IgYKmVAkSFBRkLrjgApOcnJzrq2XLln4DOTIy0nz77bd+25g/f76JiIgwS5YsyTHwy5Ur53cZY2hoqPnyyy990998842pWLGib3rt2rWmbt265pFHHjGZmZm++cHBwWbr1q251nDqB1z16tXNhg0b/JZv3brVVK5c2TcdFxfn6yCfe+65vg58tm3btpmIiIhct2+MMampqWb8+PEmISHBd2nn7NmzfcvDw8PNrl27fNMhISHmiy++8E1///33ftsvV66c2blzp18OwcHBZt++fcYYYz766CO/Y2SMMWFhYearr77K9XgYY8z27dtNWFiYb7q432fbxpExBR9LJX0cZR+ngoylkj6OCrOP4v5Msm0cGVP8n0m2jSM39lHSx9Ff92GM9+e24h5HxpT8c1txjyNj7Du3Ffc4Msa+c5tt48gY+85tfNf2fhwZY9+5je/afzrTd6RAQVOqBKlfv7558cUX81z+6aef+g3kqlWrmk2bNuWIW7BggQkPDzezZs3yi69cubJfZ7pGjRp+f/zffPONKV++vN+20tPTTa9evcxFF13k+yM704ny22+/Nenp6aZOnTrm008/9Vv+zTffmPDwcN/0Qw89ZNq0aWMOHz5shg8fbrp3726OHj1qjDHm2LFj5rrrrjOdO3f2xZ96qelfrVy50vTr18/vA6t58+Zm+vTpxhhj3n77bVOhQgXz+OOP+5bPmjXLNGnSxDfdsGFD89///tc3vXnzZlO2bFlz8uRJX/5//UC88MILzdChQ3PNyRhjhg4d6vd/M4r7fbZxHBlTsLFU0seRMQUfSyV9HBVmH8X9mWTbODKm+D+TbBtHbuyjpI8jY+w7txX3ODKmdJzbinMcGWPfua24x5Ex9p3bbBtHxth3buO7dtFvn+/agfFdO5DQlCpB+vTpY4YMGZLn8i1bthjHcXzTnTp1MpMnT841dv78+SYkJMRv4Ldt29bvEse/euONN3L8sWebM2eOiYmJMU8//bQJCQk57YkyKOjPh8Y5jpPjEsnFixebhIQE33RGRoa56qqrTKVKlUynTp1MWFiYCQ8PNwkJCSYiIsLUqlXL7Nixw2/7eX3AZUtPT/f990svvWTKlCljzjvvPBMWFmZee+01Exsba6677jrTq1cvU7ZsWd8HoDHGTJ8+3URFRZlhw4aZRx55xMTGxpqbbrrJb3t/vbd/1apVJiIiwjRq1MgMGTLETJgwwTz22GNmyJAhpnHjxqZ8+fJmzZo1vvjifp9tHkfG5G8slfRxZEzBx1JJH0eF2UdxfybZNo6MKf7PJNvGkRv7KOnjKHsfNp3binscGVO6zm3FMY6Mse/cVtzjyBj7zm22jSNj7Du38V276LfPd+3A+K4dSGhKlSCpqal+XfAzWbhw4Wn/UObPn2+Sk5N902vXrs3RBT/VjBkzzJNPPpnn8q+//tq0atXKOI6T54ly1apVfq+/fvGeOnWqmTRpUo713nnnHTN48GBz+eWXm86dO5sBAwaYZ555xvz6669+cQMHDjRHjhzJM8fcfPDBB+Zf//qXWbdunTHmz8tR+/fvb6655hrz/PPP54ifOXOmSUxMNBdeeKF56KGHzPHjx/2Owfbt23Oss3PnTjNs2DDTrl07U69ePVOvXj3Trl078+CDD+a4rLS432fbx5ExZx5LpWEcGVPwsVSSx1Fh9lHcn0k2jiNjiv8zyaZx5MY+Svo4MsbOc1txjiNjSt+5rbjGkTF2nduKcxwZY9+5zbZxZIyd5za+axft9vmuHRjftQOJY4wxAopIVlaWjh49qsjISDmO43U6KMEYSygKjCMUBcYRigLjCEWBcYSiwliCLYK9TgCFs3v3bqWlpclxHEVHRysuLs6q+KioqBJfQ1HH28jGY1TQsWRbDUVRc0kbS27UXNyfSaVhXJT0cSTZ9z7YNo7c2AfjqOSNIy9qYBx5H8937ZLBxmNk27nNxnMnXOb1pVoomClTppgaNWr47u/Nvt+3Ro0a5oknnrA+3sac3Kj5dLZs2VKg+4eLIt7GY2RbTjbWfDolYRwVZp1Ai3drH3nxYhwZY9/7YFu8rTnlpbjHUV7r2HaMAnFceDmOCrNOaRxHNubEd207jpFtOdlY8+l49R0pENCUKkHGjBljIiMjzWOPPWY+/fRTs3//frNv3z7z6aefmscee8xERUWZRx991Np4G3Nyo+Yzye0Br8UZb+Mxsi0nG2su6Ptc3PG8z3aMi5I+jgpTQ6DF25pTQd/noozPbR3bjlEgjguvx1Fh1ilt48jGnLweR7m9z8Udb+Mxsi0nG2su6PvsdXxpQlOqBKlRo4ZZtGhRnssXLlxoYmNjrY23MSc3au7Ro8dpX+3bt/frihd3vI3HyLacbKy5pI+jwqwTaPFu7MO2ceRGzSU93sac3BgXJf3cFojjorjHUWHWCbRxZGNOfNe24xjZlpONNds2jgIJz5QqQQ4dOqT69evnubxevXo6fPiwtfE25uRGzW+88YY6deqk6OjoXNfJzMx0Nd7GY2RbTjbWXNLHUWHWCbR4N/Zh2zgqTA2BFm9jTm6Mi5J+bgvEcVHc46gw6wTaOLIxJ75rn7kGxkXRxxdmHdvGUUDxuiuG/EtKSjJ9+/Y1J06cyLHsxIkTpk+fPiYpKcnaeBtzcqPmpk2bmueeey5HfLZPP/3Uryte3PE2HiPbcrKx5pI+jgqzTqDFu7EP28ZRYWoItHgbc3JjXJT0c1sgjoviHkeFWSfQxpGNOblRs23nNhuPkW052VizbeMokHClVAny5JNPqnPnzqpWrZqSkpIUHR0tx3GUlpamNWvWKDQ0VMuXL7c23sac3Kj5wgsv1CeffKKbbrop1/c1NDRUtWrVci3exmNkW0421lzSx5GN74Nt8W7sw7ZxZOP7YFu8jTm5MS5K+rktEMdFcY+jwqwTaOPIxpz4rm3HMbItJxtrtm0cBRLHGGO8TgL5d/ToUb300kvasGGD0tLSJEkxMTFq06aN+vTpo8jISKvjbcypuOMzMjKUmZmp8PDwHMciN8UdX5gaGBfe11waxlFh1gm0+OLeh43jqKA1BGK8bTm5MS5Kw7kt0MZFQePdGBeBOI5szInv2t4fIxtzsq1mG8dRoKApBQAAAAAAANcFeZ0Azk63bt2UmppaYuNtzImavY+3MSdqLvp4G3OyLd7GnKjZ+3gbc6Jm7+NtzImavY+3MSdq9j7expyoOYB5+0grnK3y5cub7777rsTG25gTNXsfb2NO1Fz08TbmZFu8jTlRs/fxNuZEzd7H25gTNXsfb2NO1Ox9vI05UXPg4kopAAAAAAAAuI6mVAkXFxenkJCQEhtvY07U7H28jTlRc9HH25iTbfE25kTN3sfbmBM1ex9vY07U7H28jTlRs/fxNuZEzYGLB50DAAAAAADAdcFeJ4CC+/7777V27VqlpqaqTJkyql27tjp16pTrz2DaGG9jToFY8zfffKN169YpLS1NjuMoOjpaiYmJSkhI8CTexpyoOX815ObYsWPavHmz2rVrV2zrBFq8jTm5UTNQEJmZmSpTpoxv+qOPPlJGRobatGmT6/+NLunxNuZUGmr+qxtvvFHjxo1TbGzsGWNtjLcxp0Cr+fDhw/r2229VvXp11ahRw/N4G3Oi5vzVUCp5/VAr5N+vv/5qrr32WuM4jnEcxwQFBZmYmBhTpkwZU758eTN9+nSr423MKRBr/uWXX8xVV11lHMcxFStWNPXq1TMJCQmmYsWKJigoyFx99dUmPT3dtXgbc6Lm/NVwOlu2bDFBQUH5ji/MOoEWb2NORRH/xx9/mH/84x+mbt26plWrVmbOnDl+y9PS0vzWCbR4G3Oyseb9+/ebtm3bmjJlyph27dqZn3/+2XTr1s13bqxXr57Zv39/qYm3MafSUPNnn32W6yskJMQsWrTIN21rvI05BWLNI0aMMMeOHTPG/PlZdsstt5igoCDf9/QePXqY48ePuxZvY07UnL8aAgVNqRLk1ltvNW3btjVbtmwxX331lbnmmmvMsGHDzLFjx8zs2bNNeHi4efnll62NtzGnQKy5f//+pmnTpmbDhg05xtiGDRtMs2bNzA033OBavI05UXP+ajgdGi5FH29jTkURP2rUKBMdHW0mT55sRo4caaKiosytt97qW56WlmYcxwnYeBtzsrHm/v37m8TERLNkyRJz/fXXm8TERHPppZeaH374wezZs8dceuml5s477yw18TbmVBpqzv7HYXbT6tTXqf94tDXexpwCseagoCBz4MABY4wx48aNM1WrVjWvv/662bdvn3njjTfMueeea8aMGeNavI05UXP+aggUNKVKkCpVqphNmzb5pn/++WcTFhbm67hOnz7dtGjRwtp4G3MKxJqjoqJybTxkW79+vYmKinIt3sacqPnM8ZUqVTrtKzIyMseXxoKuE2jxNubkRs3nnXeeeeONN3zT3377rUlISDADBw40WVlZOa6ICbR4G3Oysebq1aub9evXG2OMOXTokHEcx7z33nu+5e+//76pU6dOqYm3MafSUHPz5s1Nt27dzPbt282uXbvMrl27zM6dO01wcLBZvny5b56t8TbmFIg1O47jaz60aNHCzJ492y/fV1991TRs2NC1eBtzoub81RAoeKZUCXLy5Em/5wOVL19eJ0+e1LFjxxQeHq7OnTvrgQcesDbexpwCsWZJchxHecltWXHH25gTNZ9+WUZGhu644w41bdo01/jdu3dr9OjRZ7VOoMXbmJMbNe/bt09NmjTxTdetW1erVq1S+/bt1b9/f02aNCmg423MycaaDx8+rHPPPVeSdM455yg8PFxxcXF+66emppaaeBtzKg01f/zxxxo2bJiuueYavfTSSzr//PN9y2JjY/3WtTHexpwCsWbp/39v2rt3ry666CK/ZRdddJF2797taryNOVFz/moICF53xZB/nTp18rvEePLkyaZ69eq+6U8++cRUqVLF2ngbcwrEmvv162eaNWtmNm7caP5q48aNpkWLFqZ///6uxduYEzWfOT4xMdFMnTo1R2y23G7TKug6gRZvY05u1Fy7dm2/Kxey7du3z9SrV8907NjRb51Ai7cxJxtrrlWrlvnoo4980w8++KA5dOiQb3rLli1+58KSHm9jTqWh5mxvv/22qVGjhhk/frzJzMw0wcHBZuvWrTnibI23MadAqtlxHDNu3Dgzbdo0Exsba9asWeO3fMuWLaZSpUquxduYEzXnr4ZAQVOqBNm8ebM555xzTExMjKlVq5YpW7aseeWVV3zLp0+f7vfMF9vibcwpEGs+fPiwufzyy43jOKZSpUqmfv36pkGDBqZSpUomKCjIdO3a1Rw+fNi1eBtzouYzx48bN86kpKSYvOzZs8cMHDjQb15B1wm0eBtzcqPmm266yQwaNCjX+B9++MGcd955fs2HQIu3MScba77qqqtO2xCdPn26ad++famJtzGn0lDzqdLS0kzXrl3NJZdckq+GiG3xNuYUKDXHxcWZ+Ph43+uvY/CJJ54wF198sWvxNuZEzfmrIVA4xhjj9dVayL/U1FS9+eabysjIUPv27dWoUaMSFW9jToFYsyR99dVXWr9+vdLS0iRJMTExatOmjRo0aOBJvI05UXP+agDOxu7du/XVV1+pS5cuuS5PTU3VsmXLNGDAgICMtzEnG2s+k40bN6pcuXJ+twSW5ngbcyqpNf/73//WypUr9eSTT+brJ9tti7cxp0Cs+VQbNmxQaGio322AXsbbmBM1BxaaUgAAAAAAAHBdkNcJoOgcPnxYL7zwQomNtzGn0lxzVlZWrvFZWVnas2eP6/E25kTNRR9vY062xduYEzV7H29jTtTsfbyNOVGz9/E25kTN3sfbmBM1nzk+IHh79yCKUm4PkC1J8TbmVBprTk9PN3//+99NWFiYqVatmnnkkUfMyZMnfcv/+jPbxR1vY07UTM0cI2q2Jd7GnKjZ+3gbc6Jm7+NtzImavY+3MSdqzl8NgYKmVAmSnp5+2tcHH3yQY+DbFG9jToFY8z333GPq1atn/vvf/5pnn33WxMXFmW7dupmMjAxjzJ8fiI7juBZvY07UTM0cI2q2Jd7GnKjZ+3gbc6Jm7+NtzImavY+3MSdqzl8NgYKmVAniOI4JCgrK85W93NZ4G3MKxJpr1aplVq5c6Zv+6aefTOvWrU3nzp3N77//nqNLX9zxNuZEzdTMMaJmW+JtzImavY+3MSdq9j7expyo2ft4G3Oi5vzVEChoSpUgkZGRZuLEiWbVqlW5vp599lm/gWxbvI05BWLN4eHh5vvvv/fL8ciRI6ZNmzamffv25vvvv3c13sacqJmaOUbUbEu8jTlRs/fxNuZEzd7H25gTNXsfb2NO1Jy/GgIFTakSJDk52UycODHP5Vu2bPG75M+2eBtzCsSa69evb956660ccUePHjVt2rQxzZs39/tALO54G3OiZmrmGFGzLfE25kTN3sfbmBM1ex9vY07U7H28jTlRc/5qCBT8+l4J0qdPH4WFheW5PCYmRqNGjbI23sacArHmzp07a+7cuTniypcvr3fffTfHtoo73sacqLno423MybZ4G3OiZu/jbcyJmr2PtzEnavY+3sacqNn7eBtzoub81RAwvO6KAXDXzz//bL788ss8lx89etSsWrXKtXgbc6Lmoo+3MSfb4m3MiZq9j7cxJ2r2Pt7GnKjZ+3gbc6Jm7+NtzIma81dDoHCMMcbrxhgAAAAAAAACS7DXCaBgjh07pvnz52vdunVKS0uT4ziKjo5W27Zt1bt3b0VERFgdb2NO1Ox9vI05UTM1c4yo2ZZ4G3OiZu/jbcyJmr2PtzEnavY+3sacqDl/NQQCrpQqQbZt26ZOnTrpt99+U1JSkqKjo2WM0cGDB7V69WpFRERo2bJlatSokZXxNuZEzd7H25gTNVMzx4iabYm3MSdq9j7expyo2ft4G3OiZu/jbcyJmvNXQ8DI/51+8FpycrLp1auXycjIyLEsIyPD9O7d2yQnJ1sbb2NO1Ox9vI05UTM1c4yo2ZZ4G3OiZu/jbcyJmr2PtzEnavY+3sacqDl/NQQKmlIlSLly5czWrVvzXP7FF1+YcuXKWRtvY07U7H28jTlRc9HH25iTbfE25kTN3sfbmBM1ex9vY07U7H28jTlRs/fxNuZEzfmrIVAEeX2lFvKvUqVK+uabb/Jc/u2336pSpUrWxtuYEzV7H29jTtRc9PE25mRbvI05UbP38TbmRM3ex9uYEzV7H29jTtTsfbyNOVFz/moIGF53xZB/o0aNMlFRUWby5Mlmy5YtJjU11aSlpZktW7aYyZMnm0qVKpnRo0dbG29jTtTsfbyNOVEzNXOMqNmWeBtzombv423MiZq9j7cxJ2r2Pt7GnKg5fzUECppSJcxjjz1mqlevbhzHMUFBQSYoKMg4jmOqV69uJk6caH28jTlRs/fxNuZEzdTMMaJmW+JtzImavY+3MSdq9j7expyo2ft4G3Oi5vzVEAj49b0SaufOnUpLS5MkxcTEqHbt2iUq3sacqNn7eBtzomZq9iLexpyo2ft4G3OiZu/jbcyJmr2PtzEnavY+3sacqDl/NZRmNKUAAAAAAADgOh50XsIcP35ca9eu1bZt23Is+/333/XCCy9YHW9jTtTsfbyNOVFz0cfbmJNt8TbmRM3ex9uYEzV7H29jTtTsfbyNOVGz9/E25kTN+ashIHh79yAKYseOHSYuLs53D2pSUpLZv3+/b3laWpoJCgqyNt7GnKjZ+3gbc6JmauYYUbMt8TbmRM3ex9uYEzV7H29jTtTsfbyNOVFz/moIFFwpVYI8+OCDatq0qQ4ePKgdO3YoMjJSbdu21Z49e0pEvI05UbP38TbmRM3UzDGiZlvibcyJmr2PtzEnavY+3sacqNn7eBtzoub81RAwvO6KIf+qVatmPv/8c795gwcPNrVq1TLfffddju6qbfE25kTN3sfbmBM1UzPHiJptibcxJ2r2Pt7GnKjZ+3gbc6Jm7+NtzIma81dDoKApVYJUqFDBbNu2Lcf8u+66y9SoUcOsWbPGbyDbFm9jTtTsfbyNOVEzNXOMqNmWeBtzombv423MiZq9j7cxJ2r2Pt7GnKg5fzUECppSJUirVq3MCy+8kOuyO++801SsWNFvINsWb2NO1Ox9vI05UXPRx9uYk23xNuZEzd7H25gTNXsfb2NO1Ox9vI05UbP38TbmRM35qyFQ0JQqQcaPH2+6du2a5/I77rjDOI5jbbyNOVGz9/E25kTNRR9vY062xduYEzV7H29jTtTsfbyNOVGz9/E25kTN3sfbmBM156+GQOEYY4zXz7UCAAAAAABAYOHX9wAAAAAAAOA6mlIAAAAAAABwHU0pAAAAAAAAuI6mFAAAAAAAAFxHUwoAAAAAAACuoykFAACAM1q1apUcx9Evv/zidSoAAKCUoCkFAAACijFGHTt2VJcuXXIsmzlzpqKiorRnzx5Xc6pdu7aWLl2a67L4+Hg5jqMFCxbkWNa4cWM5jqPnn3++SPNJTk7WkCFDinSbAAAAf0VTCgAABBTHcTR37lx99NFHevrpp33zd+7cqQcffFDTpk1TrVq1inSfJ06cyHPZ559/rkOHDumyyy7LM6ZmzZqaO3eu37wNGzYoLS1NERERRZYnAACAm2hKAQCAgFOzZk1NmzZNDzzwgHbu3CljjG666SZ16NBBF110ka644gqVL19e0dHR6t+/v3766SffukuXLtUll1yiihUrqnLlyrryyiv13Xff+Zbv2rVLjuPoP//5j5KTkxUWFqaXXnopz1z+97//qUuXLgoNDc0zpm/fvlq9erX27t3rmzdnzhz17dtXwcHBfrF79uzR1VdfrfLlyysyMlLXXXedDhw44FuekpKiFi1a6MUXX1R8fLyioqLUq1cvHT16VJI0cOBArV69WtOmTZPjOHIcR7t27fKtv3nzZrVs2VLh4eFKTEzUjh07znzAAQAAckFTCgAABKQBAwaoQ4cOuvHGGzV9+nR9+eWXmjZtmpKSktSiRQtt2rRJS5cu1YEDB3Tdddf51jt27JiGDh2qjRs3asWKFQoKClKPHj2UlZXlt/0HH3xQ99xzj7Zv357rrYLZlixZoquvvvq0uUZHR6tLly6aN2+eJOm3337Tq6++qkGDBvnFGWP0t7/9TT///LNWr16t5cuX67vvvtP111/vF/fdd99p8eLFevPNN/Xmm29q9erVeuyxxyRJ06ZNU5s2bXTLLbcoNTVVqampqlmzpm/dkSNH6vHHH9emTZsUHBycIwcAAID8Cj5zCAAAQOn0zDPPqEmTJvrggw/02muvafbs2brgggs0fvx4X8ycOXNUs2ZNff3116pXr56uueYav23Mnj1b1apV07Zt29SkSRPf/CFDhqhnz56n3f++ffv02Wef6YorrjhjroMGDdL999+vkSNH6rXXXlPdunXVokULv5j33ntPn3/+uXbu3OlrJL344otq3LixNm7cqFatWkmSsrKy9Pzzz6tChQqSpP79+2vFihUaN26coqKiVLZsWYWHhysmJiZHHuPGjVNSUpIkafjw4erWrZt+//13hYWFnbEGAACAU3GlFAAACFjVqlXTrbfeqoYNG6pHjx7avHmzVq5cqfLly/teDRo0kCTfLXrfffed+vTpozp16igyMlK1a9eWpBwPR2/ZsuUZ979kyRK1bdtW55xzzhlju3Xrpl9//VVr1qzRnDlzcr1Cafv27apZs6bflU2NGjVSxYoVtX37dt+8+Ph4X0NKkqpXr66DBw+eMQdJatasmd96kvK9LgAAwKm4UgoAAAS04OBg33OZsrKy1L17d02cODFHXHYDpnv37qpZs6aeffZZxcbGKisrS02aNNEff/zhF5+fB5Dn59a9U/Ps37+/Ro0apY8++kiLFi3KEWOMkeM4Z5wfEhLit9xxnBy3H+bl1HWzt5nfdQEAAE5FUwoAAOD/XHDBBXr99dcVHx+f4wHiknTo0CFt375dTz/9tC699FJJ0tq1awu1r19//VUrV67UjBkz8r3OoEGD9K9//UvXX3+9KlWqlGN5o0aNtGfPHu3du9d3tdS2bduUnp6uhg0b5ns/ZcuWVWZmZr7jAQAACoPb9wAAAP7PnXfeqZ9//lm9e/fWxx9/rO+//17Lli3ToEGDlJmZqUqVKqly5cp65pln9O233+r999/X0KFDC7WvpUuXKiEhQXXq1Mn3Og0bNtRPP/2kuXPn5rq8Y8eOatasmfr27atPPvlEH3/8sW644QYlJSXl63bCbPHx8froo4+0a9cu/fTTT1wJBQAAigVNKQAAgP8TGxurDz/8UJmZmerSpYuaNGmie++9V1FRUQoKClJQUJAWLFigzZs3q0mTJrrvvvs0efLkQu3rf//7X75v3TtV5cqVVa5cuVyXOY6jxYsXq1KlSmrXrp06duyoOnXq6NVXXy3QPh544AGVKVNGjRo1UtWqVXM8LwsAAKAoOMYY43USAAAAgSQzM1PVqlXTO++8o4suusjrdAAAADzBlVIAAAAuO3TokO677z61atXK61QAAAA8w5VSAAAAAAAAcB1XSgEAAAAAAMB1NKUAAAAAAADgOppSAAAAAAAAcB1NKQAAAAAAALiOphQAAAAAAABcR1MKAAAAAAAArqMpBQAAAAAAANfRlAIAAAAAAIDraEoBAAAAAADAdf8PyoaW2V/jIhQAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plot_date_histogram(items, group_by=(\"year\", \"month\"))" + ] + }, + { + "cell_type": "markdown", + "id": "31dd76bc-39d9-457e-bb2e-b916b438928e", + "metadata": {}, + "source": [ + "## Bounding Boxes via flatgeobuf\n", + "Printing boudning boxes of the items that are returnded from a STAC search gives an overview of the spatial coverage.\n", + "\n", + "This function writes a flatgeobuf file to further be used for interacting with the bbox geometries. Flatgeobuf is great for fast visualization of many polygons." + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "9b6dc0f1-8775-42b2-a859-b50d33ede675", + "metadata": {}, + "outputs": [], + "source": [ + "import geopandas as gpd\n", + "from shapely.geometry import box\n", + "\n", + "def create_flatgeobuf_from_items(items, out_path, properties_to_keep=None):\n", + " \"\"\"\n", + " Create a FlatGeobuf (.fgb) file with bounding-box polygons from pystac.Item objects.\n", + "\n", + " Parameters\n", + " ----------\n", + " items : list[pystac.Item]\n", + " List of STAC items. Each must have `.bbox` and `.properties`.\n", + " out_path : str\n", + " Output FlatGeobuf path.\n", + " properties_to_keep : list[str], optional\n", + " If not None, only include these properties in the output.\n", + " \"\"\"\n", + " \n", + " records = []\n", + "\n", + " for item in items:\n", + " # Ensure bbox exists\n", + " if item.bbox is None or len(item.bbox) != 4:\n", + " continue\n", + "\n", + " minx, miny, maxx, maxy = item.bbox\n", + " poly = box(minx, miny, maxx, maxy)\n", + "\n", + " # Select properties\n", + " if properties_to_keep:\n", + " props = {k: item.properties.get(k) for k in properties_to_keep}\n", + " else:\n", + " # keep all properties, but ensure everything is serializable\n", + " props = {}\n", + " for k, v in item.properties.items():\n", + " if isinstance(v, (str, int, float, bool)) or v is None:\n", + " props[k] = v\n", + " else:\n", + " props[k] = str(v)\n", + "\n", + " props[\"geometry\"] = poly\n", + " records.append(props)\n", + "\n", + " # GeoDataFrame\n", + " gdf = gpd.GeoDataFrame(records, geometry=\"geometry\", crs=\"EPSG:4326\")\n", + "\n", + " # Write FlatGeobuf\n", + " gdf.to_file(out_path, driver=\"FlatGeobuf\")\n", + "\n", + " print(f\"Created {len(gdf)} bounding-box polygons in {out_path}\")" + ] + }, + { + "cell_type": "markdown", + "id": "6eff3537-f484-4759-86e5-8c1d786e0c90", + "metadata": {}, + "source": [ + "This function simply maps the bounding boxes to give a first indication." + ] + }, + { + "cell_type": "code", + "execution_count": 49, + "id": "ce749314-e935-4a21-9291-e5403ae02081", + "metadata": {}, + "outputs": [], + "source": [ + "from typing import List, Optional, Sequence, Union, Dict, Any\n", + "import geopandas as gpd\n", + "from shapely.geometry import box, mapping\n", + "import folium\n", + "import math\n", + "\n", + "def render_map_from_flatgeobuf(\n", + " fgb_path: str,\n", + " center: Optional[Sequence[float]] = None,\n", + " zoom_start: int = 2,\n", + " simplify_tolerance: Optional[float] = None,\n", + " style_function: Optional[callable] = None\n", + ") -> folium.Map:\n", + " \"\"\"\n", + " Render interactive Folium map from a FlatGeobuf of polygons.\n", + "\n", + " Parameters:\n", + " fgb_path: str\n", + " Path to the FlatGeobuf file created above.\n", + " center: tuple(float, float) or None\n", + " (lat, lon) center for the initial map view. If None, center is computed from dataset centroid.\n", + " zoom_start: int\n", + " Initial zoom level.\n", + " simplify_tolerance: float or None\n", + " If provided, geometries are simplified (units in degrees) to improve performance while rendering.\n", + " Try values like 0.001 (small) to 0.01 (stronger).\n", + " style_function: callable or None\n", + " Optional function that receives a feature and returns a style dict for folium.GeoJson (e.g. color, weight).\n", + "\n", + " Returns:\n", + " folium.Map\n", + " The interactive folium Map containing the polygon layer.\n", + " \"\"\"\n", + "\n", + " # Read the flatgeobuf\n", + " gdf = gpd.read_file(fgb_path)\n", + "\n", + " if gdf.empty:\n", + " raise ValueError(\"FlatGeobuf file contains no features.\")\n", + "\n", + " # Simplify geometry to speed up rendering\n", + " if simplify_tolerance is not None and simplify_tolerance > 0:\n", + " # preserve topology for polygons where possible\n", + " gdf[\"geometry\"] = gdf[\"geometry\"].simplify(simplify_tolerance, preserve_topology=True)\n", + "\n", + " # Compute center if not provided (folium wants [lat, lon])\n", + " if center is None:\n", + " bounds = gdf.total_bounds # [minx, miny, maxx, maxy]\n", + " minx, miny, maxx, maxy = bounds\n", + " center = [(miny + maxy) / 2.0, (minx + maxx) / 2.0]\n", + "\n", + " m = folium.Map(location=center, zoom_start=zoom_start, tiles=\"CartoDB positron\")\n", + "\n", + " # Convert GeoDataFrame to GeoJSON-like mapping in memory and add to map\n", + " # Use folium.GeoJson; for many features use a FeatureGroup to keep it organized\n", + " fg = folium.FeatureGroup(name=\"Bounding boxes\", overlay=True, control=True)\n", + "\n", + " # default style function if none provided\n", + " if style_function is None:\n", + " def style_function_default(feature):\n", + " return {\n", + " \"fillColor\": \"#1f9e89\",\n", + " \"color\": \"#000000\",\n", + " \"weight\": 1,\n", + " \"fillOpacity\": 0.4,\n", + " }\n", + " style_function = style_function_default\n", + "\n", + " # Add polygons to the FeatureGroup\n", + " folium.GeoJson(\n", + " gdf,\n", + " style_function=style_function\n", + " ).add_to(fg)\n", + "\n", + " m.add_child(fg)\n", + " folium.LayerControl().add_to(m)\n", + " return m" + ] + }, + { + "cell_type": "markdown", + "id": "0c8dee6f-76f9-4bef-bc0d-75a905b91063", + "metadata": {}, + "source": [ + "Here the two functions are applied in sequence to create a flateobuf file and then map the polygons on an interactive map. \n", + "\n", + "*Note: In cases where a long timeseries is requested there are usually many overlaps. In this case the next visualization might help.* " + ] + }, + { + "cell_type": "code", + "execution_count": 50, + "id": "71a43c2e-5076-462c-ae1c-e6ca018b6946", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Created 2483 bounding-box polygons in ./test_s1slc.fgb\n" + ] + }, + { + "data": { + "text/html": [ + "
Make this Notebook Trusted to load map: File -> Trust Notebook
" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "create_flatgeobuf_from_items(items = items, out_path = \"./stacsearch_bboxes.fgb\", properties_to_keep=None)\n", + "m = render_map_from_flatgeobuf(\"stacsearch_bboxes.fgb\", simplify_tolerance=0.001)\n", + "#m.save(\"bboxes_map.html\")\n", + "display(m)" + ] + }, + { + "cell_type": "markdown", + "id": "a4e9c678-2d2f-41a9-9ab8-30b6661748ae", + "metadata": {}, + "source": [ + "## Heatmap via flatgeobuf\n", + "If there are many overlaps between the bounding boxes a heatmap can show where the hotspots of the counts of overlapping bounding boxes are.\n", + "\n", + "This function computes a heatmap from a flatgeobuf file created with the function `create_flatgeobuf_from_items()` " + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "id": "4da291e6-6700-4d5f-b9f9-5740ec6f170a", + "metadata": {}, + "outputs": [], + "source": [ + "import geopandas as gpd\n", + "import numpy as np\n", + "\n", + "def compute_bbox_heatmap(\n", + " fgb_path,\n", + " grid_size=50\n", + "):\n", + " \"\"\"\n", + " Computes a fast heatmap (bbox density grid) from a FlatGeobuf file of bounding boxes.\n", + "\n", + " Parameters: \n", + " fgb_path: path to the flatgeobuf file\n", + " grid_size: number of the grids to create in x and y direction (rasterization)\n", + " Returns:\n", + " heatmap: 2D numpy array (counts)\n", + " extent: (minx, miny, maxx, maxy)\n", + " \"\"\"\n", + " gdf = gpd.read_file(fgb_path)\n", + "\n", + " # dataset bounds\n", + " minx, miny, maxx, maxy = gdf.total_bounds\n", + "\n", + " # prepare grid\n", + " x_bins = np.linspace(minx, maxx, grid_size + 1)\n", + " y_bins = np.linspace(miny, maxy, grid_size + 1)\n", + "\n", + " heatmap = np.zeros((grid_size, grid_size), dtype=np.uint32)\n", + "\n", + " # fast loop through bboxes\n", + " for geom in gdf.geometry:\n", + " bxmin, bymin, bxmax, bymax = geom.bounds\n", + "\n", + " ix1 = np.searchsorted(x_bins, bxmin) - 1\n", + " ix2 = np.searchsorted(x_bins, bxmax) - 1\n", + " iy1 = np.searchsorted(y_bins, bymin) - 1\n", + " iy2 = np.searchsorted(y_bins, bymax) - 1\n", + "\n", + " # increment all overlapping grid cells\n", + " heatmap[iy1:iy2+1, ix1:ix2+1] += 1\n", + "\n", + " return heatmap, (minx, miny, maxx, maxy)" + ] + }, + { + "cell_type": "markdown", + "id": "1cf78157-a075-418c-bb0f-ff31b9dda49f", + "metadata": {}, + "source": [ + "This function renders the heatmap created by `compute_bbox_heatmap()`." + ] + }, + { + "cell_type": "code", + "execution_count": 45, + "id": "6871af86-7813-4077-9c3c-6f9688ae1c43", + "metadata": {}, + "outputs": [], + "source": [ + "import folium\n", + "import matplotlib.pyplot as plt\n", + "from folium.raster_layers import ImageOverlay\n", + "\n", + "def map_bbox_heatmap(heatmap, extent, colormap=\"viridis\", opacity=0.6):\n", + " \"\"\"\n", + " Renders the heatmap as an ImageOverlay on a folium interactive map.\n", + "\n", + " Parameters: \n", + " heatmap: heatmap created by compute_bbox_heatmap\n", + " extent: extent created by compute_bbox_heatmap\n", + " colormap: a colormap known to matplotlib\n", + " opacity: opacity between 0 and 1\n", + " \"\"\"\n", + " minx, miny, maxx, maxy = extent\n", + "\n", + " fig, ax = plt.subplots(figsize=(6,6))\n", + " ax.axis(\"off\")\n", + "\n", + " # generate image array\n", + " img = plt.cm.get_cmap(colormap)(heatmap / heatmap.max())\n", + "\n", + " # convert to uint8 for overlay\n", + " img_uint8 = (img[:, :, :3] * 255).astype(\"uint8\")\n", + "\n", + " m = folium.Map(\n", + " location=[(miny + maxy) / 2, (minx + maxx) / 2],\n", + " zoom_start=3,\n", + " tiles=\"CartoDB positron\"\n", + " )\n", + "\n", + " ImageOverlay(\n", + " image=img_uint8,\n", + " bounds=[[miny, minx], [maxy, maxx]],\n", + " opacity=opacity,\n", + " ).add_to(m)\n", + "\n", + " folium.LayerControl().add_to(m)\n", + " return m" + ] + }, + { + "cell_type": "markdown", + "id": "e57054d0-5310-4d95-92d1-228345c6a329", + "metadata": {}, + "source": [ + "Execute the function `compute_bbox_heatmap()` and `map_bbox_heatmap()`." + ] + }, + { + "cell_type": "code", + "execution_count": 32, + "id": "49246dac-8044-4405-8926-20db9cc2dc91", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/dss/dsstbyfs02/scratch/06/di38laz/di38laz/tmp.rkjo7BaLxG/ipykernel_161795/3598274246.py:15: MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed in 3.11. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap()`` or ``pyplot.get_cmap()`` instead.\n", + " img = plt.cm.get_cmap(colormap)(heatmap / heatmap.max())\n" + ] + }, + { + "data": { + "text/html": [ + "
Make this Notebook Trusted to load map: File -> Trust Notebook
" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeQAAAHiCAYAAAA597/kAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAB6RJREFUeJzt1TEBACAMwDDAv+chYz0SBf16Z2YOALDqbQcAAIYMAAmGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAYYMAAGGDAABhgwAAR/JoAfASfVs8wAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "heatmap, extent = compute_bbox_heatmap(\"stacsearch_bboxes.fgb\", grid_size=200)\n", + "m = map_bbox_heatmap(heatmap, extent)\n", + "#m.save(\"bbox_heatmap.html\")\n", + "display(m)" + ] + } + ], + "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.11.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}