{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8dd35554-a9dd-49cf-b9fa-24fa8ae6cecf",
   "metadata": {},
   "source": [
    "# Burn scar analysis using embeddings from partial inputs\n",
    "This notebook contains a complete example for how to run Clay. It\n",
    "combines the following three different aspects\n",
    "\n",
    "1. Create single-chip datacubes with time series data for a location and a date range\n",
    "2. Run the model with partial inputs, in this case RGB + NIR\n",
    "3. Study burn scares through the embeddings generated for that datacube\n",
    "\n",
    "## Let's start with importing and creating constants"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b7bcff1e-bdb5-47f8-aa0e-d68d6fdd3476",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ensure working directory is the repo home\n",
    "import os\n",
    "\n",
    "os.chdir(\"..\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "15d65ec9-86aa-4275-89ba-ec79fdbad361",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "from pathlib import Path\n",
    "\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy\n",
    "import pandas as pd\n",
    "import pystac_client\n",
    "import rasterio\n",
    "import rioxarray  # noqa: F401\n",
    "import stackstac\n",
    "import torch\n",
    "from rasterio.enums import Resampling\n",
    "from shapely import Point\n",
    "from sklearn import decomposition\n",
    "\n",
    "from src.datamodule import ClayDataModule\n",
    "from src.model_clay import CLAYModule\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "BAND_GROUPS = {\n",
    "    \"rgb\": [\"red\", \"green\", \"blue\"],\n",
    "    \"rededge\": [\"rededge1\", \"rededge2\", \"rededge3\", \"nir08\"],\n",
    "    \"nir\": [\n",
    "        \"nir\",\n",
    "    ],\n",
    "    \"swir\": [\"swir16\", \"swir22\"],\n",
    "    \"sar\": [\"vv\", \"vh\"],\n",
    "}\n",
    "\n",
    "STAC_API = \"https://earth-search.aws.element84.com/v1\"\n",
    "COLLECTION = \"sentinel-2-l2a\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6341305-9c44-4a1e-847c-80d77b01c0bf",
   "metadata": {},
   "source": [
    "## Search for imagery over an area of interest\n",
    "In this example we use a location and date range to visualize a forest fire that happened in [Monchique in 2018](https://pt.wikipedia.org/wiki/Inc%C3%AAndio_de_Monchique_de_2018)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "a1886f5a-8669-40e7-8fae-e45619570e3c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found 12 items\n"
     ]
    }
   ],
   "source": [
    "# Point over Monchique Portugal\n",
    "poi = 37.30939, -8.57207\n",
    "\n",
    "# Dates of a large forest fire\n",
    "start = \"2018-07-01\"\n",
    "end = \"2018-09-01\"\n",
    "\n",
    "catalog = pystac_client.Client.open(STAC_API)\n",
    "\n",
    "search = catalog.search(\n",
    "    collections=[COLLECTION],\n",
    "    datetime=f\"{start}/{end}\",\n",
    "    bbox=(poi[1] - 1e-5, poi[0] - 1e-5, poi[1] + 1e-5, poi[0] + 1e-5),\n",
    "    max_items=100,\n",
    "    query={\"eo:cloud_cover\": {\"lt\": 80}},\n",
    ")\n",
    "\n",
    "items = search.get_all_items()\n",
    "\n",
    "print(f\"Found {len(items)} items\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4ba5c36-90a6-427c-80c5-2a83ad11a1b0",
   "metadata": {},
   "source": [
    "## Download the data\n",
    "Get the data into a numpy array and visualize the imagery. The burn scar is visible in the last five images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c371501c-3ef0-4507-9073-0521a1c733be",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract coordinate system from first item\n",
    "epsg = items[0].properties[\"proj:epsg\"]\n",
    "\n",
    "# Convert point into the image projection\n",
    "poidf = gpd.GeoDataFrame(\n",
    "    pd.DataFrame(),\n",
    "    crs=\"EPSG:4326\",\n",
    "    geometry=[Point(poi[1], poi[0])],\n",
    ").to_crs(epsg)\n",
    "\n",
    "coords = poidf.iloc[0].geometry.coords[0]\n",
    "\n",
    "# Create bounds of the correct size, the model\n",
    "# requires 512x512 pixels at 10m resolution.\n",
    "bounds = (\n",
    "    coords[0] - 2560,\n",
    "    coords[1] - 2560,\n",
    "    coords[0] + 2560,\n",
    "    coords[1] + 2560,\n",
    ")\n",
    "\n",
    "# Retrieve the pixel values, for the bounding box in\n",
    "# the target projection. In this example we use only\n",
    "# the RGB and NIR band groups.\n",
    "stack = stackstac.stack(\n",
    "    items,\n",
    "    bounds=bounds,\n",
    "    snap_bounds=False,\n",
    "    epsg=epsg,\n",
    "    resolution=10,\n",
    "    dtype=\"float32\",\n",
    "    rescale=False,\n",
    "    fill_value=0,\n",
    "    assets=BAND_GROUPS[\"rgb\"] + BAND_GROUPS[\"nir\"],\n",
    "    resampling=Resampling.nearest,\n",
    ")\n",
    "\n",
    "stack = stack.compute()\n",
    "\n",
    "stack.sel(band=[\"red\", \"green\", \"blue\"]).plot.imshow(\n",
    "    row=\"time\", rgb=\"band\", vmin=0, vmax=2000, col_wrap=6\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce633fb1-fc82-4c88-8204-cda47aa9c874",
   "metadata": {},
   "source": [
    "![Minicube visualization](https://github.com/Clay-foundation/model/assets/901647/c6e924e5-6ba1-4924-b99a-df8b90731a5f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77e7c22c-1bfd-4281-bb12-8330c3eedc25",
   "metadata": {},
   "source": [
    "## Write data to tif files\n",
    "To use the mini datacube in the Clay dataloader, we need to write the\n",
    "images to tif files on disk. These tif files are then used by the Clay\n",
    "data loader for creating embeddings below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6509c3b2-a67c-447d-a7a1-e5fbcc1e35b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "outdir = Path(\"data/minicubes\")\n",
    "outdir.mkdir(exist_ok=True, parents=True)\n",
    "\n",
    "# Write tile to output dir\n",
    "for tile in stack:\n",
    "    # Grid code like MGRS-29SNB\n",
    "    mgrs = str(tile.coords[\"grid:code\"].values).split(\"-\")[1]\n",
    "    date = str(tile.time.values)[:10]\n",
    "\n",
    "    name = \"{dir}/claytile_{mgrs}_{date}.tif\".format(\n",
    "        dir=outdir,\n",
    "        mgrs=mgrs,\n",
    "        date=date.replace(\"-\", \"\"),\n",
    "    )\n",
    "    tile.rio.to_raster(name, compress=\"deflate\")\n",
    "\n",
    "    with rasterio.open(name, \"r+\") as rst:\n",
    "        rst.update_tags(date=date)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebc4b6ee-db58-4005-9689-a7d0acdc6a79",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "## Create embeddings\n",
    "Now switch gears and load the tiles to create embeddings and analyze them. \n",
    "\n",
    "The model checkpoint can be loaded directly from huggingface, and the data\n",
    "directory points to the directory we created in the steps above.\n",
    "\n",
    "Note that the normalization parameters for the data module need to be \n",
    "adapted based on the band groups that were selected as partial input. The\n",
    "full set of normalization parameters can be found [here](https://github.com/Clay-foundation/model/blob/main/src/datamodule.py#L108)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d89e0135-9473-4f76-9f09-e4e295dd51c9",
   "metadata": {},
   "source": [
    "### Load the model and set up the data module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "301ee2db-c5fc-4628-b837-12e6ea477415",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total number of chips: 12\n"
     ]
    }
   ],
   "source": [
    "DATA_DIR = \"data/minicubes\"\n",
    "CKPT_PATH = \"https://huggingface.co/made-with-clay/Clay/resolve/main/Clay_v0.1_epoch-24_val-loss-0.46.ckpt\"\n",
    "\n",
    "# Load model\n",
    "rgb_model = CLAYModule.load_from_checkpoint(\n",
    "    CKPT_PATH,\n",
    "    mask_ratio=0.0,\n",
    "    band_groups={\"rgb\": (2, 1, 0), \"nir\": (3,)},\n",
    "    bands=4,\n",
    "    strict=False,  # ignore the extra parameters in the checkpoint\n",
    ")\n",
    "# Set the model to evaluation mode\n",
    "rgb_model.eval()\n",
    "\n",
    "\n",
    "# Load the datamodule, with the reduced set of\n",
    "class ClayDataModuleRGB(ClayDataModule):\n",
    "    MEAN = [\n",
    "        1369.03,  # red\n",
    "        1597.68,  # green\n",
    "        1741.10,  # blue\n",
    "        2858.43,  # nir\n",
    "    ]\n",
    "    STD = [\n",
    "        2026.96,  # red\n",
    "        2011.88,  # green\n",
    "        2146.35,  # blue\n",
    "        2016.38,  # nir\n",
    "    ]\n",
    "\n",
    "\n",
    "data_dir = Path(DATA_DIR)\n",
    "\n",
    "dm = ClayDataModuleRGB(data_dir=str(data_dir.absolute()), batch_size=20)\n",
    "dm.setup(stage=\"predict\")\n",
    "trn_dl = iter(dm.predict_dataloader())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db3f3e5e-8668-4830-9c77-cc1d8cb35234",
   "metadata": {},
   "source": [
    "### Create the embeddings for the images over the forest fire\n",
    "This will loop through the images returned by the data loader\n",
    "and evaluate the model for each one of the images. The raw\n",
    "embeddings are reduced to mean values to simplify the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c5762240-9d22-4ebd-8e39-83fc6594a459",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average embeddings have shape (12, 768)\n"
     ]
    }
   ],
   "source": [
    "embeddings = []\n",
    "\n",
    "for batch in trn_dl:\n",
    "    with torch.inference_mode():\n",
    "        # Move data from to the device of model\n",
    "        batch[\"pixels\"] = batch[\"pixels\"].to(rgb_model.device)\n",
    "        # Pass just the specific band through the model\n",
    "        batch[\"timestep\"] = batch[\"timestep\"].to(rgb_model.device)\n",
    "        batch[\"latlon\"] = batch[\"latlon\"].to(rgb_model.device)\n",
    "\n",
    "        # Pass pixels, latlon, timestep through the encoder to create encoded patches\n",
    "        (\n",
    "            unmasked_patches,\n",
    "            unmasked_indices,\n",
    "            masked_indices,\n",
    "            masked_matrix,\n",
    "        ) = rgb_model.model.encoder(batch)\n",
    "\n",
    "        embeddings.append(unmasked_patches.detach().cpu().numpy())\n",
    "\n",
    "embeddings = numpy.vstack(embeddings)\n",
    "\n",
    "embeddings_mean = embeddings[:, :-2, :].mean(axis=1)\n",
    "\n",
    "print(f\"Average embeddings have shape {embeddings_mean.shape}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72db5745-21c6-4b8e-b8f7-cb48c0f9c9ef",
   "metadata": {},
   "source": [
    "## Analyze embeddings\n",
    "Now we can make a simple analysis of the embeddings. We reduce all the\n",
    "embeddings to a single number using Principle Component Analysis. Then\n",
    "we can plot the principal components. The effect of the fire on the\n",
    "embeddings is clearly visible. We use the following color code in the graph:\n",
    "\n",
    "| Color | Interpretation |\n",
    "|---|---|\n",
    "| Green | Cloudy Images |\n",
    "| Blue | Before the fire |\n",
    "| Red | After the fire |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "88f3b2dc-8f2a-447b-a6af-b04e0d1ff61c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7f9948d29890>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAHDCAYAAAApyGCxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA7oUlEQVR4nO3deXgUZbr+8bsJJAE0jYqEAJHFg+wjEIQEjYJKZFVEhuCSGRVxOI6jyM8NmXFcJy5HB1xwGR05OgyikqgzLIoKgkMQgQR3REDDkoAgdANKCJ3n90dOWpokkEA63V39/VxXXdpvv115n4Ki71S9VeUyMxMAAICDNAj1AAAAAOoaAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADhOw1APIBTKysq0detWnXjiiXK5XKEeDgAAqAEz0549e9SqVSs1aHDkYzRRGXC2bt2q5OTkUA8DAAAcg02bNqlNmzZH7BOVAefEE0+UVL6BEhISQjwaAABQE16vV8nJyf7v8SOJyoBTcVoqISGBgAMAQISpyfQSJhkDAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHico7GSO8HSj1afrcpVq/rUinJybphmHpim0UE+phAQAiCAEHYeX2l3L0+Jc3y3fC5vKGYunW/7TRpK7T9Mg1o0I7OABAxOAUFcLG7S/l6NHvR8vXdHNAu6/pFj36/Wjd/lJOiEYGAIg0QQ04S5Ys0YgRI9SqVSu5XC69+eabR+yfk5OjQYMG6dRTT1VCQoLS0tL0zjvvBPSZMWOGXC5XpWX//v1BrATBdqDUp8e/vFmSSYc/Q81lkqTHv5yoA6W+eh8bACDyBDXg7Nu3T2eeeaaeeuqpGvVfsmSJBg0apHnz5mnVqlUaOHCgRowYofz8/IB+CQkJKioqClji4+ODUQLqyfS5S8tPS1X3gFiXyXfCJk2fu7RexwUAiExBnYMzZMgQDRkypMb9p06dGvD6L3/5i9566y3961//Uq9evfztLpdLLVu2rKthIgys31ZUp/0AANEtrCcZl5WVac+ePTr55JMD2vfu3au2bdvK5/OpZ8+euv/++wMC0OFKSkpUUlLif+31eoMyXq7+OXanJyZJxTXsBwDAUYT1JOPHHntM+/bt05gxY/xtnTt31owZM/T2229r1qxZio+P19lnn61169ZVu57s7Gy53W7/kpycXOdjvf2lHDW5q51uWTNQTxVfoVvWDFSTu9oxMbaGbhiWrpi9bSSr5hyVuRSzN1k3DEuv34EBACKSy8ysXn6Qy6Xc3FyNHDmyRv1nzZql6667Tm+99ZYuvPDCavuVlZWpd+/eOvfcc/XEE09U2aeqIzjJycnyeDxKSEioVR1Vqbj6p9IE2f/7sr6t7Rtc4lwDv2xH+ScWS2I7AgAklX9/u93uGn1/h+URnNmzZ2vcuHF67bXXjhhuJKlBgwY666yzjngEJy4uTgkJCQFLXeHqn7rzyDWjdFvbNxSzr3VAe8y+NoQbAECthN0cnFmzZunaa6/VrFmzNGzYsKP2NzMVFBSoR48e9TC6yvxX/1TnkKt/Jo4cUG/jilSPXDNKD5RewlwmAMBxCWrA2bt3r7799lv/640bN6qgoEAnn3yyTjvtNE2ePFlbtmzRyy+/LKk83PzmN7/RtGnTlJqaquLi8lmnjRs3ltvtliTde++9Sk1NVceOHeX1evXEE0+ooKBATz/9dDBLqRZX/9S92EYxhEEAwHEJ6imqlStXqlevXv4rnCZNmqRevXrp7rvvliQVFRWpsLDQ3/+5557TwYMH9fvf/15JSUn+5eabb/b32b17t66//np16dJFGRkZ2rJli5YsWaK+ffsGs5Rq1fSqHq7+AQCg/tTbJONwUptJSkdzoNSnJne1k6/plsCJsRXMpZh9bfTTXzZymgUAgOMQ8ZOMI0lsoxhN6jqt/MXhlzj/3+tJXacSbgAAqEcEnDrA1T8AAIQXTlHV8SXjXP0DAEBw1Ob7O+wuE49kXP0DAEB44BQVAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwnKAGnCVLlmjEiBFq1aqVXC6X3nzzzaN+5sMPP1RKSori4+PVoUMHPfvss5X6zJkzR127dlVcXJy6du2q3NzcIIweAABEqqAGnH379unMM8/UU089VaP+Gzdu1NChQ5Wenq78/HzddddduummmzRnzhx/n7y8PGVmZiorK0tr1qxRVlaWxowZo48//jhYZQAAgAjjMjOrlx/kcik3N1cjR46sts8dd9yht99+W1999ZW/bcKECVqzZo3y8vIkSZmZmfJ6vZo/f76/z+DBg3XSSSdp1qxZNRqL1+uV2+2Wx+NRQkLCsRUEAADqVW2+v8NqDk5eXp4yMjIC2i666CKtXLlSpaWlR+yzbNmyatdbUlIir9cbsAAAAOcKq4BTXFysxMTEgLbExEQdPHhQO3bsOGKf4uLiatebnZ0tt9vtX5KTk+t+8AAAIGyEVcCRyk9lHariDNqh7VX1ObztUJMnT5bH4/EvmzZtqsMRAwCAcNMw1AM4VMuWLSsdidm+fbsaNmyoU0455Yh9Dj+qc6i4uDjFxcXV/YABAEBYCqsjOGlpaVq4cGFA27vvvqs+ffqoUaNGR+zTv3//ehsnAAAIb0E9grN37159++23/tcbN25UQUGBTj75ZJ122mmaPHmytmzZopdffllS+RVTTz31lCZNmqTx48crLy9PL774YsDVUTfffLPOPfdcPfzww7rkkkv01ltv6b333tNHH30UzFIAAEAECeoRnJUrV6pXr17q1auXJGnSpEnq1auX7r77bklSUVGRCgsL/f3bt2+vefPmafHixerZs6fuv/9+PfHEE7rsssv8ffr3769XX31VL730kn71q19pxowZmj17tvr16xfMUgAAQASpt/vghBPugwMAQOSJ2PvgAAAA1AUCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcBwCDgAAcJx6CTjTp09X+/btFR8fr5SUFC1durTavldffbVcLlelpVu3bv4+M2bMqLLP/v3766McAAAQ5oIecGbPnq2JEydqypQpys/PV3p6uoYMGaLCwsIq+0+bNk1FRUX+ZdOmTTr55JP161//OqBfQkJCQL+ioiLFx8cHuxwAABABGgb7Bzz++OMaN26crrvuOknS1KlT9c477+iZZ55RdnZ2pf5ut1tut9v/+s0339SuXbt0zTXXBPRzuVxq2bJljcZQUlKikpIS/2uv13sspQAAgAgR1CM4Bw4c0KpVq5SRkRHQnpGRoWXLltVoHS+++KIuvPBCtW3bNqB97969atu2rdq0aaPhw4crPz+/2nVkZ2f7g5Pb7VZycnLtiwEAABEjqAFnx44d8vl8SkxMDGhPTExUcXHxUT9fVFSk+fPn+4/+VOjcubNmzJiht99+W7NmzVJ8fLzOPvtsrVu3rsr1TJ48WR6Px79s2rTp2IsCAABhL+inqKTy00mHMrNKbVWZMWOGmjVrppEjRwa0p6amKjU11f/67LPPVu/evfXkk0/qiSeeqLSeuLg4xcXFHdvgAQBAxAnqEZzmzZsrJiam0tGa7du3Vzqqczgz09///ndlZWUpNjb2iH0bNGigs846q9ojOAAAILoENeDExsYqJSVFCxcuDGhfuHCh+vfvf8TPfvjhh/r22281bty4o/4cM1NBQYGSkpKOa7wAAMAZgn6KatKkScrKylKfPn2Ulpam559/XoWFhZowYYKk8vkxW7Zs0csvvxzwuRdffFH9+vVT9+7dK63z3nvvVWpqqjp27Civ16snnnhCBQUFevrpp4NdDgAAiABBDziZmZnauXOn7rvvPhUVFal79+6aN2+e/6qooqKiSvfE8Xg8mjNnjqZNm1blOnfv3q3rr79excXFcrvd6tWrl5YsWaK+ffsGuxwAABABXGZmoR5EffN6vXK73fJ4PEpISAj1cAAAQA3U5vubZ1EBAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHIeAAAADHqZeAM336dLVv317x8fFKSUnR0qVLq+27ePFiuVyuSsvXX38d0G/OnDnq2rWr4uLi1LVrV+Xm5ga7DAAAECGCHnBmz56tiRMnasqUKcrPz1d6erqGDBmiwsLCI35u7dq1Kioq8i8dO3b0v5eXl6fMzExlZWVpzZo1ysrK0pgxY/Txxx8HuxwAABABXGZmwfwB/fr1U+/evfXMM8/427p06aKRI0cqOzu7Uv/Fixdr4MCB2rVrl5o1a1blOjMzM+X1ejV//nx/2+DBg3XSSSdp1qxZRx2T1+uV2+2Wx+NRQkJC7YsCAAD1rjbf30E9gnPgwAGtWrVKGRkZAe0ZGRlatmzZET/bq1cvJSUl6YILLtCiRYsC3svLy6u0zosuuqjadZaUlMjr9QYsAADAuYIacHbs2CGfz6fExMSA9sTERBUXF1f5maSkJD3//POaM2eOcnJy1KlTJ11wwQVasmSJv09xcXGt1pmdnS232+1fkpOTj7MyAAAQzhrWxw9xuVwBr82sUluFTp06qVOnTv7XaWlp2rRpk/7nf/5H55577jGtc/LkyZo0aZL/tdfrJeQAAOBgQT2C07x5c8XExFQ6srJ9+/ZKR2COJDU1VevWrfO/btmyZa3WGRcXp4SEhIAFAAA4V1ADTmxsrFJSUrRw4cKA9oULF6p///41Xk9+fr6SkpL8r9PS0iqt8913363VOgEAgHMF/RTVpEmTlJWVpT59+igtLU3PP/+8CgsLNWHCBEnlp4+2bNmil19+WZI0depUtWvXTt26ddOBAwf0j3/8Q3PmzNGcOXP867z55pt17rnn6uGHH9Yll1yit956S++9954++uijYJcDAAAiQNADTmZmpnbu3Kn77rtPRUVF6t69u+bNm6e2bdtKkoqKigLuiXPgwAHdeuut2rJlixo3bqxu3bpp7ty5Gjp0qL9P//799eqrr+qPf/yj/vSnP+n000/X7Nmz1a9fv2CXAwAAIkDQ74MTjrgPDgAAkSds7oMDAAAQCgQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOAQcAADgOPUScKZPn6727dsrPj5eKSkpWrp0abV9c3JyNGjQIJ166qlKSEhQWlqa3nnnnYA+M2bMkMvlqrTs378/2KUAAIAIEPSAM3v2bE2cOFFTpkxRfn6+0tPTNWTIEBUWFlbZf8mSJRo0aJDmzZunVatWaeDAgRoxYoTy8/MD+iUkJKioqChgiY+PD3Y5AAAgArjMzIL5A/r166fevXvrmWee8bd16dJFI0eOVHZ2do3W0a1bN2VmZuruu++WVH4EZ+LEidq9e/cxjcnr9crtdsvj8SghIeGY1gEAAOpXbb6/g3oE58CBA1q1apUyMjIC2jMyMrRs2bIaraOsrEx79uzRySefHNC+d+9etW3bVm3atNHw4cMrHeE5VElJibxeb8ACAACcK6gBZ8eOHfL5fEpMTAxoT0xMVHFxcY3W8dhjj2nfvn0aM2aMv61z586aMWOG3n77bc2aNUvx8fE6++yztW7duirXkZ2dLbfb7V+Sk5OPvSgAABD26mWSscvlCnhtZpXaqjJr1izdc889mj17tlq0aOFvT01N1VVXXaUzzzxT6enpeu2113TGGWfoySefrHI9kydPlsfj8S+bNm06voIAAEBYaxjMlTdv3lwxMTGVjtZs37690lGdw82ePVvjxo3T66+/rgsvvPCIfRs0aKCzzjqr2iM4cXFxiouLq93gAQBAxArqEZzY2FilpKRo4cKFAe0LFy5U//79q/3crFmzdPXVV+uf//ynhg0bdtSfY2YqKChQUlLScY8ZAABEvqAewZGkSZMmKSsrS3369FFaWpqef/55FRYWasKECZLKTx9t2bJFL7/8sqTycPOb3/xG06ZNU2pqqv/oT+PGjeV2uyVJ9957r1JTU9WxY0d5vV498cQTKigo0NNPPx3scgAAQAQIesDJzMzUzp07dd9996moqEjdu3fXvHnz1LZtW0lSUVFRwD1xnnvuOR08eFC///3v9fvf/97f/tvf/lYzZsyQJO3evVvXX3+9iouL5Xa71atXLy1ZskR9+/YNdjkAACACBP0+OOGI++AAABB5wuY+OAAAAKFAwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5DwAEAAI5TLwFn+vTpat++veLj45WSkqKlS5cesf+HH36olJQUxcfHq0OHDnr22Wcr9ZkzZ466du2quLg4de3aVbm5ucEavqP4fNLixdKsWeX/9flCPSIAAOpe0APO7NmzNXHiRE2ZMkX5+flKT0/XkCFDVFhYWGX/jRs3aujQoUpPT1d+fr7uuusu3XTTTZozZ46/T15enjIzM5WVlaU1a9YoKytLY8aM0ccffxzsciJaTo7Urp00cKB0xRXl/23XrrwdAAAncZmZBfMH9OvXT71799Yzzzzjb+vSpYtGjhyp7OzsSv3vuOMOvf322/rqq6/8bRMmTNCaNWuUl5cnScrMzJTX69X8+fP9fQYPHqyTTjpJs2bNqrTOkpISlZSU+F97vV4lJyfL4/EoISGhTuoMdzk50ujR0uF/2i5X+X/feEMaNar+xwUAQE15vV653e4afX8H9QjOgQMHtGrVKmVkZAS0Z2RkaNmyZVV+Ji8vr1L/iy66SCtXrlRpaekR+1S3zuzsbLndbv+SnJx8rCVFJJ9PuvnmyuFG+qVt4kROVwEAnCOoAWfHjh3y+XxKTEwMaE9MTFRxcXGVnykuLq6y/8GDB7Vjx44j9qlunZMnT5bH4/EvmzZtOtaSItLSpdLmzdW/byZt2lTeDwAAJ2hYHz/EVXEe5P+YWaW2o/U/vL0264yLi1NcXFytxuwkRUV12w8AgHAX1IDTvHlzxcTEVDqysn379kpHYCq0bNmyyv4NGzbUKaeccsQ+1a0z2iUl1W0/lJ/OW7q0PBQmJUnp6VJMTKhHFXnYjgCCJainqGJjY5WSkqKFCxcGtC9cuFD9+/ev8jNpaWmV+r/77rvq06ePGjVqdMQ+1a0z2qWnS23a/DKh+HAul5ScXN4PR8fVaHWD7QggqCzIXn31VWvUqJG9+OKL9uWXX9rEiROtadOm9t1335mZ2Z133mlZWVn+/hs2bLAmTZrYLbfcYl9++aW9+OKL1qhRI3vjjTf8ff7zn/9YTEyMPfTQQ/bVV1/ZQw89ZA0bNrTly5fXaEwej8ckmcfjqdtiw9icOWYuV/lSPuumfKlomzMn1COMDBXb8dBtGK7b8eBBs0WLzP75z/L/HjwY6hH9IpK2I4DwUZvv76AHHDOzp59+2tq2bWuxsbHWu3dv+/DDD/3v/fa3v7XzzjsvoP/ixYutV69eFhsba+3atbNnnnmm0jpff/1169SpkzVq1Mg6d+5sc2rxL2I0Bhyz8i+NNm0Cv1CSk/kyqamDBytvv8O/nJOTwyNIVPVn3aZNePxZR9J2BBBeavP9HfT74ISj2lxH7zTMeTh2ixeXn0Y5mkWLpAEDgj2a6oX7PY8iZTsCCD+1+f6ul6uoED5iYvjSOFaRcDXa0e555HKV3/PokktCF2wjYTsCiHw8bBOooUi4Gi0S7nkUCdsRQOQj4AA1FAlXo0XC0ZFI2I4AIh8BB6ihmBhp2rTy/z/8y7ni9dSpoZ3TFAlHRyJhOwKIfAQcoBZGjSqfpNu6dWB7mzahn7wrRc7RkXDfjgAiH1dRRdlVVKgb4Xw1WsVVVFLgZONwuYrqUOG8HQGEn9p8fxNwCDhwoJyc8qupDp1wnJxcfuonXMINANQWl4kDUW7UqPJLwTk6AiBaEXAAh+KeRwCiGZOMAQCA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4xBwAACA4/A0cQA4Ap9PWrpUKiqSkpKk9PTyJ7UDCG8EHACoRk6OdPPN0ubNv7S1aSNNmyaNGhW6cQE4Ok5RAUAVcnKk0aMDw40kbdlS3p6TE5pxAagZAg4AHMbnKz9yY1b5vYq2iRPL+wEITwQcADjM0qWVj9wcykzatKm8H4DwRMABgMMUFdVtPwD1j4ADAIdJSqrbfgDqHwEHAA6Tnl5+tZTLVfX7LpeUnFzeD0B4IuAAwGFiYsovBZcqh5yK11Oncj8cIJwRcACgCqNGSW+8IbVuHdjepk15O/fBAcIbN/oDgGqMGiVdcgl3MgYiEQEHAI4gJkYaMCDUowBQW5yiAgAAjkPAAQAAjhPUgLNr1y5lZWXJ7XbL7XYrKytLu3fvrrZ/aWmp7rjjDvXo0UNNmzZVq1at9Jvf/EZbt24N6DdgwAC5XK6AZezYscEsBQAARJCgBpwrrrhCBQUFWrBggRYsWKCCggJlZWVV2/+nn37S6tWr9ac//UmrV69WTk6OvvnmG1188cWV+o4fP15FRUX+5bnnngtmKQAAIIIEbZLxV199pQULFmj58uXq16+fJOlvf/ub0tLStHbtWnXq1KnSZ9xutxYuXBjQ9uSTT6pv374qLCzUaaed5m9v0qSJWrZsGazhAwCACBa0Izh5eXlyu93+cCNJqampcrvdWrZsWY3X4/F45HK51KxZs4D2mTNnqnnz5urWrZtuvfVW7dmzp9p1lJSUyOv1BiwAAMC5gnYEp7i4WC1atKjU3qJFCxUXF9doHfv379edd96pK664QgkJCf72K6+8Uu3bt1fLli31+eefa/LkyVqzZk2loz8VsrOzde+99x5bIQAAIOLU+gjOPffcU2mC7+HLypUrJUmuKh7kYmZVth+utLRUY8eOVVlZmaZPnx7w3vjx43XhhReqe/fuGjt2rN544w299957Wr16dZXrmjx5sjwej3/ZtGlTbcsGAAARpNZHcG688cajXrHUrl07ffrpp9q2bVul93744QclJiYe8fOlpaUaM2aMNm7cqA8++CDg6E1VevfurUaNGmndunXq3bt3pffj4uIUFxd3xHUAAADnqHXAad68uZo3b37UfmlpafJ4PFqxYoX69u0rSfr444/l8XjUv3//aj9XEW7WrVunRYsW6ZRTTjnqz/riiy9UWlqqpKSkmhcCAAAcK2iTjLt06aLBgwdr/PjxWr58uZYvX67x48dr+PDhAVdQde7cWbm5uZKkgwcPavTo0Vq5cqVmzpwpn8+n4uJiFRcX68CBA5Kk9evX67777tPKlSv13Xffad68efr1r3+tXr166eyzzw5WOQAAIIIE9T44M2fOVI8ePZSRkaGMjAz96le/0iuvvBLQZ+3atfJ4PJKkzZs36+2339bmzZvVs2dPJSUl+ZeKK69iY2P1/vvv66KLLlKnTp100003KSMjQ++9955ieAIeAACQ5DIzC/Ug6pvX65Xb7ZbH4znq/B4AABAeavP9zbOoAACA4xBwAACA4wTtRn8AgPrhO+DTZ9OX6qf1RWpyepJ63JCumFjmJCK6EXAAIIItvz1Hpz1+s3r6Nvvbtt7aRoWTpin1kVEhHBkQWpyiAoAItfz2HPV9dLRaHhJuJKmlb4v6Pjpay2/PCdHIgNAj4ABABPId8Om0x2+WZJX+IW+g8otjkx+fKN8BX72PDQgHBBwAiECfTV+qVr7N1f4j3kCm1r5N+mz60nodFxAuCDgAEIF+Wl9Up/0ApyHgAEAEanJ6zZ69V9N+gNMQcAAgAvW4IV1bY9qoTK4q3y+TS1tiktXjhvR6HhkQHgg4ABCBYmJjVDhpmiRVCjkVrzdNmsr9cBC1CDgAEKFSHxmlFbe9oeKY1gHtRTFttOK2N7gPDqIaD9vkYZsAIhx3Mka0qM33N3cyBoAIFxMbo54TB4R6GEBY4RQVAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHAIOAABwHO6DAwAA6ky43HiSgAMAAOrE8ttzdNrjN6unb7O/beutbVQ4aVq9PzqEU1QAAOC4Lb89R30fHa2Wh4QbSWrp26K+j47W8ttz6nU8BBwAAHBcfAd8Ou3xmyVZpWDRQOWPvEx+fKJ8B3z1NiYCDgAg6HwHfCqYuljL/jBLBVMX1+sXHYLvs+lL1cq3udpQ0UCm1r5N+mz60nobE3NwAABBFU7zMhAcP60vqtN+dYEjOACAoAm3eRkIjianJ9Vpv7rgMjOrt58WJrxer9xutzwejxISEkI9HABwJN8Bn7Y1aaeW1Zy6KJNLRTFt1PKnjSG5jBh155c/6y3+OTeHqqs/69p8f3MEBwAQFOE4LwPBERMbo8JJ0ySVh5lDVbzeNGlqvQZZAg4AICjCcV4Ggif1kVFacdsbKo5pHdBeFNNGK257o97nWzHJGAAQFOE4LwPBlfrIKPkeuEQFh93JuHUITkEyB4c5OAAQFPU1L6OuhMsjBlA95uAAAEIuHOdlVGf57Tna1qSdet4yUP2fukI9bxmobU3acZVXBCPgAACCJtzmZVSFS9mdiVNUnKICgKAL19M/XMoeWWrz/c0kYwBA0MXExqjnxAGhHkYln01fGnCH5cNVXMpeMH1pWIw/XINiOArqKapdu3YpKytLbrdbbrdbWVlZ2r179xE/c/XVV8vlcgUsqampAX1KSkr0hz/8Qc2bN1fTpk118cUXa/Pm6v+CAgBQlUi6lJ15QrUT1IBzxRVXqKCgQAsWLNCCBQtUUFCgrKyso35u8ODBKioq8i/z5s0LeH/ixInKzc3Vq6++qo8++kh79+7V8OHD5fPx8DYAQM1FyqXszBOqvaDNwfnqq6/UtWtXLV++XP369ZMkLV++XGlpafr666/VqVOnKj939dVXa/fu3XrzzTerfN/j8ejUU0/VK6+8oszMTEnS1q1blZycrHnz5umiiy466tiYgwMAkCLjUnbmCf0iLC4Tz8vLk9vt9ocbSUpNTZXb7dayZcuO+NnFixerRYsWOuOMMzR+/Hht377d/96qVatUWlqqjIwMf1urVq3UvXv3atdbUlIir9cbsAAAEAmXsvPIi2MTtIBTXFysFi1aVGpv0aKFiouLq/3ckCFDNHPmTH3wwQd67LHH9Mknn+j8889XSUmJf72xsbE66aSTAj6XmJhY7Xqzs7P984DcbreSk5OPozIAgJOE+6XskTRPKJzUOuDcc889lSYBH76sXLlSkuRyuSp93syqbK+QmZmpYcOGqXv37hoxYoTmz5+vb775RnPnzj3iuI603smTJ8vj8fiXTZs21aJiAIDTpT4ySok/faeCvy7Sshv/qYK/LlLLnzaGPNxIkTNPKNzU+jLxG2+8UWPHjj1in3bt2unTTz/Vtm3bKr33ww8/KDExscY/LykpSW3bttW6deskSS1bttSBAwe0a9eugKM427dvV//+/atcR1xcnOLi4mr8MwEA0SdcL2XvcUO6tt7a5qjzhHrckB6C0YWvWgec5s2bq3nz5kftl5aWJo/HoxUrVqhv376SpI8//lgej6faIFKVnTt3atOmTUpKKk+mKSkpatSokRYuXKgxY8ZIkoqKivT555/rkUceqW05AACEtYp5Qi0fHa0yuQJCzqHzhELxQMtwFrQ5OF26dNHgwYM1fvx4LV++XMuXL9f48eM1fPjwgCuoOnfurNzcXEnS3r17deuttyovL0/fffedFi9erBEjRqh58+a69NJLJUlut1vjxo3T//t//0/vv/++8vPzddVVV6lHjx668MILg1UOAAAhE+7zhMJRUO9kPHPmTN10003+K54uvvhiPfXUUwF91q5dK4/HI0mKiYnRZ599ppdfflm7d+9WUlKSBg4cqNmzZ+vEE0/0f+avf/2rGjZsqDFjxujnn3/WBRdcoBkzZigmhvQKAHCm1EdGyffAJSo47E7GHLmpGs+i4j44AABEhLC4Dw4AAECoEHAAAIDjEHAAAIDjEHAAAIDjEHAAAIDjEHAAAIDjEHAAAIDjEHAAAIDjEHAAAIDjBPVRDeGq4ubNXq83xCMBAAA1VfG9XZOHMERlwNmzZ48kKTk5OcQjAQAAtbVnzx653e4j9onKZ1GVlZVp69atOvHEE+VyuUI9nOPm9XqVnJysTZs2OfbZWtFQoxQ9dQZbNGxHanSOaKmzLpiZ9uzZo1atWqlBgyPPsonKIzgNGjRQmzZtQj2MOpeQkOD4nSMaapSip85gi4btSI3OES11Hq+jHbmpwCRjAADgOAQcAADgOAQcB4iLi9Of//xnxcXFhXooQRMNNUrRU2ewRcN2pEbniJY661tUTjIGAADOxhEcAADgOAQcAADgOAQcAADgOAQcAADgOAScOrRv374aPR8jkkVDjVL01Bls0bAdqdE5oqXOaEHAqSM+n0/jxo3T0KFD9fzzz4d6OEERDTVK0VNnsEXDdqRG54iWOqMJAaeO7N69Wz179tSgQYP00EMP6ZxzztG6detCPaw6FQ01StFTZ7BFw3akRueIljqjCffBCYK9e/dq9OjRSk1N1T333BPq4QRFNNQoRU+dwRYN25EanSNa6nQ6juDUkYqcePDgQZ1wwglKT0/XSy+9pL179wa8H8mioUYpeuoMtmjYjtTojBql6KkzmhBw6oCZyeVySZIaNix/QPvy5cvVqVMnnXDCCZLkf3/Hjh2hGeRxioYapeipM9iiYTtSozNqlKKnzmjTMNQDcIKKv/hbt27V+++/r1dffVUfffSR5s+fr7KyMjVo0EBz587Va6+9pvXr10uSHn/8cfXt2zeUw66VaKhRip46gy0atiM1OqNGKXrqjDbMwTkOXq9Xn376qRYsWKBPPvlEX3zxhWJjYzV69GiNGDFC6enpkqSFCxdqzJgxSklJUWZmpr7++mu98sorevXVV3X++eeHuIoji4YapeipM9iiYTtSozNqlKKnzqhlOCYlJSXWrFkza9asmY0dO9Yee+wx+/DDD83n8wX027BhgyUnJ9u1115ru3fvNjOzgwcP2vDhw+3WW281M7OysrJ6H39NREONZtFTZ7BFw3akxl9Eco1m0VNnNCPgHKMffvjB+vfvby6Xy1588cWA9w7dQW6++WZr3bq1eTweM/tlRxg7dqwNHjw44HPbtm0L8qhrJxpqNIueOoMtGrYjNTqjRrPoqTOaEXCO0+uvv26nnnqqde7c2ebOnRvw3rZt26xJkyb+nefgwYNmZlZUVGSpqak2efJk/460dOlSGz9+vI0fP96+/PLL+i3iKKKhRrPoqTPYomE7UqMzajSLnjqjEVdRHafRo0dr+/btuuyyy3TppZcqJydHZWVlkqR//OMfatGihUaOHCkzU0xMjCRp7ty5MjN17dpVDRqU/xGceOKJ6tSpkwoLC9WtWze9+eaboSqpkmioUYqeOoMtGrYjNTqjRil66oxKoUhVTrVnzx4rLCw0s/LDmPfee6+lpqaa1+v19/niiy/skksusWHDhtnevXv9fc3Kzwnfc8891rhx47A91BkNNZpFT53BFg3bkRrLRXqNZtFTZ7TgCE4dOuGEE5ScnCyp/LLD7t2768cff1RJSYm/z0MPPaSdO3dq/Pjxatq0qXw+n1wul8xMGzZs0JNPPqm7775bLVq08P8WEU6ioUYpeuoMtmjYjtRYLtJrlKKnzqgRwnDleIWFhZaammo9e/a0yZMn21lnnWVJSUn2wgsv+PtUJP/9+/fbuHHjrH379gHrqDjn6/V6bePGjWH3W0E01GgWPXUGWzRsR2osF+k1mkVPnU5FwKkHU6ZMseHDh9stt9xia9as8bcfOlN/0aJF1rBhQ8vNzfW/V7HjbN261Xr06GG9e/e20047zW688cawuywxGmo0i546gy0atiM1lov0Gs2ip06nIeDUk9LS0mrf2717tw0cONDOP//8gPaKHWDUqFGWkpJib7/9tuXl5dnZZ59t559/vv3www9BHXNtRUONZtFTZ7BFw3akRmfUaBY9dToJASdEVqxYYbm5uebxeOyVV16xBg0a2BdffGFm5cm/4jeDsrIyGzNmjN15553+z3755Zd2xhln2Lx580Iy9pqKhhrNoqfOYIuG7UiNzqjRLHrqjGQ8iypEDh48qHHjxqlZs2bau3evJkyYoK5du/qfe1Jh+fLl6ty5sxYsWCCPxyO3262WLVtq3bp12r9/fwgrOLpoqFGKnjqDLRq2IzU6o0YpeuqMZFxFFSJpaWnavn27Lr30Uv3www/67rvv9PXXX/t3jN27d+umm27S6NGjtWbNGjVu3FgtWrTQJZdcogEDBqhnz55h/6C3aKhRip46gy0atiM1OqNGKXrqjGihPoQEs40bN1pKSor17dvXNmzYYGblhz8bNmxoixYt8vebPn26nXnmmbZq1aqA+zJEgmio0Sx66gy2aNiO1LjI3y+SazSLnjojDQEnjKxfv95/46i//e1v1rhx44AbSf3888+WmJhor7/+esDnDh48WGlGvtfrtQ0bNti3335bP4OvoWio0Sx66gy2aNiO1OiMGs2ip85IwSmqMNKhQwc1bdpUknTeeefpjDPO8N/u2+VyafPmzYqPj1eTJk0kSX/961/l8XgUExMjl8slSfL5fFq+fLn69euniy++WGeddZbGjh2rffv2haSmw0VDjVL01Bls0bAdqdEZNUrRU2fECHXCQtXKysrsvvvus2bNmtn48ePtoYcesjZt2tjw4cNt5cqVVlpaasOHD7dmzZrZAw884P/czJkzLSUlxQYNGmT/+c9/bNWqVXbeeefZxIkTw+6+C9FQo1n01Bls0bAdqdEZNZpFT53hjIAT5j777DMbOXKkjRo1yn77299WOlz55ptvWp8+fWzWrFn2/fff23nnnWdjx44NuL/Cgw8+aF26dLE9e/bU9/BrJBpqNIueOoOtLrdjxemDcEONzqjRLHrqDEcEnAixb9++at/bv3+/mZndfffdlpaWZv/+97/N7Je7bL722mvWsWNH/+S3/Px8u+666+zzzz8P8qhrJxpqNIueOoOtLrbj+vXrAz73888/B2m0x4YanVGjWfTUGU6YgxMhKs7ZViUuLk7r1q1TXl6ezjrrLA0ZMkSSZGaSpHfffVdNmzZV+/btJUmdO3dWkyZN1KNHD1199dVhcy+GaKhRip46g+14t2OTJk3UoUMHSdJnn32ma6+9VldccYVuuummgIcrhhI1OqNGKXrqDCshDFeoQ1988YW1aNHC5s+fb2a/PODt22+/tdjYWHv55ZcD2s3KD50OGDDA4uLibPr06fU/6FqKhhrNoqfOYDvadvzf//1fMzNbu3atXXTRRXb66afbk08+aeeff761bt3a/7lwRo3OqNEseuqsTwQch1i8eLHFxMRYUVGRv83n89nw4cOtb9++tnPnzoD+h3453nbbbeZyuezKK6+st/Eei2io0Sx66gy2o23HHTt2mFn5F4bL5bJPPvnE3+93v/udZWZm1vuYa4sanVGjWfTUWZ8IOA7h8XgsLS3Npk6damblqf+BBx6wRo0a2X/+85+AvhXndfft22evvPKKnXLKKXbllVfahx9+GPB+uImGGs2ip85gq247NmzY0D766CN/v++++84uv/xye+edd/xtU6dOtQ4dOoT9PUio0Rk1mkVPnfWJgOMgU6dOtbi4OMvIyLDWrVtbt27d7Jlnnqm2/+9+9zs744wz7KqrrqrHUR6faKjRLHrqDLaqtuOTTz5pZmZfffWV3XXXXZaUlGQpKSnWoEEDu/HGG23atGl29tln23nnnRfawdcQNTqjRrPoqbO+EHAcZvPmzfbAAw/Y3//+d/vmm2/87WVlZf7f5levXm1XX321NWrUyF544QX/pYeR8tt+NNRoFj11Btuh23Ht2rX+9htvvNHOO+88mzlzpu3cudPef/99O/300+3iiy+2F154wd/30FOA4YoanVGjWfTUWR8IOFGg4i98WVmZffrpp9ahQwc766yzLDc3t8q+FTeTOvymUuF875VjrbHCvn37bNmyZbZu3br6GO4xO9Y6KwLP+vXr7aWXXrJPP/203sYcbsrKymzPnj02cuRIGzNmTMB7f/7zn23YsGEBbZG4T1DjkWusEAn7/fHUGe37PZeJR4GYmBhJ0mOPPaarrrpKycnJeuWVVzRy5Eh/nx9//DHgluGlpaX+W4d/+OGHysrK0qhRo5SZmakffvghFGUc0bHU6PP5/O998803+vvf/66RI0fq8ssv1+7du+u5gpo51jornnA8depUzZw5U0OGDNFll12mnTt3hqKMkHK5XDrhhBOUmJion376yX8priT17dtXq1at0q5du+T1eiN2n6DG6muMtP3+eOqM+v0+tPkK9SU/P99iY2NtwoQJtmvXrkrvL1iwwGJjY23KlCkB7c8++6x17tzZUlJS7Pnnn7eLL77Y2rVrZ19//XU9jbzmjrVGs/LfelavXm2vvvqquVwue/DBB62kpKQeRl17x1Onz+ezoqIiW7Nmjf3617+2c845x7Zu3VoPow4/H3zwgbVr186uueYa++KLL2zWrFnWvXt3GzlypJmZzZ07N+L3CWqsvkazyNrvj6fOaN3vCThR5IcffjCv12tmVc/R+Ne//mVdunSxAQMG2I8//mgbN260uLg4u/POO/2XKHq9XktNTfVPfAs3Na1x4MCBtm3btkrv33jjjdapUyf77LPPgj7W41GbOouLi6vst3btWjvjjDNszpw5wR9wmFq3bp1dcMEF9l//9V/WqVMnS01NDbgM3wn7BDU6Z7+vTZ3s9wQcHGb//v2Wl5dnZWVlNnr0aOvbt69t3749oE+vXr3stttuC9EIj9/+/fttxYoV5vF4zOyXeS3r1683l8tl06dPt9LS0oD3du3aZZs3bw7NgI9RRZ0VQcis8gTE5s2b27Rp0/yvPR5PwDNwosXnn39uhYWF/m116HY6ln3iwIEDYTd3pa5rDEdHq9Ep+31N6mS/J+DgMBVpf82aNeZyuWzu3LkB73/++efWv39/e/DBB/1tK1eutBkzZtjjjz8esFNFiooJeRdeeKGdc845tmXLlkp9srOz7ZRTTrHLL7/c/w9kJDl8cuWGDRvs1ltvtYSEhIDnWP3ud7+z9PR0mzBhQkT+WQbDsewT7777rv33f/+3/fGPf4yI7ch+z37vxP2eScYIUDEp7bnnnlN6err69OkT8P6yZcv0448/qlevXpKk+fPn68orr1R2drb+/e9/q3Xr1po6dWp9D/uY+Xw+uVwuLViwQB988IHuuOMOtWrVyv+eJH366afKy8tTamqqioqKdOqpp+qee+4J4ahrr7S0VGvXrtULL7ygrKwspaen69///rdeeeUVdevWTZJ08OBBnXnmmbrhhhu0YcMGnXHGGZozZ06IRx56Nd0nKtp9Pp8aN26s5s2ba8GCBXK73VqzZk29j7s22O/Z7x2534c6YSH8lJWV2U033WQXXXRRQPsnn3xiGRkZNmLECH/bhAkTbPDgwbZmzRozM5s9e7b96le/soULF9brmGurrKzM/wRfM7MOHTrYVVdd5Z+0e+hvPrfffrt16dLFf4llbm6unX766XbppZeG7T0nKsa1YsUKu+WWW+yCCy6wZs2a2Zlnnml33HGH5eTk+O+Zc/hveRWmTJliQ4cODdsa61Nt9olDXXPNNdahQ4eImNDJfs9+b+as/Z6AgyplZ2dbamqq/y/5/v377brrrrPu3bsHPAMlNzfXBg8eHPDZjh072l133VWv462tb7/91kaMGGHz58+3Rx991E466SRbsWJFpX4bNmyw22+/3Vq3bm0jRoywDRs2+N+ruHdGON9U79xzzzWXy2VXXXWVfffdd0fsW1ZWFvCP3nvvvWcnnniiffzxx8EeZkQ40j6xevVqMzMrLS31b8NPP/3UGjZsaDNnzgzZmGuL/b4c+70z9nsCDqq0Zs0aa9eunQ0bNsweffRR69Kli/Xu3dv/pOpt27bZww8/bEOHDjW3223jx4+3Xbt2WW5urnXp0sUef/zxan9DCBe33HKLNWrUyFwul1177bUBv9kdavv27fbFF1/YuHHjrE+fPlZQUFDPIz0+jz76qMXExNill15q33//fbX9Dv1y/vnnn+2xxx6zZs2a2b59++prqGGtun3iueeeM7PKvxEPGjTIBgwYUO3fq3DEfv8L9vvI3+8JOKjW9u3b7fLLL7fhw4fbNddcY99++63/KoO0tDRLTU21Bx980N566y3r06ePNW7c2AYNGmSTJk2KmIe+7dixw6655hpzuVz2hz/8odKTug+1ceNG69Gjhz388MP1OMK6UVxcbCNGjDCXy2X33ntvpTudVtixY4fl5ubakCFD7JRTTrHnn38+FMMNW1XtExXb8NBtmZOTYzExMQEPSYwU7PeB2O8jFwEHR/Xzzz8HvM7Pz7dGjRpZXl6ev2379u2Wmppq8+fPr+/h1YmlS5fapZdeah988IGZlR+iPvSqiop/EDIzMy0jIyMkY6wLixYtssmTJ1thYaGZld9PZ8WKFfbHP/7RzjvvPGvdurV1797drr32WluwYEGIRxu+Dt8nDrV//37r2LGjjRs3rh5HVPfY79nvIx0BB7WWn59vrVu3tn/+858B7eecc4795S9/CdGo6tbVV19t9913X0Dbvn37rGPHjvanP/0pRKOqW16v17p06WIul8suv/xyu//++23ZsmWhHlZEWr16tS1evNjMzB5++GE76aSTAuZtOAH7Pft9pOEycdRa586dlZaWptdff11btmyRz+dTbm6uXC6Xfvrpp1AP77jY/z3npVevXnriiSd0xRVXqKCgQB988IGuueYaNWrUSOecc06IR1k3TjzxRA0dOlQxMTGKj4/Xf//3fystLU1S+eWlqLnCwkINHDhQ559/vu6//35NmTJF7du3D/Ww6hT7Pft9xAl1wkJkKi4utqFDh9pJJ51kQ4YMsSZNmtjo0aP9h0Cd4Pvvv7dhw4ZZ586drWPHjpaYmGhz586NqEmjNVFYWGjp6enWpEkTu+eee8J+kmi42rhxo1122WXmcrns+uuvt6KiolAPqc6x3ztHNOz3LrNDHk0K1NJHH32k1atX6/TTT1dKSopatmwZ6iHVuW+++Ubx8fFq2rSpTjnllFAPJ2jeeecdZWZm6g9/+IP+/Oc/q2HDhqEeUkT64IMPdOWVV6pfv3569tlnHblPsN87h5P3ewIOgADbtm1TYmJiqIcR8VavXq2ePXv67xIMhDMn7vcEHAAA4Dj8agEAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAAByHgAMAABzn/wOZYdIZzKl5qgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pca = decomposition.PCA(n_components=1)\n",
    "pca_result = pca.fit_transform(embeddings_mean)\n",
    "\n",
    "plt.xticks(rotation=-30)\n",
    "# All points\n",
    "plt.scatter(stack.time, pca_result, color=\"blue\")\n",
    "\n",
    "# Cloudy images\n",
    "plt.scatter(stack.time[0], pca_result[0], color=\"green\")\n",
    "plt.scatter(stack.time[2], pca_result[2], color=\"green\")\n",
    "\n",
    "# After fire\n",
    "plt.scatter(stack.time[-5:], pca_result[-5:], color=\"red\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a16fbdb8-1c2d-4c84-8526-283fa14faa53",
   "metadata": {},
   "source": [
    "In the plot above, each image embedding is one point. One can clearly \n",
    "distinguish the two cloudy images and the values after the fire are\n",
    "consistently low."
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}