MMW #3638 Error Demo

3638-1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pystac_client import Client\n",
    "from rasterio.coords import BoundingBox\n",
    "from rasterio.enums import Resampling\n",
    "from rasterio.mask import mask\n",
    "from rasterio.merge import merge\n",
    "from rasterio.warp import (\n",
    "    calculate_default_transform,\n",
    "    transform_geom,\n",
    "    reproject,\n",
    ")\n",
    "from shapely.geometry import mapping, shape, box\n",
    "import json\n",
    "import numpy as np\n",
    "import os\n",
    "import rasterio as rio\n",
    "import tempfile\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_albers_crs_for_aoi(aoi):\n",
    "    \"\"\"\n",
    "    Return the Albers Equal Area Projection for the given AoI\n",
    "\n",
    "    Since we want to calculate area, we need to use an Equal Area projection,\n",
    "    but this differs based on where you are in the globe. We use rough bounding\n",
    "    boxes to see if the shape neatly fits within one of the continents. If not,\n",
    "    we fall back to a global approximation.\n",
    "    \"\"\"\n",
    "\n",
    "    if aoi.within(box(-170, 15, -50, 75)):  # North America\n",
    "        return 'EPSG:5070'\n",
    "    elif aoi.within(box(-10, 34, 40, 72)):  # Europe\n",
    "        return 'EPSG:3035'\n",
    "    elif aoi.within(box(25, -10, 180, 60)):  # Asia\n",
    "        return 'ESRI:102025'\n",
    "    elif aoi.within(box(-20, -35, 55, 38)):  # Africa\n",
    "        return 'ESRI:102022'\n",
    "    elif aoi.within(box(-90, -60, -30, 15)):  # South America\n",
    "        return 'ESRI:102033'\n",
    "    elif aoi.within(box(112, -45, 155, -10)):  # Australia\n",
    "        return 'ESRI:102034'\n",
    "    else:  # Global\n",
    "        return 'ESRI:54017'  # Behrmann Equal Area Cylindrical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "def clip_and_reproject_tile(tiff, aoi, dst_crs, resampling=Resampling.nearest):\n",
    "    \"\"\"\n",
    "    Clip a tiff to the given AoI and reproject to destination CRS\n",
    "\n",
    "    We clip the tiff first, and then reproject only the clipped bits to\n",
    "    minimize transformation. Uses nearest neighbor resampling by default.\n",
    "\n",
    "    Returns the data, transform, and meta after clipping and reprojecting.\n",
    "    \"\"\"\n",
    "    aoi_crs = aoi.crs if hasattr(aoi, 'crs') else 4326\n",
    "\n",
    "    with rio.open(tiff) as src:\n",
    "        tiff_bounds = box(*src.bounds)  # Create a Shapely box for GeoTIFF bounds\n",
    "        tiff_crs = src.crs\n",
    "\n",
    "        # Reproject AoI to GeoTIFF CRS\n",
    "        reprojected_aoi = transform_geom(aoi_crs, src.crs, aoi)\n",
    "        aoi_geom = shape(reprojected_aoi)\n",
    "        \n",
    "        # Create GeoDataFrames for plotting\n",
    "        gdf_tiff = gpd.GeoDataFrame({\"geometry\": [tiff_bounds]}, crs=tiff_crs)\n",
    "        gdf_aoi = gpd.GeoDataFrame({\"geometry\": [aoi_geom]}, crs=tiff_crs)\n",
    "\n",
    "        # Plotting\n",
    "        fig, ax = plt.subplots(figsize=(10, 10))\n",
    "\n",
    "        # Plot GeoTIFF bounds\n",
    "        gdf_tiff.plot(ax=ax, edgecolor=\"blue\", facecolor=\"none\", linewidth=2, label=\"GeoTIFF Bounds\")\n",
    "\n",
    "        # Plot reprojected AoI\n",
    "        gdf_aoi.plot(ax=ax, edgecolor=\"green\", facecolor=\"none\", linewidth=2, label=\"Reprojected AoI\")\n",
    "\n",
    "        # Customize the plot\n",
    "        ax.set_title(\"GeoTIFF Bounds and Reprojected AoI\")\n",
    "        ax.set_xlabel(\"Longitude\")\n",
    "        ax.set_ylabel(\"Latitude\")\n",
    "        plt.legend()\n",
    "        plt.grid(True)\n",
    "        plt.show()\n",
    "\n",
    "        reprojected_aoi_bbox = BoundingBox(*shape(reprojected_aoi).bounds)\n",
    "        clip_data, clip_transform = mask(src, [reprojected_aoi], crop=True)\n",
    "\n",
    "        # Define the output metadata for the reprojected clip\n",
    "        dst_transform, width, height = calculate_default_transform(\n",
    "            src.crs,\n",
    "            dst_crs,\n",
    "            clip_data.shape[2],\n",
    "            clip_data.shape[1],\n",
    "            *reprojected_aoi_bbox,\n",
    "        )\n",
    "        dst_meta = src.meta.copy()\n",
    "        dst_meta.update(\n",
    "            {\n",
    "                'crs': dst_crs,\n",
    "                'transform': dst_transform,\n",
    "                'width': width,\n",
    "                'height': height,\n",
    "            }\n",
    "        )\n",
    "\n",
    "        # Reproject the clipped data to the destination CRS\n",
    "        reprojected_clip_data = np.empty(\n",
    "            (src.count, height, width), dtype=src.dtypes[0]\n",
    "        )\n",
    "        reproject(\n",
    "            source=clip_data,\n",
    "            destination=reprojected_clip_data,\n",
    "            src_transform=clip_transform,\n",
    "            src_crs=src.crs,\n",
    "            dst_transform=dst_transform,\n",
    "            dst_crs=dst_crs,\n",
    "            resampling=resampling,\n",
    "        )\n",
    "\n",
    "        return reprojected_clip_data, dst_transform, dst_meta\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "def query_histogram(geojson, url, collection, asset, filter, cachekey=''):\n",
    "    aoi = shape(json.loads(geojson))\n",
    "\n",
    "    # Get list of intersecting tiffs from catalog\n",
    "    client = Client.open(url)\n",
    "    search = client.search(\n",
    "        collections=collection,\n",
    "        intersects=mapping(aoi),\n",
    "        filter=filter,\n",
    "    )\n",
    "    tiffs = [item.assets[asset].href for item in search.items()]\n",
    "\n",
    "    # Sort the tiffs, since the merging is order sensitive. We use the \"first\"\n",
    "    # method, using the first image for overlapping pixels.\n",
    "    tiffs = sorted(tiffs)\n",
    "\n",
    "    # Raise error if no overlapping tiffs are found\n",
    "    if not tiffs:\n",
    "        return {\n",
    "            'error': (\n",
    "                f'No overlapping tiffs found in collection {collection} '\n",
    "                f'with filter {filter} '\n",
    "                f'with AoI {geojson[:255]} ...'\n",
    "            )\n",
    "        }\n",
    "\n",
    "    # Find the Albers Equal Area CRS for this AoI\n",
    "    dst_crs = get_albers_crs_for_aoi(aoi)\n",
    "\n",
    "    # Reproject the tiffs and clip to the AoI to make tiles\n",
    "    clips = []\n",
    "    for tiff in tiffs:\n",
    "        clip_data, clip_transform, clip_meta = clip_and_reproject_tile(\n",
    "            tiff, aoi, dst_crs\n",
    "        )\n",
    "        temp_file = tempfile.NamedTemporaryFile(delete=False, suffix='.tif')\n",
    "        with rio.open(temp_file.name, 'w', **clip_meta) as dst:\n",
    "            dst.write(clip_data)\n",
    "        clips.append(temp_file.name)\n",
    "\n",
    "    # Merge the clipped rasters\n",
    "    datasets = [rio.open(clip_path) for clip_path in clips]\n",
    "    merged_data, merged_transform = merge(datasets, method='first')\n",
    "\n",
    "    # Count how many of each class type there are\n",
    "    values, counts = np.unique(merged_data, return_counts=True)\n",
    "    histogram = list(zip(values.tolist(), counts.tolist()))\n",
    "\n",
    "    # Ensure the result is JSON serializable\n",
    "    histogram = {int(value): int(count) for value, count in histogram}\n",
    "\n",
    "    # Calculate pixel size from dataset resolution\n",
    "    pixel_size = datasets[0].res[0] * datasets[0].res[1]\n",
    "\n",
    "    # Close datasets\n",
    "    for ds in datasets:\n",
    "        ds.close()\n",
    "\n",
    "    # Clean up temporary files\n",
    "    for temp_file in clips:\n",
    "        os.remove(temp_file)\n",
    "\n",
    "    result = {\n",
    "        'result': histogram,\n",
    "        'pixel_size': pixel_size,\n",
    "    }\n",
    "\n",
    "    return result\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "URL = \"https://api.impactobservatory.com/stac-aws\"\n",
    "COLLECTION = \"io-10m-annual-lulc\"\n",
    "ASSET = \"supercell\"\n",
    "FILTER={\n",
    "    \"op\": \"like\",\n",
    "    \"args\": [{\"property\": \"id\"}, f\"%-2023\"],\n",
    "}\n",
    "geojson='{\"type\":\"MultiPolygon\",\"coordinates\":[[[[-75.19678870333927,39.96553758512753],[-75.18505772457821,39.96553758512753],[-75.18505772457821,39.97452797890017],[-75.19678870333927,39.97452797890017],[-75.19678870333927,39.96553758512753]]]]}'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_452544/634836532.py:37: UserWarning: Legend does not support handles for PatchCollection instances.\n",
      "See: https://matplotlib.org/stable/tutorials/intermediate/legend_guide.html#implementing-a-custom-legend-handler\n",
      "  plt.legend()\n",
      "/tmp/ipykernel_452544/634836532.py:37: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.\n",
      "  plt.legend()\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkEAAANXCAYAAAAy/IufAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABQsklEQVR4nO3deVhV5f7//9dmBhFBcUBEcZ6nLA2H1BTNbKCyQTlpWnkq+2p6KvWUqZVpgx3rVGaDQ6V50uaOZaSi6VEzS1MzS8MhxSlTQGQQ7t8f/NgftwzCFtjq/XxcFxfute611nu/WbhfrGFvhzHGCAAAwDJeni4AAADAEwhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAysTcuXPlcDi0e/duT5dSYomJiXI4HEpMTPR0KeVq9+7dcjgcmjt3rqdLKbVJkybJ4XB4ugxcoghBqFBJSUl68MEH1aRJEwUFBSkoKEgtWrTQiBEj9NNPP5Xbdnv06CGHw3HOr0mTJkmSoqOjdd1117mso6hlatWq5RyT/x92YV+vv/56qWr08/NT/fr1NXz4cO3bt6/Me4KSyw94+V8+Pj6KjIzUXXfdpf3793u6PI87cOCAJk2apE2bNnm6FEnSbbfdJofDobFjx57XevJ/7t9//30ZVYYLjY+nC4A9vvjiC91+++3y8fFRfHy82rZtKy8vL/3yyy/66KOPNHPmTCUlJalevXplvu3HHntM99xzj/Pxhg0b9PLLL+uf//ynmjdv7pzepk2bYtcTGxurwYMHu0wLDAwsMG7mzJkKDg52mdapU6dz1lmnTh1NnTpVkpSVlaWff/5Zr7/+upYuXart27crKCjonOtA+XnyySdVv359ZWRkaN26dZo7d65Wr16trVu3KiAgwNPlFalevXo6deqUfH19y2X9Bw4c0OTJkxUdHa127dqVyzZKKiUlRZ9//rmio6P1/vvva9q0aRxJQpEIQagQu3bt0h133KF69epp2bJlioiIcJn/7LPP6rXXXpOXV/kcnIyNjXV5HBAQoJdfflmxsbHq0aNHidfTpEkT/e1vfzvnuAEDBig8PLy0ZapKlSoF1l+/fn09+OCDWrNmTYHngYrVr18/XX755ZKke+65R+Hh4Xr22Wf12Wef6bbbbquwOowxysjIKDSAF8bhcFzQIa0sffjhh8rJydHs2bN19dVXa9WqVerevbuny8IFitNhqBDPPfecTp48qTlz5hQIQJLk4+OjkSNHKioqymX6L7/8ogEDBqhq1aoKCAjQ5Zdfrs8++6zA8r///rtuvfVWVa1aVUFBQbryyiv13//+t9yeT0XKP93m4+P6N8uPP/6ofv36KSQkRMHBwerVq5fWrVvnMqao6ykKu34n/xTg6tWr1bFjRwUEBKhBgwZ65513Ciy/bds2XX311QoMDFSdOnX09NNPKzc3t8C477//Xn379lV4eLgCAwNVv359DRs27JzP+dNPP1X//v1Vu3Zt+fv7q2HDhnrqqaeUk5PjMq5Hjx5q1aqVfv75Z/Xs2VNBQUGKjIzUc889V2Cdf/zxh+Li4lSpUiXVqFFDo0ePVmZm5jlrKU63bt0k5YX8M5Vkv83/GaxatUp///vfVa1aNYWEhGjw4MH666+/XMbm/2yWLl2qyy+/XIGBgZo1a5akku37RV0TVNLfr+PHj2v06NGKjo6Wv7+/6tSpo8GDB+vo0aNKTEzUFVdcIUkaOnSo85Thmdtav369rrnmGlWpUkVBQUHq3r271qxZU2A7q1ev1hVXXKGAgAA1bNjQ+RxLY/78+YqNjVXPnj3VvHlzzZ8/v9Bxy5cvV7du3VSpUiWFhobqxhtv1Pbt20u9PVzcOBKECvHFF1+oUaNGJTollG/btm3q0qWLIiMjNW7cOFWqVEkffPCB4uLi9OGHH+qmm26SJB06dEidO3dWenq6Ro4cqWrVqmnevHm64YYbtHjxYue4spCRkaGjR4+6TKtcubL8/f1dph07dszlsbe3t8LCws65/pycHOf6s7OztX37dk2cOFGNGjVSly5dnOO2bdumbt26KSQkRI8++qh8fX01a9Ys9ejRQytXrixVn8+0c+dODRgwQHfffbeGDBmi2bNn66677lKHDh3UsmVLSdLBgwfVs2dPnT592vlzeeONNwoclTh8+LD69Omj6tWra9y4cQoNDdXu3bv10UcfnbOOuXPnKjg4WGPGjFFwcLCWL1+uJ554QikpKXr++eddxv7111+65pprdPPNN+u2227T4sWLNXbsWLVu3Vr9+vWTJJ06dUq9evXS3r17NXLkSNWuXVvvvvuuli9f7laf8uWHyDN/tiXdb/M9+OCDCg0N1aRJk7Rjxw7NnDlTe/bscV60nW/Hjh0aOHCg/v73v+vee+9V06ZNz2vfL2mdaWlp6tatm7Zv365hw4bpsssu09GjR/XZZ5/pjz/+UPPmzfXkk0/qiSee0PDhw53BsHPnzpLywka/fv3UoUMHTZw4UV5eXpozZ46uvvpqffvtt+rYsaMkacuWLc79ZdKkSTp9+rQmTpyomjVrlvjnceDAAa1YsULz5s2TJA0cOFD/+te/9Morr8jPz8857ptvvlG/fv3UoEEDTZo0SadOndK///1vdenSRT/88IOio6NLvE1c5AxQzk6cOGEkmbi4uALz/vrrL3PkyBHnV3p6unNer169TOvWrU1GRoZzWm5uruncubNp3Lixc9pDDz1kJJlvv/3WOS01NdXUr1/fREdHm5ycnALbXbRokZFkVqxYUWjN9erVM/3793eZJqnQrzlz5jjHTJw4sdAx9erVO1ebTPfu3Qtdtnnz5ub33393GRsXF2f8/PzMrl27nNMOHDhgKleubK666qoC9Zxtzpw5RpJJSkpyec6SzKpVq5zTDh8+bPz9/c0//vEP57T8fq9fv95lXJUqVVzW+fHHHxtJZsOGDed87mc7cz/I9/e//90EBQW57A/5PXvnnXec0zIzM02tWrXMLbfc4pw2Y8YMI8l88MEHzmknT540jRo1KnY/yJffr2+++cYcOXLE7Nu3zyxevNhUr17d+Pv7m3379jnHlnS/zV9nhw4dTFZWlnP6c889ZySZTz/91Dkt/2fz1VdfudRV0n0/KSmpwL5a0jqfeOIJI8l89NFHBfqSm5trjDFmw4YNBdafP79x48amb9++zrHG5P1869evb2JjY53T4uLiTEBAgNmzZ49z2s8//2y8vb0L3YcL88ILL5jAwECTkpJijDHm119/NZLMxx9/7DKuXbt2pkaNGubPP/90Ttu8ebPx8vIygwcPdk7L/xm5sw/j4sDpMJS7lJQUSSpwobCUdzqjevXqzq9XX31VUt6RlOXLl+u2225Tamqqjh49qqNHj+rPP/9U37599dtvvznvylmyZIk6duyorl27OtcbHBys4cOHa/fu3fr555/L7LnceOONSkhIcPnq27dvgXEffvihy5iiDsmfLTo62rnMl19+qRkzZujEiRPq16+fjhw5IinvaNHXX3+tuLg4NWjQwLlsRESEBg0apNWrVzt7XlotWrRw/iUvSdWrV1fTpk31+++/O6ctWbJEV155pfMv+Pxx8fHxLusKDQ2VlHcUMDs7u1R1nHlUKf/n361bN6Wnp+uXX35xGRscHOxyHZWfn586duxYoOaIiAgNGDDAOS0oKEjDhw8vVV29e/dW9erVFRUVpQEDBqhSpUr67LPPVKdOHUml22/zDR8+3OWC5fvvv18+Pj5asmSJy7j69esX2Nfc3fdLU+eHH36otm3bFnpU6VwXHG/atEm//fabBg0apD///NO5nZMnT6pXr15atWqVcnNzlZOTo6VLlyouLk5169Z1Lt+8efNCf7+KMn/+fPXv31+VK1eWJDVu3FgdOnRw+f1LTk7Wpk2bdNddd6lq1arO6W3atFFsbGyBvuPSZvXpsFWrVun555/Xxo0blZycrI8//lhxcXGlWocxRtOnT9cbb7yhPXv2KDw8XA888IAee+yx8in6IpT/H1JaWlqBebNmzVJqaqoOHTrk8kK2c+dOGWM0YcIETZgwodD1Hj58WJGRkdqzZ0+hp3/y7/ras2ePWrVqVRZPRXXq1FHv3r3POe6qq65y68LoSpUquaz/mmuuUdeuXXX55Zdr2rRpmj59uo4cOaL09HQ1bdq0wPLNmzdXbm6u9u3b5zx9VRpnvgDlCwsLc7lGpah+n11P9+7ddcstt2jy5Mn617/+pR49eiguLk6DBg0qcPrwbNu2bdPjjz+u5cuXFwh0J06ccHlcp06dAi/GYWFhLm+5sGfPHjVq1KjAuMJ6WJxXX31VTZo00YkTJzR79mytWrXK5bmUZr/N17hxY5f5wcHBioiIKPB+S/Xr1y+wLnf3/dLUuWvXLt1yyy2FjjmX3377TZI0ZMiQIsecOHFCmZmZOnXqVIFeSHk/o5IEk+3bt+vHH3/U4MGDtXPnTuf0Hj166NVXX1VKSopCQkK0Z88e53rP1rx5cy1dulQnT55UpUqVzrlNXPysDkEnT55U27ZtNWzYMN18881urWPUqFH6+uuv9cILL6h169Y6duxYgetBbFelShVFRERo69atBebl/wd+9n/4+RfZPvzww0X+JdioUaOyLfQC1aFDB1WpUkWrVq0q9bJF/aV+9gXG+by9vQudboxxa9uLFy/WunXr9Pnnn2vp0qUaNmyYpk+frnXr1hV6ZFDKuwi3e/fuCgkJ0ZNPPqmGDRsqICBAP/zwg8aOHVvgAuyyrPlcOnbs6Lw7LC4uTl27dtWgQYO0Y8cOBQcHl+t+W9I7wUqion6/8rfz/PPPF3nrfHBw8HlfoC5J7733niRp9OjRGj16dIH5H374oYYOHXre28GlxeoQ1K9fP+eFk4XJzMzUY489pvfff1/Hjx9Xq1at9Oyzzzpvqd6+fbtmzpyprVu3Ov+qKOyvNUj9+/fXW2+9pe+++87lNEpR8k/z+Pr6nvPIS7169bRjx44C0/NPm5TH+w5VtJycHOeRtOrVqysoKKjI5+zl5eW8yy7/gt3jx487T09Jcv417I569eo5/8I/U2H1SNKVV16pK6+8UlOmTNGCBQsUHx+vhQsXurxv05kSExP1559/6qOPPtJVV13lnJ6UlHReNW/dulXGmAIXG7vL29tbU6dOVc+ePfXKK69o3Lhxpdpv8/3222/q2bOn83FaWpqSk5N17bXXnnNZd/f90tTZsGHDQv+AOVNRYbthw4aSpJCQkGK3U716dQUGBpZqvzqTMUYLFixQz5499cADDxSY/9RTT2n+/PkaOnSosydF9S08PJyjQBbhmqBiPPjgg1q7dq0WLlyon376SbfeequuueYa5y/q559/rgYNGuiLL75Q/fr1FR0drXvuuYcjQYV49NFHFRQUpGHDhunQoUMF5p/9V3uNGjXUo0cPzZo1S8nJyQXG518fI0nXXnutvvvuO61du9Y57eTJk3rjjTcUHR2tFi1alOEzqXgrVqxQWlqa2rZtKynvxbdPnz769NNPXY6gHTp0SAsWLFDXrl0VEhIi6f9ehM48inTy5Enn3TPuuPbaa7Vu3Tp99913zmlHjhwpcN3TX3/9VeDnmn80oLi//POP7Jy5bFZWll577bXzqvnAgQNavHixc1p6erreeOMNt9cp5Z1q6dixo2bMmKGMjIxS7bf53njjDZdrpmbOnKnTp08X+wdaPnf3/dLUecstt2jz5s36+OOPC4zL/xnlh4bjx4+7zO/QoYMaNmyoF154odDT4fnb8fb2Vt++ffXJJ59o7969zvnbt2/X0qVLi3r6TmvWrNHu3bs1dOhQDRgwoMDX7bffrhUrVujAgQOKiIhQu3btNG/ePJd6t27dqq+//rpE4ROXDquPBBVn7969mjNnjvbu3avatWtLyjt0/NVXX2nOnDl65pln9Pvvv2vPnj1atGiR3nnnHeXk5Gj06NEaMGDAed96e6lp3LixFixYoIEDB6pp06bOd4w2xigpKUkLFiyQl5eX8wJTKe/6i65du6p169a699571aBBAx06dEhr167VH3/8oc2bN0uSxo0bp/fff1/9+vXTyJEjVbVqVc2bN09JSUn68MMPy+0NGMvDiRMnnIf1T58+7bxlOjAwUOPGjXOOe/rpp5WQkKCuXbvqgQcekI+Pj2bNmqXMzEyX98jp06eP6tatq7vvvluPPPKIvL29NXv2bFWvXt3lxaY0Hn30Ub377ru65pprNGrUKOct8vXq1XO5DmfevHl67bXXdNNNN6lhw4ZKTU3Vm2++qZCQkGJfaDp37qywsDANGTJEI0eOlMPh0Lvvvntep7fuvfdevfLKKxo8eLA2btyoiIgIvfvuu2XyDtyPPPKIbr31Vs2dO1f33XdfiffbfFlZWerVq5duu+027dixQ6+99pq6du2qG2644ZzbPp99v6R1PvLII1q8eLFuvfVWDRs2TB06dNCxY8f02Wef6fXXX1fbtm3VsGFDhYaG6vXXX1flypVVqVIlderUSfXr19dbb72lfv36qWXLlho6dKgiIyO1f/9+rVixQiEhIfr8888lSZMnT9ZXX32lbt266YEHHtDp06f173//Wy1btjznR+rMnz9f3t7e6t+/f6Hzb7jhBj322GNauHChxowZo+eff179+vVTTEyM7r77buct8lWqVHF+dA4s4Ylb0i5EOus2yi+++MJIMpUqVXL58vHxMbfddpsxxph7773XSDI7duxwLrdx40Yjyfzyyy8V/RQuCjt37jT333+/adSokQkICDCBgYGmWbNm5r777jObNm0qMH7Xrl1m8ODBplatWsbX19dERkaa6667zixevLjAuAEDBpjQ0FATEBBgOnbsaL744osi63D3FvkRI0YU+/zyb0k/cuRIseMKc/Yt8g6Hw1StWtXccMMNZuPGjQXG//DDD6Zv374mODjYBAUFmZ49e5r//e9/BcZt3LjRdOrUyfj5+Zm6deuaF198schb5M9+zvl1de/e3WXaTz/9ZLp3724CAgJMZGSkeeqpp8zbb7/tss4ffvjBDBw40NStW9f4+/ubGjVqmOuuu858//335+zFmjVrzJVXXmkCAwNN7dq1zaOPPmqWLl1a4GfWvXt307JlywLLDxkypMDbEuzZs8fccMMNJigoyISHh5tRo0aZr776qlS3yBd2q3ROTo5p2LChadiwoTl9+rQxpmT7bf46V65caYYPH27CwsJMcHCwiY+Pd7l125iifzb52zrXvl/YLfIlrdMYY/7880/z4IMPmsjISOPn52fq1KljhgwZYo4ePeoc8+mnn5oWLVoYHx+fAtv68ccfzc0332yqVatm/P39Tb169cxtt91mli1b5rKdlStXmg4dOhg/Pz/ToEED8/rrrxf5Ng/5srKyTLVq1Uy3bt2KHGOMMfXr1zft27d3Pv7mm29Mly5dTGBgoAkJCTHXX3+9+fnnn12W4Rb5S5/DmHK4evAi5HA4XO4O+89//qP4+Hht27atwIWXwcHBqlWrliZOnKhnnnnG5VD2qVOnFBQUpK+//pqPOABQpLlz52ro0KHasGGD82Lr8rJr1y41atRI7777bok+9gWwBafDitC+fXvl5OTo8OHDLu+bcqYuXbro9OnT2rVrl/Pai19//VXSpXExLoBLQ/51P+68bQNwKbM6BKWlpbm8n0RSUpI2bdqkqlWrqkmTJoqPj9fgwYM1ffp0tW/fXkeOHNGyZcvUpk0b9e/fX71799Zll12mYcOGacaMGcrNzdWIESMUGxurJk2aePCZAUCe2bNna/bs2c7PFQPwfy6eK0bLwffff6/27durffv2kqQxY8aoffv2euKJJyRJc+bM0eDBg/WPf/xDTZs2VVxcnDZs2OB8QzkvLy99/vnnCg8P11VXXaX+/furefPmWrhwoceeEwCcafjw4Tp27JgWLVrk8jYJACSuCQIAAFay+kgQAACwFyEIAABYyboLo3Nzc3XgwAFVrlz5nJ+ADAAAKpYxRqmpqapdu3a5v9mtdSHowIEDzs9VAgAAF6Z9+/a5fIpAebAuBFWuXFlSXnPzP1+pMNnZ2fr666/Vp08f+fr6VlR5Fz365j565x765j565x765r6S9C4lJUVRUVHO1+vyZF0Iyj8FFhIScs4QFBQUpJCQEHbyUqBv7qN37qFv7qN37qFv7itN7yrikhUujAYAAFYiBAEAACsRggAAgJWsuyYIAACUH2OMTp8+rZycnALzTp8+LS8vL10oH1ZBCAIAAGUiKytLycnJSk9PL3S+MUYRERHav3+/IiMj5efnV8EVuiIEAQCA85abm6ukpCR5e3urdu3a8vPzK3CHV05Ojk6cOKGTJ08qKSlJjRs3Lvc3RCwOIQgAAJy3rKws5ebmKioqSkFBQYWOyc3NVXZ2tkJCQrRv3z5lZWUpICCggiv9P1wYDQAAykxJjux48ujPmS6MKgAAACoYIQgAAFiJEAQAAKxECAIAAFYiBAEAgDJTkjdCvFDeLJEQBAAAzlv+p8IX9UaJZ8ofc65Pki9vvE8QAAA4b97e3goNDdXhw4clSUFBQYW+WWJqaqpSU1MVFhYmb29vT5TqRAgCAABlolatWpLkDEJnM8bo5MmTioiIcI71JEIQAAAoEw6HQxEREapRo4ays7MLzD99+rSWL1+udu3aFThK5AmEIAAAUKa8vb0LPdWVnZ19wVwULXFhNAAAsBQhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAAr8SnyZWjRIumJJ6TUVE9X4jkBAdL06VLz5lJGhqerubjQO/fQN/fRO/fQt/9TubL01FPSgAGersQ9hKAy9MQT0i+/eLoKzwoMzPt+4IB06pRna7nY0Dv30Df30Tv30DdXEyYQgqD/OwLk5SVFRHi2Fk8JCMj7Xrs2fyGVFr1zD31zH71zD33Lk5ws5eZe3Gc/CEHlICJC+uMPT1fhGdnZ0pIl0vbtkq+vp6u5uNA799A399E799C3PHXqSPv3e7qK88OF0QAuGFk5WTqecdzTZQCwBCEIwAXhx+QfVWVaFYU9G6b7v7jf0+UAsAAhCMAF4YNtHyjjdN4FFq9vfF2nc097uCIAlzpCEIALwj+7/VPf3PmNJKl+aH15O7w9XBGASx0XRgO4IFT2r6y2tdrq+djnNaDFADkcDk+XBOASRwgCcMEIDwrXw50f9nQZACzB6TAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArHTBhKBp06bJ4XDooYceKtH4hQsXyuFwKC4urlzrAgAAl6YLIgRt2LBBs2bNUps2bUo0fvfu3Xr44YfVrVu3cq4MAABcqjwegtLS0hQfH68333xTYWFh5xyfk5Oj+Ph4TZ48WQ0aNKiACgEAwKXIx9MFjBgxQv3791fv3r319NNPn3P8k08+qRo1aujuu+/Wt99+e87xmZmZyszMdD5OSUmRJGVnZys7O7vI5fLnFTfmbAEBUmBg3vdSLHZJcadvyEPv3EPf3Efv3EPf8rjzmleS3lVkXx3GGFNhWzvLwoULNWXKFG3YsEEBAQHq0aOH2rVrpxkzZhQ6fvXq1brjjju0adMmhYeH66677tLx48f1ySefFLmNSZMmafLkyQWmL1iwQEFBQWX0TAAAQFlIT0/XoEGDdOLECYWEhJTrtjx2JGjfvn0aNWqUEhISFBAQcM7xqampuvPOO/Xmm28qPDy8xNsZP368xowZ43yckpKiqKgo9enTp9jmZmdnKyEhQbGxsfL19S3Rtpo3lw4ckGrXlrZvL3GJlxR3+oY89M499M199M499C2PO695Jeld/hmbiuCxELRx40YdPnxYl112mXNaTk6OVq1apVdeeUWZmZny9vZ2ztu1a5d2796t66+/3jktNzdXkuTj46MdO3aoYcOGBbbj7+8vf3//AtN9fX1LtPOWdJwkZWRIp07lfbf490JS6foGV/TOPfTNffTOPbb37Xxe84rrXUX21GMhqFevXtqyZYvLtKFDh6pZs2YaO3asSwCSpGbNmhUY//jjjys1NVUvvfSSoqKiyr1mAABw6fBYCKpcubJatWrlMq1SpUqqVq2ac/rgwYMVGRmpqVOnKiAgoMD40NBQSSowHQAA4Fw8fndYcfbu3SsvL4/fxQ8AAC5BF1QISkxMLPbx2ebOnVtutQAAgEsbh1kAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKF0wImjZtmhwOhx566KEix7z55pvq1q2bwsLCFBYWpt69e+u7776ruCIBAMAl44IIQRs2bNCsWbPUpk2bYsclJiZq4MCBWrFihdauXauoqCj16dNH+/fvr6BKAQDApcLjISgtLU3x8fF68803FRYWVuzY+fPn64EHHlC7du3UrFkzvfXWW8rNzdWyZcsqqFoAAHCp8PF0ASNGjFD//v3Vu3dvPf3006VaNj09XdnZ2apatWqRYzIzM5WZmel8nJKSIknKzs5WdnZ2kcvlzytuzNkCAqTAwLzvpVjskuJO35CH3rmHvrmP3rmHvuVx5zWvJL2ryL46jDGmwrZ2loULF2rKlCnasGGDAgIC1KNHD7Vr104zZswo0fIPPPCAli5dqm3btikgIKDQMZMmTdLkyZMLTF+wYIGCgoLOp3wAAFDG0tPTNWjQIJ04cUIhISHlui2PHQnat2+fRo0apYSEhCIDTHGmTZumhQsXKjExsdjlx48frzFjxjgfp6SkOK8lKq652dnZSkhIUGxsrHx9fUtUU/Pm0oEDUu3a0vbtJX8ulxJ3+oY89M499M199M499C2PO695Jeld/hmbiuCxELRx40YdPnxYl112mXNaTk6OVq1apVdeeUWZmZny9vYudNkXXnhB06ZN0zfffHPOi6n9/f3l7+9fYLqvr2+Jdt6SjpOkjAzp1Km87xb/XkgqXd/git65h765j965x/a+nc9rXnG9q8ieeiwE9erVS1u2bHGZNnToUDVr1kxjx44tMgA999xzmjJlipYuXarLL7+8IkoFAACXII+FoMqVK6tVq1Yu0ypVqqRq1ao5pw8ePFiRkZGaOnWqJOnZZ5/VE088oQULFig6OloHDx6UJAUHBys4OLhinwAAALioefwW+eLs3btXycnJzsczZ85UVlaWBgwYoIiICOfXCy+84MEqAQDAxcjjt8ifKTExsdjHu3fvrrBaAADApe2CPhIEAABQXghBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACu5HYJ27dqlxx9/XAMHDtThw4clSV9++aW2bdtWZsUBAACUF7dC0MqVK9W6dWutX79eH330kdLS0iRJmzdv1sSJE8u0QAAAgPLgVggaN26cnn76aSUkJMjPz885/eqrr9a6devKrDgAAIDy4lYI2rJli2666aYC02vUqKGjR4+ed1EAAADlza0QFBoaquTk5ALTf/zxR0VGRp53UQAAAOXNrRB0xx13aOzYsTp48KAcDodyc3O1Zs0aPfzwwxo8eHBZ1wgAAFDm3ApBzzzzjJo1a6aoqCilpaWpRYsWuuqqq9S5c2c9/vjjZV0jAABAmfNxZyE/Pz+9+eabmjBhgrZu3aq0tDS1b99ejRs3Luv6AAAAyoVbIShf3bp1Vbdu3bKqBQAAoMKUOASNGTOmxCt98cUX3SoGAACgopQ4BP34448uj3/44QedPn1aTZs2lST9+uuv8vb2VocOHdwqZNq0aRo/frxGjRqlGTNmFDlu0aJFmjBhgnbv3q3GjRvr2Wef1bXXXuvWNgEAgL1KHIJWrFjh/PeLL76oypUra968eQoLC5Mk/fXXXxo6dKi6detW6iI2bNigWbNmqU2bNsWO+9///qeBAwdq6tSpuu6667RgwQLFxcXphx9+UKtWrUq9XQAAYC+37g6bPn26pk6d6gxAkhQWFqann35a06dPL9W60tLSFB8frzfffNNlfYV56aWXdM011+iRRx5R8+bN9dRTT+myyy7TK6+84s7TAAAAFnPrwuiUlBQdOXKkwPQjR44oNTW1VOsaMWKE+vfvr969e+vpp58uduzatWsLXJvUt29fffLJJ0Uuk5mZqczMTJfaJSk7O1vZ2dlFLpc/r7gxZwsIkAID876XYrFLijt9Qx565x765j565x76lsed17yS9K4i++pWCLrppps0dOhQTZ8+XR07dpQkrV+/Xo888ohuvvnmEq9n4cKF+uGHH7Rhw4YSjT948KBq1qzpMq1mzZo6ePBgkctMnTpVkydPLjD966+/VlBQ0Dm3mZCQUKLaJOnMg2BLlpR4sUtSafoGV/TOPfTNffTOPbb37Xxe84rrXXp6upsVlZ5bIej111/Xww8/rEGDBjkTm4+Pj+6++249//zzJVrHvn37NGrUKCUkJCggIMCdMkpk/PjxLkePUlJSFBUVpT59+igkJKTI5bKzs5WQkKDY2Fj5+vqWaFvNm0sHDki1a0vbt5936Rcld/qGPPTOPfTNffTOPfQtjzuveSXpXf4Zm4rgVggKCgrSa6+9pueff167du2SJDVs2FCVKlUq8To2btyow4cP67LLLnNOy8nJ0apVq/TKK68oMzNT3t7eLsvUqlVLhw4dcpl26NAh1apVq8jt+Pv7y9/fv8B0X1/fEu28JR0nSRkZ0qlTed8t/r2QVLq+wRW9cw99cx+9c4/tfTuf17zieleRPT2vN0usVKnSOe/oKkqvXr20ZcsWl2lDhw5Vs2bNNHbs2AIBSJJiYmK0bNkyPfTQQ85pCQkJiomJcasGAABgL7dCUM+ePeVwOIqcv3z58nOuo3LlygVua69UqZKqVavmnD548GBFRkZq6tSpkqRRo0ape/fumj59uvr376+FCxfq+++/1xtvvOHO0wAAABZzKwS1a9fO5XF2drY2bdqkrVu3asiQIWVRlyRp79698vL6v7v4O3furAULFujxxx/XP//5TzVu3FiffPIJ7xEEAABKza0Q9K9//avQ6ZMmTVJaWprbxSQmJhb7WJJuvfVW3XrrrW5vAwAAQHLzzRKL8re//U2zZ88uy1UCAACUizINQWvXri3X290BAADKilunw85+Q0RjjJKTk/X9999rwoQJZVIYAABAeXIrBIWEhLjcHebl5aWmTZvqySefVJ8+fcqsOAAAgPLiVgiaO3duGZcBAABQsdy6JqhBgwb6888/C0w/fvy4GjRocN5FAQAAlDe3QtDu3buVk5NTYHpmZqb2799/3kUBAACUt1KdDvvss8+c/166dKmqVKnifJyTk6Nly5YpOjq6zIoDAAAoL6UKQXFxcZIkh8NR4J2hfX19FR0drenTp5dZcQAAAOWlVCEoNzdXklS/fn1t2LBB4eHh5VIUAABAeXPr7rCkpKSyrgMAAKBClTgEvfzyyxo+fLgCAgL08ssvFzt25MiR510YAABAeSpxCPrXv/6l+Ph4BQQEFPkBqlLe9UKEIAAAcKErcQg68xQYp8MAAMDFzq33CXryySeVnp5eYPqpU6f05JNPnndRAAAA5c2tEDR58mSlpaUVmJ6enq7Jkyefd1EAAADlza0QZIxx+QDVfJs3b1bVqlXPuygAAIDyVqpb5MPCwuRwOORwONSkSROXIJSTk6O0tDTdd999ZV4kAABAWStVCJoxY4aMMRo2bJgmT57s8rEZfn5+io6OVkxMTJkXCQAAUNZKFYLyPyqjfv366ty5s3x9fculKAAAgPLm1jtGd+/e3fnvjIwMZWVlucwPCQk5v6oAAADKmVsXRqenp+vBBx9UjRo1VKlSJYWFhbl8AQAAXOjcCkGPPPKIli9frpkzZ8rf319vvfWWJk+erNq1a+udd94p6xoBAADKnFunwz7//HO988476tGjh4YOHapu3bqpUaNGqlevnubPn6/4+PiyrhMAAKBMuXUk6NixY2rQoIGkvOt/jh07Jknq2rWrVq1aVXbVAQAAlBO3QlCDBg2cnx/WrFkzffDBB5LyjhCdeds8AADAhcqtEDR06FBt3rxZkjRu3Di9+uqrCggI0OjRo/Xoo4+WaYEAAADlwa1rgkaPHu38d+/evfXLL79o48aNCg8P13vvvVdmxQEAAJQXt44Ena1evXq6+eabVaVKFb399ttlsUoAAIByVSYhCAAA4GJDCAIAAFYiBAEAACuV6sLom2++udj5x48fP59aAAAAKkypQtC53gOoSpUqGjx48HkVBAAAUBFKFYLmzJlTXnUAAABUKK4JAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFjJoyFo5syZatOmjUJCQhQSEqKYmBh9+eWXxS4zY8YMNW3aVIGBgYqKitLo0aOVkZFRQRUDAIBLhY8nN16nTh1NmzZNjRs3ljFG8+bN04033qgff/xRLVu2LDB+wYIFGjdunGbPnq3OnTvr119/1V133SWHw6EXX3zRA88AAABcrDwagq6//nqXx1OmTNHMmTO1bt26QkPQ//73P3Xp0kWDBg2SJEVHR2vgwIFav359hdQLAAAuHR4NQWfKycnRokWLdPLkScXExBQ6pnPnznrvvff03XffqWPHjvr999+1ZMkS3XnnnUWuNzMzU5mZmc7HKSkpkqTs7GxlZ2cXuVz+vOLGnC0gQAoMzPteisUuKe70DXnonXvom/vonXvoWx53XvNK0ruK7KvDGGMqbGuF2LJli2JiYpSRkaHg4GAtWLBA1157bZHjX375ZT388MMyxuj06dO67777NHPmzCLHT5o0SZMnTy4wfcGCBQoKCiqT5wAAAMpGenq6Bg0apBMnTigkJKRct+XxEJSVlaW9e/fqxIkTWrx4sd566y2tXLlSLVq0KDA2MTFRd9xxh55++ml16tRJO3fu1KhRo3TvvfdqwoQJha6/sCNBUVFROnr0aLHNzc7OVkJCgmJjY+Xr61ui59K8uXTggFS7trR9e4kWueS40zfkoXfuoW/uo3fuoW953HnNK0nvUlJSFB4eXiEhyOOnw/z8/NSoUSNJUocOHbRhwwa99NJLmjVrVoGxEyZM0J133ql77rlHktS6dWudPHlSw4cP12OPPSYvr4I3u/n7+8vf37/AdF9f3xLtvCUdJ0kZGdKpU3nfLf69kFS6vsEVvXMPfXMfvXOP7X07n9e84npXkT294N4nKDc31+XIzZnS09MLBB1vb29JkocPaAEAgIuMR48EjR8/Xv369VPdunWVmpqqBQsWKDExUUuXLpUkDR48WJGRkZo6daqkvLvJXnzxRbVv3955OmzChAm6/vrrnWEIAACgJDwagg4fPqzBgwcrOTlZVapUUZs2bbR06VLFxsZKkvbu3ety5Ofxxx+Xw+HQ448/rv3796t69eq6/vrrNWXKFE89BQAAcJHyaAh6++23i52fmJjo8tjHx0cTJ07UxIkTy7EqAABggwvumiAAAICKQAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABW8mgImjlzptq0aaOQkBCFhIQoJiZGX375ZbHLHD9+XCNGjFBERIT8/f3VpEkTLVmypIIqBgAAlwofT268Tp06mjZtmho3bixjjObNm6cbb7xRP/74o1q2bFlgfFZWlmJjY1WjRg0tXrxYkZGR2rNnj0JDQyu+eAAAcFHzaAi6/vrrXR5PmTJFM2fO1Lp16woNQbNnz9axY8f0v//9T76+vpKk6OjoiigVAABcYjwags6Uk5OjRYsW6eTJk4qJiSl0zGeffaaYmBiNGDFCn376qapXr65BgwZp7Nix8vb2LnSZzMxMZWZmOh+npKRIkrKzs5WdnV1kPfnzihtztoAAKTAw73spFrukuNM35KF37qFv7qN37qFvedx5zStJ7yqyrw5jjKmwrRViy5YtiomJUUZGhoKDg7VgwQJde+21hY5t1qyZdu/erfj4eD3wwAPauXOnHnjgAY0cOVITJ04sdJlJkyZp8uTJBaYvWLBAQUFBZfpcAADA+UlPT9egQYN04sQJhYSElOu2PB6CsrKytHfvXp04cUKLFy/WW2+9pZUrV6pFixYFxjZp0kQZGRlKSkpyHvl58cUX9fzzzys5ObnQ9Rd2JCgqKkpHjx4ttrnZ2dlKSEhQbGys89TbuTRvLh04INWuLW3fXqJFLjnu9A156J176Jv76J176Fsed17zStK7lJQUhYeHV0gI8vjpMD8/PzVq1EiS1KFDB23YsEEvvfSSZs2aVWBsRESEfH19XU59NW/eXAcPHlRWVpb8/PwKLOPv7y9/f/8C0319fUu085Z0nCRlZEinTuV9t/j3QlLp+gZX9M499M199M49tvftfF7ziutdRfb0gnufoNzcXJcjN2fq0qWLdu7cqdzcXOe0X3/9VREREYUGIAAAgKJ4NASNHz9eq1at0u7du7VlyxaNHz9eiYmJio+PlyQNHjxY48ePd46///77dezYMY0aNUq//vqr/vvf/+qZZ57RiBEjPPUUAADARcqjp8MOHz6swYMHKzk5WVWqVFGbNm20dOlSxcbGSpL27t0rL6//y2lRUVFaunSpRo8erTZt2igyMlKjRo3S2LFjPfUUAADARcqjIejtt98udn5iYmKBaTExMVq3bl05VQQAAGxxwV0TBAAAUBEIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABW8vF0AZei5GSpTh1PV+EZAQHS9OlS8+ZSRoanq7m40Dv30Df30Tv30Lc8+/d7uoLzRwgqQ5Ur533Pzb00dg53BAbmfT9wQDp1yrO1XGzonXvom/vonXvom6sxYzxdgfsIQWXoqaekCROk1FRPV+I5AQF532vXtvsvJHfQO/fQN/fRO/fQt/8zZgwhCP+/AQPyvmyWnS0tWSJt3y75+nq6mosLvXMPfXMfvXMPfbt0cGE0AACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKhCAAAGAlQhAAALASIQgAAFiJEAQAAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgAABgJUIQAACwEiEIAABYiRAEAACsRAgCAABWIgQBAAArEYIAAICVCEEAAMBKPp4uoKIZYyRJKSkpxY7Lzs5Wenq6UlJS5OvrWxGlXRLom/vonXvom/vonXvom/tK0rv81+f81+vyZF0ISk1NlSRFRUV5uBIAAFCU1NRUValSpVy34TAVEbUuILm5uTpw4IAqV64sh8NR5LiUlBRFRUVp3759CgkJqcAKL270zX30zj30zX30zj30zX0l6Z0xRqmpqapdu7a8vMr3qh3rjgR5eXmpTp06JR4fEhLCTu4G+uY+euce+uY+euce+ua+c/WuvI8A5ePCaAAAYCVCEAAAsBIhqAj+/v6aOHGi/P39PV3KRYW+uY/euYe+uY/euYe+ue9C6511F0YDAABIHAkCAACWIgQBAAArEYIAAICVCEEAAMBKF2UImjp1qq644gpVrlxZNWrUUFxcnHbs2OEyJiMjQyNGjFC1atUUHBysW265RYcOHXIZs3fvXvXv319BQUGqUaOGHnnkEZ0+fdplTGJioi677DL5+/urUaNGmjt3boF6Xn31VUVHRysgIECdOnXSd999V+paKsrMmTPVpk0b5xtVxcTE6MsvvyxVrTb27WzTpk2Tw+HQQw895JxG7wo3adIkORwOl69mzZqVqlYb+yZJ+/fv19/+9jdVq1ZNgYGBat26tb7//nvnfGOMnnjiCUVERCgwMFC9e/fWb7/95rKOY8eOKT4+XiEhIQoNDdXdd9+ttLQ0lzE//fSTunXrpoCAAEVFRem5554rUMuiRYvUrFkzBQQEqHXr1lqyZInL/JLUUlGio6ML7HMOh0MjRoyQxD5XlJycHE2YMEH169dXYGCgGjZsqKeeesrlM7wuuX3OXIT69u1r5syZY7Zu3Wo2bdpkrr32WlO3bl2TlpbmHHPfffeZqKgos2zZMvP999+bK6+80nTu3Nk5//Tp06ZVq1amd+/e5scffzRLliwx4eHhZvz48c4xv//+uwkKCjJjxowxP//8s/n3v/9tvL29zVdffeUcs3DhQuPn52dmz55ttm3bZu69914TGhpqDh06VOJaKtJnn31m/vvf/5pff/3V7Nixw/zzn/80vr6+ZuvWrSWq1da+nem7774z0dHRpk2bNmbUqFHO6fSucBMnTjQtW7Y0ycnJzq8jR46UuFZb+3bs2DFTr149c9ddd5n169eb33//3SxdutTs3LnTOWbatGmmSpUq5pNPPjGbN282N9xwg6lfv745deqUc8w111xj2rZta9atW2e+/fZb06hRIzNw4EDn/BMnTpiaNWua+Ph4s3XrVvP++++bwMBAM2vWLOeYNWvWGG9vb/Pcc8+Zn3/+2Tz++OPG19fXbNmypVS1VJTDhw+77G8JCQlGklmxYoUxhn2uKFOmTDHVqlUzX3zxhUlKSjKLFi0ywcHB5qWXXnKOudT2uYsyBJ3t8OHDRpJZuXKlMcaY48ePG19fX7No0SLnmO3btxtJZu3atcYYY5YsWWK8vLzMwYMHnWNmzpxpQkJCTGZmpjHGmEcffdS0bNnSZVu333676du3r/Nxx44dzYgRI5yPc3JyTO3atc3UqVNLXIunhYWFmbfeeou+lUBqaqpp3LixSUhIMN27d3eGIHpXtIkTJ5q2bdsWOo++FW3s2LGma9euRc7Pzc01tWrVMs8//7xz2vHjx42/v795//33jTHG/Pzzz0aS2bBhg3PMl19+aRwOh9m/f78xxpjXXnvNhIWFOXuZv+2mTZs6H992222mf//+Ltvv1KmT+fvf/17iWjxp1KhRpmHDhiY3N5d9rhj9+/c3w4YNc5l28803m/j4eGPMpbnPXZSnw8524sQJSVLVqlUlSRs3blR2drZ69+7tHNOsWTPVrVtXa9eulSStXbtWrVu3Vs2aNZ1j+vbtq5SUFG3bts055sx15I/JX0dWVpY2btzoMsbLy0u9e/d2jilJLZ6Sk5OjhQsX6uTJk4qJiaFvJTBixAj179+/wPOjd8X77bffVLt2bTVo0EDx8fHau3dviWu1tW+fffaZLr/8ct16662qUaOG2rdvrzfffNM5PykpSQcPHnSpt0qVKurUqZNL70JDQ3X55Zc7x/Tu3VteXl5av369c8xVV10lPz8/55i+fftqx44d+uuvv5xjiutvSWrxlKysLL333nsaNmyYHA4H+1wxOnfurGXLlunXX3+VJG3evFmrV69Wv379JF2a+9xFH4Jyc3P10EMPqUuXLmrVqpUk6eDBg/Lz81NoaKjL2Jo1a+rgwYPOMWfu4Pnz8+cVNyYlJUWnTp3S0aNHlZOTU+iYM9dxrloq2pYtWxQcHCx/f3/dd999+vjjj9WiRQv6dg4LFy7UDz/8oKlTpxaYR++K1qlTJ82dO1dfffWVZs6cqaSkJHXr1k2pqan0rRi///67Zs6cqcaNG2vp0qW6//77NXLkSM2bN89Zb359RdV78OBB1ahRw2W+j4+PqlatWib9PXP+uWrxlE8++UTHjx/XXXfdJYnf1eKMGzdOd9xxh5o1ayZfX1+1b99eDz30kOLj45315tdXVL0X2z530X+K/IgRI7R161atXr3a06VcNJo2bapNmzbpxIkTWrx4sYYMGaKVK1d6uqwL2r59+zRq1CglJCQoICDA0+VcVPL/ipSkNm3aqFOnTqpXr54++OADBQYGerCyC1tubq4uv/xyPfPMM5Kk9u3ba+vWrXr99dc1ZMgQD1d38Xj77bfVr18/1a5d29OlXPA++OADzZ8/XwsWLFDLli21adMmPfTQQ6pdu/Ylu89d1EeCHnzwQX3xxRdasWKF6tSp45xeq1YtZWVl6fjx4y7jDx06pFq1ajnHnH0Ffv7jc40JCQlRYGCgwsPD5e3tXeiYM9dxrloqmp+fnxo1aqQOHTpo6tSpatu2rV566SX6VoyNGzfq8OHDuuyyy+Tj4yMfHx+tXLlSL7/8snx8fFSzZk16V0KhoaFq0qSJdu7cyT5XjIiICLVo0cJlWvPmzZ2nEvNrOtdzOnz4sMv806dP69ixY2XS3zPnn6sWT9izZ4+++eYb3XPPPc5p7HNFe+SRR5xHg1q3bq0777xTo0ePdh79vhT3uYsyBBlj9OCDD+rjjz/W8uXLVb9+fZf5HTp0kK+vr5YtW+actmPHDu3du1cxMTGSpJiYGG3ZssXlh5WQkKCQkBDnfzwxMTEu68gfk78OPz8/dejQwWVMbm6uli1b5hxTklo8LTc3V5mZmfStGL169dKWLVu0adMm59fll1+u+Ph457/pXcmkpaVp165dioiIYJ8rRpcuXQq89cevv/6qevXqSZLq16+vWrVqudSbkpKi9evXu/Tu+PHj2rhxo3PM8uXLlZubq06dOjnHrFq1StnZ2c4xCQkJatq0qcLCwpxjiutvSWrxhDlz5qhGjRrq37+/cxr7XNHS09Pl5eUaC7y9vZWbmyvpEt3nSnwJ9QXk/vvvN1WqVDGJiYkut0Gmp6c7x9x3332mbt26Zvny5eb77783MTExJiYmxjk//xbIPn36mE2bNpmvvvrKVK9evdBbIB955BGzfft28+qrrxZ6C6S/v7+ZO3eu+fnnn83w4cNNaGioy10F56qlIo0bN86sXLnSJCUlmZ9++smMGzfOOBwO8/XXX5eoVlv7Vpgz7w4zht4V5R//+IdJTEw0SUlJZs2aNaZ3794mPDzcHD58uES12tq37777zvj4+JgpU6aY3377zcyfP98EBQWZ9957zzlm2rRpJjQ01Hz66afmp59+MjfeeGOhtyu3b9/erF+/3qxevdo0btzY5Xbl48ePm5o1a5o777zTbN261SxcuNAEBQUVuF3Zx8fHvPDCC2b79u1m4sSJhd6ufK5aKlJOTo6pW7euGTt2bIF57HOFGzJkiImMjHTeIv/RRx+Z8PBw8+ijjzrHXGr73EUZgiQV+jVnzhznmFOnTpkHHnjAhIWFmaCgIHPTTTeZ5ORkl/Xs3r3b9OvXzwQGBprw8HDzj3/8w2RnZ7uMWbFihWnXrp3x8/MzDRo0cNlGvn//+9+mbt26xs/Pz3Ts2NGsW7fOZX5Jaqkow4YNM/Xq1TN+fn6mevXqplevXs4AVNJabexbYc4OQfSucLfffruJiIgwfn5+JjIy0tx+++0u73VD34r2+eefm1atWhl/f3/TrFkz88Ybb7jMz83NNRMmTDA1a9Y0/v7+plevXmbHjh0uY/78808zcOBAExwcbEJCQszQoUNNamqqy5jNmzebrl27Gn9/fxMZGWmmTZtWoJYPPvjANGnSxPj5+ZmWLVua//73v6WupSItXbrUSCq0Bva5wqWkpJhRo0aZunXrmoCAANOgQQPz2GOPudzKfqntcw5jzngrSAAAAEtclNcEAQAAnC9CEAAAsBIhCAAAWIkQBAAArEQIAgAAViIEAQAAKxGCAACAlQhBAADASoQgABe13bt3y+FwaNOmTeWyfofDoU8++aRc1g3AswhBAM7LXXfdpbi4OI9tPyoqSsnJyWrVqpUkKTExUQ6Ho8AncwPA2Xw8XQAAnA9vb2/VqlXL02UAuAhxJAhAuVm5cqU6duwof39/RUREaNy4cTp9+rRzfo8ePTRy5Eg9+uijqlq1qmrVqqVJkya5rOOXX35R165dFRAQoBYtWuibb75xOUV15umw3bt3q2fPnpKksLAwORwO3XXXXZKk6OhozZgxw2Xd7dq1c9neb7/9pquuusq5rYSEhALPad++fbrtttsUGhqqqlWr6sYbb9Tu3bvPt1UAPIAQBKBc7N+/X9dee62uuOIKbd68WTNnztTbb7+tp59+2mXcvHnzVKlSJa1fv17PPfecnnzySWf4yMnJUVxcnIKCgrR+/Xq98cYbeuyxx4rcZlRUlD788ENJ0o4dO5ScnKyXXnqpRPXm5ubq5ptvlp+fn9avX6/XX39dY8eOdRmTnZ2tvn37qnLlyvr222+1Zs0aBQcH65prrlFWVlZp2gPgAsDpMADl4rXXXlNUVJReeeUVORwONWvWTAcOHNDYsWP1xBNPyMsr72+wNm3aaOLEiZKkxo0b65VXXtGyZcsUGxurhIQE7dq1S4mJic5TXlOmTFFsbGyh2/T29lbVqlUlSTVq1FBoaGiJ6/3mm2/0yy+/aOnSpapdu7Yk6ZlnnlG/fv2cY/7zn/8oNzdXb731lhwOhyRpzpw5Cg0NVWJiovr06VO6JgHwKEIQgHKxfft2xcTEOMOCJHXp0kVpaWn6448/VLduXUl5IehMEREROnz4sKS8ozlRUVEu1/x07Nix3OqNiopyBiBJiomJcRmzefNm7dy5U5UrV3aZnpGRoV27dpVLXQDKDyEIgEf5+vq6PHY4HMrNzS3z7Xh5eckY4zItOzu7VOtIS0tThw4dNH/+/ALzqlevfl71Aah4hCAA5aJ58+b68MMPZYxxHg1as2aNKleurDp16pRoHU2bNtW+fft06NAh1axZU5K0YcOGYpfx8/OTlHc90ZmqV6+u5ORk5+OUlBQlJSW51Ltv3z4lJycrIiJCkrRu3TqXdVx22WX6z3/+oxo1aigkJKREzwHAhYsLowGctxMnTmjTpk0uX8OHD9e+ffv0//7f/9Mvv/yiTz/9VBMnTtSYMWOc1wOdS2xsrBo2bKghQ4bop59+0po1a/T4449LkstptjPVq1dPDodDX3zxhY4cOaK0tDRJ0tVXX613331X3377rbZs2aIhQ4bI29vbuVzv3r3VpEkTDRkyRJs3b9a3335b4CLs+Ph4hYeH68Ybb9S3336rpKQkJSYmauTIkfrjjz/caR0ADyIEAThviYmJat++vcvXU089pSVLlui7775T27Ztdd999+nuu+92hpiS8Pb21ieffKK0tDRdccUVuueee5zBJCAgoNBlIiMjNXnyZI0bN041a9bUgw8+KEkaP368unfvruuuu079+/dXXFycGjZs6FzOy8tLH3/8sU6dOqWOHTvqnnvu0ZQpU1zWHRQUpFWrVqlu3bq6+eab1bx5c919993KyMjgyBBwEXKYs0+SA8AFbM2aNeratat27tzpEmIAoLQIQQAuaB9//LGCg4PVuHFj7dy5U6NGjVJYWJhWr17t6dIAXOS4MBrABS01NVVjx47V3r17FR4ert69e2v69OmeLgvAJYAjQQAAwEpcGA0AAKxECAIAAFYiBAEAACsRggAAgJUIQQAAwEqEIAAAYCVCEAAAsBIhCAAAWOn/A2+pW3iVA+L1AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1000x1000 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_452544/634836532.py:37: UserWarning: Legend does not support handles for PatchCollection instances.\n",
      "See: https://matplotlib.org/stable/tutorials/intermediate/legend_guide.html#implementing-a-custom-legend-handler\n",
      "  plt.legend()\n",
      "/tmp/ipykernel_452544/634836532.py:37: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.\n",
      "  plt.legend()\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf8AAANXCAYAAAA7IVYUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABMHElEQVR4nO3deVzU1eL/8fewg4qgoIIiuGtuuVzNLTW3vN6Kyuyq92pqebvZT9NbabdMLUtb7Hq7ltniUmne0hbrmkYqbrmX5p4a7riniAsinN8fPpivEwg4DIx4Xs/Hg4fOZ8585sxh9MXMfGZwGGOMAACANXy8PQEAAFC0iD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD9ggenTp8vhcGjv3r3enkq+JSYmyuFwKDEx0dtTKVR79+6Vw+HQ9OnTvT2V6zZ69Gg5HA5vTwNuIP5wSkpK0uOPP66aNWsqJCREISEhuuWWWzRo0CD9/PPPhXa97dq1k8PhyPNr9OjRkqS4uDj96U9/ctnHtS5ToUIF55is/6hy+nrnnXeua44BAQGqUqWKBg4cqAMHDnh8TZB/WT/YZH35+fmpYsWKeuihh3To0CFvT8/rDh8+rNGjR2vjxo3enookqUePHnI4HBo+fHiB9pP1fV+/fr2HZmYXP29PADeGb775Rg8++KD8/PzUu3dvNWzYUD4+PtqxY4c+//xzTZ48WUlJSYqNjfX4dT/77LN6+OGHnafXrVunN998U//85z9Vp04d5/YGDRrkup9OnTqpT58+LtuCg4OzjZs8ebJKlizpsq158+Z5zrNSpUoaN26cJOnSpUvatm2b3nnnHS1cuFDbt29XSEhInvtA4XnhhRdUpUoVXbx4UatXr9b06dO1YsUKbdmyRUFBQd6e3jXFxsbqwoUL8vf3L5T9Hz58WGPGjFFcXJxuvfXWQrmO/EpJSdHXX3+tuLg4ffLJJxo/fjzPHHgJ8Yf27NmjP//5z4qNjdWiRYsUFRXlcv4rr7yit99+Wz4+hfNEUadOnVxOBwUF6c0331SnTp3Url27fO+nZs2a+stf/pLnuO7duysiIuJ6p6nSpUtn23+VKlX0+OOPa+XKldluB4pW165d1bRpU0nSww8/rIiICL3yyiuaN2+eevToUWTzMMbo4sWLOf7gmROHw3FD/3DiSXPnzlVGRoamTp2qO+64Q8uWLVPbtm29PS0r8bQ/9Oqrr+rcuXOaNm1atvBLkp+fnwYPHqyYmBiX7Tt27FD37t1VpkwZBQUFqWnTppo3b162y//666964IEHVKZMGYWEhOi2227T//73v0K7PUUp62UFPz/Xn6N/+uknde3aVaGhoSpZsqQ6dOig1atXu4y51uulOb0+n/VSx4oVK9SsWTMFBQWpatWq+vDDD7NdfuvWrbrjjjsUHBysSpUqaezYscrMzMw2bv369erSpYsiIiIUHBysKlWqqH///nne5q+++krdunVTdHS0AgMDVa1aNb344ovKyMhwGdeuXTvVq1dP27ZtU/v27RUSEqKKFSvq1VdfzbbPgwcPKj4+XiVKlFC5cuU0dOhQpaWl5TmX3LRp00bSlR9ur5af+23W92DZsmX629/+prJlyyo0NFR9+vTRb7/95jI263uzcOFCNW3aVMHBwZoyZYqk/N33r/Waf37/fZ0+fVpDhw5VXFycAgMDValSJfXp00cnTpxQYmKi/vCHP0iS+vXr53xp5OrrWrNmje68806VLl1aISEhatu2rVauXJntelasWKE//OEPCgoKUrVq1Zy38XrMnDlTnTp1Uvv27VWnTh3NnDkzx3GLFy9WmzZtVKJECYWFhemee+7R9u3br/v6cG088oe++eYbVa9ePV9PfWfZunWrWrVqpYoVK2rEiBEqUaKEPv30U8XHx2vu3Lm69957JUlHjx5Vy5Ytdf78eQ0ePFhly5bVjBkzdPfdd2vOnDnOcZ5w8eJFnThxwmVbqVKlFBgY6LLt1KlTLqd9fX0VHh6e5/4zMjKc+09PT9f27ds1atQoVa9eXa1atXKO27p1q9q0aaPQ0FA9/fTT8vf315QpU9SuXTstXbr0utb5art371b37t01YMAA9e3bV1OnTtVDDz2kJk2aqG7dupKkI0eOqH379rp8+bLz+/Luu+9mexR67Ngxde7cWZGRkRoxYoTCwsK0d+9eff7553nOY/r06SpZsqSGDRumkiVLavHixXr++eeVkpKi1157zWXsb7/9pjvvvFP33XefevTooTlz5mj48OGqX7++unbtKkm6cOGCOnTooP3792vw4MGKjo7WRx99pMWLF7u1Tlmyfni6+nub3/ttlscff1xhYWEaPXq0du7cqcmTJ2vfvn3OgxGz7Ny5Uz179tTf/vY3PfLII6pVq1aB7vv5nWdqaqratGmj7du3q3///mrcuLFOnDihefPm6eDBg6pTp45eeOEFPf/88xo4cKDzB6KWLVtKuhLZrl27qkmTJho1apR8fHw0bdo03XHHHVq+fLmaNWsmSdq8ebPz/jJ69GhdvnxZo0aNUvny5fP9/Th8+LCWLFmiGTNmSJJ69uypf/3rX5o0aZICAgKc477//nt17dpVVatW1ejRo3XhwgX95z//UatWrfTjjz8qLi4u39eJXBhY7cyZM0aSiY+Pz3beb7/9Zo4fP+78On/+vPO8Dh06mPr165uLFy86t2VmZpqWLVuaGjVqOLc98cQTRpJZvny5c9vZs2dNlSpVTFxcnMnIyMh2vZ999pmRZJYsWZLjnGNjY023bt1ctknK8WvatGnOMaNGjcpxTGxsbF7LZNq2bZvjZevUqWN+/fVXl7Hx8fEmICDA7Nmzx7nt8OHDplSpUub222/PNp/fmzZtmpFkkpKSXG6zJLNs2TLntmPHjpnAwEDzj3/8w7kta73XrFnjMq506dIu+/ziiy+MJLNu3bo8b/vvXX0/yPK3v/3NhISEuNwfstbsww8/dG5LS0szFSpUMPfff79z28SJE40k8+mnnzq3nTt3zlSvXj3X+0GWrPX6/vvvzfHjx82BAwfMnDlzTGRkpAkMDDQHDhxwjs3v/TZrn02aNDGXLl1ybn/11VeNJPPVV185t2V9bxYsWOAyr/ze95OSkrLdV/M7z+eff95IMp9//nm2dcnMzDTGGLNu3bps+886v0aNGqZLly7OscZc+f5WqVLFdOrUybktPj7eBAUFmX379jm3bdu2zfj6+uZ4H87J66+/boKDg01KSooxxphffvnFSDJffPGFy7hbb73VlCtXzpw8edK5bdOmTcbHx8f06dPHuS3re+TOfRjG8LS/5VJSUiQp2wFw0pWnbSMjI51fb731lqQrj5wXL16sHj166OzZszpx4oROnDihkydPqkuXLtq1a5fzKOv58+erWbNmat26tXO/JUuW1MCBA7V3715t27bNY7flnnvuUUJCgstXly5dso2bO3euy5hrPfX4e3Fxcc7LfPvtt5o4caLOnDmjrl276vjx45KuPDvw3XffKT4+XlWrVnVeNioqSr169dKKFSuca369brnlFucjN0mKjIxUrVq19Ouvvzq3zZ8/X7fddpvzEVvWuN69e7vsKywsTNKVZ33S09Ovax5XP4uQ9f1v06aNzp8/rx07driMLVmypMtxEgEBAWrWrFm2OUdFRal79+7ObSEhIRo4cOB1zatjx46KjIxUTEyMunfvrhIlSmjevHmqVKmSpOu732YZOHCgy4F4f//73+Xn56f58+e7jKtSpUq2+5q79/3rmefcuXPVsGHDHJ9FyOtAuo0bN2rXrl3q1auXTp486byec+fOqUOHDlq2bJkyMzOVkZGhhQsXKj4+XpUrV3Zevk6dOjn++7qWmTNnqlu3bipVqpQkqUaNGmrSpInLv7/k5GRt3LhRDz30kMqUKePc3qBBA3Xq1CnbusN9Vj/tv2zZMr322mvasGGDkpOT9cUXXyg+Pv669mGM0YQJE/Tuu+9q3759ioiI0GOPPaZnn322cCbtYVn/EFNTU7OdN2XKFJ09e1ZHjx51+Q989+7dMsZo5MiRGjlyZI77PXbsmCpWrKh9+/bl+DR31lH8+/btU7169TxxU1SpUiV17Ngxz3G33367Wwf8lShRwmX/d955p1q3bq2mTZtq/PjxmjBhgo4fP67z58+rVq1a2S5fp04dZWZm6sCBA86n6a/H1f/xZgkPD3d5Dfpa6/37+bRt21b333+/xowZo3/9619q166d4uPj1atXr2wvk/ze1q1b9dxzz2nx4sXZfpA5c+aMy+lKlSpli1B4eLjLW0f37dun6tWrZxuX0xrm5q233lLNmjV15swZTZ06VcuWLXO5Lddzv81So0YNl/NLliypqKiobJ+XUKVKlWz7cve+fz3z3LNnj+6///4cx+Rl165dkqS+fftec8yZM2eUlpamCxcuZFsL6cr3KD9B3r59u3766Sf16dNHu3fvdm5v166d3nrrLaWkpCg0NFT79u1z7vf36tSpo4ULF+rcuXMqUaJEnteJ3Fkd/3Pnzqlhw4bq37+/7rvvPrf2MWTIEH333Xd6/fXXVb9+fZ06dSrba8o3stKlSysqKkpbtmzJdl7Wf1y//48u6+CxJ5988po/+VevXt2zE71BNWnSRKVLl9ayZcuu+7LXemT2+wPnsvj6+ua43Rjj1nXPmTNHq1ev1tdff62FCxeqf//+mjBhglavXp3jM0HSlYPL2rZtq9DQUL3wwguqVq2agoKC9OOPP2r48OHZDiz05Jzz0qxZM+fR/vHx8WrdurV69eqlnTt3qmTJkoV6v83vkf35UVT/vrKu57XXXrvmWwBLlixZ4AMvJenjjz+WJA0dOlRDhw7Ndv7cuXPVr1+/Al8P8s/q+Hft2tV50FFO0tLS9Oyzz+qTTz7R6dOnVa9ePb3yyivOt59t375dkydP1pYtW5w/qeb0COBG161bN73//vtau3aty9PF15L1dLa/v3+ej7RjY2O1c+fObNuznh4ujM8NKGoZGRnOZ04iIyMVEhJyzdvs4+PjfNdE1oFop0+fdj4NL8n56McdsbGxzkd0V8tpPpJ022236bbbbtNLL72kWbNmqXfv3po9e7bL5y5cLTExUSdPntTnn3+u22+/3bk9KSmpQHPesmWLjDHZDqJzl6+vr8aNG6f27dtr0qRJGjFixHXdb7Ps2rVL7du3d55OTU1VcnKy/vjHP+Z5WXfv+9czz2rVquX4g/vVrvVDZrVq1SRJoaGhuV5PZGSkgoODr+t+dTVjjGbNmqX27dvrsccey3b+iy++qJkzZ6pfv37ONbnWukVERPCo30N4zT8Xjz/+uFatWqXZs2fr559/1gMPPKA777zT+Y/g66+/VtWqVfXNN9+oSpUqiouL08MPP1ysHvlL0tNPP62QkBD1799fR48ezXb+7x+llStXTu3atdOUKVOUnJycbXzW69+S9Mc//lFr167VqlWrnNvOnTund999V3Fxcbrllls8eEuK3pIlS5SamqqGDRtKuhKdzp0766uvvnJ5xuTo0aOaNWuWWrdurdDQUEn/95/v1c8anDt3znk0tDv++Mc/avXq1Vq7dq1z2/Hjx7Md1/Dbb79l+75mPfrL7ZFe1iP5qy976dIlvf322wWa8+HDhzVnzhzntvPnz+vdd991e5/SlaeUmzVrpokTJ+rixYvXdb/N8u6777ocEzF58mRdvnw51wcNWdy971/PPO+//35t2rRJX3zxRbZxWd+jrFiePn3a5fwmTZqoWrVqev3113N82S/renx9fdWlSxd9+eWX2r9/v/P87du3a+HChde6+U4rV67U3r171a9fP3Xv3j3b14MPPqglS5bo8OHDioqK0q233qoZM2a4zHfLli367rvv8vVDF/LH6kf+udm/f7+mTZum/fv3Kzo6WtKVp+EWLFigadOm6eWXX9avv/6qffv26bPPPtOHH36ojIwMDR06VN27dy/w25SKUo0aNTRr1iz17NlTtWrVcn7CnzFGSUlJmjVrlnx8fJwHTklXXl9t3bq16tevr0ceeURVq1bV0aNHtWrVKh08eFCbNm2SJI0YMUKffPKJunbtqsGDB6tMmTKaMWOGkpKSNHfu3EL74KDCcObMGefTl5cvX3a+9Ss4OFgjRoxwjhs7dqwSEhLUunVrPfbYY/Lz89OUKVOUlpbm8h73zp07q3LlyhowYICeeuop+fr6aurUqYqMjHT5T/Z6PP300/roo4905513asiQIc63+sXGxrq8zj5jxgy9/fbbuvfee1WtWjWdPXtW7733nkJDQ3P9D7Zly5YKDw9X3759NXjwYDkcDn300UcFehr/kUce0aRJk9SnTx9t2LBBUVFR+uijjzzyiYlPPfWUHnjgAU2fPl2PPvpovu+3WS5duqQOHTqoR48e2rlzp95++221bt1ad999d57XXZD7fn7n+dRTT2nOnDl64IEH1L9/fzVp0kSnTp3SvHnz9M4776hhw4aqVq2awsLC9M4776hUqVIqUaKEmjdvripVquj9999X165dVbduXfXr108VK1bUoUOHtGTJEoWGhurrr7+WJI0ZM0YLFixQmzZt9Nhjj+ny5cv6z3/+o7p16+b50d8zZ86Ur6+vunXrluP5d999t5599lnNnj1bw4YN02uvvaauXbuqRYsWGjBggPOtfqVLl3Z+xDc8wBtvMbgR6XdvOfnmm2+MJFOiRAmXLz8/P9OjRw9jjDGPPPKIkWR27tzpvNyGDRuMJLNjx46ivgkFtnv3bvP3v//dVK9e3QQFBZng4GBTu3Zt8+ijj5qNGzdmG79nzx7Tp08fU6FCBePv728qVqxo/vSnP5k5c+ZkG9e9e3cTFhZmgoKCTLNmzcw333xzzXm4+1a/QYMG5Xr7st5ad/z48VzH5eT3b/VzOBymTJky5u677zYbNmzINv7HH380Xbp0MSVLljQhISGmffv25ocffsg2bsOGDaZ58+YmICDAVK5c2bzxxhvXfKvf729z1rzatm3rsu3nn382bdu2NUFBQaZixYrmxRdfNB988IHLPn/88UfTs2dPU7lyZRMYGGjKlStn/vSnP5n169fnuRYrV640t912mwkODjbR0dHm6aefNgsXLsz2PWvbtq2pW7dutsv37ds329sr9+3bZ+6++24TEhJiIiIizJAhQ8yCBQuu661+Ob3lKyMjw1SrVs1Uq1bNXL582RiTv/tt1j6XLl1qBg4caMLDw03JkiVN7969Xd6CZsy1vzdZ15XXfT+nt/rld57GGHPy5Enz+OOPm4oVK5qAgABTqVIl07dvX3PixAnnmK+++srccsstxs/PL9t1/fTTT+a+++4zZcuWNYGBgSY2Ntb06NHDLFq0yOV6li5dapo0aWICAgJM1apVzTvvvHPNt6tmuXTpkilbtqxp06bNNccYY0yVKlVMo0aNnKe///5706pVKxMcHGxCQ0PNXXfdZbZt2+ZyGd7qVzAOYwrhyJtiyOFwuBzt/9///le9e/fW1q1bsx20VLJkSVWoUEGjRo3Syy+/7PK04IULFxQSEqLvvvuOj3sFiqnp06erX79+WrdunfMgwsKyZ88eVa9eXR999FG+Pp4a8ASe9r+GRo0aKSMjQ8eOHXN5b/XVWrVqpcuXL2vPnj3O129/+eUXSTfHgWwACl/W6/ruvP0UcJfV8U9NTXV5z2lSUpI2btyoMmXKqGbNmurdu7f69OmjCRMmqFGjRjp+/LgWLVqkBg0aqFu3burYsaMaN26s/v37a+LEicrMzNSgQYPUqVMn1axZ04u3DEBxMHXqVE2dOtX5uf9AUSk+R1sVgvXr16tRo0Zq1KiRJGnYsGFq1KiRnn/+eUnStGnT1KdPH/3jH/9QrVq1FB8fr3Xr1jk/bMXHx0dff/21IiIidPvtt6tbt26qU6eOZs+e7bXbBKD4GDhwoE6dOqXPPvvM5e2eQGHjNX8AACxj9SN/AABsRPwBALCMdQf8ZWZm6vDhwypVqlSev/UKAIDiwhijs2fPKjo6Os8PULMu/ocPH3Z+tjoAADebAwcOuHwia06si3/Wr7A9cOCA8zPWbZKenq7vvvtOnTt3dvk95cgd6+Ye1s09rJt7bF+3lJQUxcTEODuXG+vin/VUf2hoqLXxDwkJUWhoqJX/ONzFurmHdXMP6+Ye1u2K/LykzQF/AABYhvgDAGAZ4g8AgGWse80fAIAbkTFGly9fVkZGRo7n+/r6ys/PzyNvUyf+AAB42aVLl5ScnKzz58/nOi4kJERRUVEKCAgo0PURfwAAvCgzM1NJSUny9fVVdHS0AgICsj26N8bo0qVLOn78uJKSklSjRo08P8gnN8QfAAAvunTpkjIzMxUTE6OQkJBrjgsODpa/v7/27dunS5cuKSgoyO3r5IA/AABuAPl5JF+QR/su+/HIXgAAQLFB/AEAsAzxBwDAMsQfAADLEH8AAG4AxhiPjMkP4g8AgBdl/QbCvD7g5+oxBf2thbzPHwAAL/L19VVYWJiOHTsm6cqn+OX0IT/nz5/XsWPHFBYWJl9f3wJdJ/EHAMDLKlSoIEnOHwCuJSwszDm2IIg/AABe5nA4FBUVpXLlyik9PT3HMf7+/gV+xJ+F+AMAcIPw9fX1WOBzwwF/AABYhvgDAGAZ4g8AgGWIPwAAliH+AABYhvgDAGAZ4g8AgGWIPwAAliH+AABYhvgDAGAZ4g8AgGWIPwAAliH+AABYht/qV0CffSY9/7x09qy3Z5I/QUHShAlSnTrSxYvenk3xwbq5h3VzD+vmnuK4bqVKSS++KHXvXrTXS/wL6PnnpR07vD2L/AsOvvLn4cPShQvenUtxwrq5h3VzD+vmnuK6biNHEv9iJ+sRv4+PFBXl3bnkR1DQlT+jo4vPT8Y3AtbNPaybe1g39xS3dUtOljIzvfPMMfH3kKgo6eBBb88ib+np0vz50vbtkr+/t2dTfLBu7mHd3MO6uae4rVulStKhQ965bg74AwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACzj1fiPHj1aDofD5at27drXHP/ee++pTZs2Cg8PV3h4uDp27Ki1a9cW4YwBACj+vP7Iv27dukpOTnZ+rVix4ppjExMT1bNnTy1ZskSrVq1STEyMOnfurEOHDhXhjAEAKN78vD4BPz9VqFAhX2Nnzpzpcvr999/X3LlztWjRIvXp06cwpgcAwE3H6/HftWuXoqOjFRQUpBYtWmjcuHGqXLlyvi57/vx5paenq0yZMtcck5aWprS0NOfplJQUSVJ6errS09MLNnlJQUFScPCVPz2wu0KXdZs9cdttwrq5h3VzD+vmnuK2bp7ux/XcbocxxhT8Kt3z7bffKjU1VbVq1VJycrLGjBmjQ4cOacuWLSpVqlSel3/ssce0cOFCbd26VUFBQTmOGT16tMaMGZNt+6xZsxQSElLg2wAAwI3g/Pnz6tWrl86cOaPQ0NBcx3o1/r93+vRpxcbG6o033tCAAQNyHTt+/Hi9+uqrSkxMVIMGDa45LqdH/jExMTpx4kSei5MfdepIhw9L0dHS9u0F3l2hS09PV0JCgjp16iR/f39vT6fYYN3cw7q5h3VzT3FbN0/3IyUlRREREfmKv9ef9r9aWFiYatasqd27d+c67vXXX9f48eP1/fff5xp+SQoMDFRgYGC27f7+/h65c1y8KF24cOXPYnBfc/LU7bcN6+Ye1s09rJt7isu6ebof13ObvX60/9VSU1O1Z88eRUVFXXPMq6++qhdffFELFixQ06ZNi3B2AADcHLwa/yeffFJLly7V3r179cMPP+jee++Vr6+vevbsKUnq06ePnnnmGef4V155RSNHjtTUqVMVFxenI0eO6MiRI0pNTfXWTQAAoNjx6tP+Bw8eVM+ePXXy5ElFRkaqdevWWr16tSIjIyVJ+/fvl4/P//18MnnyZF26dEndu3d32c+oUaM0evToopw6AADFllfjP3v27FzPT0xMdDm9d+/ewpsMAACWuKFe8wcAAIWP+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFjGq/EfPXq0HA6Hy1ft2rVzvcxnn32m2rVrKygoSPXr19f8+fOLaLYAANwcvP7Iv27dukpOTnZ+rVix4ppjf/jhB/Xs2VMDBgzQTz/9pPj4eMXHx2vLli1FOGMAAIo3r8ffz89PFSpUcH5FRERcc+y///1v3XnnnXrqqadUp04dvfjii2rcuLEmTZpUhDMGAKB48/P2BHbt2qXo6GgFBQWpRYsWGjdunCpXrpzj2FWrVmnYsGEu27p06aIvv/zymvtPS0tTWlqa83RKSookKT09Xenp6QWef1CQFBx85U8P7K7QZd1mT9x2m7Bu7mHd3MO6uae4rZun+3E9t9thjDEFv0r3fPvtt0pNTVWtWrWUnJysMWPG6NChQ9qyZYtKlSqVbXxAQIBmzJihnj17Ore9/fbbGjNmjI4ePZrjdYwePVpjxozJtn3WrFkKCQnx3I0BAMCLzp8/r169eunMmTMKDQ3NdaxXH/l37drV+fcGDRqoefPmio2N1aeffqoBAwZ45DqeeeYZl2cLUlJSFBMTo86dO+e5OPlRp450+LAUHS1t317g3RW69PR0JSQkqFOnTvL39/f2dIoN1s09rJt7WDf3FLd183Q/sp7Zzg+vP+1/tbCwMNWsWVO7d+/O8fwKFSpke4R/9OhRVahQ4Zr7DAwMVGBgYLbt/v7+HrlzXLwoXbhw5c9icF9z8tTttw3r5h7WzT2sm3uKy7p5uh/Xc5u9fsDf1VJTU7Vnzx5FRUXleH6LFi20aNEil20JCQlq0aJFUUwPAICbglfj/+STT2rp0qXau3evfvjhB917773y9fV1vqbfp08fPfPMM87xQ4YM0YIFCzRhwgTt2LFDo0eP1vr16/X444976yYAAFDsePVp/4MHD6pnz546efKkIiMj1bp1a61evVqRkZGSpP3798vH5/9+PmnZsqVmzZql5557Tv/85z9Vo0YNffnll6pXr563bgIAAMWOV+M/e/bsXM9PTEzMtu2BBx7QAw88UEgzAgDg5ndDveYPAAAKH/EHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALON2/Pfs2aPnnntOPXv21LFjxyRJ3377rbZu3eqxyQEAAM9zK/5Lly5V/fr1tWbNGn3++edKTU2VJG3atEmjRo3y6AQBAIBnuRX/ESNGaOzYsUpISFBAQIBz+x133KHVq1d7bHIAAMDz3Ir/5s2bde+992bbXq5cOZ04caLAkwIAAIXHrfiHhYUpOTk52/affvpJFStWLPCkAABA4XEr/n/+8581fPhwHTlyRA6HQ5mZmVq5cqWefPJJ9enTx9NzBAAAHuRW/F9++WXVrl1bMTExSk1N1S233KLbb79dLVu21HPPPefpOQIAAA/yc+dCAQEBeu+99zRy5Eht2bJFqampatSokWrUqOHp+QEAAA9zK/5ZKleurMqVK3tqLgAAoAjkO/7Dhg3L907feOMNtyYDAAAKX77j/9NPP7mc/vHHH3X58mXVqlVLkvTLL7/I19dXTZo08ewMAQCAR+U7/kuWLHH+/Y033lCpUqU0Y8YMhYeHS5J+++039evXT23atPH8LAEAgMe4dbT/hAkTNG7cOGf4JSk8PFxjx47VhAkTPDY5AADgeW7FPyUlRcePH8+2/fjx4zp79myBJwUAAAqPW/G/99571a9fP33++ec6ePCgDh48qLlz52rAgAG67777PD1HAADgQW691e+dd97Rk08+qV69eik9Pf3Kjvz8NGDAAL322msenSAAAPAst+IfEhKit99+W6+99pr27NkjSapWrZpKlCjh0ckBAADPK9CH/JQoUUINGjTw1FwAAEARcCv+7du3l8PhuOb5ixcvdntCAACgcLkV/1tvvdXldHp6ujZu3KgtW7aob9++npgXAAAoJG7F/1//+leO20ePHq3U1NQCTQgAABQut97qdy1/+ctfNHXqVE/uEgAAeJhH479q1SoFBQV5cpcAAMDD3Hra//cf5GOMUXJystavX6+RI0d6ZGIAAKBwuBX/0NBQl6P9fXx8VKtWLb3wwgvq3LmzxyYHAAA8z634T58+3cPTAAAARcWt1/yrVq2qkydPZtt++vRpVa1atcCTAgAAhcet+O/du1cZGRnZtqelpenQoUMFnhQAACg81/W0/7x585x/X7hwoUqXLu08nZGRoUWLFikuLs5jkwMAAJ53XfGPj4+XJDkcjmyf5Ofv76+4uDhNmDDBY5MDAACed13xz8zMlCRVqVJF69atU0RERKFMCgAAFB63jvZPSkry9DwAAEARyXf833zzTQ0cOFBBQUF68803cx07ePDgAk8MAAAUjnzH/1//+pd69+6toKCga/5iH+nK8QDEHwCAG1e+43/1U/087Q8AQPHl1vv8X3jhBZ0/fz7b9gsXLuiFF14o8KQAAEDhcSv+Y8aMUWpqarbt58+f15gxYwo8KQAAUHjcir8xxuUX+2TZtGmTypQpU+BJAQCAwnNdb/ULDw+Xw+GQw+FQzZo1XX4AyMjIUGpqqh599FGPTxIAAHjOdcV/4sSJMsaof//+GjNmjMvH+wYEBCguLk4tWrTw+CQBAIDnXFf8sz7St0qVKmrZsqX8/f0LZVIAAKDwuPUJf23btnX+/eLFi7p06ZLL+aGhoQWbFQAAKDRuHfB3/vx5Pf744ypXrpxKlCih8PBwly8AAHDjciv+Tz31lBYvXqzJkycrMDBQ77//vsaMGaPo6Gh9+OGHnp4jAADwILee9v/666/14Ycfql27durXr5/atGmj6tWrKzY2VjNnzlTv3r09PU8AAOAhbj3yP3XqlKpWrSrpyuv7p06dkiS1bt1ay5Yt89zsAACAx7kV/6pVqzo/37927dr69NNPJV15RuDqt/8BAIAbj1vx79evnzZt2iRJGjFihN566y0FBQVp6NChevrppz06QQAA4FluveY/dOhQ5987duyoHTt2aMOGDYqIiNDHH3/ssckBAADPc+uR/+/FxsbqvvvuU+nSpfXBBx94YpcAAKCQeCT+AACg+CD+AABYhvgDAGCZ6zrg77777sv1/NOnTxdkLgAAoAhcV/zzeg9/6dKl1adPnwJNCAAAFK7riv+0adMKax4AAKCI8Jo/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWuWHiP378eDkcDj3xxBO5jps4caJq1aql4OBgxcTEaOjQobp48WLRTBIAgJuAn7cnIEnr1q3TlClT1KBBg1zHzZo1SyNGjNDUqVPVsmVL/fLLL3rooYfkcDj0xhtvFNFsAQAo3rz+yD81NVW9e/fWe++9p/Dw8FzH/vDDD2rVqpV69eqluLg4de7cWT179tTatWuLaLYAABR/Xn/kP2jQIHXr1k0dO3bU2LFjcx3bsmVLffzxx1q7dq2aNWumX3/9VfPnz9df//rXa14mLS1NaWlpztMpKSmSpPT0dKWnpxd4/kFBUnDwlT89sLtCl3WbPXHbbcK6uYd1cw/r5p7itm6e7sf13G6HMcYU/CrdM3v2bL300ktat26dgoKC1K5dO916662aOHHiNS/z5ptv6sknn5QxRpcvX9ajjz6qyZMnX3P86NGjNWbMmGzbZ82apZCQEE/cDAAAvO78+fPq1auXzpw5o9DQ0FzHei3+Bw4cUNOmTZWQkOB8rT+v+CcmJurPf/6zxo4dq+bNm2v37t0aMmSIHnnkEY0cOTLHy+T0yD8mJkYnTpzIc3Hyo04d6fBhKTpa2r69wLsrdOnp6UpISFCnTp3k7+/v7ekUG6ybe1g397Bu7ilu6+bpfqSkpCgiIiJf8ffa0/4bNmzQsWPH1LhxY+e2jIwMLVu2TJMmTVJaWpp8fX1dLjNy5Ej99a9/1cMPPyxJql+/vs6dO6eBAwfq2WeflY9P9kMYAgMDFRgYmG27v7+/R+4cFy9KFy5c+bMY3NecPHX7bcO6uYd1cw/r5p7ism6e7sf13Gavxb9Dhw7avHmzy7Z+/fqpdu3aGj58eLbwS1ee0vh94LPGefHVCwAAihWvxb9UqVKqV6+ey7YSJUqobNmyzu19+vRRxYoVNW7cOEnSXXfdpTfeeEONGjVyPu0/cuRI3XXXXTn+sAAAALLz+tH+udm/f7/LI/3nnntODodDzz33nA4dOqTIyEjdddddeumll7w4SwAAipcbKv6JiYm5nvbz89OoUaM0atSoopsUAAA3Ga9/yA8AAChaxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyxB8AAMsQfwAALEP8AQCwDPEHAMAyN0z8x48fL4fDoSeeeCLXcadPn9agQYMUFRWlwMBA1axZU/Pnzy+aSQIAcBPw8/YEJGndunWaMmWKGjRokOu4S5cuqVOnTipXrpzmzJmjihUrat++fQoLCyuaiQIAcBPwevxTU1PVu3dvvffeexo7dmyuY6dOnapTp07phx9+kL+/vyQpLi6uCGYJAMDNw+vxHzRokLp166aOHTvmGf958+apRYsWGjRokL766itFRkaqV69eGj58uHx9fXO8TFpamtLS0pynU1JSJEnp6elKT08v8PyDgqTg4Ct/emB3hS7rNnvittuEdXMP6+Ye1s09xW3dPN2P67ndDmOMKfhVumf27Nl66aWXtG7dOgUFBaldu3a69dZbNXHixBzH165dW3v37lXv3r312GOPaffu3Xrsscc0ePBgjRo1KsfLjB49WmPGjMm2fdasWQoJCfHkzQEAwGvOnz+vXr166cyZMwoNDc11rNfif+DAATVt2lQJCQnO1/rzin/NmjV18eJFJSUlOR/pv/HGG3rttdeUnJyc42VyeuQfExOjEydO5Lk4+VGnjnT4sBQdLW3fXuDdFbr09HQlJCSoU6dOzpdOkDfWzT2sm3tYN/cUt3XzdD9SUlIUERGRr/h77Wn/DRs26NixY2rcuLFzW0ZGhpYtW6ZJkyYpLS0t21P5UVFR8vf3d9lep04dHTlyRJcuXVJAQEC26wkMDFRgYGC27f7+/h65c1y8KF24cOXPYnBfc/LU7bcN6+Ye1s09rJt7isu6ebof13ObvRb/Dh06aPPmzS7b+vXrp9q1a1/zNfxWrVpp1qxZyszMlI/PlXcp/vLLL4qKisox/AAAIDuvvc+/VKlSqlevnstXiRIlVLZsWdWrV0+S1KdPHz3zzDPOy/z973/XqVOnNGTIEP3yyy/63//+p5dfflmDBg3y1s0AAKDY8frR/rnZv3+/8xG+JMXExGjhwoUaOnSoGjRooIoVK2rIkCEaPny4F2cJAEDxckPFPzExMdfTktSiRQutXr26aCYEAMBN6Ib5eF8AAFA0iD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWMbP2xO4WSQnS5UqeXsWeQsKkiZMkOrUkS5e9PZsig/WzT2sm3tYN/cUt3U7dMh71038C6hUqSt/ZmZ69xuZX8HBV/48fFi6cMG7cylOWDdJFddIlVdKezpJx+rn6yKsm3tYN/cU13UbNqzor5P4F9CLL0ojR0pnz3p7JvkTFHTlz+jo4vGT8Y3C9nVLL/2LjvW47cqJjABVmHVQvhcj87yc7evmLtbNPcVx3YYNI/7FUvfuV76Ki/R0af58aft2yd/f27MpPmxftxX7j+nhebW069QuZfpe0heJe3Rbpbzjb/u6uYt1cw/rln8c8AcgT60rt9aOx3co02Tq1gq3qklUE29PCUAB8MgfQL4dGnZIUSWj5HA4vD0VAAVA/AHkW3SpaG9PAYAH8LQ/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZYg/AACWIf4AAFiG+AMAYBniDwCAZW6Y+I8fP14Oh0NPPPFEvsbPnj1bDodD8fHxhTovAABuNjdE/NetW6cpU6aoQYMG+Rq/d+9ePfnkk2rTpk0hzwwAgJuP13+lb2pqqnr37q333ntPY8eOzXN8RkaGevfurTFjxmj58uU6ffp0ruPT0tKUlpbmPJ2SkiJJSk9PV3p6eoHmXhxl3WYbb3tBsG7uYd3cw7q5x/Z1u57b7fX4Dxo0SN26dVPHjh3zFf8XXnhB5cqV04ABA7R8+fI8x48bN05jxozJtv27775TSEiIW3O+GSQkJHh7CsUS6+Ye1s09rJt7bF238+fP53usV+M/e/Zs/fjjj1q3bl2+xq9YsUIffPCBNm7cmO/reOaZZzRs2DDn6ZSUFMXExKhz584KDQ293ikXe+np6UpISFCnTp3k7+/v7ekUG6ybe1g397Bu7rF93bKe2c4Pr8X/wIEDGjJkiBISEhQUFJTn+LNnz+qvf/2r3nvvPUVEROT7egIDAxUYGJhtu7+/v5V3jiy23353sW7uYd3cw7q5x9Z1u57b7LX4b9iwQceOHVPjxo2d2zIyMrRs2TJNmjRJaWlp8vX1dZ63Z88e7d27V3fddZdzW2ZmpiTJz89PO3fuVLVq1YruBgAAUEx5Lf4dOnTQ5s2bXbb169dPtWvX1vDhw13CL0m1a9fONv65557T2bNn9e9//1sxMTGFPmcAAG4GXot/qVKlVK9ePZdtJUqUUNmyZZ3b+/Tpo4oVK2rcuHEKCgrKNj4sLEySsm0HAADX5vWj/XOzf/9++fjcEB9FAADATeOGin9iYmKup39v+vTphTYXAABuVjysBgDAMsQfAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsMwN9SE/RcEYI+n6fvXhzSQ9PV3nz59XSkqKlb/1yl2sm3tYN/ewbu6xfd2yupbVudxYF/+zZ89KEr8ICABwUzp79qxKly6d6xiHyc+PCDeRzMxMHT58WKVKlZLD4fD2dIpcSkqKYmJidODAAYWGhnp7OsUG6+Ye1s09rJt7bF83Y4zOnj2r6OjoPH8vjnWP/H18fFSpUiVvT8PrQkNDrfzHUVCsm3tYN/ewbu6xed3yesSfhQP+AACwDPEHAMAyxN8ygYGBGjVqlAIDA709lWKFdXMP6+Ye1s09rFv+WXfAHwAAtuORPwAAliH+AABYhvgDAGAZ4g8AgGWI/w1m8uTJatCggfNDKlq0aKFvv/3Wef7Fixc1aNAglS1bViVLltT999+vo0ePuuxj//796tatm0JCQlSuXDk99dRTunz5ssuYxMRENW7cWIGBgapevbqmT5+ebS5vvfWW4uLiFBQUpObNm2vt2rUu5+dnLt4yfvx4ORwOPfHEE85trF12o0ePlsPhcPmqXbv2dc3TtjXLcujQIf3lL39R2bJlFRwcrPr162v9+vXO840xev755xUVFaXg4GB17NhRu3btctnHqVOn1Lt3b4WGhiosLEwDBgxQamqqy5iff/5Zbdq0UVBQkGJiYvTqq69mm8tnn32m2rVrKygoSPXr19f8+fNdzs/PXIpCXFxctvubw+HQoEGDJHF/K1IGN5R58+aZ//3vf+aXX34xO3fuNP/85z+Nv7+/2bJlizHGmEcffdTExMSYRYsWmfXr15vbbrvNtGzZ0nn5y5cvm3r16pmOHTuan376ycyfP99ERESYZ555xjnm119/NSEhIWbYsGFm27Zt5j//+Y/x9fU1CxYscI6ZPXu2CQgIMFOnTjVbt241jzzyiAkLCzNHjx51jslrLt6ydu1aExcXZxo0aGCGDBni3M7aZTdq1ChTt25dk5yc7Pw6fvx4vudp45oZY8ypU6dMbGyseeihh8yaNWvMr7/+ahYuXGh2797tHDN+/HhTunRp8+WXX5pNmzaZu+++21SpUsVcuHDBOebOO+80DRs2NKtXrzbLly831atXNz179nSef+bMGVO+fHnTu3dvs2XLFvPJJ5+Y4OBgM2XKFOeYlStXGl9fX/Pqq6+abdu2meeee874+/ubzZs3X9dcisKxY8dc7msJCQlGklmyZIkxhvtbUSL+xUB4eLh5//33zenTp42/v7/57LPPnOdt377dSDKrVq0yxhgzf/584+PjY44cOeIcM3nyZBMaGmrS0tKMMcY8/fTTpm7dui7X8eCDD5ouXbo4Tzdr1swMGjTIeTojI8NER0ebcePGGWNMvubiDWfPnjU1atQwCQkJpm3bts74s3Y5GzVqlGnYsGGO57Fm1zZ8+HDTunXra56fmZlpKlSoYF577TXnttOnT5vAwEDzySefGGOM2bZtm5Fk1q1b5xzz7bffGofDYQ4dOmSMMebtt9824eHhzrXMuu5atWo5T/fo0cN069bN5fqbN29u/va3v+V7Lt4yZMgQU61aNZOZmcn9rYjxtP8NLCMjQ7Nnz9a5c+fUokULbdiwQenp6erYsaNzTO3atVW5cmWtWrVKkrRq1SrVr19f5cuXd47p0qWLUlJStHXrVueYq/eRNSZrH5cuXdKGDRtcxvj4+Khjx47OMfmZizcMGjRI3bp1y3b7WLtr27Vrl6Kjo1W1alX17t1b+/fvz/c8bV2zefPmqWnTpnrggQdUrlw5NWrUSO+9957z/KSkJB05csRlvqVLl1bz5s1d1i4sLExNmzZ1junYsaN8fHy0Zs0a55jbb79dAQEBzjFdunTRzp079dtvvznH5La++ZmLN1y6dEkff/yx+vfvL4fDwf2tiBH/G9DmzZtVsmRJBQYG6tFHH9UXX3yhW265RUeOHFFAQIDCwsJcxpcvX15HjhyRJB05csTlH0bW+Vnn5TYmJSVFFy5c0IkTJ5SRkZHjmKv3kddcitrs2bP1448/aty4cdnOY+1y1rx5c02fPl0LFizQ5MmTlZSUpDZt2ujs2bOsWS5+/fVXTZ48WTVq1NDChQv197//XYMHD9aMGTOc882a37Xme+TIEZUrV87lfD8/P5UpU8Yj63v1+XnNxRu+/PJLnT59Wg899JAk/o0WNet+q19xUKtWLW3cuFFnzpzRnDlz1LdvXy1dutTb07qhHThwQEOGDFFCQoKCgoK8PZ1io2vXrs6/N2jQQM2bN1dsbKw+/fRTBQcHe3FmN7bMzEw1bdpUL7/8siSpUaNG2rJli9555x317dvXy7MrHj744AN17dpV0dHR3p6KlXjkfwMKCAhQ9erV1aRJE40bN04NGzbUv//9b1WoUEGXLl3S6dOnXcYfPXpUFSpUkCRVqFAh2xGpWafzGhMaGqrg4GBFRETI19c3xzFX7yOvuRSlDRs26NixY2rcuLH8/Pzk5+enpUuX6s0335Sfn5/Kly/P2uVDWFiYatasqd27d3N/y0VUVJRuueUWl2116tRxvmSSNae8btOxY8dczr98+bJOnTrlkfW9+vy85lLU9u3bp++//14PP/ywcxv3t6JF/IuBzMxMpaWlqUmTJvL399eiRYuc5+3cuVP79+9XixYtJEktWrTQ5s2bXf5TSUhIUGhoqPM/qxYtWrjsI2tM1j4CAgLUpEkTlzGZmZlatGiRc0x+5lKUOnTooM2bN2vjxo3Or6ZNm6p3797Ov7N2eUtNTdWePXsUFRXF/S0XrVq10s6dO122/fLLL4qNjZUkValSRRUqVHCZb0pKitasWeOydqdPn9aGDRucYxYvXqzMzEw1b97cOWbZsmVKT093jklISFCtWrUUHh7uHJPb+uZnLkVt2rRpKleunLp16+bcxv2tiHn7iEO4GjFihFm6dKlJSkoyP//8sxkxYoRxOBzmu+++M8ZceftJ5cqVzeLFi8369etNixYtTIsWLZyXz3orTOfOnc3GjRvNggULTGRkZI5vhXnqqafM9u3bzVtvvZXjW2ECAwPN9OnTzbZt28zAgQNNWFiYy1G2ec3F264+2t8Y1i4n//jHP0xiYqJJSkoyK1euNB07djQRERHm2LFj+ZqnjWtmzJW3k/r5+ZmXXnrJ7Nq1y8ycOdOEhISYjz/+2Dlm/PjxJiwszHz11Vfm559/Nvfcc0+Ob/Vr1KiRWbNmjVmxYoWpUaOGy1v9Tp8+bcqXL2/++te/mi1btpjZs2ebkJCQbG/18/PzM6+//rrZvn27GTVqVI5v9ctrLkUlIyPDVK5c2QwfPjzbedzfig7xv8H079/fxMbGmoCAABMZGWk6dOjgDL8xxly4cME89thjJjw83ISEhJh7773XJCcnu+xj7969pmvXriY4ONhERESYf/zjHyY9Pd1lzJIlS8ytt95qAgICTNWqVc20adOyzeU///mPqVy5sgkICDDNmjUzq1evdjk/P3Pxpt/Hn7XL7sEHHzRRUVEmICDAVKxY0Tz44IMu71Vnza7t66+/NvXq1TOBgYGmdu3a5t1333U5PzMz04wcOdKUL1/eBAYGmg4dOpidO3e6jDl58qTp2bOnKVmypAkNDTX9+vUzZ8+edRmzadMm07p1axMYGGgqVqxoxo8fn20un376qalZs6YJCAgwdevWNf/73/+uey5FZeHChUZSjtfP/a3o8Ct9AQCwDK/5AwBgGeIPAIBliD8AAJYh/gAAWIb4AwBgGeIPAIBliD8AAJYh/gAAWIb4Ayg0e/fulcPh0MaNGwtl/w6HQ19++WWh7Bu4mRF/4Cb20EMPKT4+3mvXHxMTo+TkZNWrV0+SlJiYKIfDke23pQEoWn7engCAm5evr+9N8ytQgZsJj/wBSy1dulTNmjVTYGCgoqKiNGLECF2+fNl5frt27TR48GA9/fTTKlOmjCpUqKDRo0e77GPHjh1q3bq1goKCdMstt+j77793eSr+6qf99+7dq/bt20uSwsPD5XA49NBDD0mS4uLiNHHiRJd933rrrS7Xt2vXLt1+++3O60pISMh2mw4cOKAePXooLCxMZcqU0T333KO9e/cWdKmAmw7xByx06NAh/fGPf9Qf/vAHbdq0SZMnT9YHH3ygsWPHuoybMWOGSpQooTVr1ujVV1/VCy+84IxuRkaG4uPjFRISojVr1ujdd9/Vs88+e83rjImJ0dy5cyVd+d3oycnJ+ve//52v+WZmZuq+++5TQECA1qxZo3feeUfDhw93GZOenq4uXbqoVKlSWr58uVauXKmSJUvqzjvv1KVLl65neYCbHk/7AxZ6++23FRMTo0mTJsnhcKh27do6fPiwhg8frueff14+PlceFzRo0ECjRo2SJNWoUUOTJk3SokWL1KlTJyUkJGjPnj1KTEx0PrX/0ksvqVOnTjlep6+vr8qUKSNJKleunMLCwvI93++//147duzQwoULFR0dLUl6+eWX1bVrV+eY//73v8rMzNT7778vh8MhSZo2bZrCwsKUmJiozp07X98iATcx4g9YaPv27WrRooUzkpLUqlUrpaam6uDBg6pcubKkK/G/WlRUlI4dOybpyqP3mJgYl9f0mzVrVmjzjYmJcYZfklq0aOEyZtOmTdq9e7dKlSrlsv3ixYvas2dPocwLKK6IP4Br8vf3dzntcDiUmZnp8evx8fGRMcZlW3p6+nXtIzU1VU2aNNHMmTOznRcZGVmg+QE3G+IPWKhOnTqaO3eujDHOR/8rV65UqVKlVKlSpXzto1atWjpw4ICOHj2q8uXLS5LWrVuX62UCAgIkXTle4GqRkZFKTk52nk5JSVFSUpLLfA8cOKDk5GRFRUVJklavXu2yj8aNG+u///2vypUrp9DQ0HzdBsBWHPAH3OTOnDmjjRs3unwNHDhQBw4c0P/7f/9PO3bs0FdffaVRo0Zp2LBhztf789KpUydVq1ZNffv21c8//6yVK1fqueeekySXlxOuFhsbK4fDoW+++UbHjx9XamqqJOmOO+7QRx99pOXLl2vz5s3q27evfH19nZfr2LGjatasqb59+2rTpk1avnx5toMLe/furYiICN1zzz1avny5kpKSlJiYqMGDB+vgwYPuLB1w0yL+wE0uMTFRjRo1cvl68cUXNX/+fK1du1YNGzbUo48+qgEDBjjjnR++vr768ssvlZqaqj/84Q96+OGHnUEOCgrK8TIVK1bUmDFjNGLECJUvX16PP/64JOmZZ55R27Zt9ac//UndunVTfHy8qlWr5rycj4+PvvjiC124cEHNmjXTww8/rJdeesll3yEhIVq2bJkqV66s++67T3Xq1NGAAQN08eJFngkAfsdhfv9CGwC4aeXKlWrdurV2797tEm8ANxbiD8BtX3zxhUqWLKkaNWpo9+7dGjJkiMLDw7VixQpvTw1ALjjgD4Dbzp49q+HDh2v//v2KiIhQx44dNWHCBG9PC0AeeOQPAIBlOOAPAADLEH8AACxD/AEAsAzxBwDAMsQfAADLEH8AACxD/AEAsAzxBwDAMv8fH6o0JVp5dl8AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1000x1000 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "ename": "ValueError",
     "evalue": "Input shapes do not overlap raster.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mWindowError\u001b[0m                               Traceback (most recent call last)",
      "File \u001b[0;32m~/dev/model-my-watershed/.venv/lib/python3.10/site-packages/rasterio/mask.py:80\u001b[0m, in \u001b[0;36mraster_geometry_mask\u001b[0;34m(dataset, shapes, all_touched, invert, crop, pad, pad_width)\u001b[0m\n\u001b[1;32m     79\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[0;32m---> 80\u001b[0m     window \u001b[38;5;241m=\u001b[39m \u001b[43mgeometry_window\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdataset\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mshapes\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpad_x\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mpad_x\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpad_y\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mpad_y\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     82\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m WindowError:\n\u001b[1;32m     83\u001b[0m     \u001b[38;5;66;03m# If shapes do not overlap raster, raise Exception or UserWarning\u001b[39;00m\n\u001b[1;32m     84\u001b[0m     \u001b[38;5;66;03m# depending on value of crop\u001b[39;00m\n",
      "File \u001b[0;32m~/dev/model-my-watershed/.venv/lib/python3.10/site-packages/rasterio/features.py:477\u001b[0m, in \u001b[0;36mgeometry_window\u001b[0;34m(dataset, shapes, pad_x, pad_y, north_up, rotated, pixel_precision, boundless)\u001b[0m\n\u001b[1;32m    476\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m boundless:\n\u001b[0;32m--> 477\u001b[0m     window \u001b[38;5;241m=\u001b[39m \u001b[43mwindow\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mintersection\u001b[49m\u001b[43m(\u001b[49m\u001b[43mraster_window\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    479\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m window\n",
      "File \u001b[0;32m~/dev/model-my-watershed/.venv/lib/python3.10/site-packages/rasterio/windows.py:775\u001b[0m, in \u001b[0;36mWindow.intersection\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m    763\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"Return the intersection of this window and another\u001b[39;00m\n\u001b[1;32m    764\u001b[0m \n\u001b[1;32m    765\u001b[0m \u001b[38;5;124;03mParameters\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    773\u001b[0m \u001b[38;5;124;03mWindow\u001b[39;00m\n\u001b[1;32m    774\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[0;32m--> 775\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mintersection\u001b[49m\u001b[43m(\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mother\u001b[49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/dev/model-my-watershed/.venv/lib/python3.10/site-packages/rasterio/windows.py:125\u001b[0m, in \u001b[0;36miter_args.<locals>.wrapper\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    124\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(args) \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m1\u001b[39m \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(args[\u001b[38;5;241m0\u001b[39m], Iterable):\n\u001b[0;32m--> 125\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mfunction\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    126\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n",
      "File \u001b[0;32m~/dev/model-my-watershed/.venv/lib/python3.10/site-packages/rasterio/windows.py:239\u001b[0m, in \u001b[0;36mintersection\u001b[0;34m(*windows)\u001b[0m\n\u001b[1;32m    226\u001b[0m \u001b[38;5;250m\u001b[39m\u001b[38;5;124;03m\"\"\"Innermost extent of window intersections.\u001b[39;00m\n\u001b[1;32m    227\u001b[0m \n\u001b[1;32m    228\u001b[0m \u001b[38;5;124;03mWill raise WindowError if windows do not intersect.\u001b[39;00m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m    237\u001b[0m \u001b[38;5;124;03mWindow\u001b[39;00m\n\u001b[1;32m    238\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[0;32m--> 239\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mfunctools\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mreduce\u001b[49m\u001b[43m(\u001b[49m\u001b[43m_intersection\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mwindows\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/dev/model-my-watershed/.venv/lib/python3.10/site-packages/rasterio/windows.py:257\u001b[0m, in \u001b[0;36m_intersection\u001b[0;34m(w1, w2)\u001b[0m\n\u001b[1;32m    256\u001b[0m \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[0;32m--> 257\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m WindowError(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mIntersection is empty \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mw1\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mw2\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n",
      "\u001b[0;31mWindowError\u001b[0m: Intersection is empty Window(col_off=23928, row_off=89570, width=101, height=101) Window(col_off=0, row_off=0, width=51218, height=89289)",
      "\nDuring handling of the above exception, another exception occurred:\n",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[38], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43mquery_histogram\u001b[49m\u001b[43m(\u001b[49m\u001b[43mgeojson\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mURL\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mCOLLECTION\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mASSET\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mFILTER\u001b[49m\u001b[43m)\u001b[49m\n",
      "Cell \u001b[0;32mIn[36], line 33\u001b[0m, in \u001b[0;36mquery_histogram\u001b[0;34m(geojson, url, collection, asset, filter, cachekey)\u001b[0m\n\u001b[1;32m     31\u001b[0m clips \u001b[38;5;241m=\u001b[39m []\n\u001b[1;32m     32\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m tiff \u001b[38;5;129;01min\u001b[39;00m tiffs:\n\u001b[0;32m---> 33\u001b[0m     clip_data, clip_transform, clip_meta \u001b[38;5;241m=\u001b[39m \u001b[43mclip_and_reproject_tile\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m     34\u001b[0m \u001b[43m        \u001b[49m\u001b[43mtiff\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43maoi\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdst_crs\u001b[49m\n\u001b[1;32m     35\u001b[0m \u001b[43m    \u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m     36\u001b[0m     temp_file \u001b[38;5;241m=\u001b[39m tempfile\u001b[38;5;241m.\u001b[39mNamedTemporaryFile(delete\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mFalse\u001b[39;00m, suffix\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m.tif\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m     37\u001b[0m     \u001b[38;5;28;01mwith\u001b[39;00m rio\u001b[38;5;241m.\u001b[39mopen(temp_file\u001b[38;5;241m.\u001b[39mname, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mw\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mclip_meta) \u001b[38;5;28;01mas\u001b[39;00m dst:\n",
      "Cell \u001b[0;32mIn[35], line 42\u001b[0m, in \u001b[0;36mclip_and_reproject_tile\u001b[0;34m(tiff, aoi, dst_crs, resampling)\u001b[0m\n\u001b[1;32m     39\u001b[0m plt\u001b[38;5;241m.\u001b[39mshow()\n\u001b[1;32m     41\u001b[0m reprojected_aoi_bbox \u001b[38;5;241m=\u001b[39m BoundingBox(\u001b[38;5;241m*\u001b[39mshape(reprojected_aoi)\u001b[38;5;241m.\u001b[39mbounds)\n\u001b[0;32m---> 42\u001b[0m clip_data, clip_transform \u001b[38;5;241m=\u001b[39m \u001b[43mmask\u001b[49m\u001b[43m(\u001b[49m\u001b[43msrc\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43m[\u001b[49m\u001b[43mreprojected_aoi\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcrop\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m)\u001b[49m\n\u001b[1;32m     44\u001b[0m \u001b[38;5;66;03m# Define the output metadata for the reprojected clip\u001b[39;00m\n\u001b[1;32m     45\u001b[0m dst_transform, width, height \u001b[38;5;241m=\u001b[39m calculate_default_transform(\n\u001b[1;32m     46\u001b[0m     src\u001b[38;5;241m.\u001b[39mcrs,\n\u001b[1;32m     47\u001b[0m     dst_crs,\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m     50\u001b[0m     \u001b[38;5;241m*\u001b[39mreprojected_aoi_bbox,\n\u001b[1;32m     51\u001b[0m )\n",
      "File \u001b[0;32m~/dev/model-my-watershed/.venv/lib/python3.10/site-packages/rasterio/mask.py:178\u001b[0m, in \u001b[0;36mmask\u001b[0;34m(dataset, shapes, all_touched, invert, nodata, filled, crop, pad, pad_width, indexes)\u001b[0m\n\u001b[1;32m    175\u001b[0m     \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m    176\u001b[0m         nodata \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m0\u001b[39m\n\u001b[0;32m--> 178\u001b[0m shape_mask, transform, window \u001b[38;5;241m=\u001b[39m \u001b[43mraster_geometry_mask\u001b[49m\u001b[43m(\u001b[49m\n\u001b[1;32m    179\u001b[0m \u001b[43m    \u001b[49m\u001b[43mdataset\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mshapes\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mall_touched\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mall_touched\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43minvert\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43minvert\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcrop\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mcrop\u001b[49m\u001b[43m,\u001b[49m\n\u001b[1;32m    180\u001b[0m \u001b[43m    \u001b[49m\u001b[43mpad\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mpad\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpad_width\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mpad_width\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    182\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m indexes \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m    183\u001b[0m     out_shape \u001b[38;5;241m=\u001b[39m (dataset\u001b[38;5;241m.\u001b[39mcount, ) \u001b[38;5;241m+\u001b[39m shape_mask\u001b[38;5;241m.\u001b[39mshape\n",
      "File \u001b[0;32m~/dev/model-my-watershed/.venv/lib/python3.10/site-packages/rasterio/mask.py:86\u001b[0m, in \u001b[0;36mraster_geometry_mask\u001b[0;34m(dataset, shapes, all_touched, invert, crop, pad, pad_width)\u001b[0m\n\u001b[1;32m     82\u001b[0m \u001b[38;5;28;01mexcept\u001b[39;00m WindowError:\n\u001b[1;32m     83\u001b[0m     \u001b[38;5;66;03m# If shapes do not overlap raster, raise Exception or UserWarning\u001b[39;00m\n\u001b[1;32m     84\u001b[0m     \u001b[38;5;66;03m# depending on value of crop\u001b[39;00m\n\u001b[1;32m     85\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m crop:\n\u001b[0;32m---> 86\u001b[0m         \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mInput shapes do not overlap raster.\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m     87\u001b[0m     \u001b[38;5;28;01melse\u001b[39;00m:\n\u001b[1;32m     88\u001b[0m         warnings\u001b[38;5;241m.\u001b[39mwarn(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mshapes are outside bounds of raster. \u001b[39m\u001b[38;5;124m'\u001b[39m\n\u001b[1;32m     89\u001b[0m                       \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mAre they in different coordinate reference systems?\u001b[39m\u001b[38;5;124m'\u001b[39m)\n",
      "\u001b[0;31mValueError\u001b[0m: Input shapes do not overlap raster."
     ]
    }
   ],
   "source": [
    "query_histogram(geojson, URL, COLLECTION, ASSET, FILTER)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.10.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}