From 22dd6fcf9e4671f1799073d6758ee1b0bdf397e2 Mon Sep 17 00:00:00 2001 From: Daniel Wiesmann Date: Fri, 24 May 2024 10:01:40 +0100 Subject: [PATCH 1/2] Add notebook showing how to run v1 --- docs/clay-v1-wall-to-wall.ipynb | 1648 +++++++++++++++++++++++++++++++ 1 file changed, 1648 insertions(+) create mode 100644 docs/clay-v1-wall-to-wall.ipynb diff --git a/docs/clay-v1-wall-to-wall.ipynb b/docs/clay-v1-wall-to-wall.ipynb new file mode 100644 index 00000000..4dcac9b3 --- /dev/null +++ b/docs/clay-v1-wall-to-wall.ipynb @@ -0,0 +1,1648 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0cc5e729-9116-4ec9-bf1e-8346cbccdf7b", + "metadata": {}, + "source": [ + "## Run Clay v1\n", + "\n", + "This notebook shows how to run Clay v1 wall-to-wall, from downloading imagery\n", + "to training a tiny fine tuning head. This will include the following steps:\n", + "\n", + "1. Set a location and date range of interest\n", + "2. Download Sentinel-2 imagery for this specification\n", + "3. Load the model checkpoint\n", + "4. Prepare data into a format for the model\n", + "5. Run the model on the imagery\n", + "6. Analyise the model embeddings output using PCA\n", + "7. Train a Support Vector Machines fine tuning head" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "add63cd9", + "metadata": {}, + "outputs": [], + "source": [ + "# Add the repo root to the sys path for the model import below\n", + "import sys\n", + "\n", + "sys.path.append(\"..\")" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "6a17b8a8-a9c6-4053-833e-de97287fae49", + "metadata": {}, + "outputs": [], + "source": [ + "import math\n", + "\n", + "import geopandas as gpd\n", + "import numpy as np\n", + "import pandas as pd\n", + "import pystac_client\n", + "import stackstac\n", + "import torch\n", + "import yaml\n", + "from box import Box\n", + "from matplotlib import pyplot as plt\n", + "from rasterio.enums import Resampling\n", + "from shapely import Point\n", + "from sklearn import decomposition, svm\n", + "from stacchip.processors.prechip import normalize_timestamp\n", + "from torchvision.transforms import v2\n", + "\n", + "from src.model import ClayMAEModule" + ] + }, + { + "cell_type": "markdown", + "id": "beac6394-9762-422b-9f5d-82d226018c0c", + "metadata": {}, + "source": [ + "### Specify location and date of interest\n", + "In this example we will use a location in Portugal where a forest fire happened. We will run the model over the time period of the fire and analyse the model embeddings." + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "08d7787d-1506-4de7-89dc-c1054910acf7", + "metadata": {}, + "outputs": [], + "source": [ + "# Point over Monchique Portugal\n", + "lat, lon = 37.30939, -8.57207\n", + "\n", + "# Dates of a large forest fire\n", + "start = \"2018-07-01\"\n", + "end = \"2018-09-01\"" + ] + }, + { + "cell_type": "markdown", + "id": "2bd226c9-003b-4867-a64a-8ae887e7e20a", + "metadata": {}, + "source": [ + "### Get data from STAC catalog\n", + "\n", + "Based on the location and date we can obtain a stack of imagery using stackstac. Let's start with finding the STAC items we want to analyse." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "2e80743c-7c77-459b-9984-f6c26cdff549", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/tam/apps/miniforge3/envs/claymodel/lib/python3.11/site-packages/pystac_client/item_search.py:850: FutureWarning: get_all_items() is deprecated, use item_collection() instead.\n", + " warnings.warn(\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found 12 items\n" + ] + } + ], + "source": [ + "STAC_API = \"https://earth-search.aws.element84.com/v1\"\n", + "COLLECTION = \"sentinel-2-l2a\"\n", + "\n", + "# Search the catalogue\n", + "catalog = pystac_client.Client.open(STAC_API)\n", + "search = catalog.search(\n", + " collections=[COLLECTION],\n", + " datetime=f\"{start}/{end}\",\n", + " bbox=(lon - 1e-5, lat - 1e-5, lon + 1e-5, lat + 1e-5),\n", + " max_items=100,\n", + " query={\"eo:cloud_cover\": {\"lt\": 80}},\n", + ")\n", + "\n", + "all_items = search.get_all_items()\n", + "\n", + "# Reduce to one per date (there might be some duplicates\n", + "# based on the location)\n", + "items = []\n", + "dates = []\n", + "for item in all_items:\n", + " if item.datetime.date() not in dates:\n", + " items.append(item)\n", + " dates.append(item.datetime.date())\n", + "\n", + "print(f\"Found {len(items)} items\")" + ] + }, + { + "cell_type": "markdown", + "id": "5b7c68ae-7c8a-446a-8bc7-5afba70183c2", + "metadata": {}, + "source": [ + "### Create a bounding box around the point of interest\n", + "\n", + "This is needed in the projection of the data so that we can generate image chips of the right size." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "0f3573b5-5a00-47d9-a648-5c4d7cd2c996", + "metadata": {}, + "outputs": [], + "source": [ + "# Extract coordinate system from first item\n", + "epsg = items[0].properties[\"proj:epsg\"]\n", + "\n", + "# Convert point of interest into the image projection\n", + "# (assumes all images are in the same projection)\n", + "poidf = gpd.GeoDataFrame(\n", + " pd.DataFrame(),\n", + " crs=\"EPSG:4326\",\n", + " geometry=[Point(lon, lat)],\n", + ").to_crs(epsg)\n", + "\n", + "coords = poidf.iloc[0].geometry.coords[0]\n", + "\n", + "# Create bounds in projection\n", + "size = 256\n", + "gsd = 10\n", + "bounds = (\n", + " coords[0] - (size * gsd) // 2,\n", + " coords[1] - (size * gsd) // 2,\n", + " coords[0] + (size * gsd) // 2,\n", + " coords[1] + (size * gsd) // 2,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "bbbd3f67-5f2c-46dc-9ee1-2ef1f50fa032", + "metadata": {}, + "source": [ + "### Retrieve the imagery data." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "8b8d3824-e48c-4f9d-9c7b-181c0800f96f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Working with stack of size (12, 4, 256, 256)\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'stackstac-7cbad7c129d678be53c9b6676bee564b' (time: 12,\n",
+       "                                                                band: 4,\n",
+       "                                                                y: 256, x: 256)> Size: 13MB\n",
+       "array([[[[ 9136.,  9232.,  9544., ...,  1258.,  1120.,   930.],\n",
+       "         [ 9616.,  9768.,  9840., ...,  1230.,  1208.,  1030.],\n",
+       "         [ 9992., 10008., 10000., ...,  1418.,  1336.,  1242.],\n",
+       "         ...,\n",
+       "         [  811.,   655.,   688., ...,   385.,   362.,   461.],\n",
+       "         [  798.,   675.,   727., ...,   394.,   415.,   402.],\n",
+       "         [  888.,   673.,   642., ...,   403.,   454.,   393.]],\n",
+       "\n",
+       "        [[ 8656.,  8656.,  8864., ...,  1500.,  1428.,  1220.],\n",
+       "         [ 9016.,  9160.,  9224., ...,  1546.,  1522.,  1360.],\n",
+       "         [ 9248.,  9328.,  9384., ...,  1620.,  1542.,  1482.],\n",
+       "         ...,\n",
+       "         [ 1010.,   831.,   853., ...,   277.,   276.,   336.],\n",
+       "         [ 1016.,   930.,   927., ...,   276.,   317.,   293.],\n",
+       "         [ 1112.,   885.,   827., ...,   299.,   369.,   293.]],\n",
+       "\n",
+       "        [[ 8416.,  8416.,  8640., ...,  1598.,  1466.,  1138.],\n",
+       "         [ 8744.,  8880.,  8928., ...,  1498.,  1522.,  1284.],\n",
+       "         [ 8952.,  8944.,  8960., ...,  1542.,  1478.,  1448.],\n",
+       "         ...,\n",
+       "...\n",
+       "         [  652.,   640.,   638., ...,   590.,   821.,  1008.],\n",
+       "         [  622.,   676.,   630., ...,   606.,  1092.,   726.],\n",
+       "         [  864.,   786.,   569., ...,   766.,  1068.,   630.]],\n",
+       "\n",
+       "        [[  201.,   213.,   195., ...,  1138.,  1058.,   749.],\n",
+       "         [  196.,   198.,   169., ...,   861.,   784.,   768.],\n",
+       "         [  216.,   178.,   191., ...,   870.,   806.,   820.],\n",
+       "         ...,\n",
+       "         [  857.,   838.,   846., ...,   622.,   800.,  1332.],\n",
+       "         [  922.,   848.,   771., ...,   786.,  1046.,   912.],\n",
+       "         [ 1118.,  1010.,   735., ...,   755.,   977.,   686.]],\n",
+       "\n",
+       "        [[ 3264.,  3352.,  3304., ...,  3160.,  3296.,  3376.],\n",
+       "         [ 3356.,  3300.,  3212., ...,  3188.,  3272.,  3064.],\n",
+       "         [ 3288.,  3372.,  3344., ...,  3136.,  3200.,  2932.],\n",
+       "         ...,\n",
+       "         [ 1320.,  1468.,  1298., ...,  2492.,  2556.,  3018.],\n",
+       "         [ 1630.,  1694.,  1250., ...,  2318.,  2684.,  2894.],\n",
+       "         [ 2190.,  2072.,  1288., ...,  2544.,  2942.,  2928.]]]],\n",
+       "      dtype=float32)\n",
+       "Coordinates: (12/53)\n",
+       "  * time                                     (time) datetime64[ns] 96B 2018-0...\n",
+       "    id                                       (time) <U24 1kB 'S2B_29SNB_20180...\n",
+       "  * band                                     (band) <U5 80B 'blue' ... 'nir'\n",
+       "  * x                                        (x) float64 2kB 5.366e+05 ... 5....\n",
+       "  * y                                        (y) float64 2kB 4.131e+06 ... 4....\n",
+       "    platform                                 (time) <U11 528B 'sentinel-2b' ....\n",
+       "    ...                                       ...\n",
+       "    gsd                                      int64 8B 10\n",
+       "    title                                    (band) <U20 320B 'Blue (band 2) ...\n",
+       "    common_name                              (band) <U5 80B 'blue' ... 'nir'\n",
+       "    center_wavelength                        (band) float64 32B 0.49 ... 0.842\n",
+       "    full_width_half_max                      (band) float64 32B 0.098 ... 0.145\n",
+       "    epsg                                     int64 8B 32629\n",
+       "Attributes:\n",
+       "    spec:        RasterSpec(epsg=32629, bounds=(536640.79691545, 4128000.7407...\n",
+       "    crs:         epsg:32629\n",
+       "    transform:   | 10.00, 0.00, 536640.80|\\n| 0.00,-10.00, 4130560.74|\\n| 0.0...\n",
+       "    resolution:  10
" + ], + "text/plain": [ + " Size: 13MB\n", + "array([[[[ 9136., 9232., 9544., ..., 1258., 1120., 930.],\n", + " [ 9616., 9768., 9840., ..., 1230., 1208., 1030.],\n", + " [ 9992., 10008., 10000., ..., 1418., 1336., 1242.],\n", + " ...,\n", + " [ 811., 655., 688., ..., 385., 362., 461.],\n", + " [ 798., 675., 727., ..., 394., 415., 402.],\n", + " [ 888., 673., 642., ..., 403., 454., 393.]],\n", + "\n", + " [[ 8656., 8656., 8864., ..., 1500., 1428., 1220.],\n", + " [ 9016., 9160., 9224., ..., 1546., 1522., 1360.],\n", + " [ 9248., 9328., 9384., ..., 1620., 1542., 1482.],\n", + " ...,\n", + " [ 1010., 831., 853., ..., 277., 276., 336.],\n", + " [ 1016., 930., 927., ..., 276., 317., 293.],\n", + " [ 1112., 885., 827., ..., 299., 369., 293.]],\n", + "\n", + " [[ 8416., 8416., 8640., ..., 1598., 1466., 1138.],\n", + " [ 8744., 8880., 8928., ..., 1498., 1522., 1284.],\n", + " [ 8952., 8944., 8960., ..., 1542., 1478., 1448.],\n", + " ...,\n", + "...\n", + " [ 652., 640., 638., ..., 590., 821., 1008.],\n", + " [ 622., 676., 630., ..., 606., 1092., 726.],\n", + " [ 864., 786., 569., ..., 766., 1068., 630.]],\n", + "\n", + " [[ 201., 213., 195., ..., 1138., 1058., 749.],\n", + " [ 196., 198., 169., ..., 861., 784., 768.],\n", + " [ 216., 178., 191., ..., 870., 806., 820.],\n", + " ...,\n", + " [ 857., 838., 846., ..., 622., 800., 1332.],\n", + " [ 922., 848., 771., ..., 786., 1046., 912.],\n", + " [ 1118., 1010., 735., ..., 755., 977., 686.]],\n", + "\n", + " [[ 3264., 3352., 3304., ..., 3160., 3296., 3376.],\n", + " [ 3356., 3300., 3212., ..., 3188., 3272., 3064.],\n", + " [ 3288., 3372., 3344., ..., 3136., 3200., 2932.],\n", + " ...,\n", + " [ 1320., 1468., 1298., ..., 2492., 2556., 3018.],\n", + " [ 1630., 1694., 1250., ..., 2318., 2684., 2894.],\n", + " [ 2190., 2072., 1288., ..., 2544., 2942., 2928.]]]],\n", + " dtype=float32)\n", + "Coordinates: (12/53)\n", + " * time (time) datetime64[ns] 96B 2018-0...\n", + " id (time) " + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAHSCAYAAAAkMCseAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA4IUlEQVR4nO3de3hU5bn+8XtIYAiBhGMOmHAQ5SyigBAqiLKJglLZgIVNq9AqlYqiZAdtsFalP5tW1OIRNy0gFi1eJVBQsAUrAVqgEiSonItIQkg4iGYANYfJ+/sjzciQAAlmDu/M93Ndc8mseVfyPK9ZWXfWrLXGYYwxAgAAsESDQBcAAABQF4QXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVIgNdQH2rqKjQkSNH1KxZMzkcjkCXAwAAasEYo1OnTqlt27Zq0ODCx1ZCLrwcOXJEycnJgS4DAABcgvz8fCUlJV1wTMiFl2bNmkmqbD4mJibA1QAAgNpwuVxKTk727McvJOTCS9VbRTExMYQXAAAsU5tTPjhhFwAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwSsjdpM5XSsvcemXVRh04WqhO8Ym679ZBatQwItBlAQAQdggvtfDwwmV6bteDcjc9XLmgSEr/Z5LSuj+vp388OrDFAQAQZnjb6CIeXrhMsw+NlTv6sNdyd3SBZh8aq4cXLgtQZQAAhCfCywWUlrn13K4HJRnp3I9acBhJ0nO7HlJpmdvvtQEAEK58Gl7mzp2rXr16eT4kMSUlRe++++4F11m/fr369Omjxo0b6/LLL9err77qyxIv6JVVGyvfKjrfZ0Q5jNxN8/XKqo1+rQsAgHDm0/CSlJSk3/zmN8rJyVFOTo5uuukm3X777dq5c2eN4w8ePKgRI0Zo0KBB2r59u2bOnKlp06YpKyvLl2We14GjhfU6DgAAfHc+PWF35MiRXs+feuopzZ07V1u2bFGPHj2qjX/11VfVrl07zZkzR5LUrVs35eTk6JlnntGYMWN8WWqNOsUnSkW1HAcAAPzCb+e8uN1uLVmyRGfOnFFKSkqNYzZv3qzU1FSvZTfffLNycnJUVlZW4zolJSVyuVxej/py362DFHE6STLned/IOBRxOln33Tqo3r4nAAC4MJ+Hl48//lhNmzaV0+nUlClTtHz5cnXv3r3GsUVFRYqPj/daFh8fr/Lycp04caLGdTIzMxUbG+t5JCcn11vtjRpGKK3785VPzg0w/3me1n0O93sBAMCPfB5eunTpotzcXG3ZskU/+9nPNHHiRO3ateu84x0O75BgjKlxeZWMjAwVFxd7Hvn5+fVXvKSnfzxaM9ovVcSZy7yWR5xJ0oz2S7nPCwAAfubzm9Q1atRIV1xxhSSpb9++2rp1q55//nn93//9X7WxCQkJKiryPsnk2LFjioyMVKtWrWr8+k6nU06ns/4LP8vTPx6t/1d2O3fYBQAgCPj9DrvGGJWUlNT4WkpKit5++22vZWvWrFHfvn3VsGFDf5R3Xo0aRuihUUMCWgMAAPDx20YzZ87Uxo0b9dlnn+njjz/Wo48+quzsbP3whz+UVPmWz1133eUZP2XKFB06dEhpaWnavXu3FixYoPnz5ys9Pd2XZQIAAIv49MjL0aNHdeedd6qwsFCxsbHq1auX/vrXv2rYsGGSpMLCQuXl5XnGd+zYUatXr9b06dP18ssvq23btnrhhRcCcpk0AAAITg5TdUZsiHC5XIqNjVVxcbFiYmICXQ4AAKiFuuy/+WwjAABgFcILAACwit+vNkJ4Ky1zc8k5AOA7IbzAbx5euEzP7Xqw8pO6JalISv9nktK6P8/N/gAAtcbbRvCLhxcu0+xDY+WOPuy13B1doNmHxurhhcsCVBkAwDaEF/hcaZlbz+16UJKRzv2UB0flxW7P7XpIpWVuv9cGALAP4QU+98qqjZVvFZ3nw7nlMHI3zdcrqzb6tS4AgJ0IL/C5A0cL63UcACC8EV7gc53iE+t1HAAgvBFe4HP33TpIEaeTJHOe942MQxGnk3XfrYP8WxgAwEqEF/hco4YRSuv+fOWTcwPMf56ndZ/D/V4AALVCeIFfPP3j0ZrRfqkizlzmtTziTJJmtF/KfV4AALXGBzPCr7jDLgCgJnXZf3OHXfhVo4YRemjUkECXAQCwGG8bAQAAqxBeAACAVQgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVCC8AAMAqhBcAAGAVwgsAALAK4QUAAFiF8AIAAKxCeAEAAFYhvAAAAKsQXgAAgFUILwAAwCo+DS+ZmZnq16+fmjVrpri4OI0aNUp79+694DrZ2dlyOBzVHnv27PFlqQAAwBI+DS/r16/X1KlTtWXLFq1du1bl5eVKTU3VmTNnLrru3r17VVhY6HlceeWVviwVAABYItKXX/yvf/2r1/OFCxcqLi5O27Zt0+DBgy+4blxcnJo3b+7D6gAAgI38es5LcXGxJKlly5YXHXvNNdcoMTFRQ4cO1bp16847rqSkRC6Xy+sBAABCl9/CizFGaWlpuv7669WzZ8/zjktMTNS8efOUlZWlZcuWqUuXLho6dKg2bNhQ4/jMzEzFxsZ6HsnJyb5qAQAABAGHMcb44xtNnTpVq1at0j/+8Q8lJSXVad2RI0fK4XBo5cqV1V4rKSlRSUmJ57nL5VJycrKKi4sVExPznesGAAC+53K5FBsbW6v9t1+OvDzwwANauXKl1q1bV+fgIkkDBgzQ/v37a3zN6XQqJibG6wEAAEKXT0/YNcbogQce0PLly5Wdna2OHTte0tfZvn27EhMT67k6AABgI5+Gl6lTp+rNN9/UihUr1KxZMxUVFUmSYmNjFRUVJUnKyMhQQUGBXn/9dUnSnDlz1KFDB/Xo0UOlpaVavHixsrKylJWV5ctSAQCAJXwaXubOnStJGjJkiNfyhQsXatKkSZKkwsJC5eXleV4rLS1Venq6CgoKFBUVpR49emjVqlUaMWKEL0sFAACW8NsJu/5SlxN+AABAcAi6E3YBAADqC+EFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVCC8AAMAqhBcAAGAVwgsAALAK4QUAAFiF8AIAAKxCeAEAAFYhvAAAAKsQXgAAgFUILwAAwCqEFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVCC8AAMAqhBcAAGAVn4aXzMxM9evXT82aNVNcXJxGjRqlvXv3XnS99evXq0+fPmrcuLEuv/xyvfrqq74sEwAAWMSn4WX9+vWaOnWqtmzZorVr16q8vFypqak6c+bMedc5ePCgRowYoUGDBmn79u2aOXOmpk2bpqysLF+WCgAALOEwxhh/fbPjx48rLi5O69ev1+DBg2sc88gjj2jlypXavXu3Z9mUKVO0Y8cObd68+aLfw+VyKTY2VsXFxYqJiam32gEAgO/UZf/t13NeiouLJUktW7Y875jNmzcrNTXVa9nNN9+snJwclZWVVRtfUlIil8vl9QAAAKHLb+HFGKO0tDRdf/316tmz53nHFRUVKT4+3mtZfHy8ysvLdeLEiWrjMzMzFRsb63kkJyfXe+0AACB4+C283H///froo4/0pz/96aJjHQ6H1/Oqd7bOXS5JGRkZKi4u9jzy8/Prp2AAABCUIv3xTR544AGtXLlSGzZsUFJS0gXHJiQkqKioyGvZsWPHFBkZqVatWlUb73Q65XQ667VeAAAQvHx65MUYo/vvv1/Lli3T+++/r44dO150nZSUFK1du9Zr2Zo1a9S3b181bNjQV6UCAABL+DS8TJ06VYsXL9abb76pZs2aqaioSEVFRfr66689YzIyMnTXXXd5nk+ZMkWHDh1SWlqadu/erQULFmj+/PlKT0/3ZakAAMASPg0vc+fOVXFxsYYMGaLExETP46233vKMKSwsVF5enud5x44dtXr1amVnZ6t379761a9+pRdeeEFjxozxZakAAMASfr3Piz9wnxcAAOwTtPd5AQAA+K4ILwAAwCqEFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVCC8AAMAqhBcAAGAVwgsAALAK4QUAAFiF8AIAAKxCeAEAAFYhvAAAAKsQXgAAgFUILwAAwCqEFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVXwaXjZs2KCRI0eqbdu2cjgc+stf/nLB8dnZ2XI4HNUee/bs8WWZAADAIpG+/OJnzpzR1VdfrR//+McaM2ZMrdfbu3evYmJiPM/btGnji/IAAICFfBpehg8fruHDh9d5vbi4ODVv3rz+CwIAANYLynNerrnmGiUmJmro0KFat27dBceWlJTI5XJ5PQAAQOgKqvCSmJioefPmKSsrS8uWLVOXLl00dOhQbdiw4bzrZGZmKjY21vNITk72Y8UAAMDfHMYY45dv5HBo+fLlGjVqVJ3WGzlypBwOh1auXFnj6yUlJSopKfE8d7lcSk5OVnFxsdd5MwAAIHi5XC7FxsbWav8dVEdeajJgwADt37//vK87nU7FxMR4PQAAQOgK+vCyfft2JSYmBroMAAAQJHx6tdHp06f173//2/P84MGDys3NVcuWLdWuXTtlZGSooKBAr7/+uiRpzpw56tChg3r06KHS0lItXrxYWVlZysrK8mWZAADAIj4NLzk5Obrxxhs9z9PS0iRJEydO1GuvvabCwkLl5eV5Xi8tLVV6eroKCgoUFRWlHj16aNWqVRoxYoQvywQAABbx2wm7/lKXE34AAEBwCKkTdgEAAM5GeAEAAFYhvAAAAKsQXgAAgFUILwAAwCqEFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVCC8AAMAqhBcAAGAVwgsAALAK4QUAAFiF8AIAAKxCeAEAAFYhvAAAAKsQXgAAgFUILwAAwCqEFwAAYBXCCwAAsArhBQAAWMWn4WXDhg0aOXKk2rZtK4fDob/85S8XXWf9+vXq06ePGjdurMsvv1yvvvqqL0sEAACW8Wl4OXPmjK6++mq99NJLtRp/8OBBjRgxQoMGDdL27ds1c+ZMTZs2TVlZWb4sEwAAWCTSl198+PDhGj58eK3Hv/rqq2rXrp3mzJkjSerWrZtycnL0zDPPaMyYMT6qEgAA2CSoznnZvHmzUlNTvZbdfPPNysnJUVlZWYCqAgAAwcSnR17qqqioSPHx8V7L4uPjVV5erhMnTigxMbHaOiUlJSopKfE8d7lcPq8TAAAETlAdeZEkh8Ph9dwYU+PyKpmZmYqNjfU8kpOTfV4jAAAInKAKLwkJCSoqKvJaduzYMUVGRqpVq1Y1rpORkaHi4mLPIz8/3x+lAgCAAAmqt41SUlL09ttvey1bs2aN+vbtq4YNG9a4jtPplNPp9Ed5AAAgCPj0yMvp06eVm5ur3NxcSZWXQufm5iovL09S5VGTu+66yzN+ypQpOnTokNLS0rR7924tWLBA8+fPV3p6ui/LBAAAFvHpkZecnBzdeOONnudpaWmSpIkTJ+q1115TYWGhJ8hIUseOHbV69WpNnz5dL7/8stq2basXXniBy6QBAICHw1SdERsiXC6XYmNjVVxcrJiYmECXAwAAaqEu+++gOmEXAADgYggvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVCC8AAMAqhBcAAGAVwgsAALAK4QUAAFiF8AIAAKxCeAEAAFYhvAAAAKsQXgAAgFUILwAAwCqEFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBV/BJeXnnlFXXs2FGNGzdWnz59tHHjxvOOzc7OlsPhqPbYs2ePP0oFAABBzufh5a233tJDDz2kRx99VNu3b9egQYM0fPhw5eXlXXC9vXv3qrCw0PO48sorfV0qAACwgM/Dy3PPPae7775b99xzj7p166Y5c+YoOTlZc+fOveB6cXFxSkhI8DwiIiJ8XSoAALCAT8NLaWmptm3bptTUVK/lqamp2rRp0wXXveaaa5SYmKihQ4dq3bp15x1XUlIil8vl9QAAAKHLp+HlxIkTcrvdio+P91oeHx+voqKiGtdJTEzUvHnzlJWVpWXLlqlLly4aOnSoNmzYUOP4zMxMxcbGeh7Jycn13gcAAAgekf74Jg6Hw+u5MabasipdunRRly5dPM9TUlKUn5+vZ555RoMHD642PiMjQ2lpaZ7nLpeLAAMAQAjz6ZGX1q1bKyIiotpRlmPHjlU7GnMhAwYM0P79+2t8zel0KiYmxusBAABCl0/DS6NGjdSnTx+tXbvWa/natWs1cODAWn+d7du3KzExsb7LAwAAFvL520ZpaWm688471bdvX6WkpGjevHnKy8vTlClTJFW+7VNQUKDXX39dkjRnzhx16NBBPXr0UGlpqRYvXqysrCxlZWX5ulQAAGABn4eXcePG6fPPP9esWbNUWFionj17avXq1Wrfvr0kqbCw0OueL6WlpUpPT1dBQYGioqLUo0cPrVq1SiNGjPB1qQAAwAIOY4wJdBH1yeVyKTY2VsXFxZz/AgCAJeqy/+azjQAAgFUILwAAwCqEFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVCC8AAMAqhBcAAGAVwgsAALAK4QUAAFiF8AIAAKxCeAEAAFYhvAAAAKsQXgAAgFUILwAAwCqEFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVfwSXl555RV17NhRjRs3Vp8+fbRx48YLjl+/fr369Omjxo0b6/LLL9err77qjzIBAIAFfB5e3nrrLT300EN69NFHtX37dg0aNEjDhw9XXl5ejeMPHjyoESNGaNCgQdq+fbtmzpypadOmKSsry9elAgAACziMMcaX36B///669tprNXfuXM+ybt26adSoUcrMzKw2/pFHHtHKlSu1e/duz7IpU6Zox44d2rx580W/n8vlUmxsrIqLixUTE1M/TQAAAJ+qy/7bp0deSktLtW3bNqWmpnotT01N1aZNm2pcZ/PmzdXG33zzzcrJyVFZWVm18SUlJXK5XF4PAAAQunwaXk6cOCG32634+Hiv5fHx8SoqKqpxnaKiohrHl5eX68SJE9XGZ2ZmKjY21vNITk6uvwYAAEDQ8csJuw6Hw+u5MabasouNr2m5JGVkZKi4uNjzyM/Pr4eKAQBAsIr05Rdv3bq1IiIiqh1lOXbsWLWjK1USEhJqHB8ZGalWrVpVG+90OuV0OuuvaAAAENR8euSlUaNG6tOnj9auXeu1fO3atRo4cGCN66SkpFQbv2bNGvXt21cNGzb0Wa0AAMAOPn/bKC0tTX/4wx+0YMEC7d69W9OnT1deXp6mTJkiqfJtn7vuusszfsqUKTp06JDS0tK0e/duLViwQPPnz1d6erqvSwUAABbw6dtGkjRu3Dh9/vnnmjVrlgoLC9WzZ0+tXr1a7du3lyQVFhZ63fOlY8eOWr16taZPn66XX35Zbdu21QsvvKAxY8b4ulQAAGABn9/nxd+4zwsAAPYJmvu8AAAA1DfCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwCuEFAABYhfACAACsQngBAABWIbwAAACrEF4AAIBVCC8AAMAqhBcAAGAVwgsAALAK4QUAAFiF8AIAAKxCeAEAAFaJDHQBAOrO7ZY2bpQKC6XERGnQICkiItBVAYB/EF5CCDu0+hHs87hsmfTgg9Lhw98uS0qSnn9eGj06cHWdK9jnUbKjRgDV8bZRiFi2TOrQQbrxRmnChMr/duhQuRy1F+zzuGyZNHasd3CRpIKCyuXBVGcwz6NkR40AauYwxphAF1GfXC6XYmNjVVxcrJiYmECX4xdVO7Rz/086HJX/Xbo0uP4iD1bBPo9ud+XO9dzgUsXhqDwCc/BgYI8eBPs8SnbUCISbuuy/CS+Ws2WHFuxsmMfs7MqjAxezbp00ZIivq6mZDfNoQ41AOKrL/pu3jSy3ceP5fwlLlX9Z5udXjsP52TCPhYX1O84XbJhHG2oEcGGEF8vZsEOzgQ3zmJhYv+N8wYZ5tKFGABdGeLGcDTs0G9gwj4MGVb6dUXVexrkcDik5uXJcoNgwjzbUCODCCC+Ws2GHZgMb5jEiovJy6Kp6zlb1fM6cwJ6nYcM82lAjgAsjvFjOhh2aDWyZx9GjK6+Euewy7+VJScFxhYwN82hDjQAujPASAoJ9h2YLW+Zx9Gjps88qryp6883K/x48GFz1Bfs82lAjgPPz6aXSX3zxhaZNm6aVK1dKkr7//e/rxRdfVPPmzc+7zqRJk7Ro0SKvZf3799eWLVtq9T3D7VLps3G30PrBPNYPG+bRhhqBcBE093kZPny4Dh8+rHnz5kmSfvrTn6pDhw56++23z7vOpEmTdPToUS1cuNCzrFGjRmrZsmWtvmc4hxcAAGxVl/23zz7baPfu3frrX/+qLVu2qH///pKk3//+90pJSdHevXvVpUuX867rdDqVkJDgq9IAAIDFfHbOy+bNmxUbG+sJLpI0YMAAxcbGatOmTRdcNzs7W3FxcercubMmT56sY8eOnXdsSUmJXC6X1wMAAIQun4WXoqIixcXFVVseFxenoqKi8643fPhwvfHGG3r//ff17LPPauvWrbrppptUUlJS4/jMzEzFxsZ6HsnJyfXWAwAACD51Di9PPPGEHA7HBR85OTmSJEcNN1IwxtS4vMq4ceN06623qmfPnho5cqTeffdd7du3T6tWrapxfEZGhoqLiz2P/Pz8urYEAAAsUudzXu6//36NHz/+gmM6dOigjz76SEePHq322vHjxxUfH1/r75eYmKj27dtr//79Nb7udDrldDpr/fUAAIDd6hxeWrdurdatW190XEpKioqLi/XBBx/ouuuukyT961//UnFxsQYOHFjr7/f5558rPz9fidyrGwAAyIfnvHTr1k233HKLJk+erC1btmjLli2aPHmybrvtNq8rjbp27arly5dLkk6fPq309HRt3rxZn332mbKzszVy5Ei1bt1a//3f/+2rUgEAgEV8eofdN954Q1dddZVSU1OVmpqqXr166Y9//KPXmL1796q4uFiSFBERoY8//li33367OnfurIkTJ6pz587avHmzmjVr5stSAQCAJXx6k7pA4CZ1AADYpy77bz7bCAAAWIXwAgAArEJ4AQAAViG8AAAAqxBeAACAVQgvAADAKoQXAABgFcILAACwSp0/2wgA4D/uUrc+fmWjvjpQqCadEnXVfYMU0Sgi0GUBAUV4AYAgteXhZWr33IPq7T7sWXYkPUl5ac9rwNOjA1gZEFi8bQQAQWjLw8t03eyxSjgruEhSgrtA180eqy0PLwtQZUDgEV4AIMi4S91q99yDkky1X9INVPlxdMnPPSR3qdvvtQHBgPACAEHm41c2qq378Hl/QTeQ0WXufH38yka/1gXfcpe6lTsnW5se+JNy52QTTi+Ac14AIMh8daCwXsch+HF+U90QXgAgyDTplFiv43yNK6K+m6rzm/SftwSrJLgLlDB7rLZoKQHmHA5jjLn4MHu4XC7FxsaquLhYMTExgS4HAOrMXerW0SYdlOAu8JzjcrYKOVQYkaSErw4GPCRUHTFoe/YRgwiOGNTWt/+va36bMJj+X/taXfbfnPMCAEEmolGE8tKel1S58zpb1fP8tDkB35nZdEVUsJ5PYtv5TcEyj4QXAAhCA54erQ9mLFVRxGVeywsjkvTBjMC/jWDTFVFbHl6mo006qPf0GzXwpQnqPf1GHW3SISjClU3nNwXTPBJeACBIDXh6tOK/+ky5v1unTfe/qdzfrVPCVwcDHlwke44YBPvRIVvObwq2eeScFwBAnW164E8a+NKEi4+7/00NfPF//FBRdTacT2LD+U3+mkfOeQEA+JQNRwxsODpkw/lNwTiPhBcAQJ1ddd8gHYlIqrbDrVIhhwoiknXVfYP8XNm3bDmfJNjPbwrGeeQ+LwCAOqs6YpAwe6wq5PB6y+PsIwaXBfCIgQ1Hh6oMeHq03P/vduWec7+cQM5flWCcR855AQBcspru81IQkaz8tDkBP2Jgw/kkNvDXPHLOCwDAL4L5iigbziexQTDOI0deAAAhLZiPDtnE1/NYl/034QUAEPL4/KX64ct5JLwQXgAAsArnvAAAgJBFeAEAAFYhvAAAAKsQXgAAgFUILwAAwCo+DS9PPfWUBg4cqCZNmqh58+a1WscYoyeeeEJt27ZVVFSUhgwZop07d/qyTAAAYBGfhpfS0lLdcccd+tnPflbrdZ5++mk999xzeumll7R161YlJCRo2LBhOnXqlA8rBQAAtvBpeHnyySc1ffp0XXXVVbUab4zRnDlz9Oijj2r06NHq2bOnFi1apK+++kpvvvmmL0sFAACWCKpzXg4ePKiioiKlpqZ6ljmdTt1www3atGlTjeuUlJTI5XJ5PQAAQOiKDHQBZysqKpIkxcfHey2Pj4/XoUOHalwnMzNTTz75ZLXlhBgAAOxRtd+uzY3/6xxennjiiRrDwtm2bt2qvn371vVLezgc3p9aaYyptqxKRkaG0tLSPM8LCgrUvXt3JScnX/L3BwAAgXHq1CnFxsZecEydw8v999+v8ePHX3BMhw4d6vplJUkJCQmSKo/AJCYmepYfO3as2tGYKk6nU06n0/O8adOmys/PV7Nmzc4beGzicrmUnJys/Pz8kP2spnDoUQqfPn0tHOaRHkNHuPRZH4wxOnXqlNq2bXvRsXUOL61bt1br1q0vqbCL6dixoxISErR27Vpdc801kiqvWFq/fr1++9vf1uprNGjQQElJST6pL5BiYmJC/gc/HHqUwqdPXwuHeaTH0BEufX5XFzviUsWnJ+zm5eUpNzdXeXl5crvdys3NVW5urk6fPu0Z07VrVy1fvlxS5dtFDz30kH79619r+fLl+uSTTzRp0iQ1adJEEyZM8GWpAADAEj49YfeXv/ylFi1a5HledTRl3bp1GjJkiCRp7969Ki4u9ox5+OGH9fXXX+u+++7TF198of79+2vNmjVq1qyZL0sFAACW8Gl4ee211/Taa69dcMy5ZxU7HA498cQTeuKJJ3xXmEWcTqcef/xxr/N6Qk049CiFT5++Fg7zSI+hI1z69DeHqc01SQAAAEEiqG5SBwAAcDGEFwAAYBXCCwAAsArhBQAAWIXwUk+OHj2qkydPBroMn6JH1Fa4zGM49EmPCEaEl3qwf/9+ff/739f//M//aNeuXYEuxyfoEbUVLvMYDn3SI4IV4aUeXHnllXrmmWfUtWtXDR48OCQ3AHpEbYXLPIZDn/SIYMV9Xr6jiooKNWjwbQZMS0vTzp07tWzZMkVHRwewsvpDj6HRoz+EyzyGQ5/0GBo9hiqOvFyiiooKGWPUoEEDlZWVqaKiQpLUr18/HT9+3PPcZvQYGj36Q7jMYzj0SY+h0WOo8+nHA4Sys9N6w4YNPf/OycnR0aNHq33sgY3oMTR69Idwmcdw6JMeQ6PHUEd4qaOcnBz97W9/U1lZmdxut4wxioiIUEFBgT799FP94x//0MKFC63+6HN6DI0e/SFc5jEc+qTH0OgxXHDOSx2cOXPG8+nWkydPVn5+vtxut+Lj41VRUaFWrVpp0qRJuvrqqz3JvqCgQG63WzExMWrevHkAq68degyNHv0hXOYxHPqkx9DoMawY1MnOnTtNbGysmT179gXHHT9+3IwcOdJ07drVtG/f3vTp08fs2LHDT1V+N/T4LZt79Idwmcdw6JMev2Vzj+GC8HIJtm7dapxOp3nssce8lldUVBhjjNm8ebPp2LGj6d27t1m4cKH529/+ZmbMmGHi4+PNoUOHAlFyndFjaPToD+Eyj+HQJz2GRo/hgPByibZt22YcDof59a9/bYwxpry83BhjzMcff2ySk5PNsGHDzBdffOEZf/LkSXPNNdeYd955JxDlXhJ6DI0e/SFc5jEc+qTH0Ogx1HGp9CW69tprtX37djmdTp06dUoRERGSpAkTJujyyy/XO++8o+bNm6u8vFyS1KJFCx08eFCFhYWBLLtO6DE0evSHcJnHcOiTHkOjx5AX6PRku7KyMs+/X3rpJdOpUyezZ8+eaq89+eSTJjk52Rw5csTvNX5X9BgaPfpDuMxjOPRJj6HRY6jiyMt3FBn57dXmhw8fVosWLdS+fXuv115//XW99957GjdunGJjYyVJX375pVasWKGioiL/F11H9BgaPfpDuMxjOPRJj6HRY6givNQD85+rzaOiopSYmKjGjRt7Xps3b57mz5+vFi1a6L777lOTJk0kSVu2bNGSJUs0duxY7d69OyB11wU9hkaP/hAu8xgOfdJjaPQYkgJ41Cfk7NmzxzRr1sxMnjzZ/P73vzdjx441vXr1MhMmTDD79+83xhjjdrs94zds2GCGDBliHA6HKSoqClTZdUKPodGjP4TLPIZDn/QYGj2GEsJLPcvNzTXDhg0zAwYMMCkpKWbRokWe90nP/sE3xph169aZhIQEc9tttwWi1EtGj6HRoz+EyzyGQ5/0GBo9hgrCiw988803pqSkxGtZ1T0Eqv79/vvvm/j4eHP77bfXOCbY0WNo9OgP4TKP4dAnPYZGj6GA8OJDNf0wu91us27dOhMXF+f1g+92u6384afH0OjRH8JlHsOhT3oMjR5tRnjxsVOnTpns7GzPIcfz/eDbjB4r2d6jP4TLPIZDn/RYyfYebcUHM/rY8ePH9V//9V/q37+/7rjjDv3oRz/SwIEDtXz5ckmVZ7o7HA7P+AULFmjjxo267LLLNGnSJF1xxRWBKr3W6DE0evSHcJnHcOiTHkOjR2sFLDaFka1bt5ro6GjjcDjM6NGjPcvLy8u9DjXOnz/fNGjQwHz/+983t912m4mPjzcfffRRIEquM3oMjR79IVzmMRz6pMfQ6NFG3OfFD/r27asPPvhALVq0UGJioiTJ7XYrIiLCK7Xn5+drwoQJWrp0qd5++23dc889Gj9+vI4fPx6o0muNHkOjR38Il3kMhz7pMTR6tFHkxYegPnTv3l1r1qxRv379dNVVV+mnP/2ppMof+DfeeEP79u1TcXGxoqKi1LBhQ0nSj370I61YsUKFhYVq06ZNIMuvFXoMjR79IVzmMRz6pMfQ6NE2HHnxoz59+ujf//63evXq5UnsjzzyiGbNmqV27dopISFBf/7zn/XrX/9ap0+f1pYtW7Rz584AV1039BgaPfpDuMxjOPRJj6HRo1UC/b5VONuwYYNxOBxm06ZNnmV/+ctfTEREhBkxYoRxOBzm8ccfD1yB9YAeQ6NHfwiXeQyHPukxNHoMZoSXAMrNzTXdunUz27ZtM2VlZcbtdpvi4mLTr18/s23bNpOfnx/oEr8zegyNHv0hXOYxHPqkx9DoMZjxtlEAJSYmKjIyUitWrFBkZKQaNGigTz/9VHv27NHp06eVlJQkqfLkMFvRY2j06A/hMo/h0Cc9hkaPwYz7vARYbm6uhg0bpmHDhsnpdGrTpk3q3r27Fi9erNOnT8sYo4SEBM/Z7VL1ewsEO3oMjR79IVzmMRz6pMfQ6DFYceQlwHr37q0NGzaoVatWOnTokMaPH6+MjAxFR0frjTfeUO/evXXo0CFFRESotLRUkjw/+BUVFVakenoMjR79IVzmMRz6pMfQ6DFo+f+dKtTE7XabsrIyr2VlZWVmypQp5rLLLjOfffaZMcaYr7/+2qxatcpMmDDBDB061IwfP97s3LkzECXXGT2GRo/+EC7zGA590mNo9BhsuM9LkGjQoIEaNPA+EBYZGalXXnlFTZs21datW9W6dWtlZGToH//4hxo0aKDbbrtN+/bt08CBA/Wvf/1LXbp0CVD1tUOPodGjP4TLPIZDn/QYGj0GnUCnJ9ROSUmJmTJliunXr5+ZNWuW122px40bZ9LT063/gDB6DI0e/SFc5jEc+qTH0OjR3zjyYgG3263nnntOH374ocaNG6dp06bJ4XDI/Odc6xMnTqhp06bVkr9N6DE0evSHcJnHcOiTHkOjx0Bgtizw9ddf66233tLVV1+tqVOnqmHDhnK73XI4HDp27JhOnTqldu3aqaKiQpJUVlYW4Irrjh5Do0d/qM08Jicnyxhj9QmR4dAnPYZGjwERoCM+qIPly5ebiIgI8+WXXxpjjNeJYT/5yU9Mq1atTFFRkTHGmBMnTpjJkydb92mm9BgaPfpDbeaxsLDQa52zD9PbIhz6pMfQ6DEQOPJigdatW6t79+46c+aMpMoTwSoqKjR16lQtWbJEixcvVnx8vNxut1q1aqXk5GT17t1bO3bsCHDltUePodGjP9RmHhMSElRUVKSlS5dq0qRJmjVrlvbv3x/gyusmHPqkx9DoMSACnZ5wcXl5eSYxMdHMnDnTfPLJJ2bVqlVm1KhRpmnTpmbp0qU1pvSZM2eaJk2amO3bt/u/4EtAj6HRoz+cbx6jo6PN0qVLTXl5uTGm8i/ejh07mpSUFPODH/zAtG7d2qojWeHQJz2GRo+BQHixxLZt20znzp1Nhw4dTIcOHcy1115rPvzww2rjzj5jffbs2SY6Otrs3r3bn6VeMnr8ls09+sO589i7d2+Tk5PjNeYXv/iFGTp0qOf5Y489Znr16mWOHj3q73IvWTj0SY+VbO/R3/h4AIscP35chYWFcjqdateunaKiojy3mjbn3HL6+PHj2rVrl2688UbFx8drx44diouLC2D1tUOPodGjP5w9j0lJSYqOjtaOHTuUnZ2tf/7zn4qOjlaDBg00f/58SdLWrVv1ox/9SH/+85/Vq1evAFdfe+HQJz2GRo/+xKXSFmnTpo3atGnjeV5eXu55/7RBgwYqLy9XeXm5nnnmGa1bt07Z2dm68847de2111qzw6PH0OjRH86dxyNHjmjWrFn68ssvdfnllys+Pl7PPPOMrrzySk2cOFEfffSRXC6XTp06FcCq6y4c+qTH0OjRrwJ30Af1pby83Hz66afmvvvuMz169DCxsbFmypQp5m9/+1u1cbaiR+9xqC4nJ8e0adPGvP76655l77zzjmnSpIm54YYbTLdu3cyPf/zjAFZYP8KhT3oMjR59iSMvlnv99de1bds2/eEPf1DPnj1100036bHHHlPz5s3VsGFDr7FVn2pqG3oMjR59rbS0VC1atFDfvn0lVX5671VXXaWOHTvqqaeeUlJSktq3by9JXp/ya5tw6JMeQ6NHX+JSaYsdP35cM2bM0I4dO/SrX/1Kf//73zVnzhy1adPGa4dXddMzG9FjaPToD506dZLT6dSiRYskVX56b8OGDVVcXKyvvvrKsyMoKyuzekcQDn3SY2j06EucsGu5kydPqqysTPHx8Z5l5j8nfebl5Sk6OlqtWrWqltzNOSeGBrNL7fHsccHuu/TIX2Xfys3N1bBhw3TzzTcrOjpa2dnZatKkiXJycnT8+HE1bNjQ+u1BCo8+L7VHyZ4+v0uP4b7dc+TFci1btvTs8KpyaNVG+/7776tTp046dOiQIiIiVFpa6rVuRUWFFbegr0uPZ99eu6ioSO+9956Kior8X3Qd1aXHs4/ALF++XLNnz9ahQ4f8X3QQ6t27tzZu3KioqCgdOXJEKSkp2rBhgyIiIrRq1aqQ2B6k8OizLj3aut3XpUe2+3P4/Swb+E1FRYW59957TfPmzc3hw4c9y5csWWJ+9rOfme9973vmzjvvNDt37gxgld9NVY8tWrQwBQUFXq9t2LDBjB071gwYMCBke3zkkUfMmDFjTKdOnbih1VlqOqnZ7XaH3PYQDn1eqMdQ2e7r2iPbPTepC3lut9s89thjZvny5cYYYzIyMkxUVJS54447zNSpU83dd99tmjdvbvbs2RPYQr8Dt9ttnnjiCbNs2TKv5RUVFebAgQPme9/7nmnZsqU5fvx4gCr87s7t8ey78Z48edL8/Oc/N+3atTMHDx4MUIXB6+wb/oXy9hAOfZ7bYyhu9xfqke3+W4SXEHf2D/vDDz9snE6n+cMf/mBOnjzpWT5y5Ejz85//PBDl1auSkhLPv0tLS40xxvz97383jRs3Nj/5yU9Mfn5+oEqrN2f3eO7yIUOGmJdfftnPFdklXLaHcOnTGLb7cN3uuVQ6xFWdNzF79mw9++yz+uMf/6jx48dXO5nts88+C0B19atRo0YyxqiiokINGzbUe++9p+HDh3s+6CwxMTHQJV6SqpvXGWPUqFEjSd/e2E6SZ7nL5dKnn34ayFKDXrhsD+HSp8R2H7bbfSCTE/xjz549plu3buapp54y33zzjddrubm5plOnTmbevHkBqq5+VR1yfe+990xkZKS55557zJEjRzyvn/1x9LYpKioyZ86cMadOnTIul8ucOXPGfPXVV2b37t1m9uzZpkGDBmbNmjWBLjPohcv2EC59GsN2H47bPUdewkBBQYEKCwuVmpoqp9PpWf7ZZ5/ptddeU+vWrdWjR48AVvjdfPPNN/rkk0/Us2dPNW7cWO+//75uueWWan95ud1uz18tf/7zn3Xdddd57qUQ7ObNm6dp06apc+fOcrlcioqK0jfffKPS0lI1btxYp06d0ty5c3XDDTdIqux17dq16t27txISEgJcfXCpzfZw9dVXe5bv27dPcXFxat68eQCqvXRs92z3oYzwEgbKysrUvn17XXnllZ5le/fu1YIFC7R8+XLNmDFDAwcOlCQ9/fTTKigoUFRUlO655x5dccUVgSq71r7++muNHj1aY8eO1S233KJbb71VkyZN0pNPPun1C6zqngiLFi3Sb37zGxUWFmr9+vVeO6pg1b9/f0VGRqpFixZatGiRIiMjderUKVVUVOiyyy5TdHS01+ceLViwQHPnzlVkZKQWLFignj17BrD64HKx7eGRRx7RddddJ0navXu3Fi5cqP379+u3v/2tOnfuHKiy64ztnu0+pLf7QB/6ge8dPnzYJCQkmPT0dLNr1y6zYsUKM3bsWNO1a1fz7LPPesb99re/NV27djWpqalmzJgxpk2bNtZchvevf/3LxMbGGofDUe2Q8dkWLFhgBg8ebCZPnmwmTpxooqOjTW5urp+rvTQ7duwwjRo1Mr/85S89y869MuFse/bsMf/7v/9rWrdubd1VJb50oe3hd7/7XbXxK1euNBMmTDCNGjUyBw4c8H/Bl4jt/lts96GH8BImPvzwQ9OpUyfTo0cPExkZaX7wgx+YxYsXe43JzMw0o0aNMp9//rkxxpgZM2aYHj16WHOp4a5du0zz5s3ND3/4wxrvm5Cbm2smTJhgWrRo4dmo09PTTePGja35ZZ2Tk2Oio6PNI488ct4x5/4yu/vuu83EiROtft+/vtW0PSxZssTzenl5udc83nXXXcbhcJj3338/EOVeMrZ7tvtQ3e4JL2Hk2LFjZt++fWbfvn3mzJkznuXp6enmwQcfNGPGjPH6FNMDBw6Yrl27mg8++CAQ5V6Sbdu2maZNm5oXX3yx2i+yb775xmzdutUMHz7cDBs2zLNRP/7442bGjBmBKPeSfPDBB6Z58+bm73//u9fy0tJSzy+wqktGjam8VHbYsGF+rdEGZ28PZ1+Keu7Pzb333muio6PN6tWr/V1ivWC7Z7sPRYSXMHfDDTeYpKQkM2PGDDN9+nQTGRlpnnzySeN2u82sWbNMmzZtrLlLZZV9+/aZtWvXGmOq/zVijDFHjx41gwcPNk899ZS/S6s3R44cMYcOHbrgmIqKCnP8+HGTkpJihg8f7qfK7Hb2/VGM+Ta4vPPOOwGqyDfY7u3Edv8tTtgNY2vXrtWHH36o1atX6/rrr5ckfe9739Mdd9yhTZs2KTc3V/fee6/atm0b4Err5sorr/ScpPj444+rf//+uu222yRV3juhTZs2iouL04EDB7zWM5Z8mJskzwmJJ0+e1PTp0+V2u3X69GkZY9SgQQNVVFSosLBQJ0+elDFGS5YskWRXj/525MgRxcXFea5Muffee/XGG2/orbfe0q233hrg6uoP2z3bfSggvISxVq1aKT4+3nMDJEkaOnSounfvrl/84hfq3LmzmjVrpqioqABWeek+//xzvfvuuzpz5oxuu+02zwYuSSdOnKj24XQ2btynTp3SH//4RzkcDj377LM6fvy4SktL1bJlSzmdTrVp00bjxo3z/D+2sUd/+Oabb/TUU09Jkn73u99p2rRpeuONN7RkyZKQCi4S2z3bfWggvISxyy67TE2bNtWqVas8l4bu2LFDpaWlatKkidcleDZq1aqVFi5cqOuvv17R0dGaOHGirrjiCmVkZGj//v1asGCBJLv/Kmnfvr0++eQTXXfddTLGeHbA56q6Wydq1rhxY919993q16+f3n33XX3xxRchd8SlCts9230oILyEsfj4eC1cuFA33XSTDh8+rNOnT+vLL79UmzZt1KZNm0CXVy+uuuoqrV+/XpMmTdKKFSt08OBBOZ1Opaenq2/fvpLs/6uke/fuys7O1uDBg/XFF19o1qxZ1caE6i+w+nTttddqx44dGjJkiK677jrdcsstgS7JJ9ju2e5DgcMYYwJdBAJr165devnll3XgwAH16tVL48eP17XXXhvosurVkSNHtHfvXu3atUs333yzEhIS1LRp00CXVa8+/PBD9e3bVy+++KKmTp0a6HKstWPHDl1zzTXKyMjQr371q5DdAbDdh4Zw3e4JL5D07Qd+nf3BX7DPrl27dOLECQ0aNMj6vywD6eOPP9aKFSs0bdo0xcTEBLocn2G7Dw3huN0TXuDF5veBgfoUTjt0tnvYhvACAACsEppv5gIAgJBFeAEAAFYhvAAAAKsQXgAAgFUILwAAwCqEFwAAYBXCCwAAsArhBQAAWIXwAgAArEJ4AQAAViG8AAAAq/x/0DbJkIHviEUAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Run PCA\n", + "pca = decomposition.PCA(n_components=1)\n", + "pca_result = pca.fit_transform(embeddings)\n", + "\n", + "plt.xticks(rotation=-45)\n", + "\n", + "# Plot all points in blue first\n", + "plt.scatter(stack.time, pca_result, color=\"blue\")\n", + "\n", + "# Re-plot cloudy images in green\n", + "plt.scatter(stack.time[0], pca_result[0], color=\"green\")\n", + "plt.scatter(stack.time[2], pca_result[2], color=\"green\")\n", + "\n", + "# Color all images after fire in red\n", + "plt.scatter(stack.time[-5:], pca_result[-5:], color=\"red\")" + ] + }, + { + "cell_type": "markdown", + "id": "b38b70a6-2156-41f8-967e-a490cc8e2778", + "metadata": {}, + "source": [ + "### And finally, some finetuning\n", + "\n", + "We are going to train a classifier head on the embeddings and use it to detect fires." + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "1da07de0-b8f2-46c9-bd2a-58b15ca2224f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Matched 5 out of 5 correctly\n" + ] + } + ], + "source": [ + "# Label the images we downloaded\n", + "# 0 = Cloud\n", + "# 1 = Forest\n", + "# 2 = Fire\n", + "labels = np.array([0, 1, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2])\n", + "\n", + "# Split into fit and test manually, ensuring we have all 3 classes in both sets\n", + "fit = [0, 1, 3, 4, 7, 8, 9]\n", + "test = [2, 5, 6, 10, 11]\n", + "\n", + "# Train a support vector machine model\n", + "clf = svm.SVC()\n", + "clf.fit(embeddings[fit] + 100, labels[fit])\n", + "\n", + "# Predict classes on test set\n", + "prediction = clf.predict(embeddings[test] + 100)\n", + "\n", + "# Perfect match for SVM\n", + "match = np.sum(labels[test] == prediction)\n", + "print(f\"Matched {match} out of {len(test)} correctly\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "claymodel", + "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.8" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From 514215bfc20fa8af3ddcf47a1207489f06641251 Mon Sep 17 00:00:00 2001 From: Daniel Wiesmann Date: Fri, 24 May 2024 10:04:06 +0100 Subject: [PATCH 2/2] Add notebook to toc --- docs/_config.yml | 1 + docs/_toc.yml | 2 ++ 2 files changed, 3 insertions(+) diff --git a/docs/_config.yml b/docs/_config.yml index 10147fa0..2e3abae7 100644 --- a/docs/_config.yml +++ b/docs/_config.yml @@ -15,6 +15,7 @@ execute: - partial-inputs.ipynb - tutorial_digital_earth_pacific_patch_level.ipynb - patch_level_cloud_cover.ipynb + - clay-v1-wall-to-wall.ipynb # Define the name of the latex output file for PDF builds latex: diff --git a/docs/_toc.yml b/docs/_toc.yml index 1704f2ec..f6fb351a 100644 --- a/docs/_toc.yml +++ b/docs/_toc.yml @@ -26,6 +26,8 @@ parts: file: data_sampling - caption: Running the model chapters: + - title: Clay v1 wall-to-wall example + file: clay-v1-wall-to-wall - title: Run over a region file: run_region - title: Generating embeddings